Code Documentation

Phantom Designer

class mri_distortion_toolkit.PhantomDesigner.PhantomSlice(slice_shape='rectangle', slice_thickness=30, HVL_x=250, HVL_Y=250, z_pos=0, hole_depth=17, hole_spacing=25, hole_radius=4.35, hole_stop=None, GridOffset=0, HoleCentroids='cartesian', DSV=150, LoadRegion=None, ReferenceCrosshairRadius=None, GuideRods=None, hole_start=None, bottom_cut=0)

Draws a single slice of the phantom All dimensions in mm All inputs are optional; calling without defining will just use default value.

Parameters:

  • slice_shape (string) – Shape of base slice. Currently can only be ‘rectangle’ or ‘ellipse’ but others could be added

  • slice_thickness (float) – Thickness of slice in Z

  • HVL_x (float) – Half value of slice in X (Left/Right). For ellipse this is equivalent to major radius

  • HVL_Y (float) – Half value of slice in Y (Up/Down). For ellipse this is equivalent to major radius

  • hole_depth (float) – Depth of holes. If value is > slice_thickness will cut completely through

  • hole_spacing (float) – distance between hole centroids - only approximate, as there are integers involved!

  • hole_radius (float) – Radius of holes

  • hole_start (None or float) – x/y coordinate to start holes at. Default is None, which will try and choose a sensible value.

  • hole_stop (None or float) – x/y coordinate to stop holes at. Default is None, which will try and choose a sensible value.

  • GridOffset (float) – move the entire hole distrubution by this much. Allows for offset distributions in differenct slices

  • HoleCentroids (string ('cartesian' or 'ROI_polar') or 2D list/numpy array) – ‘cartesian’ or ‘ROI_polar’ will produce either cartesian or polar grids based on hole_spacing. User can instead specify the location of all hole centroids as a 2D list of X/Y points e.g. [[x1,x2,x3],[y1,y2,y3]]. The value of GridOffset will be applied to these points. If user values are entered, hole_stop and hole_start are ignored.

  • z_pos (float) – z position of the middle of the hole cutout.

  • DSV (None or float) – Circles representing this raddi at a given z value will be etched on each slice; if you choose ROI_polar for marker positions, the code will attempt to ensure good coverage over the DSV surface`

  • LoadRegion (None or Dict) – Can specify a cylindrical or spherical load region which will be cut away geom each slice, e.g. | LoadRegion = {‘shape’: ‘cylinder’, ‘radius’: 100, ‘height’: 200} or, | LoadRegion = {‘shape’: ‘sphere’, ‘radius’: 100} (this is not well tested), or | LoadRegion={‘shape’: ‘rectangle’, ‘width’: 70, ‘height’: 200}, or | LoadRegion = None (default) | You can can also optionally add a key z_offset to move the load position along the z axis

  • GuideRods (None or Dict) – Adds one guide rod in each corner. Position is distance from edge. GuideRods = {‘radius’: 10, ‘position’: 20,’length’: 350} See: https://www.essentracomponents.com/en-au/p/fully-threaded-studs-rods/0060110000vr

  • ReferenceCrosshairRadius (float) – if not None, a crosshair of markers will be drawn inside ReferenceCrosshairRadius. Other markers will not be added in this region.

  • bottom_cut (float) – this cuts away the bottom of the slice. Since MRI scanners often have limits on how low the bed can go, it can be useful to lower the phantom as much as possible.

add_full_scale_drawing()

Add a basic drawing of the phantom that can be printed off and used to manually drill the holes

ToDo:: Figure out which paper size needs to be used depending on phantom size

draw_DSV()
draw_Load()

Draw the load. To assist in visualisation.

draw_slice()

This is the method which actually draws the slice in FreeCAD Prior to calling this method, the code should work without any calls to FreeCAD :return:

Marker Analysis

class mri_distortion_toolkit.MarkerAnalysis.MarkerVolume(input_data, ImExtension='dcm', r_min=None, r_max=None, iterative_segmentation=False, threshold=None, n_markers_expected=None, fat_shift_direction=None, verbose=False, gaussian_image_filter_sd=1, correct_fat_water_shift=False, marker_size_lower_tol=0.9, marker_size_upper_tol=1)

This class excepts either a path to a folder of dicom files or a numpy array of marker positions. In the former case, it will attempt to automatically extract the marker positions. In the latter, it will just use the input marker positions. Once the instance is created, it contains a number of handy plots and functions.

Parameters:

  • input_data_path (path, string, np.ndarray) – Either a path to a folder of dicoms, or a numpy array with three columns corresponding to the [x,y,z] marker positions

  • ImExtension (str, optional) – change to whatever extension your dicom files, most commonly either dcm or IMA

  • r_min (float ,optional) – any markers at radii less than this will be discounted.

  • r_max (float, optional) – any markers at radii greater than this will be discounted.

  • iterative_segmentation (bool, optional) – If set to true, a slower iterative method will be used to find the threshold value used to segment the markers. Otherwise, otsu’s method is used.

  • threshold (float, optional) – Manually set the threshold value. If this is left as None otsu’s method is used. This will replace the iterative segmentation method regardless of flag.

  • n_markers_expected (int, optional) – if you know how many markers you expect, you can enter the value here. The code will then warn you if it finds a different number.

  • gaussian_image_filter_sd (float, optional) – size of filter to blur input data with. Higher values are better at removing noise prior to thresholding, but may also blur out small/low signal markers!

  • correct_fat_water_shift (bool, optional) – if True, will attempt to automatically correct marker positions for calculated fat water shift

  • marker_size_lower_tol (float, optional) – lower acceptance window for a marker compared to the median size of markers. e.g. markermarker_size_tol=0.5 will result in markers being accepted within when (1-marker_size_lower_tol)*median_marker_size <= marker_size <= (1+markermarker_size_tol)*marker_size_upper_tol

  • marker_size_upper_tol (float, optional) – upper acceptance window for a marker compared to the median size of markers.

  • verbose (bool, optional) – if True, prints a lot of info about the marker extraction process to the screen. mostly useful during debugging

  • fat_shift_direction (int, optional) – We aren’t sure yet if we can confidently predict the direction that the fat water shift will be in. if you set correct_fat_water_shift=True, you can change the direction it shifts in here. it should take a value of -1, or 1. pull requests making this easier to use very welcome!

export_to_slicer(save_path=None, filename='slicer_centroids')

export a json file that can be read in by slicer, allowing a good way to visualise marker segmentation performance. This file will be saved at the same spot as input data if dicom was input, otherwise it will be saved to save_path.

perturb_marker_positions(random_perturbation, systemic_perturbation=0)

add random noise to the marker positions. Useful to perform sensitivity analysis. operates in place.

Parameters:

  • random_perturbation (float) – upper and lower limits of random_perturbation to each marker coordinates. uniform distribution is used

  • systemic_perturbation (float, optional) – constant offset added to all markers

plot_3D_markers(title='3D marker positions')

Just a quick plot of the marker positions; very useful sanity check!

rotate_markers(xaxis_angle=0, yaxis_angle=0, zaxis_angle=0)

rotate markers using euler angles

Parameters:

  • xaxis_angle – rotation around x angle in degrees

  • yaxis_angle – rotation around y angle in degrees

  • zaxis_angle – rotation around z angle in degrees

save_dicom_data(save_path=None, filename='dicom_data')

save the dicom data as json. This file will be saved at the same spot as input data if dicom was input, otherwise it will be saved to save_path.

translate_markers(x_shift=0, y_shift=0, z_shift=0)

translate markers

Parameters:

  • x_shift – x translation in mm

  • y_shift – y translation in mm

  • z_shift – z translation in mm

Returns:

class mri_distortion_toolkit.MarkerAnalysis.MatchedMarkerVolumes(GroundTruthData, DistortedData, reverse_gradient_data=None, WarpSearchData=True, AutomatchMarkers=True, AllowDoubleMatching=False, sorting_method='radial', n_refernce_markers=0, skip_unmatchable_markers=True)
Parameters:

  • GroundTruthData (MarkerVolume object, or path to dicom folder, or numpy array) –

  • DistortedData (MarkerVolume object, or path to dicom folder, or numpy array) –

  • reverse_gradient_data (None, or MarkerVolume object, or path to dicom folder, or numpy array) –

  • WarpSearchData (boolean) – if True, position of found markers is used to update expected positions for ground truth. Recomended is True if there is substantial distortion present.

  • AutomatchMarkers (boolean) – if True, will automatically match the distorted data to the ground truth using a search algorith. If False, will simply use input data in the order it is entered. This is useful if you have sets of markers that were independently matched through some other process.

  • AllowDoubleMatching (bool, optional) – when False, each ground truth marker is only allowed to match with one distorted marker. However, sometimes this can result in a chain reaction where one marker gets mismatched then every other marker gets mismatched too. In these cases you can set this parameter to True, and any double matched markers will simply be deleted.

  • sorting_method (str, optional) – ‘radial’ or ‘nearest’. This is only important if WarpSearchData is True; in that case, we are building a motion model on-the-fly based on previously seen markers, so the order we see them can be important. For data on a sphere, you may have more success by setting this to ‘nearest’, but this is a bit of an untested feature

  • n_refernce_markers (int, optional) – the n inner_most Reference markers are used to align the ground truth to the MR volume before matching. this can correct for setup errors. Note that we always move the reference. the best way for this to not matter either way is to not have any set up error!

  • skip_unmatchable_markers (bool, optional) – if true, distorted markers where the nearest centroid marker is above the maximum distortion threshold are removed

plot_3D_markers(add_arrows=True, title='3D marker positions')

Works very similarly to the MarkerVolume version, but plots both sets of markers and adds arrows

Parameters:

  • add_arrows (bool, optional) – if True, arrows drawn between matches

  • title (str, optional) – plot title

plot_compressed_markers(z_max=20, z_min=-20, add_arrows=True, title=None)

compresses the 3D markers in the z plane, allowing a 2D visualisation.

Parameters:

  • z_max (float, optional) – maximum z coordiante of markers to includes

  • z_min (float, optional) – minimum z coordiante of markers to includes

  • add_arrows (bool, optional) – if True, arrows drawn between matches

  • title (str, optional) – plot title

Field Calculation

class mri_distortion_toolkit.FieldCalculation.ConvertMatchedMarkersToBz(MatchedCentroids, dicom_data)

an object to calculate and store magnetic fields. Gradient fields are in T/m B0 fields are in T

Parameters:

  • MatchedCentroids – a pandas data frame containing columns [‘x_gt’, ‘y_gt’, ‘z_gt’, ‘x_gnl’, ‘y_gnl’, ‘z_gnl’]. If you want to calculate B0, is should also have [‘x_B0’, ‘y_B0’, ‘z_B0’]. Such a dataframe is normally created within the MatchedMarkerVolumes class. All coordinates in mm.

  • dicom_data (dictionary or path) –

    either a dictionary or a path to a .json file. The minimume required parameters are below. dicom_data is calculated automatically when MR data is used to create a MarkerVolume, and can be exported to json using:

    marker_volume.save_dicom_data().

    dicom_data example:

    dicom_data = {'FOV': [400, 400, 400],
'bandwidth': 203,
'gama': 42.6,
'image_size': [256, 256, 256]}

Harmonic Analysis

class mri_distortion_toolkit.Harmonics.SphericalHarmonicFit(input_Bz_data, r_outer=150, n_order=5, AssessHarmonicPk_Pk=True, QuantifyFit=True, TrimDataBy_r_outer=False, scale=1)

Uses single value decomposition to fit spherical harmonics to InputData.

Let the legendre basis be called L, the Harmonics be H and the magnetic field B. L will have size [n_coords, n_harmonics], Harmonics has size [n_harmonics, 1] and B has size [n_coords, 1]

Our proposition is that there exists some set of harmonics such that

\[L.H = B\]

We can find H by inverting L:

\[L^{-1}.L.H = L^{-1}.B \ H = L^{-1}.B\]

This task is performed using numpys pseudo_inverse functionality to invert L

Parameters:

  • _input_Bz_data (Pandas.DataFrame) – A pandas dataframe with columns x, y, z, Bz. Data should be in mm and T

  • r_outer (float, optional) – radius of sphere of interest. If AssessHarmonicPk_Pk=True, data is reconstucted on r_outer. if TrimDataBy_r_outer=True, data lying outside r_outer will be deleted. Finally, r_outer is used to set the limits on plotting data.

  • n_order (int, optional) – order of harmonic fit. higher order means more terms to fit, which may result in a better fit but also is slower. If the data you input doesn’t satisfy some nyquist criterion (I need to figure this out!) then you will get a spurious result. For most MRI application 8 should be a reasonable number.

  • AssessHarmonicPk_Pk (bool, optional) – if True, data is reconstructed on the surface of r_outer, and the pk-pk contribution from each harmonic calculated. This is a very useful way to interpret harmonics which are otherwise quite abstract.

  • QuantifyFit (bool, optional) – if True, reconstruct data at all input data points, and compare input data point to reconstructed data point. This is useful to check how well the fit has worked.

  • TrimDataBy_r_outer (bool, optional) – if True, any input data outside r_outer is deleted

  • scale (float, optional) – Value to scale the harmonics by. The main intention of this parameter is to normalise gradient field harmonics to gradient strength. harmonics will be multiplied by the value of scale

plot_cut_planes(resolution=2.5, AddColorBar=True, quantity='uT', vmin=None, vmax=None)

Reconstruct the Bz field at the cardinal planes. Note this is basically the same code copied three times in a row, one for each plane.

Parameters:

  • resolution (float, optional) – the resolution of reconstructed data. Default is 2.5 mm

  • AddColorBar (Boolean, optional) – Add a color bar to each plot

  • quantity (string, optional) – quantity to use in reconstruct_Bz; must match the options for that function

  • vmin (float, optional) – same as matplotlib.pyplot.imshow; vmin and vmax control the color map bounds.

  • vmax (float, optional) – same as matplotlib.pyplot.imshow; vmin and vmax control the color map bounds.

plot_harmonics_pk_pk(cut_off=0.1, title=None, return_axs=False, plot_percentage_of_dominant=False, drop_dominant_harmonic=False)

produces a barplot of harmonics.

Parameters:

  • cut_off (float, optional) – cutoff point relative to highest harmonic. e.g. cut_off=.1 means that harmonics which produce less than 10% the pk-pk perturbation of the dominant harmonic are ignored.

  • title – title of plot

  • return_axs – if True, will return the plotting axs for further manipulation

  • plot_percentage_of_dominant – if True, switches from absolute pk-pk to percentage

  • drop_dominant_harmonic – if True, drops dominant harmonic before plotting

print_key_harmonics(cut_off=0.01)

print the harmonics with value > cut_off to the terminal in pk-pk.

Parameters:

cut_off (float, optional) – cutoff point relative to highest harmonic. e.g. cut_off=.1 means that harmonics which produce less than 10% the pk-pk perturbation of the dominant harmonic are ignored.

Reports

class mri_distortion_toolkit.Reports.DefaultTestSuite

these are the tests which are run if no others are specified

test_distortion_less_than_2mm_within_r_100()

distortion test to be run by MRI_QA_Reporter

test_distortion_less_than_4mm_within_r_150()

distortion test to be run by MRI_QA_Reporter

class mri_distortion_toolkit.Reports.Elekta_Distortion_tests

Geometric tests for Elekta system as described here: https://aapm.onlinelibrary.wiley.com/doi/epdf/10.1002/mp.14764

test_distortion_less_than_1_mm_within_r_100_mm()

distortion test to be run by MRI_QA_Reporter

test_distortion_less_than_2_mm_within_r_150_mm()

distortion test to be run by MRI_QA_Reporter

test_distortion_less_than_4_mm_within_r_200_mm()

distortion test to be run by MRI_QA_Reporter

class mri_distortion_toolkit.Reports.MRI_QA_Reporter(MatchedMarkerVolume=None, gradient_harmonics=None, B0_harmonics=None, recon_coords_cartesian=None, tests_to_run=<class 'mri_distortion_toolkit.Reports.DefaultTestSuite'>, dicom_data=None, r_outer=None, show_plots=False, style='dark')

generate a report based on either some matched data, or harmonics. In the latter case, report data is reconstructed.

Parameters:

  • MatchedMarkerVolume (pandas.DataFrame, optional) – pandas dataframe containing [‘x_gt’, ‘y_gt’, ‘z_gt’, ‘x_gnl’, ‘y_gnl’, ‘z_gnl’] You must enter one of MatchedMarkerVolume or gradient_harmonics

  • gradient_harmonics (list, optional) – list of either three pandas series or three paths to csv files from SphericalHarmonicFit. List elements should correspond to [G_x_harmonics, G_y_harmonics, G_z_harmonics]. You must enter one of MatchedMarkerVolume or gradient_harmonics

  • B0_harmonics (pandas.Series or path, optional) – pandas series from SphericalHarmonicFit or path to saved csv file. If entered, is used to report and plot B0 homogeneity

  • recon_coords (pandas.DataFrame, optional) – The coordinates you want to calculate the fields at. a pandas Dataframe which can contain whatever you want, but MUST contain columns called: x, y, z. If left as None a default is generated.

  • tests_to_run (class, optional) – a class containing any tests to run. The class in this instance serves only as a name space. All methods inside the class starting with ‘test’ are identified as tests. The test methods will have access to everything inside this class.

  • dicom_data (dict or path to json, optional) – either a dictionary or path to a saved json file of the dicom_data object generated by MarkerVolume. This data is required for harmonic reconstruction. For MatchedMarkerVolume it is not required, but it does contain useful information about the acquisition so still a good idea to use it.

  • r_outer (float, optional) – for the harmonic approach, r_outer is used to determine where the data is reconstructed. It is also used to control the limits in which data is plotted

  • show_plots (bool, optional) – if True, all plots will be shown in your browser. Otherwise they are only saved

  • style (str) – control the report style. At the moment can only be ‘light’ or ‘dark’

write_html_report(output_folder=None, report_name=None)

Generates a html report encompassing available acquisition information, test results, and intercative plotly figures

Parameters:

output_folder (Path or string, optional) – folder to write report. Defaults to ~/Documents/MR_QA_Reports

class mri_distortion_toolkit.Reports.TG284_Tests

tests defined here https://aapm.onlinelibrary.wiley.com/doi/full/10.1002/mp.14695

test_distortion_less_than_1mm_within_r_100()

distortion test to be run by MRI_QA_Reporter

test_distortion_less_than_2mm_within_r_200()

distortion test to be run by MRI_QA_Reporter

distortion_correction

class mri_distortion_toolkit.DistortionCorrection.DistortionCorrectorBase(ImageDirectory, gradient_harmonics, B0_harmonics=None, dicom_data=None, ImExtension='.dcm', correct_through_plane=True, correct_B0=False, B0_direction='forward', pad=0)

This is the base class which other distortion correction methods can inherit from, and includes the core required behavior.

Parameters:

  • ImageDirectory (pathlib.Path or str) – Directory containing images to correct. Note that all images matching ImExtension will be read in; therefore only one series should be inside ImageDirectory (we dont check this!)

  • gradient_harmonics (list) – a list of three harmonics, which should either be csv files exported from an instance of SphericalHarmonicFit, or SphericalHarmonicFit.harmonics

  • B0_harmonics (SphericalHarmonicFit.harmonics or csv, optional) – Instance of SphericalHarmonicFit.harmonics describing B0, or else an equivalent csv file

  • ImExtension (str, optional) – extension of files to read in ImageDirectory, e.g. ‘dcm’, or ‘IMA’

  • correct_through_plane (bool. optional) – if True, through plane (3D) correction is carried out, which is roughly twice as slow as 2D only

  • correct_B0 (bool, optional) – if True, and if B0_harmonics is supplied, will also correct for B0 effects. Probably only works for standard (non EPI) sequences.

  • B0_direction (str, optional) – ‘forward’ or ‘back’. Determines whether the B0 harmonics should be added (forward) or subtracted (backward) from the gradient harmonics

  • pad (int, optional) – pixels to 0 pad image volume. This can be useful in the k-space domain

correct_all_images()

This loops through the image array, and calls _correct_image on each slice. If correct_through_plane is True, this process occurs twice, with the ‘slice direction’ being changed in the second loop. :return:

save_all_images(save_loc=None, DSV_radius=None, grid=True)

save corrected data as png

Parameters:

save_loc (string or path) – path to save data at.

save_all_images_as_dicom(save_loc=None)

save corrected data as dicom

Parameters:

save_loc (string or path) – path to save data at.

class mri_distortion_toolkit.DistortionCorrection.ImageDomainDistortionCorrector(ImageDirectory, gradient_harmonics, B0_harmonics=None, dicom_data=None, ImExtension='.dcm', correct_through_plane=True, correct_B0=False, B0_direction='forward', pad=0)

This performs image correction in the image domain by constructing a spline between the linear and distorted indices.

class mri_distortion_toolkit.DistortionCorrection.KspaceDistortionCorrector(pad=10, **kwds)

This performs correction using a least squares based approach in the k-space domain, and is based on this work

utilities

General purpose reusable functions

class mri_distortion_toolkit.utilities.bcolors

This is just here to enable me to print pretty colors to the linux terminal

mri_distortion_toolkit.utilities.build_dicom_affine(Dicomfiles)

Build the dicom affine matrix. The reference I used to create the affine is NIbabel (we are using the r/c defined under ‘DICOM affines again’)

Parameters:

Dicomfiles (list) – list of pydicom.dataset.FileDataset instances (created using pydicom.read_file)

Returns:

CoordinateMatrix, a 4x4 matrix which transforms the image based coordinates (i,j,k) to world coordinates (x,y,z)

mri_distortion_toolkit.utilities.combine_harmonics(harmonics_1, harmonics_2, operation='add', n_order_return=None)

This is a general purpose function to combine two harmonic series. If the series are of the same order, this is the default order of the returned series. If they are different, by default the returned series will match the higher order. In either case, you can specify the return order to be lower than the default.

Parameters:

  • harmonics_1 (array-like) – first series to combine

  • harmonics_2 (array-like) – second series to combine

  • operation (str, optional) – ‘add’ or ‘subtract’

  • n_order_return (int, optional) – order of returned series

Returns:

combined_harmonics

mri_distortion_toolkit.utilities.compare_recon_report_with_ground_truth_report(ground_truth_report, recon_report, tolerance=1)

This is to compare the distortion data contained within two Reports.MRI_QA_recon_reporter objects. The original purpose was to compare a recon_report generated from direct data, and one generated via harmonics, although in principle it can take any two reports.

Parameters:

  • ground_truth_report – Reports.MRI_QA_recon_reporter instance

  • recon_report – Reports.MRI_QA_recon_reporter instance

  • tolerance – tolerance in mm

mri_distortion_toolkit.utilities.convert_cartesian_to_spherical(InputCoords)

Converts cartesian coordinates [x,y,z] to spherical coordinates [r, azimuth, elevation]. If r, azimuth, or elevation already exist as column names they are overwritten. The convention used here matches the convention matlab uses, except that we add pi to the azimuth and pi/2 to the elevation such that

  • 0 <= azimuth <= 2*pi

  • 0 <= elevation <= pi

There was once a very good reason for doing this; I can’t remember what it was anymore but I like to stick to tradition.

Parameters:

InputCoords – a data frame containing (at least) columns x, y, z.

Returns InputCoords:

same data frame with additional columns [r, azimuth, elevation]

mri_distortion_toolkit.utilities.convert_spherical_harmonics(harmonics, input_format='full', output_format='none')

function to convert between different harmonic formalisms. At the moment only works between ‘fully normalised’ and ‘unnormalised’ but we could add more later if we need to. example:

g_x_harmonics_no_norm = convert_spherical_harmonics(g_x_harmonics, input_format='full', output_format='none')
Parameters:

  • harmonics – a pandas series or numpy array of size (n_order+1)^2, where the harmonics are ordered as A00, A10, A11, B11, A20, A21, B21, A22, B22…etc.

  • input_format – format of input harmonic normalisation; ‘none’ or ‘full’

  • output_format – format of output harmonic normalisation; ‘none’ or ‘full’

Returns:

converted harmonics as a pandas series

mri_distortion_toolkit.utilities.convert_spherical_to_cartesian(InputCoords)

Converts spherical coordinates [r, azimuth, elevation] to cartesian coordinates [x,y,z]. if the cartesian columns ‘x’, ‘y’, ‘z’ already exist, they will be overwritten

Parameters:

InputCoords – a data frame containing columns r, azimuth, elevation

Returns InputCoords:

same data frame with additional columns [x,y,z]

mri_distortion_toolkit.utilities.dicom_to_numpy(path_to_dicoms, FilesToReadIn=None, file_extension='dcm', return_XYZ=False, zero_pad=0, enforce_increasing_coords=False)

This function does two things:

  1. Loads all the dicom files in path_to_dicoms in a numpy array. data is ordered according to [row, column, slice]

  2. Builds the matrix that transforms the [r,c,s] indices into [x,y,z] coordinates

Parameters:

  • path_to_dicoms (pathlib.Path or string) – folder where all dicoms images are stored

  • FilesToReadIn (list, optional) – By default, this function will read in all files of file type file_extension. If you only want to read in a subset of files, you can specify them with FilesToReadIn

  • file_extension (string) – extension of dicom files, e.g. ‘dcm’ or ‘IMA’

  • return_XYZ (bool, optional) – if True, will build and return three matrics containing the X, Y, and Z coords at each i,j,k pixel.

  • zero_pad (int, optional) – this many zeros will be placed around the returned volume and the coordinates updated accordingly. main use case was for distortion correction. Note that all returned objects will be consistently affected by this parameter

Returns:

ImageArray: a numpy array of voxels

Returns:

dicom_affine: a matrix that transforms the [r,c,s] indices into [x,y,z] coordinates

Returns:

(X, Y, Z): Coordinate arrays the same size as ImageArray. Optional

mri_distortion_toolkit.utilities.enumerate_subfolders(data_loc)

A simple function that prints all the subfolders in data_loc as a dict, e.g:

{'1': 'folder1',
 '2': 'folder2'} etc

I find this is a useful to create a dict object that I copy to the start of my analysis scripts

mri_distortion_toolkit.utilities.generate_harmonic_names(n_order)

generate the names of each harmonic. Used to label the columns in the harmonics dataframe

Parameters:

n_order (int) – order of harmonics

Returns:

coeef_names: a list of names

mri_distortion_toolkit.utilities.generate_legendre_basis(InputCoords, n_order)

Produce spherical base at all coordinate values. This utilises the ‘fully normalised’ convention described here <https://au.mathworks.com/help/matlab/ref/legendre.html>`_ Note that we scale all r values by 1000, which is to help ensure numerical stability. This is based on the assumption that this input data is in mm. The code will throw a warning if this does not appear to the case.

