Odd/Even matrix size fix, DC fix, new padding features, test suite (#41)

* Checkpoint

* Refactor: use numpy padding, add edge as default, 50 pixels per side by default

* Update readme

* Remove local files

* Customize for ISMRM

* Update gitignore

* Incomplete checkpoint - fixing fft

* Fix

* Decide on nan replacement implementtion

* Fix bug

* Add savefig

* Create some initial integration tests - contains some failures to be resolved later (#38) (#42) (#43)

* Refactor analytical comparison function to have more explicit args  (#36)

* Refactor function args for geometry examples

* Reorgnize test folder

* Write tests for compare_to_analytical and create internal python function (non-cli)

* Add gitignore cases

* Add integration tests

* Change on push

* Add name

* Change name

* Try coveralls setup

* Try coveralls setup

* Try coveralls setup

* Add expected fails

* Add xfails

* Add coverage report

* Add coverage report

* Add coverage report

* Add badges

* Fix fftshift bug by using ifftshift

* Fix floating-point precision issues in compute_fieldmap tests

* Remove xfail test tag for test resolved by ifftshift fix

* Add boundary exclusion mask to explude Gibbs ringing when testing against analytical solution

* Fix buffer=0 edge case in compute_bz function

* Test 50 px buffer

* Fix meshgrid indexing for non-square matrices in Spherical class

* Add indexing='ij' to Cylindrical class meshgrid calls

* Add comprehensive Cylindrical geometry test coverage

* Fix incorrect DC component adjustment and add validation tests

* Cleanup

* Add tests for neighbouring voxel fields (all three directions) of a single voxel of water

* Update coverage badge in README.md

* Update coverage badge to reflect custom branch

* Update README.md to reflect custom branch status

* Fix badge link in README for test workflow

Author: mathieuboudreau

.github/workflows/run_tests.yml

@@ -5,16 +5,15 @@ name: Run tests
 
 on:
   push:
-    branches: [ "main" ]
+    branches: [ "main", '*/*' ]
   pull_request:
-    branches: [ "main" ]
+    branches: [ "main", '*/*'  ]
 
 permissions:
   contents: read
 
 jobs:
   build:
-
     runs-on: ubuntu-latest
 
     steps:
@@ -30,4 +29,7 @@ jobs:
         if [ -f requirements.txt ]; then pip install -r requirements.txt; fi
     - name: Test with pytest
       run: |
-        pytest -v
+        pytest -v --cov functions/ --cov-report=lcov
+    # Step 4: Submit to coveralls
+    - name: Submit to coveralls
+      uses: coverallsapp/github-action@v2.3.6

.gitignore

@@ -5,4 +5,13 @@
 .pytest_cache/
 .ipynb_checkpoints/
 .vscode/
-build/
+build/
+data/
+trash/
+.DS_Store
+*.nii.gz
+*.nii
+run.py
+tests.py
+spherical_analytical_vs_simulated.png
+cylindrical_analytical_vs_simulated.png

.pytest.ini

@@ -0,0 +1,3 @@
+[pytest]
+markers =
+    unit: small functions

README.md

@@ -1,4 +1,7 @@
 # susceptibility-to-fieldmap-fft
+[![Run tests](https://github.com/shimming-toolbox/susceptibility-to-fieldmap-fft/actions/workflows/run_tests.yml/badge.svg?branch=mb%2Fcustom_pad)](https://github.com/shimming-toolbox/susceptibility-to-fieldmap-fft/actions/workflows/run_tests.yml)
+[![Coverage Status](https://coveralls.io/repos/github/shimming-toolbox/susceptibility-to-fieldmap-fft/badge.svg?branch=mb/custom_pad)](https://coveralls.io/github/shimming-toolbox/susceptibility-to-fieldmap-fft?branch=mb/custom_pad)
+
 
 # Table of contents
 1. [Theory](#theory)
@@ -113,6 +116,8 @@ The `compute_fieldmap` command allows computation of a $B_0$ fieldmap based on a
 **Inputs** 
 - input_file : path to the susceptibility distribution (NIfTI file)
 - output_file : path for the fieldmap (NIfTI file)
+- -b, buffer (optional, default=50 voxels): Value-padding the edges of the volume
+- -m, method (optional, default=edge): Method used for value-padding, see np.pad options
 
 **Output** 
 The calculated fieldmap at the specified path.

functions/analytical_cases.py

@@ -3,6 +3,7 @@
 from matplotlib import pyplot as plt
 from scipy.ndimage import rotate
 import click
+import copy
 
 class Visualization:
     """
@@ -90,24 +91,26 @@ def plot_comparaison_analytical(self, Bz_analytical, simulated_Bz, geometry_type
         axes[0].plot(np.linspace(-dimensions[0]//2, dimensions[0]//2, dimensions[0]), simulated_Bz[:, dimensions[0]//2, dimensions[0]//2],'--', label='Simulated')
         axes[0].set_xlabel('x position [mm]')
         axes[0].set_ylabel('Field variation [ppm]')
+
         axes[0].set_ylim(vmin, vmax)
         axes[0].legend()
 
-        axes[1].plot(np.linspace(-dimensions[0]//2, dimensions[0]//2, dimensions[0]), Bz_analytical[dimensions[0]//2, :, dimensions[0]//2], label='Theory')
-        axes[1].plot(np.linspace(-dimensions[0]//2, dimensions[0]//2, dimensions[0]), simulated_Bz[dimensions[0]//2, :, dimensions[0]//2],'--', label='Simulated')
+        axes[1].plot(np.linspace(-dimensions[1]//2, dimensions[1]//2, dimensions[1]), Bz_analytical[dimensions[0]//2, :, dimensions[2]//2], label='Theory')
+        axes[1].plot(np.linspace(-dimensions[1]//2, dimensions[1]//2, dimensions[1]), simulated_Bz[dimensions[0]//2, :, dimensions[2]//2],'--', label='Simulated')
         axes[1].set_xlabel('y position [mm]')
         axes[1].set_ylabel('Field variation [ppm]')
         axes[1].set_ylim(vmin, vmax)
         axes[1].legend()
 
-        axes[2].plot(np.linspace(-dimensions[0]//2, dimensions[0]//2, dimensions[0]), Bz_analytical[dimensions[0]//2, dimensions[0]//2, :], label='Theory')
-        axes[2].plot(np.linspace(-dimensions[0]//2, dimensions[0]//2, dimensions[0]), simulated_Bz[dimensions[0]//2, dimensions[0]//2, :],'--', label='Simulated')
+        axes[2].plot(np.linspace(-dimensions[2]//2, dimensions[2]//2, dimensions[2]), Bz_analytical[dimensions[0]//2, dimensions[1]//2, :], label='Theory')
+        axes[2].plot(np.linspace(-dimensions[2]//2, dimensions[2]//2, dimensions[2]), simulated_Bz[dimensions[0]//2, dimensions[1]//2, :],'--', label='Simulated')
         axes[2].set_xlabel('z position [mm]')
         axes[2].set_ylabel('Field variation [ppm]')
         axes[2].set_ylim(vmin, vmax)
         axes[2].legend()
 
         plt.tight_layout()
+        plt.savefig(f'{geometry_type}_analytical_vs_simulated.png', dpi=300)
         plt.show()
 
 class Spherical(Visualization):
@@ -142,7 +145,7 @@ def mask(self):
         """
         [x, y, z] = np.meshgrid(np.linspace(-(self.matrix[0]-1)/2, (self.matrix[0]-1)/2, self.matrix[0]),
                                 np.linspace(-(self.matrix[1]-1)/2, (self.matrix[1]-1)/2, self.matrix[1]),
-                                np.linspace(-(self.matrix[2]-1)/2, (self.matrix[2]-1)/2, self.matrix[2]))
+                                np.linspace(-(self.matrix[2]-1)/2, (self.matrix[2]-1)/2, self.matrix[2]), indexing='ij')
 
         r = np.sqrt(x**2 + y**2 + z**2)
 
@@ -168,7 +171,7 @@ def analytical_sol(self):
 
         [x, y, z] = np.meshgrid(np.linspace(-(self.matrix[0]-1)/2, (self.matrix[0]-1)/2, self.matrix[0]),
                                 np.linspace(-(self.matrix[1]-1)/2, (self.matrix[1]-1)/2, self.matrix[1]),
-                                np.linspace(-(self.matrix[2]-1)/2, (self.matrix[2]-1)/2, self.matrix[2]))
+                                np.linspace(-(self.matrix[2]-1)/2, (self.matrix[2]-1)/2, self.matrix[2]), indexing='ij')
 
 
         r = np.sqrt(x**2 + y**2 + z**2)
@@ -213,9 +216,9 @@ def mask(self):
         """
         [x, y, z] = np.meshgrid(np.linspace(-(self.matrix[0]-1)/2, (self.matrix[0]-1)/2, self.matrix[0]),
                                 np.linspace(-(self.matrix[1]-1)/2, (self.matrix[1]-1)/2, self.matrix[1]),
-                                np.linspace(-(self.matrix[2]-1)/2, (self.matrix[2]-1)/2, self.matrix[2]))
+                                np.linspace(-(self.matrix[2]-1)/2, (self.matrix[2]-1)/2, self.matrix[2]), indexing='ij')
 
-        r = x**2 + y**2 
+        r = x**2 + y**2
 
         mask = r <= self.R**2
 
@@ -247,7 +250,7 @@ def analytical_sol(self):
 
         [x, y, z] = np.meshgrid(np.linspace(-(self.matrix[0]-1)/2, (self.matrix[0]-1)/2, self.matrix[0]),
                                 np.linspace(-(self.matrix[1]-1)/2, (self.matrix[1]-1)/2, self.matrix[1]),
-                                np.linspace(-(self.matrix[2]-1)/2, (self.matrix[2]-1)/2, self.matrix[2]))
+                                np.linspace(-(self.matrix[2]-1)/2, (self.matrix[2]-1)/2, self.matrix[2]), indexing='ij')
 
         r = np.sqrt(x**2 + y**2 + z**2)
 
@@ -273,19 +276,52 @@ def analytical_sol(self):
 
         return Bz_analytical
 
+def analytical_sphere_external_field(dx, dy, dz, chi_internal, chi_external, radius):
+    """
+    Calculate analytical external field for a sphere at relative position (dx, dy, dz).
+
+    Formula: Bz/B0 = 1/3 * (chi_i - chi_e) * a^3/r^3 * (3*cos^2(theta) - 1) + 1/3 * chi_e
+
+    Args:
+        dx, dy, dz: Position relative to sphere center (in voxels)
+        chi_internal: Susceptibility inside sphere (ppm)
+        chi_external: Susceptibility outside sphere (ppm)
+        radius: Sphere radius (in voxels)
+
+    Returns:
+        Field value at position (dx, dy, dz) in ppm
+    """
+    r = np.sqrt(dx**2 + dy**2 + dz**2)
+
+    if r == 0:
+        # At center
+        return chi_external / 3.0
+
+    # cos(theta) where theta is angle from z-axis
+    cos_theta = dz / r
+
+    # External field formula
+    field = (1.0/3.0) * (chi_internal - chi_external) * (radius**3 / r**3) * (3 * cos_theta**2 - 1) + chi_external / 3.0
+
+    return field
+
 @click.command(help="Compare the analytical solution to the simulated solution for a spherical or cylindrical geometry.")
 @click.option('-t', '--geometry-type',required=True, 
               type=click.Choice(['spherical', 'cylindrical']), 
               help='Type of geometry for the simulation')
-@click.option('-b', '--buffer', default=2, 
+@click.option('-b', '--buffer', default=50, 
               help='Buffer value for zero-padding.')
-def compare_to_analytical(geometry_type, buffer):
+def compare_to_analytical(geometry_type, buffer, matrix=[128,128,128], image_res=[1,1,1], radius=15, chi=9):
     """
     Main function to compare simulated fields to analytical solutions.
 
     Parameters:
     - geometry_type (str): The type of geometry to simulate ('spherical' or 'cylindrical').
     - buffer (float): The buffer size for the simulation.
+    - matrix ([int, int, int]): Volume dimensions in terms of number of pixels.
+    - image_res ([float, float, float]): Image resolution in mm ([x,y,z]).
+    - radius (float): Radius (mm) of the geometrical object.
+    - chi (float): Susceptibility difference in ppm.
 
     Returns:
     - None
@@ -301,10 +337,35 @@ def compare_to_analytical(geometry_type, buffer):
             and a susceptibility difference of 9 ppm for the spherical and cylindrical geometries.
     """
 
-    matrix = np.array([128,128,128])
-    image_res = np.array([1,1,1]) # mm
-    R = 15 # mm
-    sus_diff = 9 # ppm
+    compare_to_analytical_internal(geometry_type, buffer, matrix, image_res, radius, chi)
+
+
+def compare_to_analytical_internal(geometry_type, buffer, matrix=[128,128,128], image_res=[1,1,1], radius=15, chi=9):
+
+    # Type check the matrix argument
+    if not isinstance(matrix, (list, tuple, np.ndarray)):
+        raise TypeError("matrix must be a list, tuple, or numpy array.")
+    if len(matrix) != 3:
+        raise ValueError("matrix must have 3 elements (x, y, z dimensions).")
+    if not all(isinstance(dim, int) for dim in matrix):
+        raise TypeError("All elements of matrix must be integers.")
+    if not all(dim > 0 for dim in matrix):
+        raise ValueError("All matrix dimensions must be positive.")
+    matrix = np.array(matrix) # Convert to numpy array for easier use later
+
+    # Type check the image_res argument
+    if not isinstance(image_res, (list, tuple, np.ndarray)):
+        raise TypeError("image_res must be a list, tuple, or numpy array.")
+    if len(image_res) != 3:
+        raise ValueError("image_res must have 3 elements (x, y, z resolutions).")
+    if not all(isinstance(res, (int, float)) for res in image_res):
+        raise TypeError("All elements of image_res must be numbers (int or float).")
+    if not all(res > 0 for res in image_res):
+        raise ValueError("All image_res values must be positive.")
+    image_res = np.array(image_res)  # Convert to numpy array
+
+    R = radius # mm
+    sus_diff = chi # ppm
 
     dicto = {'spherical': Spherical(matrix, image_res, R, sus_diff),
               'cylindrical': Cylindrical(matrix, image_res, R, sus_diff)}
@@ -315,9 +376,12 @@ def compare_to_analytical(geometry_type, buffer):
 
     # compute Bz variation
     calculated_Bz = compute_bz(sus_dist, image_res, buffer)
+
     # analytical solution
     Bz_analytical = geometry.analytical_sol()
 
     # plot the results
     geometry.plot_susceptibility_and_fieldmap(sus_dist, calculated_Bz, geometry_type)
     geometry.plot_comparaison_analytical(Bz_analytical, calculated_Bz, geometry_type)
+
+    return calculated_Bz, Bz_analytical

functions/compute_fieldmap.py

@@ -2,6 +2,7 @@
 import nibabel as nib
 import click
 from time import perf_counter
+from pathlib import Path
 
 def is_nifti(filepath):
     """
@@ -39,7 +40,7 @@ def load_sus_dist(filepath):
     return susceptibility_distribution, image_resolution, affine_matrix
 
 
-def compute_bz(susceptibility_distribution, image_resolution=np.array([1,1,1]), buffer=1):
+def compute_bz(susceptibility_distribution, image_resolution=np.array([1,1,1]), buffer=50, mode='edge'):
     """
     Compute the Bz field variation in ppm based on a susceptibility distribution
     using a Fourier-based method.
@@ -49,39 +50,67 @@ def compute_bz(susceptibility_distribution, image_resolution=np.array([1,1,1]),
 
         image_resolution (numpy.ndarray, optional): The resolution of the image in each dimension. Defaults to [1, 1, 1].
 
-        buffer (int, optional): The buffer size for the k-space grid. Defaults to 1.
+        buffer (int, optional): The buffer size (voxels) for the k-space grid on each size. Defaults to 50.
+
+        mode (str, optional): The padding mode for the susceptibility distribution. Defaults to 'edge'.
 
     Returns:
         volume_without_buffer (numpy.ndarray): The computed magnetic field Bz in ppm.
 
     """
 
+    # Pad the susceptibility distribution
+
+    if mode == 'b0SimISMRM':
+            susceptibility_distribution=np.pad(susceptibility_distribution, ((buffer, buffer), (0,buffer), (buffer,buffer)), mode='edge')
+            susceptibility_distribution=np.pad(susceptibility_distribution, ((0,0),(buffer,0),(0,0)), mode='constant', constant_values=0.35)
+    else:
+        susceptibility_distribution=np.pad(susceptibility_distribution, buffer, mode=mode)
+
     # dimensions needs to be a numpy.array
     dimensions = np.array(susceptibility_distribution.shape)
 
-    # creating the k-space grid with the buffer
-    new_dimensions = buffer*np.array(dimensions)
     kmax = 1/(2*image_resolution)
 
-    [kx, ky, kz] = np.meshgrid(np.linspace(-kmax[0], kmax[0], new_dimensions[0]),
-                                np.linspace(-kmax[1], kmax[1], new_dimensions[1]),
-                                np.linspace(-kmax[2], kmax[2], new_dimensions[2]), indexing='ij')
+    interval = 2 * kmax / dimensions
+
+
+    kx_min_shift = (dimensions[0]%2)*interval[0]/2
+    ky_min_shift = (dimensions[1]%2)*interval[1]/2
+    kz_min_shift = (dimensions[2]%2)*interval[2]/2
+
+    kx_max_shift = -interval[0] + (dimensions[0]%2)*interval[0]/2
+    ky_max_shift = -interval[1] + (dimensions[1]%2)*interval[1]/2
+    kz_max_shift = -interval[2] + (dimensions[2]%2)*interval[2]/2 
+
+
+    [kx, ky, kz] = np.meshgrid(np.linspace(-kmax[0] + kx_min_shift, kmax[0] + kx_max_shift, dimensions[0]),
+                                np.linspace(-kmax[1] + ky_min_shift, kmax[1] + ky_max_shift, dimensions[1]),
+                                np.linspace(-kmax[2] + kz_min_shift, kmax[2] + kz_max_shift, dimensions[2]), indexing='ij')
 
     # FFT procedure
     # undetermined at the center of k-space
     k2 = kx**2 + ky**2 + kz**2
 
     with np.errstate(divide='ignore', invalid='ignore'):
-        kernel = np.fft.fftshift(1/3 - kz**2/k2)
-        kernel[0,0,0] = 1/3
+        x_kernel = 1/3 - kz**2/k2
+
+        x_kernel[int(dimensions[0]/2-1/2*(dimensions[0]%2)), int(dimensions[1]/2-1/2*(dimensions[1]%2)), int(dimensions[2]/2-1/2*(dimensions[2]%2))] = 1/3
+        
+        kernel = np.fft.ifftshift(x_kernel)
+
+    FFT_chi = np.fft.fftn(susceptibility_distribution, dimensions)
 
-    FFT_chi = np.fft.fftn(susceptibility_distribution, new_dimensions)
-    FFT_chi[0,0,0] = FFT_chi[0,0,0] + np.prod(new_dimensions)*susceptibility_distribution[0,0,0]
     Bz_fft = kernel*FFT_chi
 
     # retrive the inital FOV
+
     volume_with_buffer = np.real(np.fft.ifftn(Bz_fft))
-    volume_without_buffer = volume_with_buffer[0:dimensions[0], 0:dimensions[1], 0:dimensions[2]]
+
+    if buffer == 0:
+        volume_without_buffer = volume_with_buffer
+    else:
+        volume_without_buffer = volume_with_buffer[buffer:-buffer, buffer:-buffer, buffer:-buffer]
 
     return volume_without_buffer
 
@@ -97,6 +126,7 @@ def save_to_nifti(data, affine_matrix, output_path):
     Returns:
         None
     """
+    data=data.astype(np.float32)
     nifti_image = nib.Nifti1Image(data, affine_matrix)
     nib.save(nifti_image, output_path)
 
@@ -107,7 +137,11 @@ def save_to_nifti(data, affine_matrix, output_path):
               help="Input susceptibility distribution, supported extensions: .nii, .nii.gz")
 @click.option('-o', '--output', 'output_file', type=click.Path(), default='fieldmap.nii.gz',
               help="Output fieldmap, supported extensions: .nii, .nii.gz")
-def compute_fieldmap(input_file, output_file):
+@click.option('-b', '--buffer', 'buffer', type=int, default=50, required=False,
+                help="Buffer size (voxels) for the k-space grid on each side")
+@click.option('-m', '--mode', 'mode', type=str, default='edge', required=False,
+                help="Padding mode for the susceptibility distribution")
+def compute_fieldmap(input_file, output_file, buffer, mode):
     """
     Main procedure for performing the simulation.
 
@@ -123,13 +157,18 @@ def compute_fieldmap(input_file, output_file):
         print('Start')
         susceptibility_distribution, image_resolution, affine_matrix = load_sus_dist(input_file)
         print('Susceptibility distribution loaded')
-        fieldmap = compute_bz(susceptibility_distribution, image_resolution)
+        fieldmap = compute_bz(susceptibility_distribution, image_resolution, buffer, mode)
         print('Fieldmap simulated')
+        
+        # Check if all subdirectories exist and create them if not
+        output_path = Path(output_file)
+        output_path.parent.mkdir(parents=True, exist_ok=True)
+        
         save_to_nifti(fieldmap, affine_matrix, output_file)
         print('Saving to NIfTI format')
         end_time = perf_counter()
         print(f'End. Runtime: {end_time-start_time:.2f} seconds')
     else:
         print("The input file must be NIfTI.")
 
-    
+    

requirements.txt

@@ -0,0 +1 @@
+