From f9ccc87940cb349aeff1e12ed4680364ca3715f8 Mon Sep 17 00:00:00 2001 From: charles pageot Date: Fri, 5 Jul 2024 10:47:30 -0400 Subject: [PATCH 01/21] Update setup.py file to include entry points --- setup.py | 28 ++++++++++++++++++++++++++-- 1 file changed, 26 insertions(+), 2 deletions(-) diff --git a/setup.py b/setup.py index 8f3e113..bb5f6d8 100644 --- a/setup.py +++ b/setup.py @@ -1,6 +1,30 @@ from setuptools import find_packages, setup +from os import path + +# Get the directory where this current file is saved +here = path.abspath(path.dirname(__file__)) + +with open(path.join(here, "README.rst"), encoding="utf-8") as f: + long_description = f.read() setup( - name="fft_simulation", - packages=find_packages(), + name="susceptibility-to-fieldmap-fft", + python_requires=">=3.8.8", + description="Code to compute the magnetic field variation from a given susceptibility distribution placed in a constant magnetic field.", + long_description=long_description, + url="https://github.com/shimming-toolbox/susceptibility-to-fieldmap-fft", + entry_points={ + 'console_scripts': [ + "compute_fieldmap=fft_simulation.fft_simulation:main", + "analytical_cases=fft_simulation.analytical_cases:main" + ] + }, + packages=find_packages(exclude=["docs"]), + install_requires=[ + "click", + "numpy>=1.24.4", + "nibabel>=5.2.1", + "matplotlib>=3.7.5", + "scipy>=1.10.1" + ], ) From 5b4c179793461363125253b7a4c5a797bfae56e2 Mon Sep 17 00:00:00 2001 From: charles pageot Date: Fri, 5 Jul 2024 10:48:26 -0400 Subject: [PATCH 02/21] Added documentation --- fft_simulation/analytical_cases.py | 117 +++++++++++++++++++++++++---- fft_simulation/fft_simulation.py | 30 ++------ 2 files changed, 108 insertions(+), 39 deletions(-) diff --git a/fft_simulation/analytical_cases.py b/fft_simulation/analytical_cases.py index a2a1bb0..af8aa57 100644 --- a/fft_simulation/analytical_cases.py +++ b/fft_simulation/analytical_cases.py @@ -1,11 +1,26 @@ -from fft_simulation import compute_bz +from fft_simulation.fft_simulation import compute_bz import numpy as np from matplotlib import pyplot as plt from scipy.ndimage import rotate -import argparse import click class Spherical: + """ + Represents a spherical object in 3D space. + + Attributes: + - matrix (ndarray): The dimensions of the matrix representing the object (i.e. [128, 128, 128]). + - image_res (ndarray): The resolution of the image in mm (i.e. [1, 1, 1] mm). + - R (int): The radius of the sphere in mm. + - sus_diff (float): The susceptibility difference of the sphere + (sus_diff = susceptibility_in_the_sphere - susceptibility_outside). + + Methods: + - mask(): Generates a mask representing the spherical object. + - volume(): Create a volume of the sphere with the corresponding suceptibility value. + - analyticial_sol(): Calculates the analytical solution for the magnetic field inside and outside the sphere. + """ + def __init__(self, matrix, image_res, R, sus_diff): self.matrix = matrix self.image_res = image_res @@ -13,6 +28,12 @@ def __init__(self, matrix, image_res, R, sus_diff): self.sus_diff = sus_diff def mask(self): + """ + Generates a mask representing the spherical object. + + Returns: + - mask (ndarray): A boolean array representing the mask. + """ [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])) @@ -22,9 +43,21 @@ def mask(self): return r**2 < self.R**2 def volume(self): + """ + Create a volume of the sphere with the corresponding suceptibility value. + + Returns: + - volume (ndarray): A 3D array representing the distribution of suceptibility. + """ return np.where(self.mask() == True, self.sus_diff, 0) def analyticial_sol(self): + """ + Calculates the analytical solution for the magnetic field inside and outside the sphere. + + Returns: + - Bz_analytical (ndarray): A 3D array representing the analytical solution for the magnetic field. + """ mask = self.mask() [x, y, z] = np.meshgrid(np.linspace(-(self.matrix[0]-1)/2, (self.matrix[0]-1)/2, self.matrix[0]), @@ -47,6 +80,23 @@ def analyticial_sol(self): class Cylindrical: + """ + Represents a cylindrical object in 3D space. + + Parameters: + - matrix (ndarray): The dimensions of the matrix representing the object (i.e. [128, 128, 128]). + - image_res (ndarray): The resolution of the image in mm (i.e. [1, 1, 1] mm). + - R (int): The radius of the cylinder in mm. + - sus_diff (float): The susceptibility difference of the cylinder + (sus_diff = susceptibility_in_the_cylinder - susceptibility_outside). + - theta (float, optional): The rotation angle about the y-axis. Default is pi/2. + + Methods: + - mask(): Generates a mask representing the cylinder. + - volume(): Create a volume of the sphere with the corresponding suceptibility value. + - analytical_sol(): Calculates the analytical solution for the magnetic field inside and outside the cylinder. + """ + def __init__(self, matrix, image_res, R, sus_diff, theta=np.pi/2): self.matrix = matrix self.image_res = image_res @@ -55,6 +105,12 @@ def __init__(self, matrix, image_res, R, sus_diff, theta=np.pi/2): self.theta = theta # rotation angle about the y-axis def mask(self): + """ + Generates a mask representing the cylinder. + + Returns: + - mask (ndarray): The mask representing the cylinder. + """ [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])) @@ -67,10 +123,22 @@ def mask(self): return rotate(mask, self.theta*180/np.pi, axes=(0, 2), reshape=False, order=1) def volume(self): + """ + Create a volume of the sphere with the corresponding suceptibility value. + + Returns: + - volume (ndarray): A 3D array representing the distribution of suceptibility. + """ return np.where(self.mask() == True, self.sus_diff, 0) -# TODO: For now the analytical solution is only correct for theta = 90 def analytical_sol(self): + """ + Calculates the analytical solution for the magnetic field inside and outside the cylinder. + + Returns: + - Bz_analytical_x (ndarray): The analytical solution for the magnetic field along the x-axis. + - Bz_analytical_y (ndarray): The analytical solution for the magnetic field along the y-axis. + """ phi_x = 0 phi_y = np.pi/2 @@ -96,8 +164,34 @@ def analytical_sol(self): return Bz_analytical_x, Bz_analytical_y - +@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, + help='Buffer value for zero-padding.') def main(geometry_type, buffer): + """ + Main function for running the Fourier-based simulation with analytical cases. + + Parameters: + - geometry_type (str): The type of geometry to simulate ('spherical' or 'cylindrical'). + - buffer (float): The buffer size for the simulation. + + Returns: + - None + + This function performs the following steps: + 1. Initializes the necessary variables and parameters. + 2. Creates the susceptibility geometry based on the specified geometry type. + 3. Computes the Bz variation using the computed susceptibility distribution and buffer size. + 4. Plots sections of the susceptibility distribution and Bz field variation. + 5. Plots the analytical solution and simulated results for the Bz field variation. + + Note: The function uses a matrix size of [128, 128, 128], an image resolution of [1, 1, 1] mm, a radius of 15 mm + 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 @@ -146,6 +240,7 @@ def main(geometry_type, buffer): plt.tight_layout() plt.show() + # comparaison with the analytical solution fig, axes = plt.subplots(1, 3, figsize=(10, 3), dpi=120) @@ -218,6 +313,8 @@ def main(geometry_type, buffer): plt.tight_layout() plt.show() + # comparaison with the analytical solution + fig, axes = plt.subplots(1, 3, figsize=(10, 3), dpi=120) axes[0].plot(np.linspace(-64,64, 128), Bz_analytical_x[:,63,63], label='Theory') @@ -243,15 +340,3 @@ def main(geometry_type, buffer): plt.tight_layout() plt.show() - -@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, - help='Buffer value for zero-padding.') -def cli(geometry_type, buffer): - main(geometry_type, buffer) - -if __name__ == "__main__": - cli() \ No newline at end of file diff --git a/fft_simulation/fft_simulation.py b/fft_simulation/fft_simulation.py index ab5617d..b2d2b3a 100644 --- a/fft_simulation/fft_simulation.py +++ b/fft_simulation/fft_simulation.py @@ -102,31 +102,19 @@ def save_to_nifti(data, image_resolution, output_path): -@click.command(help="Compute the magnetic field Bz in ppm from a susceptibility distribution in NIfTI format.") +@click.command(help="Compute the magnetic field variation in ppm from a susceptibility distribution in NIfTI format.") @click.argument('input_file', required=True, type=click.Path(exists=True)) @click.argument('output_file', required=True, type=click.Path()) def main(input_file, output_file): - # parser = argparse.ArgumentParser() - # parser.add_argument("-i", - # dest="input_file", - # type=str, - # required=True, - # help="Path to the NIfTI file input.") - - # parser.add_argument("-o", - # dest="output_file", - # type=str, - # required=True, - # help="Path to the NIfTI file ouput.") - # args = parser.parse_args() - """ - This script processes a NIfTI file by performing a simple operation on the voxel data - and saves the modified data to a new NIfTI file. + Main procedure for performing the simulation. - INPUT_FILE: Path to the input susceptibility distribution in NIfTI format. + Args: + input_file (str): Path to the susceptibility distribution in NIfTI format. + output_file (str): Path for the computed fieldmap in NIfTI format. - OUTPUT_FILE: Path to the output fieldmap where the data will be saved in NIfTI format. + Returns: + None """ if is_nifti(input_file): sus_dist, img_res = load_sus_dist(input_file) @@ -135,8 +123,4 @@ def main(input_file, output_file): else: print("The input file must be NIfTI.") - -if __name__ == "__main__": - main() - \ No newline at end of file From 1de698f9c5235114026407ac41ef5aa176e4f071 Mon Sep 17 00:00:00 2001 From: Charles Pageot <167803616+CharlesPageot@users.noreply.github.com> Date: Fri, 5 Jul 2024 11:00:45 -0400 Subject: [PATCH 03/21] Update README.md --- README.md | 24 ++++++++++++------------ 1 file changed, 12 insertions(+), 12 deletions(-) diff --git a/README.md b/README.md index e2b7305..2691aa2 100644 --- a/README.md +++ b/README.md @@ -89,20 +89,20 @@ Navigate to the project directory cd susceptibility-to-fieldmap-fft ``` -Install the requirements +Install the package ``` -pip install -r requirements.txt +pip install . ``` ## Usage -To execute the scripts, you need to naviagte into the fft_simulation folder. +Once the package is install, the command can be run directly from the terminal. Here is the description of the two commands available. -### Analytical cases +### analytical_cases -The _analytical_case.py_ script allows for comparaison between simulated and analytical results for a spherical and cylindrical phantom. +The _analytical_cases_ command allow for comparaison between simulated and analytical results for a spherical and cylindrical phantom. -**Arguments** +**Options** - -t, type : 'spherical' or 'cylindrical' - -b, buffer (optional, default=2): Buffer value for zero-padding around the phantom @@ -111,24 +111,24 @@ Plots to visialize the results Example: ``` -python analytical_cases.py -t "spherical" +analytical_cases -t "spherical" ``` -### fft_simulation +### compute_fieldmap -The _fft_simulation.py_ script allow computation of a $B_0$ fieldmap based on a susceptibility distribution given as an input. +The _compute_fieldmap_ command allow computation of a $B_0$ fieldmap based on a susceptibility distribution given as an input. **Arguments** -- -input_file : path to the susceptibility distribution (NIfTI file) -- -output_file : path for the fieldmap (NIfTI file) +- input_file : path to the susceptibility distribution (NIfTI file) +- output_file : path for the fieldmap (NIfTI file) **Return** The calculated fieldmap at the specified path. Example: ``` -python fft_simulation.py "\your\sus_dist\path.nii" "\yout\fieldmap\path.nii" +compute_fieldmap "\your\sus_dist\path.nii" "\yout\fieldmap\path.nii" ``` ## References : From 168be37a22707263e28ad752c02f42123193a65b Mon Sep 17 00:00:00 2001 From: Charles Pageot <167803616+CharlesPageot@users.noreply.github.com> Date: Fri, 5 Jul 2024 11:06:34 -0400 Subject: [PATCH 04/21] Update README.md --- README.md | 6 +++--- 1 file changed, 3 insertions(+), 3 deletions(-) diff --git a/README.md b/README.md index 9468303..084ed13 100644 --- a/README.md +++ b/README.md @@ -100,11 +100,11 @@ pip install . ``` ## Usage -Once the package is install, the command can be run directly from the terminal. Here is the description of the two commands available. +Once the package is install, the commands can be run directly from the terminal. Here is the description of the two commands available. ### analytical_cases -The _analytical_cases_ command allow for comparaison between simulated and analytical results for a spherical and cylindrical phantom. +The _analytical_cases_ command allows for comparaison between simulated and analytical results for a spherical and cylindrical phantom. **Options** - -t, type : 'spherical' or 'cylindrical' @@ -121,7 +121,7 @@ analytical_cases -t "spherical" ### compute_fieldmap -The _compute_fieldmap_ command allow computation of a $B_0$ fieldmap based on a susceptibility distribution given as an input. +The _compute_fieldmap_ command allows computation of a $B_0$ fieldmap based on a susceptibility distribution given as an input. **Arguments** - input_file : path to the susceptibility distribution (NIfTI file) From c1ef3174acfc8af2d3703c15f902fc8daeb3acf2 Mon Sep 17 00:00:00 2001 From: charles pageot Date: Fri, 5 Jul 2024 11:23:17 -0400 Subject: [PATCH 05/21] Updated setup.py --- setup.py | 2 +- 1 file changed, 1 insertion(+), 1 deletion(-) diff --git a/setup.py b/setup.py index bb5f6d8..7052236 100644 --- a/setup.py +++ b/setup.py @@ -4,7 +4,7 @@ # Get the directory where this current file is saved here = path.abspath(path.dirname(__file__)) -with open(path.join(here, "README.rst"), encoding="utf-8") as f: +with open(path.join(here, "README.md"), encoding="utf-8") as f: long_description = f.read() setup( From eebb7f5ea16b2abf16100b396905f94e9cff9d59 Mon Sep 17 00:00:00 2001 From: Charles Pageot <167803616+CharlesPageot@users.noreply.github.com> Date: Fri, 5 Jul 2024 11:28:10 -0400 Subject: [PATCH 06/21] Update README.md --- README.md | 6 ++++++ 1 file changed, 6 insertions(+) diff --git a/README.md b/README.md index 084ed13..214a6a5 100644 --- a/README.md +++ b/README.md @@ -93,6 +93,12 @@ Navigate to the project directory cd susceptibility-to-fieldmap-fft ``` +(temporary) Go to the develop branch + +``` +git checkout develop +``` + Install the package ``` From 3472b1ad0ab2abd0e2fbdff75e6d17d6daa1b39a Mon Sep 17 00:00:00 2001 From: Charles Pageot <167803616+CharlesPageot@users.noreply.github.com> Date: Mon, 8 Jul 2024 10:51:40 -0400 Subject: [PATCH 07/21] Update README.md --- README.md | 2 +- 1 file changed, 1 insertion(+), 1 deletion(-) diff --git a/README.md b/README.md index 214a6a5..6197bd3 100644 --- a/README.md +++ b/README.md @@ -138,7 +138,7 @@ The calculated fieldmap at the specified path. Example: ``` -compute_fieldmap "\your\sus_dist\path.nii" "\yout\fieldmap\path.nii" +compute_fieldmap "inpath/susceptibility_distribution.nii.gz" "outpath/fieldmap.nii.gz" ``` ## References From 611222fc7a54bc1a6ae336e806b623147bdd3316 Mon Sep 17 00:00:00 2001 From: Charles Pageot <167803616+CharlesPageot@users.noreply.github.com> Date: Tue, 9 Jul 2024 09:43:58 -0400 Subject: [PATCH 08/21] Update README.md --- README.md | 2 +- 1 file changed, 1 insertion(+), 1 deletion(-) diff --git a/README.md b/README.md index 6197bd3..c76cadb 100644 --- a/README.md +++ b/README.md @@ -106,7 +106,7 @@ pip install . ``` ## Usage -Once the package is install, the commands can be run directly from the terminal. Here is the description of the two commands available. +Once the package is installed, the commands can be run directly from the terminal. Here is the description of the two commands available. ### analytical_cases From df7f37949f93c66d626c833b93e28f47e90a6be2 Mon Sep 17 00:00:00 2001 From: Charles Pageot <167803616+CharlesPageot@users.noreply.github.com> Date: Tue, 9 Jul 2024 09:56:45 -0400 Subject: [PATCH 09/21] Update README.md --- README.md | 2 +- 1 file changed, 1 insertion(+), 1 deletion(-) diff --git a/README.md b/README.md index c76cadb..42cce9f 100644 --- a/README.md +++ b/README.md @@ -113,7 +113,7 @@ Once the package is installed, the commands can be run directly from the termina The _analytical_cases_ command allows for comparaison between simulated and analytical results for a spherical and cylindrical phantom. **Options** -- -t, type : 'spherical' or 'cylindrical' +- -t, geometry type : 'spherical' or 'cylindrical' - -b, buffer (optional, default=2): Buffer value for zero-padding around the phantom **Return** From 044c161be7cdc27a1af9f86eb524973f64222c70 Mon Sep 17 00:00:00 2001 From: charles pageot Date: Tue, 9 Jul 2024 11:11:42 -0400 Subject: [PATCH 10/21] Added titles to figures and change to use [] operator for masking --- fft_simulation/analytical_cases.py | 51 ++++++++++++++++++------------ 1 file changed, 31 insertions(+), 20 deletions(-) diff --git a/fft_simulation/analytical_cases.py b/fft_simulation/analytical_cases.py index af8aa57..2b88c8c 100644 --- a/fft_simulation/analytical_cases.py +++ b/fft_simulation/analytical_cases.py @@ -68,11 +68,9 @@ def analyticial_sol(self): r = np.sqrt(x**2 + y**2 + z**2) Bz_out = self.sus_diff/3 * (self.R/r)**3 * (3*z**2/r**2 - 1) - Bz_in = np.zeros(self.matrix) - + Bz_out[mask] = 0 - Bz_in = np.where(mask, Bz_in, 0) - Bz_out = np.where(~mask, Bz_out, 0) + Bz_in = np.zeros(self.matrix) Bz_analytical = Bz_out + Bz_in @@ -129,7 +127,7 @@ def volume(self): Returns: - volume (ndarray): A 3D array representing the distribution of suceptibility. """ - return np.where(self.mask() == True, self.sus_diff, 0) + return np.where(self.mask(), self.sus_diff, 0) def analytical_sol(self): """ @@ -139,6 +137,9 @@ def analytical_sol(self): - Bz_analytical_x (ndarray): The analytical solution for the magnetic field along the x-axis. - Bz_analytical_y (ndarray): The analytical solution for the magnetic field along the y-axis. """ + + mask = self.mask() + phi_x = 0 phi_y = np.pi/2 @@ -150,14 +151,14 @@ def analytical_sol(self): # solution along the x-axis: phi = 0 Bz_out_x = self.sus_diff/2 * (self.R/r)**2 * np.sin(self.theta)**2 * np.cos(2*phi_x) - Bz_out_x = np.where(~self.mask(), Bz_out_x, 0) + Bz_out_x[mask] = 0 # solution along the x-axis: phi = 90 Bz_out_y = self.sus_diff/2 * (self.R/r)**2 * np.sin(self.theta)**2 * np.cos(2*phi_y) - Bz_out_y = np.where(~self.mask(), Bz_out_y, 0) + Bz_out_y[mask] = 0 Bz_in = np.zeros(self.matrix) + self.sus_diff/6 * (3*np.cos(self.theta) - 1) - Bz_in = np.where(self.mask(), Bz_in, 0) + Bz_in[~mask] = 0 Bz_analytical_x = Bz_out_x + Bz_in Bz_analytical_y = Bz_out_y + Bz_in @@ -208,17 +209,21 @@ def main(geometry_type, buffer): # plot sections of the chi distribution fig, axes = plt.subplots(2, 3, figsize=(10, 5), dpi=120) + fig.suptitle('Susceptibility distribution (top) and Bz field variation (bottom) for a spherical geometry') h = axes[0,0].imshow(sus_dist[matrix[0] // 2, :, :], origin='lower') - axes[0,0].set_title('Slice along y-z plane') + axes[0,0].set_title('Y-Z plane') + axes[0,0].axis("off") plt.colorbar(h, label='Susceptibility [ppm]') h = axes[0,1].imshow(sus_dist[:, matrix[0] // 2, :], origin='lower') - axes[0,1].set_title('Slice along x-z plane') + axes[0,1].set_title('Z-X plane') + axes[0,1].axis("off") plt.colorbar(h, label='Susceptibility [ppm]') h = axes[0,2].imshow(sus_dist[:, :, matrix[0] // 2], origin='lower') - axes[0,2].set_title('Slice along y-z plane') + axes[0,2].set_title('X-Y plane') + axes[0,2].axis("off") plt.colorbar(h, label='Susceptibility [ppm]') # plot section of the b0 field variation @@ -227,15 +232,15 @@ def main(geometry_type, buffer): vmax = np.max(calculated_Bz)*1.1 h = axes[1,0].imshow(calculated_Bz[63, :, :],vmin=vmin, vmax=vmax, origin='lower') - axes[0,0].set_title('Y-Z plane') + axes[1,0].axis("off") plt.colorbar(h, label='Field variation [ppm]') h = axes[1,1].imshow(calculated_Bz[:, 63, :],vmin=vmin, vmax=vmax, origin='lower') - axes[0,1].set_title('X-Z plane') + axes[1,1].axis("off") plt.colorbar(h, label='Field variation [ppm]') h = axes[1,2].imshow(calculated_Bz[:, :, 63],vmin=vmin, vmax=vmax, origin='lower') - axes[0,2].set_title('X-Y plane') + axes[1,2].axis("off") plt.colorbar(h, label='Field variation [ppm]') plt.tight_layout() plt.show() @@ -243,6 +248,7 @@ def main(geometry_type, buffer): # comparaison with the analytical solution fig, axes = plt.subplots(1, 3, figsize=(10, 3), dpi=120) + fig.suptitle('Analytical solution and simulated results for the Bz field variation for a spherical geometry') axes[0].plot(np.linspace(-64,64, 128), Bz_analytical[:,63,63], label='Theory') axes[0].plot(np.linspace(-64,64, 128), calculated_Bz[:,63,63],'--', label='Simulated') @@ -281,17 +287,21 @@ def main(geometry_type, buffer): # plot sections of the chi distribution fig, axes = plt.subplots(2, 3, figsize=(10, 5), dpi=120) + fig.suptitle('Susceptibility distribution (top) and Bz field variation (bottom) for a cylindrical geometry') h = axes[0,0].imshow(sus_dist[matrix[0] // 2, :, :], origin='lower') - axes[0,0].set_title('Slice along y-z plane') + axes[0,0].set_title('Y-Z plane') + axes[0,0].axis("off") plt.colorbar(h, label='Susceptibility [ppm]') h = axes[0,1].imshow(sus_dist[:, matrix[0] // 2, :], origin='lower') - axes[0,1].set_title('Slice along x-z plane') + axes[0,1].set_title('X-Z plane') + axes[0,1].axis("off") plt.colorbar(h, label='Susceptibility [ppm]') h = axes[0,2].imshow(sus_dist[:, :, matrix[0] // 2], origin='lower') - axes[0,2].set_title('Slice along y-z plane') + axes[0,2].set_title('X-Y plane') + axes[0,2].axis("off") plt.colorbar(h, label='Susceptibility [ppm]') # plot section of the b0 field variation @@ -300,15 +310,15 @@ def main(geometry_type, buffer): vmax = np.max(calculated_Bz)*1.1 h = axes[1,0].imshow(calculated_Bz[63, :, :],vmin=vmin, vmax=vmax, origin='lower') - axes[0,0].set_title('Y-Z plane') + axes[1,0].axis("off") plt.colorbar(h, label='Field variation [ppm]') h = axes[1,1].imshow(calculated_Bz[:, 63, :],vmin=vmin, vmax=vmax, origin='lower') - axes[0,1].set_title('X-Z plane') + axes[1,1].axis("off") plt.colorbar(h, label='Field variation [ppm]') h = axes[1,2].imshow(calculated_Bz[:, :, 63],vmin=vmin, vmax=vmax, origin='lower') - axes[0,2].set_title('X-Y plane') + axes[1,2].axis("off") plt.colorbar(h, label='Field variation [ppm]') plt.tight_layout() plt.show() @@ -316,6 +326,7 @@ def main(geometry_type, buffer): # comparaison with the analytical solution fig, axes = plt.subplots(1, 3, figsize=(10, 3), dpi=120) + fig.suptitle('Analytical solution and simulated results for the Bz field variation for a cylindrical geometry') axes[0].plot(np.linspace(-64,64, 128), Bz_analytical_x[:,63,63], label='Theory') axes[0].plot(np.linspace(-64,64, 128), calculated_Bz[:,63,63],'--', label='Simulated') From e2fc3df0bff19449f4629d8b7c71a498c63f9ab3 Mon Sep 17 00:00:00 2001 From: charles pageot Date: Tue, 9 Jul 2024 11:14:29 -0400 Subject: [PATCH 11/21] Small change to avoid hardcode --- fft_simulation/fft_simulation.py | 2 +- 1 file changed, 1 insertion(+), 1 deletion(-) diff --git a/fft_simulation/fft_simulation.py b/fft_simulation/fft_simulation.py index b2d2b3a..1b549e6 100644 --- a/fft_simulation/fft_simulation.py +++ b/fft_simulation/fft_simulation.py @@ -12,7 +12,7 @@ def is_nifti(filepath): Returns: bool: True if the file is a NIfTI file, False otherwise. """ - if filepath[-4:] == '.nii' or filepath[-7:] == '.nii.gz': + if filepath.endswith('.nii') or filepath.endswith('.nii.gz'): return True else: return False From c7a3d0b38dfc7c5617ece17052c7b8c3a9d4e713 Mon Sep 17 00:00:00 2001 From: Charles Pageot <167803616+CharlesPageot@users.noreply.github.com> Date: Wed, 10 Jul 2024 11:02:06 -0400 Subject: [PATCH 12/21] Update README.md --- README.md | 14 +++++++++----- 1 file changed, 9 insertions(+), 5 deletions(-) diff --git a/README.md b/README.md index 42cce9f..6119a0e 100644 --- a/README.md +++ b/README.md @@ -81,30 +81,34 @@ These final equations are the ones used in the function **compute_bz**, which ca ## Installation -First, clone the repository +- Clone the repository ``` git clone https://github.com/shimming-toolbox/susceptibility-to-fieldmap-fft.git +cd susceptibility-to-fieldmap-fft ``` -Navigate to the project directory +- Create a virtual environnement ``` -cd susceptibility-to-fieldmap-fft +conda create --name python=3.9 +conda activate ``` -(temporary) Go to the develop branch +- Go to the develop branch (temporary) ``` git checkout develop ``` -Install the package +- Install the package ``` pip install . ``` +You will need to ```conda activate ``` each time you want to use the package. + ## Usage Once the package is installed, the commands can be run directly from the terminal. Here is the description of the two commands available. From 96a552be556c6b3e4034a399b34b74eade093df8 Mon Sep 17 00:00:00 2001 From: charles pageot Date: Wed, 10 Jul 2024 16:56:31 -0400 Subject: [PATCH 13/21] Change from click.argument to click.option --- fft_simulation/fft_simulation.py | 6 ++++-- 1 file changed, 4 insertions(+), 2 deletions(-) diff --git a/fft_simulation/fft_simulation.py b/fft_simulation/fft_simulation.py index 1b549e6..dec44ee 100644 --- a/fft_simulation/fft_simulation.py +++ b/fft_simulation/fft_simulation.py @@ -103,8 +103,10 @@ def save_to_nifti(data, image_resolution, output_path): @click.command(help="Compute the magnetic field variation in ppm from a susceptibility distribution in NIfTI format.") -@click.argument('input_file', required=True, type=click.Path(exists=True)) -@click.argument('output_file', required=True, type=click.Path()) +@click.option('-i','--input','input_file', type=click.Path(exists=True), required=True, + 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 main(input_file, output_file): """ Main procedure for performing the simulation. From dc2604a70296d090765659fc1868826198a219cd Mon Sep 17 00:00:00 2001 From: charles pageot Date: Thu, 11 Jul 2024 10:53:26 -0400 Subject: [PATCH 14/21] Change the main() procedures to something more descriptive --- fft_simulation/analytical_cases.py | 4 ++-- fft_simulation/fft_simulation.py | 2 +- setup.py | 4 ++-- 3 files changed, 5 insertions(+), 5 deletions(-) diff --git a/fft_simulation/analytical_cases.py b/fft_simulation/analytical_cases.py index 2b88c8c..7aec733 100644 --- a/fft_simulation/analytical_cases.py +++ b/fft_simulation/analytical_cases.py @@ -171,9 +171,9 @@ def analytical_sol(self): help='Type of geometry for the simulation') @click.option('-b', '--buffer', default=2, help='Buffer value for zero-padding.') -def main(geometry_type, buffer): +def compare_to_analytical(geometry_type, buffer): """ - Main function for running the Fourier-based simulation with analytical cases. + Main function to compare simulated fields to analytical solutions. Parameters: - geometry_type (str): The type of geometry to simulate ('spherical' or 'cylindrical'). diff --git a/fft_simulation/fft_simulation.py b/fft_simulation/fft_simulation.py index dec44ee..a452cf9 100644 --- a/fft_simulation/fft_simulation.py +++ b/fft_simulation/fft_simulation.py @@ -107,7 +107,7 @@ def save_to_nifti(data, image_resolution, 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 main(input_file, output_file): +def compute_fieldmap(input_file, output_file): """ Main procedure for performing the simulation. diff --git a/setup.py b/setup.py index 7052236..5a14c9b 100644 --- a/setup.py +++ b/setup.py @@ -15,8 +15,8 @@ url="https://github.com/shimming-toolbox/susceptibility-to-fieldmap-fft", entry_points={ 'console_scripts': [ - "compute_fieldmap=fft_simulation.fft_simulation:main", - "analytical_cases=fft_simulation.analytical_cases:main" + "compute_fieldmap=fft_simulation.fft_simulation:compute_fieldmap", + "analytical_cases=fft_simulation.analytical_cases:compare_to_analytical" ] }, packages=find_packages(exclude=["docs"]), From fc17441d6e53ab5286ed17b0d979a4d3bd5b5ebd Mon Sep 17 00:00:00 2001 From: Charles Pageot <167803616+CharlesPageot@users.noreply.github.com> Date: Fri, 12 Jul 2024 11:57:05 -0400 Subject: [PATCH 15/21] Update README.md Co-authored-by: Julien Cohen-Adad --- README.md | 4 ---- 1 file changed, 4 deletions(-) diff --git a/README.md b/README.md index 6119a0e..0215f0c 100644 --- a/README.md +++ b/README.md @@ -95,10 +95,6 @@ conda create --name python=3.9 conda activate ``` -- Go to the develop branch (temporary) - -``` -git checkout develop ``` - Install the package From f1f0e83e910848260beb0e321fd6c160c60bd596 Mon Sep 17 00:00:00 2001 From: Charles Pageot <167803616+CharlesPageot@users.noreply.github.com> Date: Fri, 12 Jul 2024 12:13:02 -0400 Subject: [PATCH 16/21] Update README.md Co-authored-by: Julien Cohen-Adad --- README.md | 2 +- 1 file changed, 1 insertion(+), 1 deletion(-) diff --git a/README.md b/README.md index 0215f0c..2bded93 100644 --- a/README.md +++ b/README.md @@ -127,7 +127,7 @@ analytical_cases -t "spherical" ### compute_fieldmap -The _compute_fieldmap_ command allows computation of a $B_0$ fieldmap based on a susceptibility distribution given as an input. +The `compute_fieldmap` command allows computation of a $B_0$ fieldmap based on a susceptibility distribution given as an input. **Arguments** - input_file : path to the susceptibility distribution (NIfTI file) From 9b22e764b34e3a6ac5e9d482b9e3bc5f54243850 Mon Sep 17 00:00:00 2001 From: Charles Pageot <167803616+CharlesPageot@users.noreply.github.com> Date: Fri, 12 Jul 2024 12:13:29 -0400 Subject: [PATCH 17/21] Update README.md Co-authored-by: Julien Cohen-Adad --- README.md | 2 +- 1 file changed, 1 insertion(+), 1 deletion(-) diff --git a/README.md b/README.md index 2bded93..a015d6d 100644 --- a/README.md +++ b/README.md @@ -117,7 +117,7 @@ The _analytical_cases_ command allows for comparaison between simulated and anal - -b, buffer (optional, default=2): Buffer value for zero-padding around the phantom **Return** -Plots to visialize the results +Plots to visualize the results Example: ``` From a0a5e5c093b0887d325240eec5328d8c5a2b3aca Mon Sep 17 00:00:00 2001 From: charles pageot Date: Fri, 12 Jul 2024 12:27:29 -0400 Subject: [PATCH 18/21] Updated README.md --- README.md | 40 +++++++++++++++++++++------------------- 1 file changed, 21 insertions(+), 19 deletions(-) diff --git a/README.md b/README.md index a015d6d..f9ad13d 100644 --- a/README.md +++ b/README.md @@ -95,8 +95,6 @@ conda create --name python=3.9 conda activate ``` -``` - - Install the package ``` @@ -108,38 +106,42 @@ You will need to ```conda activate ``` each time you w ## Usage Once the package is installed, the commands can be run directly from the terminal. Here is the description of the two commands available. -### analytical_cases +### compute_fieldmap -The _analytical_cases_ command allows for comparaison between simulated and analytical results for a spherical and cylindrical phantom. +The `compute_fieldmap` command allows computation of a $B_0$ fieldmap based on a susceptibility distribution given as an input. -**Options** -- -t, geometry type : 'spherical' or 'cylindrical' -- -b, buffer (optional, default=2): Buffer value for zero-padding around the phantom +**Inputs** +- input_file : path to the susceptibility distribution (NIfTI file) +- output_file : path for the fieldmap (NIfTI file) -**Return** -Plots to visualize the results +**Output** +The calculated fieldmap at the specified path. Example: ``` -analytical_cases -t "spherical" +compute_fieldmap "inpath/susceptibility_distribution.nii.gz" "outpath/fieldmap.nii.gz" ``` +### analytical_cases -### compute_fieldmap - -The `compute_fieldmap` command allows computation of a $B_0$ fieldmap based on a susceptibility distribution given as an input. +The _analytical_cases_ command allows for comparaison between simulated and analytical results for a spherical and cylindrical phantom. -**Arguments** -- input_file : path to the susceptibility distribution (NIfTI file) -- output_file : path for the fieldmap (NIfTI file) +**Inputs** +- -t, geometry type : 'spherical' or 'cylindrical' +- -b, buffer (optional, default=2): Buffer value for zero-padding around the phantom -**Return** -The calculated fieldmap at the specified path. +**Outputs** +Plots to visualize the results Example: ``` -compute_fieldmap "inpath/susceptibility_distribution.nii.gz" "outpath/fieldmap.nii.gz" +analytical_cases -t "spherical" ``` +The figures generated would be + +![alt text](Figure_1.png) + +![alt text](Figure_2.png) ## References From cb2cf79090afd9889b53915c520689e4d5fb9f55 Mon Sep 17 00:00:00 2001 From: charles pageot Date: Fri, 12 Jul 2024 13:51:59 -0400 Subject: [PATCH 19/21] Created a parent class for plotting --- fft_simulation/analytical_cases.py | 267 ++++++++++++++--------------- 1 file changed, 125 insertions(+), 142 deletions(-) diff --git a/fft_simulation/analytical_cases.py b/fft_simulation/analytical_cases.py index 7aec733..7d20bd6 100644 --- a/fft_simulation/analytical_cases.py +++ b/fft_simulation/analytical_cases.py @@ -4,7 +4,113 @@ from scipy.ndimage import rotate import click -class Spherical: +class Visualization: + """ + A class that provides visualization methods for susceptibility distribution + and Bz field variation. + """ + + def plot_figure1(self, sus_dist, simulated_Bz, geometry_type): + """ + Plot the susceptibility distribution and Bz field variation for a given geometry type. + + Parameters: + sus_dist (ndarray): Array representing the susceptibility distribution. + simulated_Bz (ndarray): Array representing the Bz field variation. + geometry_type (str): Type of geometry. + + Returns: + None + """ + + dimensions_1 = np.array(sus_dist.shape) + dimensions_2 = np.array(simulated_Bz.shape) + + fig, axes = plt.subplots(2, 3, figsize=(10, 5), dpi=120) + fig.suptitle(f'Susceptibility distribution (top) and Bz field variation (bottom) for a {geometry_type} geometry') + + h = axes[0,0].imshow(sus_dist[dimensions_1[0] // 2, :, :], origin='lower') + axes[0,0].set_title('Y-Z plane') + axes[0,0].axis("off") + plt.colorbar(h, label='Susceptibility [ppm]') + + h = axes[0,1].imshow(sus_dist[:, dimensions_1[0] // 2, :], origin='lower') + axes[0,1].set_title('Z-X plane') + axes[0,1].axis("off") + plt.colorbar(h, label='Susceptibility [ppm]') + + h = axes[0,2].imshow(sus_dist[:, :, dimensions_1[0] // 2], origin='lower') + axes[0,2].set_title('X-Y plane') + axes[0,2].axis("off") + plt.colorbar(h, label='Susceptibility [ppm]') + + # plot section of the b0 field variation + + vmin = np.min(simulated_Bz)*1.1 + vmax = np.max(simulated_Bz)*1.1 + + h = axes[1,0].imshow(simulated_Bz[dimensions_2[0] // 2, :, :], vmin=vmin, vmax=vmax, origin='lower') + axes[1,0].set_title('Y-Z plane') + axes[1,0].axis("off") + plt.colorbar(h, label='Bz [T]') + + h = axes[1,1].imshow(simulated_Bz[:, dimensions_2[0] // 2, :], vmin=vmin, vmax=vmax, origin='lower') + axes[1,1].set_title('Z-X plane') + axes[1,1].axis("off") + plt.colorbar(h, label='Bz [T]') + + h = axes[1,2].imshow(simulated_Bz[:, :, dimensions_2[0] // 2], vmin=vmin, vmax=vmax, origin='lower') + axes[1,2].set_title('X-Y plane') + axes[1,2].axis("off") + plt.colorbar(h, label='Bz [T]') + + plt.show() + + def plot_figure2(self, Bz_analytical, simulated_Bz, geometry_type): + """ + Plot the analytical solution and simulated results for the Bz field variation. + + Parameters: + Bz_analytical (ndarray): Array representing the analytical solution for the Bz field variation. + simulated_Bz (ndarray): Array representing the simulated Bz field variation. + geometry_type (str): Type of geometry. + + Returns: + None + """ + vmin = np.min(simulated_Bz)*1.1 + vmax = np.max(simulated_Bz)*1.1 + + dimensions = np.array(Bz_analytical.shape) + + fig, axes = plt.subplots(1, 3, figsize=(10, 3), dpi=120) + fig.suptitle(f'Analytical solution and simulated results for the Bz field variation for a {geometry_type} geometry') + + axes[0].plot(np.linspace(-dimensions[0]//2, dimensions[0]//2, dimensions[0]), Bz_analytical[:, dimensions[0]//2, dimensions[0]//2], label='Theory') + 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].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].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.show() + +class Spherical(Visualization): """ Represents a spherical object in 3D space. @@ -67,17 +173,13 @@ def analyticial_sol(self): r = np.sqrt(x**2 + y**2 + z**2) - Bz_out = self.sus_diff/3 * (self.R/r)**3 * (3*z**2/r**2 - 1) - Bz_out[mask] = 0 - - Bz_in = np.zeros(self.matrix) - - Bz_analytical = Bz_out + Bz_in + Bz_analytical = self.sus_diff/3 * (self.R/r)**3 * (3*z**2/r**2 - 1) + Bz_analytical[mask] = 0 # set the field inside the sphere to zero return Bz_analytical -class Cylindrical: +class Cylindrical(Visualization): """ Represents a cylindrical object in 3D space. @@ -163,7 +265,13 @@ def analytical_sol(self): Bz_analytical_x = Bz_out_x + Bz_in Bz_analytical_y = Bz_out_y + Bz_in - return Bz_analytical_x, Bz_analytical_y + # Create a single 3D array for the analytical solution. Only the lines along the cartesian axes are filled + Bz_analytical = np.zeros(self.matrix) + Bz_analytical[:, self.matrix[0] //2, self.matrix[0] //2] = Bz_analytical_x[:, self.matrix[0] //2, self.matrix[0] //2] + Bz_analytical[self.matrix[0] //2, :, self.matrix[0] //2] = Bz_analytical_y[self.matrix[0] //2, :, self.matrix[0] //2] + Bz_analytical[self.matrix[0] //2, self.matrix[0] //2, :] = Bz_analytical_x[self.matrix[0] //2, self.matrix[0] //2, :] + + return Bz_analytical @click.command(help="Compare the analytical solution to the simulated solution for a spherical or cylindrical geometry.") @click.option('-t', '--geometry-type',required=True, @@ -207,75 +315,14 @@ def compare_to_analytical(geometry_type, buffer): # analytical solution Bz_analytical = sphere.analyticial_sol() - # plot sections of the chi distribution - fig, axes = plt.subplots(2, 3, figsize=(10, 5), dpi=120) - fig.suptitle('Susceptibility distribution (top) and Bz field variation (bottom) for a spherical geometry') - h = axes[0,0].imshow(sus_dist[matrix[0] // 2, :, :], origin='lower') - axes[0,0].set_title('Y-Z plane') - axes[0,0].axis("off") - plt.colorbar(h, label='Susceptibility [ppm]') - - h = axes[0,1].imshow(sus_dist[:, matrix[0] // 2, :], origin='lower') - axes[0,1].set_title('Z-X plane') - axes[0,1].axis("off") - plt.colorbar(h, label='Susceptibility [ppm]') - - h = axes[0,2].imshow(sus_dist[:, :, matrix[0] // 2], origin='lower') - axes[0,2].set_title('X-Y plane') - axes[0,2].axis("off") - plt.colorbar(h, label='Susceptibility [ppm]') + sphere.plot_figure1(sus_dist, calculated_Bz, geometry_type) + sphere.plot_figure2(Bz_analytical, calculated_Bz, geometry_type) - # plot section of the b0 field variation - - vmin = np.min(calculated_Bz)*1.1 - vmax = np.max(calculated_Bz)*1.1 - - h = axes[1,0].imshow(calculated_Bz[63, :, :],vmin=vmin, vmax=vmax, origin='lower') - axes[1,0].axis("off") - plt.colorbar(h, label='Field variation [ppm]') - - h = axes[1,1].imshow(calculated_Bz[:, 63, :],vmin=vmin, vmax=vmax, origin='lower') - axes[1,1].axis("off") - plt.colorbar(h, label='Field variation [ppm]') - - h = axes[1,2].imshow(calculated_Bz[:, :, 63],vmin=vmin, vmax=vmax, origin='lower') - axes[1,2].axis("off") - plt.colorbar(h, label='Field variation [ppm]') - plt.tight_layout() - plt.show() - - # comparaison with the analytical solution - - fig, axes = plt.subplots(1, 3, figsize=(10, 3), dpi=120) - fig.suptitle('Analytical solution and simulated results for the Bz field variation for a spherical geometry') - - axes[0].plot(np.linspace(-64,64, 128), Bz_analytical[:,63,63], label='Theory') - axes[0].plot(np.linspace(-64,64, 128), calculated_Bz[:,63,63],'--', label='Simulated') - axes[0].set_ylim(vmin, vmax) - axes[0].set_xlabel('x position [mm]') - axes[0].set_ylabel('Field variation [ppm]') - axes[0].legend() - - axes[1].plot(np.linspace(-64,64, 128), Bz_analytical[63,:,63], label='Theory') - axes[1].plot(np.linspace(-64,64, 128), calculated_Bz[63,:,63],'--', label='Simulated') - axes[1].set_ylim(vmin, vmax) - axes[1].set_xlabel('y position [mm]') - axes[1].set_ylabel('Field variation [ppm]') - axes[1].legend() - - axes[2].plot(np.linspace(-64,64, 128), Bz_analytical[63,63,:], label='Theory') - axes[2].plot(np.linspace(-64,64, 128), calculated_Bz[63,63,:],'--', label='Simulated') - axes[2].set_ylim(vmin, vmax) - axes[2].set_xlabel('z position [mm]') - axes[2].set_ylabel('Field variation [ppm]') - axes[2].legend() - - plt.tight_layout() - plt.show() - - else: # type is "cylindrical" + + else: + geometry_type = 'cylindrical' # create the susceptibility geometry cylinder = Cylindrical(matrix, image_res, R, sus_diff) sus_dist = cylinder.volume() @@ -283,71 +330,7 @@ def compare_to_analytical(geometry_type, buffer): # compute Bz variation calculated_Bz = compute_bz(sus_dist, image_res, buffer) - Bz_analytical_x, Bz_analytical_y = cylinder.analytical_sol() - - # plot sections of the chi distribution - fig, axes = plt.subplots(2, 3, figsize=(10, 5), dpi=120) - fig.suptitle('Susceptibility distribution (top) and Bz field variation (bottom) for a cylindrical geometry') - - h = axes[0,0].imshow(sus_dist[matrix[0] // 2, :, :], origin='lower') - axes[0,0].set_title('Y-Z plane') - axes[0,0].axis("off") - plt.colorbar(h, label='Susceptibility [ppm]') + Bz_analytical = cylinder.analytical_sol() - h = axes[0,1].imshow(sus_dist[:, matrix[0] // 2, :], origin='lower') - axes[0,1].set_title('X-Z plane') - axes[0,1].axis("off") - plt.colorbar(h, label='Susceptibility [ppm]') - - h = axes[0,2].imshow(sus_dist[:, :, matrix[0] // 2], origin='lower') - axes[0,2].set_title('X-Y plane') - axes[0,2].axis("off") - plt.colorbar(h, label='Susceptibility [ppm]') - - # plot section of the b0 field variation - - vmin = np.min(calculated_Bz)*1.1 - vmax = np.max(calculated_Bz)*1.1 - - h = axes[1,0].imshow(calculated_Bz[63, :, :],vmin=vmin, vmax=vmax, origin='lower') - axes[1,0].axis("off") - plt.colorbar(h, label='Field variation [ppm]') - - h = axes[1,1].imshow(calculated_Bz[:, 63, :],vmin=vmin, vmax=vmax, origin='lower') - axes[1,1].axis("off") - plt.colorbar(h, label='Field variation [ppm]') - - h = axes[1,2].imshow(calculated_Bz[:, :, 63],vmin=vmin, vmax=vmax, origin='lower') - axes[1,2].axis("off") - plt.colorbar(h, label='Field variation [ppm]') - plt.tight_layout() - plt.show() - - # comparaison with the analytical solution - - fig, axes = plt.subplots(1, 3, figsize=(10, 3), dpi=120) - fig.suptitle('Analytical solution and simulated results for the Bz field variation for a cylindrical geometry') - - axes[0].plot(np.linspace(-64,64, 128), Bz_analytical_x[:,63,63], label='Theory') - axes[0].plot(np.linspace(-64,64, 128), calculated_Bz[:,63,63],'--', label='Simulated') - axes[0].set_ylim(vmin, vmax) - axes[0].set_xlabel('x position [mm]') - axes[0].set_ylabel('Field variation [ppm]') - axes[0].legend() - - axes[1].plot(np.linspace(-64,64, 128), Bz_analytical_y[63,:,63], label='Theory') - axes[1].plot(np.linspace(-64,64, 128), calculated_Bz[63,:,63],'--', label='Simulated') - axes[1].set_ylim(vmin, vmax) - axes[1].set_xlabel('y position [mm]') - axes[1].set_ylabel('Field variation [ppm]') - axes[1].legend() - - axes[2].plot(np.linspace(-64,64, 128), Bz_analytical_x[63,63,:], label='Theory') - axes[2].plot(np.linspace(-64,64, 128), calculated_Bz[63,63,:],'--', label='Simulated') - axes[2].set_ylim(vmin, vmax) - axes[2].set_xlabel('z position [mm]') - axes[2].set_ylabel('Field variation [ppm]') - axes[2].legend() - - plt.tight_layout() - plt.show() + cylinder.plot_figure1(sus_dist, calculated_Bz, geometry_type) + cylinder.plot_figure2(Bz_analytical, calculated_Bz, geometry_type) \ No newline at end of file From c0ad4c8e505122a68f1fe8bf0acb7e34c98aa539 Mon Sep 17 00:00:00 2001 From: charles pageot Date: Fri, 12 Jul 2024 15:16:49 -0400 Subject: [PATCH 20/21] Updated analytical_cases.py --- fft_simulation/analytical_cases.py | 43 +++++++++++------------------- 1 file changed, 15 insertions(+), 28 deletions(-) diff --git a/fft_simulation/analytical_cases.py b/fft_simulation/analytical_cases.py index 7d20bd6..705301c 100644 --- a/fft_simulation/analytical_cases.py +++ b/fft_simulation/analytical_cases.py @@ -10,7 +10,7 @@ class Visualization: and Bz field variation. """ - def plot_figure1(self, sus_dist, simulated_Bz, geometry_type): + def plot_susceptibility_and_fieldmap(self, sus_dist, simulated_Bz, geometry_type): """ Plot the susceptibility distribution and Bz field variation for a given geometry type. @@ -66,7 +66,7 @@ def plot_figure1(self, sus_dist, simulated_Bz, geometry_type): plt.show() - def plot_figure2(self, Bz_analytical, simulated_Bz, geometry_type): + def plot_comparaison_analytical(self, Bz_analytical, simulated_Bz, geometry_type): """ Plot the analytical solution and simulated results for the Bz field variation. @@ -157,7 +157,7 @@ def volume(self): """ return np.where(self.mask() == True, self.sus_diff, 0) - def analyticial_sol(self): + def analytical_sol(self): """ Calculates the analytical solution for the magnetic field inside and outside the sphere. @@ -306,31 +306,18 @@ def compare_to_analytical(geometry_type, buffer): R = 15 # mm sus_diff = 9 # ppm - if geometry_type == 'spherical': - # create the susceptibility geometry - sphere = Spherical(matrix, image_res, R, sus_diff) - sus_dist = sphere.volume() - # compute Bz variation - calculated_Bz = compute_bz(sus_dist, image_res, buffer) - # analytical solution - Bz_analytical = sphere.analyticial_sol() - - - sphere.plot_figure1(sus_dist, calculated_Bz, geometry_type) - sphere.plot_figure2(Bz_analytical, calculated_Bz, geometry_type) - - - else: - - geometry_type = 'cylindrical' - # create the susceptibility geometry - cylinder = Cylindrical(matrix, image_res, R, sus_diff) - sus_dist = cylinder.volume() + dicto = {'spherical': Spherical(matrix, image_res, R, sus_diff), + 'cylindrical': Cylindrical(matrix, image_res, R, sus_diff)} - # compute Bz variation - calculated_Bz = compute_bz(sus_dist, image_res, buffer) + # create the susceptibility geometry + geometry = dicto[geometry_type] + sus_dist = geometry.volume() - Bz_analytical = cylinder.analytical_sol() + # compute Bz variation + calculated_Bz = compute_bz(sus_dist, image_res, buffer) + # analytical solution + Bz_analytical = geometry.analytical_sol() - cylinder.plot_figure1(sus_dist, calculated_Bz, geometry_type) - cylinder.plot_figure2(Bz_analytical, calculated_Bz, geometry_type) \ No newline at end of file + # plot the results + geometry.plot_susceptibility_and_fieldmap(sus_dist, calculated_Bz, geometry_type) + geometry.plot_comparaison_analytical(Bz_analytical, calculated_Bz, geometry_type) From b7e722a6dac5adf52d60b1fe5f49caa53fd777c3 Mon Sep 17 00:00:00 2001 From: Charles Pageot <167803616+CharlesPageot@users.noreply.github.com> Date: Mon, 15 Jul 2024 10:25:22 -0400 Subject: [PATCH 21/21] Update README.md --- README.md | 6 ++++-- 1 file changed, 4 insertions(+), 2 deletions(-) diff --git a/README.md b/README.md index f9ad13d..cc21b74 100644 --- a/README.md +++ b/README.md @@ -139,9 +139,11 @@ analytical_cases -t "spherical" ``` The figures generated would be -![alt text](Figure_1.png) +![Figure_1](https://github.com/user-attachments/assets/a2f49a83-dc8f-4d94-a5ba-075e7b8675aa) + + +![Figure_2](https://github.com/user-attachments/assets/c4c2d34c-7153-43d8-ad21-8af6f49427de) -![alt text](Figure_2.png) ## References