Parameters:

  • InputCoords (Dataframe) – a data frame which can contain whatever you want, but MUST contain columns called elevation, azimuth, and r

  • n_order (int) – order of harmonic fit

Returns legendre_basis:

a pandas data frame of size [Ncoords, (n_order +1)**2)]

mri_distortion_toolkit.utilities.get_dicom_data(dicom_data)

figures out whether dicom data is a dict or a path to a json file if the latter, reads into a dict and returns

mri_distortion_toolkit.utilities.get_harmonics(Gx_Harmonics, Gy_Harmonics, Gz_Harmonics, B0_Harmonics=None)

return the gradient spherical harmonics as pandas series. this function is simply a clean way to handle the different ways users can specify harmonics to other componets of this code

Parameters:

Gx_Harmonics – either a pandas series or a path to a csv. If it is already a series we do nothing, but handle this option for cleaner code elsewhere. Gy_Harmonics and Gz_Harmonics are the same.

Returns:

Gx_Harmonics, Gy_Harmonics, Gz_Harmonics as Pandase series

mri_distortion_toolkit.utilities.plot_MarkerVolume_overlay(MarkerVolumeList, legend=None)

Plot overlaid 3D scatter plots of the marker positions in each MarkerVolume

Parameters:

  • MarkerVolumeList – list of MarkerVolume instances

  • legend – legend to display on plot

