20172674721

Committed 12 Dec 2025 04:10PM UTC coverage: 81.818%. First build

Build # 20172674721

Build Type

Pull #41

github

Committed by

mathieuboudreau

Commit Message

Fix buffer=0 edge case in compute_bz function

Pull Request Pull Request #41: Odd/Even matrix size fix + new padding features

Run Details

50 of 56 new or added lines in 2 files covered. (89.29%)

180 of 220 relevant lines covered (81.82%)

1.64 hits per line

Source File
Press 'n' to go to next uncovered line, 'b' for previous

76.47

/functions/compute_fieldmap.py

import numpy as np
import nibabel as nib
import click
from time import perf_counter
from pathlib import Path

def is_nifti(filepath):
    """
    Check if the given filepath represents a NIfTI file.

    Args:
        filepath (str): The path of the file to check.

    Returns:
        bool: True if the file is a NIfTI file, False otherwise.
    """
    if filepath.endswith('.nii') or filepath.endswith('.nii.gz'):
        return True
    else:
        return False

def load_sus_dist(filepath):
    """
    Load the susceptibility distribution from a given file.

    Args:
        filepath (str): The path to the file containing the susceptibility distribution.

    Returns:
        tuple: A tuple containing the loaded susceptibility distribution as a numpy array 
            and the image resolution as a numpy array.
    """

    image = nib.load(filepath)
    susceptibility_distribution = image.get_fdata()
    header = image.header 
    image_resolution = np.array(header.get_zooms())
    affine_matrix = image.affine

    return susceptibility_distribution, image_resolution, affine_matrix


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.

    Args:
        susceptibility_distribution (numpy.ndarray): The 3D array representing the susceptibility distribution.

        image_resolution (numpy.ndarray, optional): The resolution of the image in each dimension. Defaults to [1, 1, 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)

    kmax = 1/(2*image_resolution)

    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'):
        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[0,0,0] = FFT_chi[0,0,0] + np.prod(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))

    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

def save_to_nifti(data, affine_matrix, output_path):
    """
    Save data to NIfTI format to a specified output path.

    Args:
        data (np.array): the data to save
        affine_matrix (np.array): the affine matrix of the data
        output_path (str): the output path to save the file

    Returns:
        None
    """
    data=data.astype(np.float32)
    nifti_image = nib.Nifti1Image(data, affine_matrix)
    nib.save(nifti_image, output_path)