Returns:

None

mri_distortion_toolkit.utilities.plot_compressed_MarkerVolumes(MarkerVolumeList, p_max=20, p_min=-20, title=None, legend=None, projection_direction='z')

plot overlay of compressed marker volumes (z coord is ignored)

Parameters:

  • MarkerVolumeList (list) – list of MarkerVolume instances

  • p_max (float, optional) – max data to include in projection_direction

  • p_min (float, optional) – min z data to include in projection_direction

  • title (str, optional) – title of plot

  • legend (list, optional) – legend of plot

Returns:

None

mri_distortion_toolkit.utilities.plot_distortion_xyz_hist(MatchedMarkerVolume)

plot overlaid x, y, z distortion for an input instance of [MatchedMarkerVolume](https://acrf-image-x-institute.github.io/MRI_DistortionQA/code_docs.html#MRI_DistortionQA.MarkerAnalysis.MatchedMarkerVolumes)

Parameters:

MatchedMarkerVolume – an instance of MatchedMarkerVolumes

mri_distortion_toolkit.utilities.plot_harmonics_over_sphere(harmonics, radius, title=None)

plot the reconstructed harmonics over the surface of a sphere I’ve used plotly instead of matplotlib for this as I wasnt happy with the matplotlib result.

Parameters:

  • harmonics – either an instance of SphericalHarmonicFit or the equivalent data frame loaded from csv

  • radius (float) – radius of sphere

Returns:

None

mri_distortion_toolkit.utilities.plot_matched_volume_hist(VolumeList, legend=None)

create a histogram of total distortion for all input volumes

Parameters:
Returns:

None

mri_distortion_toolkit.utilities.printProgressBar(iteration, total, prefix='', suffix='', decimals=1, length=100, fill='█')

Call in a loop to create terminal progress bar. credit here

Parameters:

  • iteration – current iteration (Int)

  • total – total iterations (Int)

  • prefix – prefix string (Str)

  • suffix – suffix string (Str)

  • decimals – positive number of decimals in percent complete (Int)

  • length – character length of bar (Int)

  • fill – bar fill character (Str)

Returns:

progress_string: the string to print

mri_distortion_toolkit.utilities.print_dict(dict)

just provides a nicer display of the input dict than the default print function, got sick of writing the same code over and over…

mri_distortion_toolkit.utilities.read_opera_table_file(data_loc)

reads a table file from opera Tosca and returns the data as a data frame that can be read into spherical harmonic fit. Expected column order is X [MM], Y [MM], Z [MM], BZ [T], with 6 header lines. No error handling is included for incorrect input

Parameters:

data_loc (str or pathlib.Path) – location of table file

Returns:

Bz_data, a pandas dataframe

mri_distortion_toolkit.utilities.reconstruct_Bz(harmonics, coords, quantity='uT', r_outer=None)

Reconstruct Bz at an arbitrary set of coordinates. n_order of reconstruction is derived from the shape of harmonics

Parameters:

  • harmonics – a pandas series of size (n_order+1)^2, such as that defined within FieldAnalysis.SphericalHarmonicFit

  • coords – a pandas Dataframe which can contain whatever you want, but MUST contain columns called: elevation, azimuth, and r

  • quantity –

    • ‘T’: return raw field in T

    • ’uT’: return field in uT, with A00 set to 0.

  • r_outer (float or None) – optional parameter; if you try to reconstruct points with r > r_outer a warning is thrown

Returns Bz_recon:

Pandas series with the value of Bz at each element in coords

mri_distortion_toolkit.utilities.sort_dicom_slices(dicom_datasets, filter_non_image_data=True)

sort slices by instance number

Parameters:

  • dicom_datasets – list of pydicom files

  • filter_non_image_data (boolean, optional) – if True, any dicom file not containing a pixel_data attribute is removed

Returns:

sorted list of pydicom files

calculate_harmonics

mri_distortion_toolkit.calculate_harmonics.calculate_harmonics(MagneticFields, n_order=5, scale=None)

This function is essentially a wrapper of convenience. Given a set of B fields and dicom_data, it will calculate and return the Gx, Gy, and Gz harmonics.

Parameters:

  • MagneticFields (pandas dataframe) – a pandas dataframe with columns [‘x’, ‘y’, ‘z’, ‘B_Gx’, ‘B_Gy’, ‘B_Gz’, ‘B0’]. B0 is optional. x, y, z are in mm and the fields are in any unit - the harmonics will reflect the units

  • n_order (int or list, optional) – order of harmonic fit

  • scale – Can pass a list of four values, the four returned harmonic lists will be scalealised accordingly. e.g. scale = [2,2,2,2] will multiply all harmonics by 2. This is useful to normalise the fields to some value

Returns:

G_x_Harmonics, G_y_Harmonics, G_z_Harmonics, B0_Harmonics