@click.command(help="Compute the magnetic field variation in ppm from a susceptibility distribution in NIfTI format.")
@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")
@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.

    Args:
        input_file (str): Path to the susceptibility distribution in NIfTI format.
        output_file (str): Path for the computed fieldmap in NIfTI format.

    Returns:
        None
    """
    if is_nifti(input_file):
        start_time = perf_counter()
        print('Start')
        susceptibility_distribution, image_resolution, affine_matrix = load_sus_dist(input_file)
        print('Susceptibility distribution loaded')
        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.")

    

1	import numpy as np	2✔
2	import nibabel as nib	2✔
3	import click	2✔
4	from time import perf_counter	2✔
5	from pathlib import Path	2✔
6
7	def is_nifti(filepath):	2✔
8	"""
9	Check if the given filepath represents a NIfTI file.
10
11	Args:
12	filepath (str): The path of the file to check.
13
14	Returns:
15	bool: True if the file is a NIfTI file, False otherwise.
16	"""
17	if filepath.endswith('.nii') or filepath.endswith('.nii.gz'):	2✔
18	return True	2✔
19	else:
20	return False	2✔
21
22	def load_sus_dist(filepath):	2✔
23	"""
24	Load the susceptibility distribution from a given file.
25
26	Args:
27	filepath (str): The path to the file containing the susceptibility distribution.
28
29	Returns:
30	tuple: A tuple containing the loaded susceptibility distribution as a numpy array
31	and the image resolution as a numpy array.
32	"""
33
34	image = nib.load(filepath)	2✔
35	susceptibility_distribution = image.get_fdata()	2✔
36	header = image.header	2✔
37	image_resolution = np.array(header.get_zooms())	2✔
38	affine_matrix = image.affine	2✔
39
40	return susceptibility_distribution, image_resolution, affine_matrix	2✔
41
42
43	def compute_bz(susceptibility_distribution, image_resolution=np.array([1,1,1]), buffer=50, mode='edge'):	2✔
44	"""
45	Compute the Bz field variation in ppm based on a susceptibility distribution
46	using a Fourier-based method.
47
48	Args:
49	susceptibility_distribution (numpy.ndarray): The 3D array representing the susceptibility distribution.
50
51	image_resolution (numpy.ndarray, optional): The resolution of the image in each dimension. Defaults to [1, 1, 1].
52
53	buffer (int, optional): The buffer size (voxels) for the k-space grid on each size. Defaults to 50.
54
55	mode (str, optional): The padding mode for the susceptibility distribution. Defaults to 'edge'.
56
57	Returns:
58	volume_without_buffer (numpy.ndarray): The computed magnetic field Bz in ppm.
59
60	"""
61
62	# Pad the susceptibility distribution
63
64	if mode == 'b0SimISMRM':	2✔
NEW 65	susceptibility_distribution=np.pad(susceptibility_distribution, ((buffer, buffer), (0,buffer), (buffer,buffer)), mode='edge')	×
NEW 66	susceptibility_distribution=np.pad(susceptibility_distribution, ((0,0),(buffer,0),(0,0)), mode='constant', constant_values=0.35)	×
67	else:
68	susceptibility_distribution=np.pad(susceptibility_distribution, buffer, mode=mode)	2✔
69
70	# dimensions needs to be a numpy.array
71	dimensions = np.array(susceptibility_distribution.shape)	2✔
72
73	kmax = 1/(2*image_resolution)	2✔
74
75	interval = 2 * kmax / dimensions	2✔
76
77
78	kx_min_shift = (dimensions[0]%2)*interval[0]/2	2✔
79	ky_min_shift = (dimensions[1]%2)*interval[1]/2	2✔
80	kz_min_shift = (dimensions[2]%2)*interval[2]/2	2✔
81
82	kx_max_shift = -interval[0] + (dimensions[0]%2)*interval[0]/2	2✔
83	ky_max_shift = -interval[1] + (dimensions[1]%2)*interval[1]/2	2✔
84	kz_max_shift = -interval[2] + (dimensions[2]%2)*interval[2]/2	2✔
85
86
87	[kx, ky, kz] = np.meshgrid(np.linspace(-kmax[0] + kx_min_shift, kmax[0] + kx_max_shift, dimensions[0]),	2✔
88	np.linspace(-kmax[1] + ky_min_shift, kmax[1] + ky_max_shift, dimensions[1]),
89	np.linspace(-kmax[2] + kz_min_shift, kmax[2] + kz_max_shift, dimensions[2]), indexing='ij')
90
91	# FFT procedure
92	# undetermined at the center of k-space
93	k2 = kx2 + ky2 + kz**2	2✔
94
95	with np.errstate(divide='ignore', invalid='ignore'):	2✔
96	x_kernel = 1/3 - kz**2/k2	2✔
97
98	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	2✔
99
100	kernel = np.fft.ifftshift(x_kernel)	2✔
101
102	FFT_chi = np.fft.fftn(susceptibility_distribution, dimensions)	2✔
103
104	FFT_chi[0,0,0] = FFT_chi[0,0,0] + np.prod(dimensions)*susceptibility_distribution[0,0,0]	2✔
105
106	Bz_fft = kernel*FFT_chi	2✔
107
108	# retrive the inital FOV
109
110	volume_with_buffer = np.real(np.fft.ifftn(Bz_fft))	2✔
111
112	if buffer == 0:	2✔
113	volume_without_buffer = volume_with_buffer	2✔
114	else:
115	volume_without_buffer = volume_with_buffer[buffer:-buffer, buffer:-buffer, buffer:-buffer]	2✔
116
117	return volume_without_buffer	2✔
118
119	def save_to_nifti(data, affine_matrix, output_path):	2✔
120	"""
121	Save data to NIfTI format to a specified output path.
122
123	Args:
124	data (np.array): the data to save
125	affine_matrix (np.array): the affine matrix of the data
126	output_path (str): the output path to save the file
127
128	Returns:
129	None
130	"""
131	data=data.astype(np.float32)	2✔
132	nifti_image = nib.Nifti1Image(data, affine_matrix)	2✔
133	nib.save(nifti_image, output_path)	2✔
134
135
136
137	@click.command(help="Compute the magnetic field variation in ppm from a susceptibility distribution in NIfTI format.")	2✔
138	@click.option('-i','--input','input_file', type=click.Path(exists=True), required=True,	2✔
139	help="Input susceptibility distribution, supported extensions: .nii, .nii.gz")
140	@click.option('-o', '--output', 'output_file', type=click.Path(), default='fieldmap.nii.gz',	2✔
141	help="Output fieldmap, supported extensions: .nii, .nii.gz")
142	@click.option('-b', '--buffer', 'buffer', type=int, default=50, required=False,	2✔
143	help="Buffer size (voxels) for the k-space grid on each side")
144	@click.option('-m', '--mode', 'mode', type=str, default='edge', required=False,	2✔
145	help="Padding mode for the susceptibility distribution")
146	def compute_fieldmap(input_file, output_file, buffer, mode):	2✔
147	"""
148	Main procedure for performing the simulation.
149
150	Args:
151	input_file (str): Path to the susceptibility distribution in NIfTI format.
152	output_file (str): Path for the computed fieldmap in NIfTI format.
153
154	Returns:
155	None
156	"""
157	if is_nifti(input_file):	×
158	start_time = perf_counter()	×
159	print('Start')	×
160	susceptibility_distribution, image_resolution, affine_matrix = load_sus_dist(input_file)	×
161	print('Susceptibility distribution loaded')	×
NEW 162	fieldmap = compute_bz(susceptibility_distribution, image_resolution, buffer, mode)	×
163	print('Fieldmap simulated')	×
164
165	# Check if all subdirectories exist and create them if not
NEW 166	output_path = Path(output_file)	×
NEW 167	output_path.parent.mkdir(parents=True, exist_ok=True)	×
168
169	save_to_nifti(fieldmap, affine_matrix, output_file)	×
170	print('Saving to NIfTI format')	×
171	end_time = perf_counter()	×
172	print(f'End. Runtime: {end_time-start_time:.2f} seconds')	×
173	else:
174	print("The input file must be NIfTI.")	×
175
176

shimming-toolbox / susceptibility-to-fieldmap-fft / 20172674721

Source File Press 'n' to go to next uncovered line, 'b' for previous

Source File
Press 'n' to go to next uncovered line, 'b' for previous