From ade21471184d6e031536ead9707614e78c234705 Mon Sep 17 00:00:00 2001 From: Alexandre D'Astous Date: Wed, 7 Oct 2020 14:22:12 -0400 Subject: [PATCH 01/51] Initial statis shim set up --- shimmingtoolbox/shim/__init__.py | 0 shimmingtoolbox/shim/realtime_zshim.py | 82 ++++++++++++++++++++++++++ 2 files changed, 82 insertions(+) create mode 100644 shimmingtoolbox/shim/__init__.py create mode 100644 shimmingtoolbox/shim/realtime_zshim.py diff --git a/shimmingtoolbox/shim/__init__.py b/shimmingtoolbox/shim/__init__.py new file mode 100644 index 000000000..e69de29bb diff --git a/shimmingtoolbox/shim/realtime_zshim.py b/shimmingtoolbox/shim/realtime_zshim.py new file mode 100644 index 000000000..b87898d23 --- /dev/null +++ b/shimmingtoolbox/shim/realtime_zshim.py @@ -0,0 +1,82 @@ +#!/usr/bin/python3 +# -*- coding: utf-8 -*- + +import numpy as np +import os +import nibabel +from matplotlib.figure import Figure + +from shimmingtoolbox.optimizer.sequential import sequential_zslice +from shimmingtoolbox.coils.siemens_basis import siemens_basis +from shimmingtoolbox.masking.shapes import shapes +from shimmingtoolbox import __dir_testing__ +from shimmingtoolbox import __dir_shimmingtoolbox__ + + +def realtime_zshim(): + fname = os.path.join(__dir_shimmingtoolbox__, __dir_testing__, 'nifti', 'sub-example', 'fmap', 'sub-example_fieldmap.nii.gz') + nii = nibabel.load(fname) + fieldmaps = nii.get_fdata() + affine = nii.affine + + nx, ny, nz, nt = fieldmaps.shape + + # Set up coils + coord_vox = np.meshgrid(np.array(range(nx)), np.array(range(ny)), np.array(range(nz)), indexing='ij') + coord_phys = [np.zeros_like(coord_vox[0]), np.zeros_like(coord_vox[1]), np.zeros_like(coord_vox[2])] + for ix in range(nx): + for iy in range(ny): + for iz in range(nz): + coord_phys_list = np.dot([coord_vox[i][ix, iy, iz] for i in range(3)], affine[0:3, 0:3]) + affine[0:3, + 3] + for i in range(3): + coord_phys[i][ix, iy, iz] = coord_phys_list[i] + + # coord_phys was checked and has the correct scanner coordinates + # TODO: Better code ^ and add as a function + + basis = siemens_basis(coord_phys[0], coord_phys[1], coord_phys[2]) + + # Set up mask + full_mask = shapes(fieldmaps[:, :, :, 0], 'cube', center_dim1=round(nx/2)-5, len_dim1=40, len_dim2=40, len_dim3=nz) + + currents = np.zeros([8, nt]) + shimmed = np.zeros_like(fieldmaps) + masked_fieldmaps = np.zeros_like(fieldmaps) + masked_shimmed = np.zeros_like(shimmed) + for i_t in range(nt): + currents[:, i_t] = sequential_zslice(fieldmaps[:, :, :, i_t], basis, full_mask, z_slices=np.array(range(nz))) + shimmed[:, :, :, i_t] = fieldmaps[:, :, :, i_t] + np.sum(currents[:, i_t] * basis, axis=3, keepdims=False) + masked_fieldmaps[:, :, :, i_t] = full_mask * fieldmaps[:, :, :, i_t] + masked_shimmed[:, :, :, i_t] = full_mask * shimmed[:, :, :, i_t] + + i_t = 0 + # Plot results + fig = Figure(figsize=(10, 10)) + # FigureCanvas(fig) + ax = fig.add_subplot(2, 2, 1) + im = ax.imshow(masked_fieldmaps[:-1, :-1, 0, i_t]) + fig.colorbar(im) + ax.set_title("Masked unshimmed") + ax = fig.add_subplot(2, 2, 2) + im = ax.imshow(masked_shimmed[:-1, :-1, 0, i_t]) + fig.colorbar(im) + ax.set_title("Masked shimmed") + ax = fig.add_subplot(2, 2, 3) + im = ax.imshow(fieldmaps[:-1, :-1, 0, i_t]) + fig.colorbar(im) + ax.set_title("Unshimmed") + ax = fig.add_subplot(2, 2, 4) + im = ax.imshow(shimmed[:-1, :-1, 0, i_t]) + fig.colorbar(im) + ax.set_title("Shimmed") + + print(f"\nThe associated current coefficients are : {currents[:, i_t]}") + + fname_figure = os.path.join(__dir_shimmingtoolbox__, 'realtime_zshim_plot.png') + fig.savefig(fname_figure) + return fname_figure + + +if __name__ == '__main__': + realtime_zshim() From 01df7f11dbd4e11583f91996b14f7cf595b16a4c Mon Sep 17 00:00:00 2001 From: Julien Cohen-Adad Date: Wed, 7 Oct 2020 16:00:09 -0400 Subject: [PATCH 02/51] Moved realtime_zshim to CLI --- shimmingtoolbox/{shim => cli}/realtime_zshim.py | 0 1 file changed, 0 insertions(+), 0 deletions(-) rename shimmingtoolbox/{shim => cli}/realtime_zshim.py (100%) diff --git a/shimmingtoolbox/shim/realtime_zshim.py b/shimmingtoolbox/cli/realtime_zshim.py similarity index 100% rename from shimmingtoolbox/shim/realtime_zshim.py rename to shimmingtoolbox/cli/realtime_zshim.py From 6c05cad531e421a5c92dd80be156bc4901016d7b Mon Sep 17 00:00:00 2001 From: Julien Cohen-Adad Date: Wed, 7 Oct 2020 16:36:18 -0400 Subject: [PATCH 03/51] Added test for realtime_zshim --- shimmingtoolbox/cli/realtime_zshim.py | 44 ++++++++++++--------------- test/cli/test_cli_realtime_zshim.py | 32 +++++++++++++++++++ 2 files changed, 52 insertions(+), 24 deletions(-) create mode 100644 test/cli/test_cli_realtime_zshim.py diff --git a/shimmingtoolbox/cli/realtime_zshim.py b/shimmingtoolbox/cli/realtime_zshim.py index b87898d23..ac35f4170 100644 --- a/shimmingtoolbox/cli/realtime_zshim.py +++ b/shimmingtoolbox/cli/realtime_zshim.py @@ -1,9 +1,10 @@ #!/usr/bin/python3 # -*- coding: utf-8 -*- +import click import numpy as np import os -import nibabel +import nibabel as nib from matplotlib.figure import Figure from shimmingtoolbox.optimizer.sequential import sequential_zslice @@ -13,32 +14,27 @@ from shimmingtoolbox import __dir_shimmingtoolbox__ -def realtime_zshim(): - fname = os.path.join(__dir_shimmingtoolbox__, __dir_testing__, 'nifti', 'sub-example', 'fmap', 'sub-example_fieldmap.nii.gz') - nii = nibabel.load(fname) - fieldmaps = nii.get_fdata() - affine = nii.affine +CONTEXT_SETTINGS = dict(help_option_names=['-h', '--help']) - nx, ny, nz, nt = fieldmaps.shape +@click.command( + context_settings=CONTEXT_SETTINGS, + help=f"Perform realtime z-shimming." +) +@click.option("-coil", help="Coil basis to use for shimming. Enter multiple files if you wish to use more than one set " + "of shim coils (eg: Siemens gradient/shim coils and external custom coils).") +@click.option("-fmap", help="B0 fieldmap. For realtime shimming, this should be a 4d file (4th dimension being time") +@click.option("-mask", help="3D nifti file with voxels between 0 and 1 used to weight the spatial region to shim.") +@click.option("-verbose", is_flag=True, help="Be more verbose.") +def realtime_zshim(coil, fmap, mask, verbose): - # Set up coils - coord_vox = np.meshgrid(np.array(range(nx)), np.array(range(ny)), np.array(range(nz)), indexing='ij') - coord_phys = [np.zeros_like(coord_vox[0]), np.zeros_like(coord_vox[1]), np.zeros_like(coord_vox[2])] - for ix in range(nx): - for iy in range(ny): - for iz in range(nz): - coord_phys_list = np.dot([coord_vox[i][ix, iy, iz] for i in range(3)], affine[0:3, 0:3]) + affine[0:3, - 3] - for i in range(3): - coord_phys[i][ix, iy, iz] = coord_phys_list[i] + nii_coil = nib.load(coil) + nii_fmap = nib.load(fmap) + # TODO: check good practice below + if mask is not None: + nii_mask = nib.load(mask) + else: + nii_mask = None - # coord_phys was checked and has the correct scanner coordinates - # TODO: Better code ^ and add as a function - - basis = siemens_basis(coord_phys[0], coord_phys[1], coord_phys[2]) - - # Set up mask - full_mask = shapes(fieldmaps[:, :, :, 0], 'cube', center_dim1=round(nx/2)-5, len_dim1=40, len_dim2=40, len_dim3=nz) currents = np.zeros([8, nt]) shimmed = np.zeros_like(fieldmaps) diff --git a/test/cli/test_cli_realtime_zshim.py b/test/cli/test_cli_realtime_zshim.py new file mode 100644 index 000000000..09d91fe4f --- /dev/null +++ b/test/cli/test_cli_realtime_zshim.py @@ -0,0 +1,32 @@ +#!usr/bin/env python3 +# coding: utf-8 + +import os +import pathlib +import tempfile + +from click.testing import CliRunner +from shimmingtoolbox.cli.realtime_zshim import realtime_zshim +from shimmingtoolbox import __dir_testing__ + + +def test_cli_realtime_zshim(): + with tempfile.TemporaryDirectory(prefix='st_' + pathlib.Path(__file__).stem) as tmp: + runner = CliRunner() + + fname_fieldmap = os.path.join(__dir_testing__, 'nifti', 'sub-example', 'fmap', 'sub-example_fieldmap.nii.gz') + + # Set up mask + # TODO + full_mask = shapes(fieldmaps[:, :, :, 0], 'cube', center_dim1=round(nx / 2) - 5, len_dim1=40, len_dim2=40, + len_dim3=nz) + + path_dicoms = os.path.join(__dir_testing__, 'dicom_unsorted') + path_nifti = os.path.join(tmp, 'nifti') + subject_id = 'sub-test' + + result = runner.invoke(realtime_zshim, ['-i', path_dicoms, '-o', path_nifti, '-s', subject_id]) + + assert result.exit_code == 0 + # Check that dicom_to_nifti was invoked, not if files were actually created (API test already looks for that) + assert os.path.isdir(path_nifti) \ No newline at end of file From 4b8294b7b9b00e43e200d045d0f35698a27c1e61 Mon Sep 17 00:00:00 2001 From: Alexandre D'Astous Date: Thu, 8 Oct 2020 19:46:09 -0400 Subject: [PATCH 04/51] Completed test set up --- shimmingtoolbox/cli/realtime_zshim.py | 49 ++++++++++++++------------- test/cli/test_cli_realtime_zshim.py | 34 +++++++++++++------ 2 files changed, 50 insertions(+), 33 deletions(-) diff --git a/shimmingtoolbox/cli/realtime_zshim.py b/shimmingtoolbox/cli/realtime_zshim.py index ac35f4170..7c9387f42 100644 --- a/shimmingtoolbox/cli/realtime_zshim.py +++ b/shimmingtoolbox/cli/realtime_zshim.py @@ -8,43 +8,50 @@ from matplotlib.figure import Figure from shimmingtoolbox.optimizer.sequential import sequential_zslice -from shimmingtoolbox.coils.siemens_basis import siemens_basis -from shimmingtoolbox.masking.shapes import shapes -from shimmingtoolbox import __dir_testing__ from shimmingtoolbox import __dir_shimmingtoolbox__ - CONTEXT_SETTINGS = dict(help_option_names=['-h', '--help']) + @click.command( context_settings=CONTEXT_SETTINGS, help=f"Perform realtime z-shimming." ) -@click.option("-coil", help="Coil basis to use for shimming. Enter multiple files if you wish to use more than one set " - "of shim coils (eg: Siemens gradient/shim coils and external custom coils).") -@click.option("-fmap", help="B0 fieldmap. For realtime shimming, this should be a 4d file (4th dimension being time") -@click.option("-mask", help="3D nifti file with voxels between 0 and 1 used to weight the spatial region to shim.") +@click.option("-coil", required=True, type=click.Path(), help="Coil basis to use for shimming. Enter multiple files if " + "you wish to use more than one set of shim coils (eg: " + "Siemens gradient/shim coils and external custom coils).") +@click.option("-fmap", required=True, type=click.Path(), help="B0 fieldmap. For realtime shimming, this should be a 4d " + "file (4th dimension being time") +@click.option("-mask", type=click.Path(), help="3D nifti file with voxels between 0 and 1 used to weight the spatial " + "region to shim.") @click.option("-verbose", is_flag=True, help="Be more verbose.") def realtime_zshim(coil, fmap, mask, verbose): - nii_coil = nib.load(coil) + coil = nib.load(coil).get_fdata() nii_fmap = nib.load(fmap) + fieldmap = nii_fmap.get_fdata() + + # TODO: Error handling might move to API + if fieldmap.ndim != 4: + raise RuntimeError('fmap must be 4d (x, y, z, t)') + # TODO: check good practice below if mask is not None: - nii_mask = nib.load(mask) + mask = nib.load(mask).get_fdata() else: - nii_mask = None + mask = np.ones_like(fieldmap) + nx, ny, nz, nt = fieldmap.shape currents = np.zeros([8, nt]) - shimmed = np.zeros_like(fieldmaps) - masked_fieldmaps = np.zeros_like(fieldmaps) + shimmed = np.zeros_like(fieldmap) + masked_fieldmaps = np.zeros_like(fieldmap) masked_shimmed = np.zeros_like(shimmed) for i_t in range(nt): - currents[:, i_t] = sequential_zslice(fieldmaps[:, :, :, i_t], basis, full_mask, z_slices=np.array(range(nz))) - shimmed[:, :, :, i_t] = fieldmaps[:, :, :, i_t] + np.sum(currents[:, i_t] * basis, axis=3, keepdims=False) - masked_fieldmaps[:, :, :, i_t] = full_mask * fieldmaps[:, :, :, i_t] - masked_shimmed[:, :, :, i_t] = full_mask * shimmed[:, :, :, i_t] + currents[:, i_t] = sequential_zslice(fieldmap[:, :, :, i_t], coil, mask, z_slices=np.array(range(nz))) + shimmed[:, :, :, i_t] = fieldmap[:, :, :, i_t] + np.sum(currents[:, i_t] * coil, axis=3, keepdims=False) + masked_fieldmaps[:, :, :, i_t] = mask * fieldmap[:, :, :, i_t] + masked_shimmed[:, :, :, i_t] = mask * shimmed[:, :, :, i_t] i_t = 0 # Plot results @@ -59,7 +66,7 @@ def realtime_zshim(coil, fmap, mask, verbose): fig.colorbar(im) ax.set_title("Masked shimmed") ax = fig.add_subplot(2, 2, 3) - im = ax.imshow(fieldmaps[:-1, :-1, 0, i_t]) + im = ax.imshow(fieldmap[:-1, :-1, 0, i_t]) fig.colorbar(im) ax.set_title("Unshimmed") ax = fig.add_subplot(2, 2, 4) @@ -67,12 +74,8 @@ def realtime_zshim(coil, fmap, mask, verbose): fig.colorbar(im) ax.set_title("Shimmed") - print(f"\nThe associated current coefficients are : {currents[:, i_t]}") + click.echo(f"\nThe associated current coefficients are : {currents[:, i_t]}") fname_figure = os.path.join(__dir_shimmingtoolbox__, 'realtime_zshim_plot.png') fig.savefig(fname_figure) return fname_figure - - -if __name__ == '__main__': - realtime_zshim() diff --git a/test/cli/test_cli_realtime_zshim.py b/test/cli/test_cli_realtime_zshim.py index 09d91fe4f..e447ac3b6 100644 --- a/test/cli/test_cli_realtime_zshim.py +++ b/test/cli/test_cli_realtime_zshim.py @@ -4,9 +4,14 @@ import os import pathlib import tempfile +import nibabel as nib +import numpy as np from click.testing import CliRunner from shimmingtoolbox.cli.realtime_zshim import realtime_zshim +from shimmingtoolbox.masking.shapes import shapes +from shimmingtoolbox.coils.coordinates import generate_meshgrid +from shimmingtoolbox.coils.siemens_basis import siemens_basis from shimmingtoolbox import __dir_testing__ @@ -14,19 +19,28 @@ def test_cli_realtime_zshim(): with tempfile.TemporaryDirectory(prefix='st_' + pathlib.Path(__file__).stem) as tmp: runner = CliRunner() - fname_fieldmap = os.path.join(__dir_testing__, 'nifti', 'sub-example', 'fmap', 'sub-example_fieldmap.nii.gz') + fname_fieldmap = os.path.join(__dir_testing__, 'sub-fieldmap', 'fmap', 'sub-example_fieldmap.nii.gz') + nii_fmap = nib.load(fname_fieldmap) + fmap = nii_fmap.get_fdata() + affine = nii_fmap.affine # Set up mask - # TODO - full_mask = shapes(fieldmaps[:, :, :, 0], 'cube', center_dim1=round(nx / 2) - 5, len_dim1=40, len_dim2=40, - len_dim3=nz) + nx, ny, nz, nt = fmap.shape + mask = shapes(fmap[:, :, :, 0], 'cube', center_dim1=round(nx / 2) - 5, len_dim1=40, len_dim2=40, + len_dim3=nz) - path_dicoms = os.path.join(__dir_testing__, 'dicom_unsorted') - path_nifti = os.path.join(tmp, 'nifti') - subject_id = 'sub-test' + nii_mask = nib.Nifti1Image(mask.astype(int), affine) + fname_mask = os.path.join(tmp, 'mask.nii.gz') + nib.save(nii_mask, fname_mask) - result = runner.invoke(realtime_zshim, ['-i', path_dicoms, '-o', path_nifti, '-s', subject_id]) + # Set up coils + coord_phys = generate_meshgrid(fmap.shape[0:3], affine) + coil_profile = siemens_basis(coord_phys[0], coord_phys[1], coord_phys[2]) + + nii_coil = nib.Nifti1Image(coil_profile, affine) + fname_coil = os.path.join(tmp, 'coil_profile.nii.gz') + nib.save(nii_coil, fname_coil) + + result = runner.invoke(realtime_zshim, ['-fmap', fname_fieldmap, '-coil', fname_coil, '-mask', fname_mask]) assert result.exit_code == 0 - # Check that dicom_to_nifti was invoked, not if files were actually created (API test already looks for that) - assert os.path.isdir(path_nifti) \ No newline at end of file From 718115e91b8e5eced630a8033e297c8d41e4cdbf Mon Sep 17 00:00:00 2001 From: Alexandre D'Astous Date: Thu, 8 Oct 2020 19:46:40 -0400 Subject: [PATCH 05/51] Added to set up as a CLI --- setup.py | 1 + 1 file changed, 1 insertion(+) diff --git a/setup.py b/setup.py index e6e0fd2bb..bc6c8319a 100644 --- a/setup.py +++ b/setup.py @@ -22,6 +22,7 @@ "st_referencemaps=shimmingtoolbox.cli.referencemaps:main", "st_b0maps=shimmingtoolbox.cli.b0map:main", "st_download_data=shimmingtoolbox.cli.download_data:download_data", + "st_realtime_zshim=shimmingtoolbox.cli.realtime_zshim:realtime_zshim", ] }, packages=find_packages(exclude=["contrib", "docs", "tests"]), From 4ab1962f7dc470009a6b0c6d2a789a9497464cf0 Mon Sep 17 00:00:00 2001 From: Alexandre D'Astous Date: Tue, 20 Oct 2020 18:21:27 -0400 Subject: [PATCH 06/51] Reduce shimming volume --- test/cli/test_cli_realtime_zshim.py | 6 ++++-- 1 file changed, 4 insertions(+), 2 deletions(-) diff --git a/test/cli/test_cli_realtime_zshim.py b/test/cli/test_cli_realtime_zshim.py index e447ac3b6..5f566ae94 100644 --- a/test/cli/test_cli_realtime_zshim.py +++ b/test/cli/test_cli_realtime_zshim.py @@ -26,8 +26,10 @@ def test_cli_realtime_zshim(): # Set up mask nx, ny, nz, nt = fmap.shape - mask = shapes(fmap[:, :, :, 0], 'cube', center_dim1=round(nx / 2) - 5, len_dim1=40, len_dim2=40, - len_dim3=nz) + mask = shapes(fmap[:, :, :, 0], 'cube', + center_dim1=int(fmap.shape[0] / 2 - 8), + center_dim2=int(fmap.shape[1] / 2 - 20), + len_dim1=15, len_dim2=5, len_dim3=nz) nii_mask = nib.Nifti1Image(mask.astype(int), affine) fname_mask = os.path.join(tmp, 'mask.nii.gz') From 40e228fd277bcc6c7a22a19d2af78f6b213ebf9c Mon Sep 17 00:00:00 2001 From: Alexandre D'Astous Date: Tue, 20 Oct 2020 18:21:49 -0400 Subject: [PATCH 07/51] Add plot of currents through time --- shimmingtoolbox/cli/realtime_zshim.py | 10 ++++++++++ 1 file changed, 10 insertions(+) diff --git a/shimmingtoolbox/cli/realtime_zshim.py b/shimmingtoolbox/cli/realtime_zshim.py index 7c9387f42..b85948c09 100644 --- a/shimmingtoolbox/cli/realtime_zshim.py +++ b/shimmingtoolbox/cli/realtime_zshim.py @@ -78,4 +78,14 @@ def realtime_zshim(coil, fmap, mask, verbose): fname_figure = os.path.join(__dir_shimmingtoolbox__, 'realtime_zshim_plot.png') fig.savefig(fname_figure) + + fig = Figure(figsize=(10, 10)) + for i_t in range(nt): + ax = fig.add_subplot(nt, 1, i_t + 1) + ax.plot(np.arange(8), currents[:, i_t]) + ax.set_title(f"Time {i_t}") + + fname_figure = os.path.join(__dir_shimmingtoolbox__, 'realtime_zshim_currents.png') + fig.savefig(fname_figure) + return fname_figure From 1058ca20a53cf4d224425e17445e65e4164715e8 Mon Sep 17 00:00:00 2001 From: Alexandre D'Astous Date: Tue, 20 Oct 2020 19:13:28 -0400 Subject: [PATCH 08/51] Adds ability to see output error messages from runner.invoke --- test/cli/test_cli_realtime_zshim.py | 3 ++- 1 file changed, 2 insertions(+), 1 deletion(-) diff --git a/test/cli/test_cli_realtime_zshim.py b/test/cli/test_cli_realtime_zshim.py index 5f566ae94..c4fc42cfc 100644 --- a/test/cli/test_cli_realtime_zshim.py +++ b/test/cli/test_cli_realtime_zshim.py @@ -43,6 +43,7 @@ def test_cli_realtime_zshim(): fname_coil = os.path.join(tmp, 'coil_profile.nii.gz') nib.save(nii_coil, fname_coil) - result = runner.invoke(realtime_zshim, ['-fmap', fname_fieldmap, '-coil', fname_coil, '-mask', fname_mask]) + result = runner.invoke(realtime_zshim, ['-fmap', fname_fieldmap, '-coil', fname_coil, '-mask', fname_mask], + catch_exceptions=False) assert result.exit_code == 0 From 404e2e25ec8b3b94f0906210c8e92ac027ba5c47 Mon Sep 17 00:00:00 2001 From: Julien Cohen-Adad Date: Wed, 21 Oct 2020 15:53:13 -0400 Subject: [PATCH 09/51] Changed variable name --- shimmingtoolbox/cli/realtime_zshim.py | 25 +++++++++++++------------ 1 file changed, 13 insertions(+), 12 deletions(-) diff --git a/shimmingtoolbox/cli/realtime_zshim.py b/shimmingtoolbox/cli/realtime_zshim.py index b85948c09..694275d90 100644 --- a/shimmingtoolbox/cli/realtime_zshim.py +++ b/shimmingtoolbox/cli/realtime_zshim.py @@ -17,18 +17,19 @@ context_settings=CONTEXT_SETTINGS, help=f"Perform realtime z-shimming." ) -@click.option("-coil", required=True, type=click.Path(), help="Coil basis to use for shimming. Enter multiple files if " - "you wish to use more than one set of shim coils (eg: " - "Siemens gradient/shim coils and external custom coils).") -@click.option("-fmap", required=True, type=click.Path(), help="B0 fieldmap. For realtime shimming, this should be a 4d " - "file (4th dimension being time") -@click.option("-mask", type=click.Path(), help="3D nifti file with voxels between 0 and 1 used to weight the spatial " - "region to shim.") +@click.option('-coil', 'fname_coil', required=True, type=click.Path(), + help="Coil basis to use for shimming. Enter multiple files if " + "you wish to use more than one set of shim coils (eg: " + "Siemens gradient/shim coils and external custom coils).") +@click.option('-fmap', 'fname_fmap', required=True, type=click.Path(), + help="B0 fieldmap. For realtime shimming, this should be a 4d file (4th dimension being time") +@click.option('-mask', 'fname_mask', type=click.Path(), + help="3D nifti file with voxels between 0 and 1 used to weight the spatial region to shim.") @click.option("-verbose", is_flag=True, help="Be more verbose.") -def realtime_zshim(coil, fmap, mask, verbose): +def realtime_zshim(fname_coil, fname_fmap, fname_mask, verbose): - coil = nib.load(coil).get_fdata() - nii_fmap = nib.load(fmap) + coil = nib.load(fname_coil).get_fdata() + nii_fmap = nib.load(fname_fmap) fieldmap = nii_fmap.get_fdata() # TODO: Error handling might move to API @@ -36,8 +37,8 @@ def realtime_zshim(coil, fmap, mask, verbose): raise RuntimeError('fmap must be 4d (x, y, z, t)') # TODO: check good practice below - if mask is not None: - mask = nib.load(mask).get_fdata() + if fname_mask is not None: + mask = nib.load(fname_mask).get_fdata() else: mask = np.ones_like(fieldmap) From 6855ae4959d90e51f11fbd49c8b199e479d7bd5a Mon Sep 17 00:00:00 2001 From: Julien Cohen-Adad Date: Wed, 21 Oct 2020 16:57:17 -0400 Subject: [PATCH 10/51] Clarified comments, added TODO, minor refactoring - Added bounds to sequential_zslice() because of bug: https://github.com/shimming-toolbox/shimming-toolbox/issues/152 - Commented out Click decorator because not working in Pycharm: https://github.com/shimming-toolbox/shimming-toolbox/issues/156 --- shimmingtoolbox/cli/realtime_zshim.py | 73 +++++++++++++++++++-------- 1 file changed, 52 insertions(+), 21 deletions(-) diff --git a/shimmingtoolbox/cli/realtime_zshim.py b/shimmingtoolbox/cli/realtime_zshim.py index 694275d90..557007f2d 100644 --- a/shimmingtoolbox/cli/realtime_zshim.py +++ b/shimmingtoolbox/cli/realtime_zshim.py @@ -5,6 +5,7 @@ import numpy as np import os import nibabel as nib +# TODO: remove matplotlib import from matplotlib.figure import Figure from shimmingtoolbox.optimizer.sequential import sequential_zslice @@ -13,21 +14,31 @@ CONTEXT_SETTINGS = dict(help_option_names=['-h', '--help']) -@click.command( - context_settings=CONTEXT_SETTINGS, - help=f"Perform realtime z-shimming." -) -@click.option('-coil', 'fname_coil', required=True, type=click.Path(), - help="Coil basis to use for shimming. Enter multiple files if " - "you wish to use more than one set of shim coils (eg: " - "Siemens gradient/shim coils and external custom coils).") -@click.option('-fmap', 'fname_fmap', required=True, type=click.Path(), - help="B0 fieldmap. For realtime shimming, this should be a 4d file (4th dimension being time") -@click.option('-mask', 'fname_mask', type=click.Path(), - help="3D nifti file with voxels between 0 and 1 used to weight the spatial region to shim.") -@click.option("-verbose", is_flag=True, help="Be more verbose.") -def realtime_zshim(fname_coil, fname_fmap, fname_mask, verbose): - +# @click.command( +# context_settings=CONTEXT_SETTINGS, +# help=f"Perform realtime z-shimming." +# ) +# @click.option('-coil', 'fname_coil', required=True, type=click.Path(), +# help="Coil basis to use for shimming. Enter multiple files if " +# "you wish to use more than one set of shim coils (eg: " +# "Siemens gradient/shim coils and external custom coils).") +# @click.option('-fmap', 'fname_fmap', required=True, type=click.Path(), +# help="B0 fieldmap. For realtime shimming, this should be a 4d file (4th dimension being time") +# @click.option('-mask', 'fname_mask', type=click.Path(), +# help="3D nifti file with voxels between 0 and 1 used to weight the spatial region to shim.") +# @click.option("-verbose", is_flag=True, help="Be more verbose.") +def realtime_zshim(fname_coil, fname_fmap, fname_mask, verbose=True): + """ + + Args: + fname_coil: Pointing to coil profile. 4-dimensional: x, y, z, coil. + fname_fmap: + fname_mask: + verbose: + + Returns: + + """ coil = nib.load(fname_coil).get_fdata() nii_fmap = nib.load(fname_fmap) fieldmap = nii_fmap.get_fdata() @@ -35,6 +46,7 @@ def realtime_zshim(fname_coil, fname_fmap, fname_mask, verbose): # TODO: Error handling might move to API if fieldmap.ndim != 4: raise RuntimeError('fmap must be 4d (x, y, z, t)') + nx, ny, nz, nt = fieldmap.shape # TODO: check good practice below if fname_mask is not None: @@ -42,17 +54,19 @@ def realtime_zshim(fname_coil, fname_fmap, fname_mask, verbose): else: mask = np.ones_like(fieldmap) - nx, ny, nz, nt = fieldmap.shape + # Setup coil + n_coils = coil.shape[-1] + currents = np.zeros([n_coils, nt]) - currents = np.zeros([8, nt]) shimmed = np.zeros_like(fieldmap) masked_fieldmaps = np.zeros_like(fieldmap) masked_shimmed = np.zeros_like(shimmed) for i_t in range(nt): - currents[:, i_t] = sequential_zslice(fieldmap[:, :, :, i_t], coil, mask, z_slices=np.array(range(nz))) - shimmed[:, :, :, i_t] = fieldmap[:, :, :, i_t] + np.sum(currents[:, i_t] * coil, axis=3, keepdims=False) - masked_fieldmaps[:, :, :, i_t] = mask * fieldmap[:, :, :, i_t] - masked_shimmed[:, :, :, i_t] = mask * shimmed[:, :, :, i_t] + currents[:, i_t] = sequential_zslice(fieldmap[..., i_t], coil, mask, z_slices=np.array(range(nz)), + bounds=[(-np.inf, np.inf)]*8) + shimmed[..., i_t] = fieldmap[..., i_t] + np.sum(currents[:, i_t] * coil, axis=3, keepdims=False) + masked_fieldmaps[..., i_t] = mask * fieldmap[..., i_t] + masked_shimmed[..., i_t] = mask * shimmed[..., i_t] i_t = 0 # Plot results @@ -81,6 +95,7 @@ def realtime_zshim(fname_coil, fname_fmap, fname_mask, verbose): fig.savefig(fname_figure) fig = Figure(figsize=(10, 10)) + # TODO: display time as X -- add label for Y for i_t in range(nt): ax = fig.add_subplot(nt, 1, i_t + 1) ax.plot(np.arange(8), currents[:, i_t]) @@ -89,4 +104,20 @@ def realtime_zshim(fname_coil, fname_fmap, fname_mask, verbose): fname_figure = os.path.join(__dir_shimmingtoolbox__, 'realtime_zshim_currents.png') fig.savefig(fname_figure) + # TODO: fetch PMU timing + get_acquisition_time() + + + # TODO: + # fit PMU and fieldmap values + # do regression to separate static componant and RIRO component + # output coefficient with proper scaling + return fname_figure + + +if __name__ == '__main__': + fname_coil='/Users/julien/code/shimming-toolbox/shimming-toolbox/test_realtime_zshim/coil_profile.nii.gz' + fname_fmap='/Users/julien/code/shimming-toolbox/shimming-toolbox/test_realtime_zshim/sub-example_fieldmap.nii.gz' + fname_mask='/Users/julien/code/shimming-toolbox/shimming-toolbox/test_realtime_zshim/mask.nii.gz' + realtime_zshim(fname_coil, fname_fmap, fname_mask) From 88ab640b13caed1e542c2aad973861fb42cc346e Mon Sep 17 00:00:00 2001 From: Julien Cohen-Adad Date: Wed, 21 Oct 2020 17:00:49 -0400 Subject: [PATCH 11/51] Added TODO --- shimmingtoolbox/cli/realtime_zshim.py | 3 +-- 1 file changed, 1 insertion(+), 2 deletions(-) diff --git a/shimmingtoolbox/cli/realtime_zshim.py b/shimmingtoolbox/cli/realtime_zshim.py index 557007f2d..8435a1430 100644 --- a/shimmingtoolbox/cli/realtime_zshim.py +++ b/shimmingtoolbox/cli/realtime_zshim.py @@ -105,8 +105,7 @@ def realtime_zshim(fname_coil, fname_fmap, fname_mask, verbose=True): fig.savefig(fname_figure) # TODO: fetch PMU timing - get_acquisition_time() - + # get_pmu_data(fname_pmu, fname_fieldmap) # TODO: # fit PMU and fieldmap values From 541e6c0e8564a557a83f110fd18829ea2591d6e7 Mon Sep 17 00:00:00 2001 From: Alexandre D'Astous Date: Thu, 22 Oct 2020 16:46:44 -0400 Subject: [PATCH 12/51] Change debug paths and plot correct currents --- shimmingtoolbox/cli/realtime_zshim.py | 26 ++++++++++++++------------ 1 file changed, 14 insertions(+), 12 deletions(-) diff --git a/shimmingtoolbox/cli/realtime_zshim.py b/shimmingtoolbox/cli/realtime_zshim.py index 8435a1430..6daa98770 100644 --- a/shimmingtoolbox/cli/realtime_zshim.py +++ b/shimmingtoolbox/cli/realtime_zshim.py @@ -5,8 +5,9 @@ import numpy as np import os import nibabel as nib -# TODO: remove matplotlib import +# TODO: remove matplotlib and dirtesting import from matplotlib.figure import Figure +from shimmingtoolbox import __dir_testing__ from shimmingtoolbox.optimizer.sequential import sequential_zslice from shimmingtoolbox import __dir_shimmingtoolbox__ @@ -63,7 +64,7 @@ def realtime_zshim(fname_coil, fname_fmap, fname_mask, verbose=True): masked_shimmed = np.zeros_like(shimmed) for i_t in range(nt): currents[:, i_t] = sequential_zslice(fieldmap[..., i_t], coil, mask, z_slices=np.array(range(nz)), - bounds=[(-np.inf, np.inf)]*8) + bounds=[(-np.inf, np.inf)]*n_coils) shimmed[..., i_t] = fieldmap[..., i_t] + np.sum(currents[:, i_t] * coil, axis=3, keepdims=False) masked_fieldmaps[..., i_t] = mask * fieldmap[..., i_t] masked_shimmed[..., i_t] = mask * shimmed[..., i_t] @@ -95,11 +96,10 @@ def realtime_zshim(fname_coil, fname_fmap, fname_mask, verbose=True): fig.savefig(fname_figure) fig = Figure(figsize=(10, 10)) - # TODO: display time as X -- add label for Y - for i_t in range(nt): - ax = fig.add_subplot(nt, 1, i_t + 1) - ax.plot(np.arange(8), currents[:, i_t]) - ax.set_title(f"Time {i_t}") + for i_coil in range(n_coils): + ax = fig.add_subplot(n_coils, 1, i_coil + 1) + ax.plot(np.arange(nt), currents[i_coil, :]) + ax.set_title(f"Channel {i_coil}") fname_figure = os.path.join(__dir_shimmingtoolbox__, 'realtime_zshim_currents.png') fig.savefig(fname_figure) @@ -115,8 +115,10 @@ def realtime_zshim(fname_coil, fname_fmap, fname_mask, verbose=True): return fname_figure -if __name__ == '__main__': - fname_coil='/Users/julien/code/shimming-toolbox/shimming-toolbox/test_realtime_zshim/coil_profile.nii.gz' - fname_fmap='/Users/julien/code/shimming-toolbox/shimming-toolbox/test_realtime_zshim/sub-example_fieldmap.nii.gz' - fname_mask='/Users/julien/code/shimming-toolbox/shimming-toolbox/test_realtime_zshim/mask.nii.gz' - realtime_zshim(fname_coil, fname_fmap, fname_mask) +fname_coil = os.path.join(__dir_testing__, 'test_realtime_zshim', 'coil_profile.nii.gz') +fname_fmap = os.path.join(__dir_testing__, 'test_realtime_zshim', 'sub-example_fieldmap.nii.gz') +fname_mask = os.path.join(__dir_testing__, 'test_realtime_zshim', 'mask.nii.gz') +# fname_coil='/Users/julien/code/shimming-toolbox/shimming-toolbox/test_realtime_zshim/coil_profile.nii.gz' +# fname_fmap='/Users/julien/code/shimming-toolbox/shimming-toolbox/test_realtime_zshim/sub-example_fieldmap.nii.gz' +# fname_mask='/Users/julien/code/shimming-toolbox/shimming-toolbox/test_realtime_zshim/mask.nii.gz' +realtime_zshim(fname_coil, fname_fmap, fname_mask) From 59c4bb451d1d7859433468d40080646ca107888e Mon Sep 17 00:00:00 2001 From: Alexandre D'Astous Date: Thu, 22 Oct 2020 17:16:14 -0400 Subject: [PATCH 13/51] Sanity check for z component --- shimmingtoolbox/cli/realtime_zshim.py | 5 +++++ 1 file changed, 5 insertions(+) diff --git a/shimmingtoolbox/cli/realtime_zshim.py b/shimmingtoolbox/cli/realtime_zshim.py index 6daa98770..f8b7cda04 100644 --- a/shimmingtoolbox/cli/realtime_zshim.py +++ b/shimmingtoolbox/cli/realtime_zshim.py @@ -40,7 +40,12 @@ def realtime_zshim(fname_coil, fname_fmap, fname_mask, verbose=True): Returns: """ + # When using only z channnel TODO:Remove + # coil = np.expand_dims(nib.load(fname_coil).get_fdata()[:, :, :, 2], -1) + + # When using all channels TODO: Keep coil = nib.load(fname_coil).get_fdata() + nii_fmap = nib.load(fname_fmap) fieldmap = nii_fmap.get_fdata() From ad3de27ba1e74b6c9beb14b2c8a51ecf20805b0b Mon Sep 17 00:00:00 2001 From: Alexandre D'Astous Date: Thu, 22 Oct 2020 17:49:30 -0400 Subject: [PATCH 14/51] Added TODO for pmu_data, Related to #149 --- shimmingtoolbox/cli/realtime_zshim.py | 11 ++++++++--- 1 file changed, 8 insertions(+), 3 deletions(-) diff --git a/shimmingtoolbox/cli/realtime_zshim.py b/shimmingtoolbox/cli/realtime_zshim.py index f8b7cda04..3f68133d6 100644 --- a/shimmingtoolbox/cli/realtime_zshim.py +++ b/shimmingtoolbox/cli/realtime_zshim.py @@ -28,7 +28,7 @@ # @click.option('-mask', 'fname_mask', type=click.Path(), # help="3D nifti file with voxels between 0 and 1 used to weight the spatial region to shim.") # @click.option("-verbose", is_flag=True, help="Be more verbose.") -def realtime_zshim(fname_coil, fname_fmap, fname_mask, verbose=True): +def realtime_zshim(fname_coil, fname_fmap, fname_mask, fname_resp, verbose=True): """ Args: @@ -110,7 +110,11 @@ def realtime_zshim(fname_coil, fname_fmap, fname_mask, verbose=True): fig.savefig(fname_figure) # TODO: fetch PMU timing - # get_pmu_data(fname_pmu, fname_fieldmap) + # get_pmu_data(fname_pmu, fname_fieldmap) # Currently needs another Json file + # acq_times = get_acquisition_times(fname_fieldmap) # Currently needs another Json file + # pmu = PmuResp(fname_resp) + # acq_pressures = pmu.interp_resp_trace(acq_times) + # return acq_pressures, acq_times # TODO: # fit PMU and fieldmap values @@ -123,7 +127,8 @@ def realtime_zshim(fname_coil, fname_fmap, fname_mask, verbose=True): fname_coil = os.path.join(__dir_testing__, 'test_realtime_zshim', 'coil_profile.nii.gz') fname_fmap = os.path.join(__dir_testing__, 'test_realtime_zshim', 'sub-example_fieldmap.nii.gz') fname_mask = os.path.join(__dir_testing__, 'test_realtime_zshim', 'mask.nii.gz') +fname_resp = os.path.join(__dir_testing__, 'realtime_zshimming_data', 'PMUresp_signal.resp') # fname_coil='/Users/julien/code/shimming-toolbox/shimming-toolbox/test_realtime_zshim/coil_profile.nii.gz' # fname_fmap='/Users/julien/code/shimming-toolbox/shimming-toolbox/test_realtime_zshim/sub-example_fieldmap.nii.gz' # fname_mask='/Users/julien/code/shimming-toolbox/shimming-toolbox/test_realtime_zshim/mask.nii.gz' -realtime_zshim(fname_coil, fname_fmap, fname_mask) +realtime_zshim(fname_coil, fname_fmap, fname_mask, fname_resp) From 6fe32ad45af03701d872799f366ff034e7ccdaf6 Mon Sep 17 00:00:00 2001 From: Alexandre D'Astous Date: Sun, 25 Oct 2020 22:57:28 -0400 Subject: [PATCH 15/51] Added pmu timestamps and pressures --- shimmingtoolbox/cli/realtime_zshim.py | 21 +++++++++++++-------- 1 file changed, 13 insertions(+), 8 deletions(-) diff --git a/shimmingtoolbox/cli/realtime_zshim.py b/shimmingtoolbox/cli/realtime_zshim.py index 3f68133d6..bb9cf72c3 100644 --- a/shimmingtoolbox/cli/realtime_zshim.py +++ b/shimmingtoolbox/cli/realtime_zshim.py @@ -5,11 +5,14 @@ import numpy as np import os import nibabel as nib +import json # TODO: remove matplotlib and dirtesting import from matplotlib.figure import Figure from shimmingtoolbox import __dir_testing__ from shimmingtoolbox.optimizer.sequential import sequential_zslice +from shimmingtoolbox.load_nifti import get_acquisition_times +from shimmingtoolbox.pmu import PmuResp from shimmingtoolbox import __dir_shimmingtoolbox__ CONTEXT_SETTINGS = dict(help_option_names=['-h', '--help']) @@ -28,7 +31,7 @@ # @click.option('-mask', 'fname_mask', type=click.Path(), # help="3D nifti file with voxels between 0 and 1 used to weight the spatial region to shim.") # @click.option("-verbose", is_flag=True, help="Be more verbose.") -def realtime_zshim(fname_coil, fname_fmap, fname_mask, fname_resp, verbose=True): +def realtime_zshim(fname_coil, fname_fmap, fname_mask, fname_resp, fname_json, verbose=True): """ Args: @@ -109,12 +112,13 @@ def realtime_zshim(fname_coil, fname_fmap, fname_mask, fname_resp, verbose=True) fname_figure = os.path.join(__dir_shimmingtoolbox__, 'realtime_zshim_currents.png') fig.savefig(fname_figure) - # TODO: fetch PMU timing - # get_pmu_data(fname_pmu, fname_fieldmap) # Currently needs another Json file - # acq_times = get_acquisition_times(fname_fieldmap) # Currently needs another Json file - # pmu = PmuResp(fname_resp) - # acq_pressures = pmu.interp_resp_trace(acq_times) - # return acq_pressures, acq_times + # Fetch PMU timing + # TODO: Add json to fieldmap instead of asking for another json file + with open(fname_json) as json_file: + json_data = json.load(json_file) + acq_timestamps = get_acquisition_times(nii_fmap, json_data) + pmu = PmuResp(fname_resp) + acq_pressures = pmu.interp_resp_trace(acq_timestamps) # TODO: # fit PMU and fieldmap values @@ -128,7 +132,8 @@ def realtime_zshim(fname_coil, fname_fmap, fname_mask, fname_resp, verbose=True) fname_fmap = os.path.join(__dir_testing__, 'test_realtime_zshim', 'sub-example_fieldmap.nii.gz') fname_mask = os.path.join(__dir_testing__, 'test_realtime_zshim', 'mask.nii.gz') fname_resp = os.path.join(__dir_testing__, 'realtime_zshimming_data', 'PMUresp_signal.resp') +fname_json = os.path.join(__dir_testing__, 'test_realtime_zshim', 'sub-example_magnitude1.json') # fname_coil='/Users/julien/code/shimming-toolbox/shimming-toolbox/test_realtime_zshim/coil_profile.nii.gz' # fname_fmap='/Users/julien/code/shimming-toolbox/shimming-toolbox/test_realtime_zshim/sub-example_fieldmap.nii.gz' # fname_mask='/Users/julien/code/shimming-toolbox/shimming-toolbox/test_realtime_zshim/mask.nii.gz' -realtime_zshim(fname_coil, fname_fmap, fname_mask, fname_resp) +realtime_zshim(fname_coil, fname_fmap, fname_mask, fname_resp, fname_json) From 9dfaed08db3f7a214c3500604760d4f075d59349 Mon Sep 17 00:00:00 2001 From: Alexandre D'Astous Date: Sun, 25 Oct 2020 23:02:48 -0400 Subject: [PATCH 16/51] Add resp path, json path and fix fieldmap path --- test/cli/test_cli_realtime_zshim.py | 13 +++++++++++-- 1 file changed, 11 insertions(+), 2 deletions(-) diff --git a/test/cli/test_cli_realtime_zshim.py b/test/cli/test_cli_realtime_zshim.py index c4fc42cfc..91f5713b8 100644 --- a/test/cli/test_cli_realtime_zshim.py +++ b/test/cli/test_cli_realtime_zshim.py @@ -19,7 +19,8 @@ def test_cli_realtime_zshim(): with tempfile.TemporaryDirectory(prefix='st_' + pathlib.Path(__file__).stem) as tmp: runner = CliRunner() - fname_fieldmap = os.path.join(__dir_testing__, 'sub-fieldmap', 'fmap', 'sub-example_fieldmap.nii.gz') + fname_fieldmap = os.path.join(__dir_testing__, 'realtime_zshimming_data', 'nifti', 'sub-example', 'fmap', + 'sub-example_fieldmap.nii.gz') nii_fmap = nib.load(fname_fieldmap) fmap = nii_fmap.get_fdata() affine = nii_fmap.affine @@ -43,7 +44,15 @@ def test_cli_realtime_zshim(): fname_coil = os.path.join(tmp, 'coil_profile.nii.gz') nib.save(nii_coil, fname_coil) - result = runner.invoke(realtime_zshim, ['-fmap', fname_fieldmap, '-coil', fname_coil, '-mask', fname_mask], + # Path for resp data + fname_resp = os.path.join(__dir_testing__, 'realtime_zshimming_data', 'PMUresp_signal.resp') + + # Path for json file + fname_json = os.path.join(__dir_testing__, 'realtime_zshimming_data', 'nifti', 'sub-example', 'fmap', + 'sub-example_magnitude1.json') + + result = runner.invoke(realtime_zshim, ['-fmap', fname_fieldmap, '-coil', fname_coil, '-mask', fname_mask, + '-resp', fname_resp, '-json', fname_json], catch_exceptions=False) assert result.exit_code == 0 From 010df17dd95e16905ea0a8d3b6148b82dcf0dff5 Mon Sep 17 00:00:00 2001 From: Alexandre D'Astous Date: Sun, 25 Oct 2020 23:14:23 -0400 Subject: [PATCH 17/51] Udated click options and commented debug to verify test --- shimmingtoolbox/cli/realtime_zshim.py | 52 +++++++++++++++------------ 1 file changed, 29 insertions(+), 23 deletions(-) diff --git a/shimmingtoolbox/cli/realtime_zshim.py b/shimmingtoolbox/cli/realtime_zshim.py index bb9cf72c3..64e028112 100644 --- a/shimmingtoolbox/cli/realtime_zshim.py +++ b/shimmingtoolbox/cli/realtime_zshim.py @@ -18,19 +18,24 @@ CONTEXT_SETTINGS = dict(help_option_names=['-h', '--help']) -# @click.command( -# context_settings=CONTEXT_SETTINGS, -# help=f"Perform realtime z-shimming." -# ) -# @click.option('-coil', 'fname_coil', required=True, type=click.Path(), -# help="Coil basis to use for shimming. Enter multiple files if " -# "you wish to use more than one set of shim coils (eg: " -# "Siemens gradient/shim coils and external custom coils).") -# @click.option('-fmap', 'fname_fmap', required=True, type=click.Path(), -# help="B0 fieldmap. For realtime shimming, this should be a 4d file (4th dimension being time") -# @click.option('-mask', 'fname_mask', type=click.Path(), -# help="3D nifti file with voxels between 0 and 1 used to weight the spatial region to shim.") -# @click.option("-verbose", is_flag=True, help="Be more verbose.") +@click.command( + context_settings=CONTEXT_SETTINGS, + help=f"Perform realtime z-shimming." +) +@click.option('-coil', 'fname_coil', required=True, type=click.Path(), + help="Coil basis to use for shimming. Enter multiple files if " + "you wish to use more than one set of shim coils (eg: " + "Siemens gradient/shim coils and external custom coils).") +@click.option('-fmap', 'fname_fmap', required=True, type=click.Path(), + help="B0 fieldmap. For realtime shimming, this should be a 4d file (4th dimension being time") +@click.option('-mask', 'fname_mask', type=click.Path(), + help="3D nifti file with voxels between 0 and 1 used to weight the spatial region to shim.") +@click.option('-resp', 'fname_resp', type=click.Path(), + help="Siemens respiratory file containing pressure data.") +# TODO: Remove json file as input +@click.option('-json', 'fname_json', type=click.Path(), + help="Filename of json corresponding BIDS sidecar.") +@click.option("-verbose", is_flag=True, help="Be more verbose.") def realtime_zshim(fname_coil, fname_fmap, fname_mask, fname_resp, fname_json, verbose=True): """ @@ -38,6 +43,7 @@ def realtime_zshim(fname_coil, fname_fmap, fname_mask, fname_resp, fname_json, v fname_coil: Pointing to coil profile. 4-dimensional: x, y, z, coil. fname_fmap: fname_mask: + fname_resp: verbose: Returns: @@ -127,13 +133,13 @@ def realtime_zshim(fname_coil, fname_fmap, fname_mask, fname_resp, fname_json, v return fname_figure - -fname_coil = os.path.join(__dir_testing__, 'test_realtime_zshim', 'coil_profile.nii.gz') -fname_fmap = os.path.join(__dir_testing__, 'test_realtime_zshim', 'sub-example_fieldmap.nii.gz') -fname_mask = os.path.join(__dir_testing__, 'test_realtime_zshim', 'mask.nii.gz') -fname_resp = os.path.join(__dir_testing__, 'realtime_zshimming_data', 'PMUresp_signal.resp') -fname_json = os.path.join(__dir_testing__, 'test_realtime_zshim', 'sub-example_magnitude1.json') -# fname_coil='/Users/julien/code/shimming-toolbox/shimming-toolbox/test_realtime_zshim/coil_profile.nii.gz' -# fname_fmap='/Users/julien/code/shimming-toolbox/shimming-toolbox/test_realtime_zshim/sub-example_fieldmap.nii.gz' -# fname_mask='/Users/julien/code/shimming-toolbox/shimming-toolbox/test_realtime_zshim/mask.nii.gz' -realtime_zshim(fname_coil, fname_fmap, fname_mask, fname_resp, fname_json) +# Debug +# fname_coil = os.path.join(__dir_testing__, 'test_realtime_zshim', 'coil_profile.nii.gz') +# fname_fmap = os.path.join(__dir_testing__, 'test_realtime_zshim', 'sub-example_fieldmap.nii.gz') +# fname_mask = os.path.join(__dir_testing__, 'test_realtime_zshim', 'mask.nii.gz') +# fname_resp = os.path.join(__dir_testing__, 'realtime_zshimming_data', 'PMUresp_signal.resp') +# fname_json = os.path.join(__dir_testing__, 'test_realtime_zshim', 'sub-example_magnitude1.json') +# # fname_coil='/Users/julien/code/shimming-toolbox/shimming-toolbox/test_realtime_zshim/coil_profile.nii.gz' +# # fname_fmap='/Users/julien/code/shimming-toolbox/shimming-toolbox/test_realtime_zshim/sub-example_fieldmap.nii.gz' +# # fname_mask='/Users/julien/code/shimming-toolbox/shimming-toolbox/test_realtime_zshim/mask.nii.gz' +# realtime_zshim(fname_coil, fname_fmap, fname_mask, fname_resp, fname_json) From cd862185e69e221080dbffa4ddd176f0ed4b971f Mon Sep 17 00:00:00 2001 From: Julien Cohen-Adad Date: Wed, 28 Oct 2020 15:07:31 -0400 Subject: [PATCH 18/51] Added TODOs --- shimmingtoolbox/cli/realtime_zshim.py | 16 ++++++++++++++-- 1 file changed, 14 insertions(+), 2 deletions(-) diff --git a/shimmingtoolbox/cli/realtime_zshim.py b/shimmingtoolbox/cli/realtime_zshim.py index 64e028112..929018e90 100644 --- a/shimmingtoolbox/cli/realtime_zshim.py +++ b/shimmingtoolbox/cli/realtime_zshim.py @@ -49,7 +49,7 @@ def realtime_zshim(fname_coil, fname_fmap, fname_mask, fname_resp, fname_json, v Returns: """ - # When using only z channnel TODO:Remove + # When using only z channnel (corresponding to the 2nd index) TODO:Remove # coil = np.expand_dims(nib.load(fname_coil).get_fdata()[:, :, :, 2], -1) # When using all channels TODO: Keep @@ -124,13 +124,25 @@ def realtime_zshim(fname_coil, fname_fmap, fname_mask, fname_resp, fname_json, v json_data = json.load(json_file) acq_timestamps = get_acquisition_times(nii_fmap, json_data) pmu = PmuResp(fname_resp) + # TODO: deal with saturation acq_pressures = pmu.interp_resp_trace(acq_timestamps) # TODO: # fit PMU and fieldmap values # do regression to separate static componant and RIRO component # output coefficient with proper scaling - + # field(i_vox) = a(i_vox) * acq_pressures + b(i_vox) + # could use: https://scikit-learn.org/stable/modules/generated/sklearn.linear_model.LinearRegression.html + # Note: strong spatial autocorrelation on the a and b coefficients. Ie: two adjacent voxels are submitted to similar + # static B0 field and RIRO component. --> we need to find a way to account for that + # solution 1: post-fitting regularization. + # pros: easy to implement + # cons: fit is less robust to noise + # solution 2: accounting for regularization during fitting + # pros: fitting more robust to noise + # cons: (from Ryan): regularized fitting took a lot of time on Matlab + + # return fname_figure # Debug From a3d7e9d7bc9330b5c4b1f45b893deaa737f4f42a Mon Sep 17 00:00:00 2001 From: Alexandre D'Astous Date: Wed, 28 Oct 2020 22:17:01 -0400 Subject: [PATCH 19/51] Add fit for static and riro linear regression between pressure and field --- setup.py | 1 + shimmingtoolbox/cli/realtime_zshim.py | 29 ++++++++++++++++++++++++++- test/cli/test_cli_realtime_zshim.py | 4 ++-- 3 files changed, 31 insertions(+), 3 deletions(-) diff --git a/setup.py b/setup.py index 6cc12dac1..202419ba1 100644 --- a/setup.py +++ b/setup.py @@ -38,6 +38,7 @@ "matplotlib~=3.1.2", "pytest~=4.6.3", "pytest-cov~=2.5.1", + "sklearn~=0.0", ], extras_require={ 'docs': ["sphinx>=1.6", "sphinx_rtd_theme>=0.2.4"], diff --git a/shimmingtoolbox/cli/realtime_zshim.py b/shimmingtoolbox/cli/realtime_zshim.py index 929018e90..c6fb4b8b1 100644 --- a/shimmingtoolbox/cli/realtime_zshim.py +++ b/shimmingtoolbox/cli/realtime_zshim.py @@ -6,6 +6,7 @@ import os import nibabel as nib import json +from sklearn.linear_model import LinearRegression # TODO: remove matplotlib and dirtesting import from matplotlib.figure import Figure from shimmingtoolbox import __dir_testing__ @@ -50,6 +51,7 @@ def realtime_zshim(fname_coil, fname_fmap, fname_mask, fname_resp, fname_json, v """ # When using only z channnel (corresponding to the 2nd index) TODO:Remove + # TODO: Z index seems to be 0 (maybe because rotation of niftis # coil = np.expand_dims(nib.load(fname_coil).get_fdata()[:, :, :, 2], -1) # When using all channels TODO: Keep @@ -142,7 +144,32 @@ def realtime_zshim(fname_coil, fname_fmap, fname_mask, fname_resp, fname_json, v # pros: fitting more robust to noise # cons: (from Ryan): regularized fitting took a lot of time on Matlab - # + riro = np.zeros_like(fieldmap[:, :, :, 0]) + static = np.zeros_like(fieldmap[:, :, :, 0]) + # TODO: This is ugly + for i_x in range(fieldmap.shape[0]): + for i_y in range(fieldmap.shape[1]): + for i_z in range(fieldmap.shape[2]): + reg = LinearRegression().fit(acq_pressures.reshape(-1, 1), masked_fieldmaps[i_x, i_y, i_z, :]) + riro[i_x, i_y, i_z] = reg.coef_ + static[i_x, i_y, i_z] = reg.intercept_ + + # Plot results + fig = Figure(figsize=(10, 10)) + # FigureCanvas(fig) + ax = fig.add_subplot(2, 1, 1) + im = ax.imshow(riro[:-1, :-1, 0]) + fig.colorbar(im) + ax.set_title("RIRO") + ax = fig.add_subplot(2, 1, 2) + im = ax.imshow(static[:-1, :-1, 0]) + fig.colorbar(im) + ax.set_title("Static") + + fname_figure = os.path.join(__dir_shimmingtoolbox__, 'realtime_zshim_Riro_Static.png') + fig.savefig(fname_figure) + + return fname_figure # Debug diff --git a/test/cli/test_cli_realtime_zshim.py b/test/cli/test_cli_realtime_zshim.py index 91f5713b8..e35a97a3e 100644 --- a/test/cli/test_cli_realtime_zshim.py +++ b/test/cli/test_cli_realtime_zshim.py @@ -29,8 +29,8 @@ def test_cli_realtime_zshim(): nx, ny, nz, nt = fmap.shape mask = shapes(fmap[:, :, :, 0], 'cube', center_dim1=int(fmap.shape[0] / 2 - 8), - center_dim2=int(fmap.shape[1] / 2 - 20), - len_dim1=15, len_dim2=5, len_dim3=nz) + center_dim2=int(fmap.shape[1] / 2 - 5), + len_dim1=15, len_dim2=25, len_dim3=nz) nii_mask = nib.Nifti1Image(mask.astype(int), affine) fname_mask = os.path.join(tmp, 'mask.nii.gz') From 40395cc6396f924348ff94cccbe6320956bdbb92 Mon Sep 17 00:00:00 2001 From: Alexandre D'Astous Date: Wed, 28 Oct 2020 23:08:00 -0400 Subject: [PATCH 20/51] Added plots and refactored other plots --- shimmingtoolbox/cli/realtime_zshim.py | 133 +++++++++++++++++--------- 1 file changed, 90 insertions(+), 43 deletions(-) diff --git a/shimmingtoolbox/cli/realtime_zshim.py b/shimmingtoolbox/cli/realtime_zshim.py index c6fb4b8b1..d8ad0052e 100644 --- a/shimmingtoolbox/cli/realtime_zshim.py +++ b/shimmingtoolbox/cli/realtime_zshim.py @@ -77,48 +77,11 @@ def realtime_zshim(fname_coil, fname_fmap, fname_mask, fname_resp, fname_json, v shimmed = np.zeros_like(fieldmap) masked_fieldmaps = np.zeros_like(fieldmap) - masked_shimmed = np.zeros_like(shimmed) for i_t in range(nt): currents[:, i_t] = sequential_zslice(fieldmap[..., i_t], coil, mask, z_slices=np.array(range(nz)), - bounds=[(-np.inf, np.inf)]*n_coils) + bounds=[(-np.inf, np.inf)] * n_coils) shimmed[..., i_t] = fieldmap[..., i_t] + np.sum(currents[:, i_t] * coil, axis=3, keepdims=False) masked_fieldmaps[..., i_t] = mask * fieldmap[..., i_t] - masked_shimmed[..., i_t] = mask * shimmed[..., i_t] - - i_t = 0 - # Plot results - fig = Figure(figsize=(10, 10)) - # FigureCanvas(fig) - ax = fig.add_subplot(2, 2, 1) - im = ax.imshow(masked_fieldmaps[:-1, :-1, 0, i_t]) - fig.colorbar(im) - ax.set_title("Masked unshimmed") - ax = fig.add_subplot(2, 2, 2) - im = ax.imshow(masked_shimmed[:-1, :-1, 0, i_t]) - fig.colorbar(im) - ax.set_title("Masked shimmed") - ax = fig.add_subplot(2, 2, 3) - im = ax.imshow(fieldmap[:-1, :-1, 0, i_t]) - fig.colorbar(im) - ax.set_title("Unshimmed") - ax = fig.add_subplot(2, 2, 4) - im = ax.imshow(shimmed[:-1, :-1, 0, i_t]) - fig.colorbar(im) - ax.set_title("Shimmed") - - click.echo(f"\nThe associated current coefficients are : {currents[:, i_t]}") - - fname_figure = os.path.join(__dir_shimmingtoolbox__, 'realtime_zshim_plot.png') - fig.savefig(fname_figure) - - fig = Figure(figsize=(10, 10)) - for i_coil in range(n_coils): - ax = fig.add_subplot(n_coils, 1, i_coil + 1) - ax.plot(np.arange(nt), currents[i_coil, :]) - ax.set_title(f"Channel {i_coil}") - - fname_figure = os.path.join(__dir_shimmingtoolbox__, 'realtime_zshim_currents.png') - fig.savefig(fname_figure) # Fetch PMU timing # TODO: Add json to fieldmap instead of asking for another json file @@ -146,17 +109,54 @@ def realtime_zshim(fname_coil, fname_fmap, fname_mask, fname_resp, fname_json, v riro = np.zeros_like(fieldmap[:, :, :, 0]) static = np.zeros_like(fieldmap[:, :, :, 0]) - # TODO: This is ugly for i_x in range(fieldmap.shape[0]): for i_y in range(fieldmap.shape[1]): for i_z in range(fieldmap.shape[2]): - reg = LinearRegression().fit(acq_pressures.reshape(-1, 1), masked_fieldmaps[i_x, i_y, i_z, :]) + # TODO: Fit for -masked_field? + reg = LinearRegression().fit(acq_pressures.reshape(-1, 1), -masked_fieldmaps[i_x, i_y, i_z, :]) riro[i_x, i_y, i_z] = reg.coef_ static[i_x, i_y, i_z] = reg.intercept_ - # Plot results + # ================ PLOTS ================ + + # Calculate masked shim for spherical harmonics plot + masked_shimmed = np.zeros_like(shimmed) + for i_t in range(nt): + masked_shimmed[..., i_t] = mask * shimmed[..., i_t] + + # Plot unshimmed vs shimmed and their mask for spherical harmonics + i_t = 0 + fig = Figure(figsize=(10, 10)) + ax = fig.add_subplot(2, 2, 1) + im = ax.imshow(masked_fieldmaps[:-1, :-1, 0, i_t]) + fig.colorbar(im) + ax.set_title("Masked unshimmed") + ax = fig.add_subplot(2, 2, 2) + im = ax.imshow(masked_shimmed[:-1, :-1, 0, i_t]) + fig.colorbar(im) + ax.set_title("Masked shimmed") + ax = fig.add_subplot(2, 2, 3) + im = ax.imshow(fieldmap[:-1, :-1, 0, i_t]) + fig.colorbar(im) + ax.set_title("Unshimmed") + ax = fig.add_subplot(2, 2, 4) + im = ax.imshow(shimmed[:-1, :-1, 0, i_t]) + fig.colorbar(im) + ax.set_title("Shimmed") + fname_figure = os.path.join(__dir_shimmingtoolbox__, 'realtime_zshim_sphharm_shimmed.png') + fig.savefig(fname_figure) + + # Plot the coil coefs through time + fig = Figure(figsize=(10, 10)) + for i_coil in range(n_coils): + ax = fig.add_subplot(n_coils, 1, i_coil + 1) + ax.plot(np.arange(nt), currents[i_coil, :]) + ax.set_title(f"Channel {i_coil}") + fname_figure = os.path.join(__dir_shimmingtoolbox__, 'realtime_zshim_sphharm_currents.png') + fig.savefig(fname_figure) + + # Plot Static and RIRO fig = Figure(figsize=(10, 10)) - # FigureCanvas(fig) ax = fig.add_subplot(2, 1, 1) im = ax.imshow(riro[:-1, :-1, 0]) fig.colorbar(im) @@ -165,10 +165,57 @@ def realtime_zshim(fname_coil, fname_fmap, fname_mask, fname_resp, fname_json, v im = ax.imshow(static[:-1, :-1, 0]) fig.colorbar(im) ax.set_title("Static") + fname_figure = os.path.join(__dir_shimmingtoolbox__, 'realtime_zshim_riro_static.png') + fig.savefig(fname_figure) + + # Calculate fitted and shimmed for pressure fitted plot + fitted_fieldmap = riro * acq_pressures + static + shimmed_pressure_fitted = np.expand_dims(fitted_fieldmap, 2) + masked_fieldmaps - fname_figure = os.path.join(__dir_shimmingtoolbox__, 'realtime_zshim_Riro_Static.png') + # Plot pressure fitted fieldmap + fig = Figure(figsize=(10, 10)) + ax = fig.add_subplot(3, 1, 1) + im = ax.imshow(masked_fieldmaps[:-1, :-1, 0, i_t]) + fig.colorbar(im) + ax.set_title("fieldmap") + ax = fig.add_subplot(3, 1, 2) + im = ax.imshow(fitted_fieldmap[:-1, :-1, i_t]) + fig.colorbar(im) + ax.set_title("Fit") + ax = fig.add_subplot(3, 1, 3) + im = ax.imshow(shimmed_pressure_fitted[:-1, :-1, 0, i_t]) + fig.colorbar(im) + ax.set_title("Shimmed (fit + fieldmap") + fname_figure = os.path.join(__dir_shimmingtoolbox__, 'realtime_zshim_pressure_fitted.png') fig.savefig(fname_figure) + # Reshape pmu datapoints to fit those of the acquisition + pmu_times = np.linspace(pmu.start_time_mdh, pmu.stop_time_mdh, len(pmu.data)) + pmu_times_within_range = pmu_times[pmu_times > acq_timestamps[0]] + pmu_data_within_range = pmu.data[pmu_times > acq_timestamps[0]] + pmu_data_within_range = pmu_data_within_range[pmu_times_within_range < acq_timestamps[fieldmap.shape[3] - 1]] + pmu_times_within_range = pmu_times_within_range[pmu_times_within_range < acq_timestamps[fieldmap.shape[3] - 1]] + + # Calc fieldmap average within mask + fieldmap_avg = np.zeros([fieldmap.shape[3]]) + for i_time in range(nt): + masked_array = np.ma.array(fieldmap[:, :, :, i_time], mask=mask == False) + fieldmap_avg[i_time] = np.ma.average(masked_array) + + # Plot pmu vs B0 in masked region + fig = Figure(figsize=(10, 10)) + ax = fig.add_subplot(211) + ax.plot(acq_timestamps / 1000, acq_pressures, label='Interpolated pressures') + # ax.plot(pmu_times / 1000, pmu.data, label='Raw pressures') + ax.plot(pmu_times_within_range / 1000, pmu_data_within_range, label='Pmu pressures') + ax.legend() + ax.set_title("Pressure [-2048, 2047] vs time (s) ") + ax = fig.add_subplot(212) + ax.plot(acq_timestamps / 1000, fieldmap_avg, label='Mean B0') + ax.legend() + ax.set_title("Fieldmap average over unmasked region (Hz) vs time (s)") + fname_figure = os.path.join(__dir_shimmingtoolbox__, 'realtime_zshim_pmu_vs_B0.png') + fig.savefig(fname_figure) return fname_figure From 94a685f588574d8d3bf16251ecfe90b4f83e493d Mon Sep 17 00:00:00 2001 From: Alexandre D'Astous Date: Wed, 28 Oct 2020 23:20:42 -0400 Subject: [PATCH 21/51] Add comments --- shimmingtoolbox/cli/realtime_zshim.py | 8 +++++--- 1 file changed, 5 insertions(+), 3 deletions(-) diff --git a/shimmingtoolbox/cli/realtime_zshim.py b/shimmingtoolbox/cli/realtime_zshim.py index d8ad0052e..d0683ba6f 100644 --- a/shimmingtoolbox/cli/realtime_zshim.py +++ b/shimmingtoolbox/cli/realtime_zshim.py @@ -50,13 +50,14 @@ def realtime_zshim(fname_coil, fname_fmap, fname_mask, fname_resp, fname_json, v Returns: """ + # Load coil # When using only z channnel (corresponding to the 2nd index) TODO:Remove # TODO: Z index seems to be 0 (maybe because rotation of niftis # coil = np.expand_dims(nib.load(fname_coil).get_fdata()[:, :, :, 2], -1) - # When using all channels TODO: Keep coil = nib.load(fname_coil).get_fdata() + # Load fieldmap nii_fmap = nib.load(fname_fmap) fieldmap = nii_fmap.get_fdata() @@ -65,16 +66,16 @@ def realtime_zshim(fname_coil, fname_fmap, fname_mask, fname_resp, fname_json, v raise RuntimeError('fmap must be 4d (x, y, z, t)') nx, ny, nz, nt = fieldmap.shape + # Load mask # TODO: check good practice below if fname_mask is not None: mask = nib.load(fname_mask).get_fdata() else: mask = np.ones_like(fieldmap) - # Setup coil + # Shim using sequencer and optimizer n_coils = coil.shape[-1] currents = np.zeros([n_coils, nt]) - shimmed = np.zeros_like(fieldmap) masked_fieldmaps = np.zeros_like(fieldmap) for i_t in range(nt): @@ -107,6 +108,7 @@ def realtime_zshim(fname_coil, fname_fmap, fname_mask, fname_resp, fname_json, v # pros: fitting more robust to noise # cons: (from Ryan): regularized fitting took a lot of time on Matlab + # Shim using PMU riro = np.zeros_like(fieldmap[:, :, :, 0]) static = np.zeros_like(fieldmap[:, :, :, 0]) for i_x in range(fieldmap.shape[0]): From ea1e7d1214e146f87ab5ec8c4d974a9de38c477f Mon Sep 17 00:00:00 2001 From: Alexandre D'Astous Date: Mon, 2 Nov 2020 21:59:57 -0500 Subject: [PATCH 22/51] Rescale PMU and correct index for z component --- shimmingtoolbox/cli/realtime_zshim.py | 7 +++---- 1 file changed, 3 insertions(+), 4 deletions(-) diff --git a/shimmingtoolbox/cli/realtime_zshim.py b/shimmingtoolbox/cli/realtime_zshim.py index d0683ba6f..2b0e80322 100644 --- a/shimmingtoolbox/cli/realtime_zshim.py +++ b/shimmingtoolbox/cli/realtime_zshim.py @@ -51,9 +51,8 @@ def realtime_zshim(fname_coil, fname_fmap, fname_mask, fname_resp, fname_json, v """ # Load coil - # When using only z channnel (corresponding to the 2nd index) TODO:Remove - # TODO: Z index seems to be 0 (maybe because rotation of niftis - # coil = np.expand_dims(nib.load(fname_coil).get_fdata()[:, :, :, 2], -1) + # When using only z channel (corresponding to index 0) TODO:Remove + # coil = np.expand_dims(nib.load(fname_coil).get_fdata()[:, :, :, 0], -1) # When using all channels TODO: Keep coil = nib.load(fname_coil).get_fdata() @@ -91,7 +90,7 @@ def realtime_zshim(fname_coil, fname_fmap, fname_mask, fname_resp, fname_json, v acq_timestamps = get_acquisition_times(nii_fmap, json_data) pmu = PmuResp(fname_resp) # TODO: deal with saturation - acq_pressures = pmu.interp_resp_trace(acq_timestamps) + acq_pressures = pmu.interp_resp_trace(acq_timestamps) + 2048 # [0, 4095] # TODO: # fit PMU and fieldmap values From b0eba3d63573d227e7dceadc10fc5f663c1d6624 Mon Sep 17 00:00:00 2001 From: Alexandre D'Astous Date: Mon, 2 Nov 2020 22:21:57 -0500 Subject: [PATCH 23/51] Add ability to mask with threshold --- test/cli/test_cli_realtime_zshim.py | 19 +++++++++++++------ 1 file changed, 13 insertions(+), 6 deletions(-) diff --git a/test/cli/test_cli_realtime_zshim.py b/test/cli/test_cli_realtime_zshim.py index e35a97a3e..1fdfc5316 100644 --- a/test/cli/test_cli_realtime_zshim.py +++ b/test/cli/test_cli_realtime_zshim.py @@ -10,6 +10,7 @@ from click.testing import CliRunner from shimmingtoolbox.cli.realtime_zshim import realtime_zshim from shimmingtoolbox.masking.shapes import shapes +from shimmingtoolbox.masking.threshold import threshold from shimmingtoolbox.coils.coordinates import generate_meshgrid from shimmingtoolbox.coils.siemens_basis import siemens_basis from shimmingtoolbox import __dir_testing__ @@ -26,13 +27,19 @@ def test_cli_realtime_zshim(): affine = nii_fmap.affine # Set up mask - nx, ny, nz, nt = fmap.shape - mask = shapes(fmap[:, :, :, 0], 'cube', - center_dim1=int(fmap.shape[0] / 2 - 8), - center_dim2=int(fmap.shape[1] / 2 - 5), - len_dim1=15, len_dim2=25, len_dim3=nz) + # Cube + # nx, ny, nz, nt = fmap.shape + # mask = shapes(fmap[:, :, :, 0], 'cube', + # center_dim1=int(fmap.shape[0] / 2 - 8), + # center_dim2=int(fmap.shape[1] / 2 - 5), + # len_dim1=15, len_dim2=25, len_dim3=nz) + # Threshold + fname_mag = os.path.join(__dir_testing__, 'realtime_zshimming_data', 'nifti', 'sub-example', 'fmap', + 'sub-example_magnitude1.nii.gz') + mag = nib.load(fname_mag).get_fdata() + mask = threshold(mag, thr=50) - nii_mask = nib.Nifti1Image(mask.astype(int), affine) + nii_mask = nib.Nifti1Image(mask.astype(int)[:, :, :, 0], affine) fname_mask = os.path.join(tmp, 'mask.nii.gz') nib.save(nii_mask, fname_mask) From 976f13e2c71e79b1977eccffd21b7848b1cd2117 Mon Sep 17 00:00:00 2001 From: Alexandre D'Astous Date: Tue, 3 Nov 2020 00:22:07 -0500 Subject: [PATCH 24/51] Add resampling of fieldmap, riro and static to target acquisition --- setup.py | 1 + shimmingtoolbox/cli/realtime_zshim.py | 40 +++++++++++++++++++++++++-- test/cli/test_cli_realtime_zshim.py | 6 +++- 3 files changed, 44 insertions(+), 3 deletions(-) diff --git a/setup.py b/setup.py index 202419ba1..aba63801e 100644 --- a/setup.py +++ b/setup.py @@ -39,6 +39,7 @@ "pytest~=4.6.3", "pytest-cov~=2.5.1", "sklearn~=0.0", + "nilearn~=0.6.2" ], extras_require={ 'docs': ["sphinx>=1.6", "sphinx_rtd_theme>=0.2.4"], diff --git a/shimmingtoolbox/cli/realtime_zshim.py b/shimmingtoolbox/cli/realtime_zshim.py index 2b0e80322..4550a5707 100644 --- a/shimmingtoolbox/cli/realtime_zshim.py +++ b/shimmingtoolbox/cli/realtime_zshim.py @@ -7,6 +7,7 @@ import nibabel as nib import json from sklearn.linear_model import LinearRegression +from nilearn.image import resample_img # TODO: remove matplotlib and dirtesting import from matplotlib.figure import Figure from shimmingtoolbox import __dir_testing__ @@ -33,11 +34,13 @@ help="3D nifti file with voxels between 0 and 1 used to weight the spatial region to shim.") @click.option('-resp', 'fname_resp', type=click.Path(), help="Siemens respiratory file containing pressure data.") +@click.option('-anat', 'fname_anat', type=click.Path(), + help="Filename of the anatomical image to apply the correction.") # TODO: Remove json file as input @click.option('-json', 'fname_json', type=click.Path(), help="Filename of json corresponding BIDS sidecar.") @click.option("-verbose", is_flag=True, help="Be more verbose.") -def realtime_zshim(fname_coil, fname_fmap, fname_mask, fname_resp, fname_json, verbose=True): +def realtime_zshim(fname_coil, fname_fmap, fname_mask, fname_resp, fname_json, fname_anat, verbose=True): """ Args: @@ -62,7 +65,7 @@ def realtime_zshim(fname_coil, fname_fmap, fname_mask, fname_resp, fname_json, v # TODO: Error handling might move to API if fieldmap.ndim != 4: - raise RuntimeError('fmap must be 4d (x, y, z, t)') + raise RuntimeError("fmap must be 4d (x, y, z, t)") nx, ny, nz, nt = fieldmap.shape # Load mask @@ -72,6 +75,12 @@ def realtime_zshim(fname_coil, fname_fmap, fname_mask, fname_resp, fname_json, v else: mask = np.ones_like(fieldmap) + # Load anat + nii_anat = nib.load(fname_anat) + anat = nii_anat.get_fdata() + if anat.ndim != 3: + raise RuntimeError("Anatomical image must be in 3d") + # Shim using sequencer and optimizer n_coils = coil.shape[-1] currents = np.zeros([n_coils, nt]) @@ -118,6 +127,20 @@ def realtime_zshim(fname_coil, fname_fmap, fname_mask, fname_resp, fname_json, v riro[i_x, i_y, i_z] = reg.coef_ static[i_x, i_y, i_z] = reg.intercept_ + # Resample masked_fieldmaps, riro and static to target anatomical image + nii_masked_fmap = nib.Nifti1Image(masked_fieldmaps, nii_fmap.affine) + nii_riro = nib.Nifti1Image(riro, nii_fmap.affine) + nii_static = nib.Nifti1Image(static, nii_fmap.affine) + + target_affine = nii_anat.affine[:-1, :-1] + nii_resampled_fmap = resample_img(nii_masked_fmap, target_affine=target_affine, interpolation='nearest') + nii_resampled_riro = resample_img(nii_riro, target_affine=target_affine, interpolation='nearest') + nii_resampled_static = resample_img(nii_static, target_affine=target_affine, interpolation='nearest') + + nib.save(nii_resampled_fmap, os.path.join(__dir_shimmingtoolbox__, 'resampled_fmap.nii.gz')) + nib.save(nii_resampled_riro, os.path.join(__dir_shimmingtoolbox__, 'resampled_riro.nii.gz')) + nib.save(nii_resampled_static, os.path.join(__dir_shimmingtoolbox__, 'resampled_static.nii.gz')) + # ================ PLOTS ================ # Calculate masked shim for spherical harmonics plot @@ -218,6 +241,19 @@ def realtime_zshim(fname_coil, fname_fmap, fname_mask, fname_resp, fname_json, v fname_figure = os.path.join(__dir_shimmingtoolbox__, 'realtime_zshim_pmu_vs_B0.png') fig.savefig(fname_figure) + # Show anatomical image + fig = Figure(figsize=(10, 10)) + ax = fig.add_subplot(2, 1, 1) + im = ax.imshow(anat[:-1, :-1, 10]) + fig.colorbar(im) + ax.set_title("Anatomical image [:-1, :-1, 10]") + ax = fig.add_subplot(2, 1, 2) + im = ax.imshow(nii_resampled_fmap.get_fdata()[0, :-1, :-1, 0]) + fig.colorbar(im) + ax.set_title("Resampled fieldmap [0, :-1, :-1, 0]") + fname_figure = os.path.join(__dir_shimmingtoolbox__, 'reatime_zshime_anat.png') + fig.savefig(fname_figure) + return fname_figure # Debug diff --git a/test/cli/test_cli_realtime_zshim.py b/test/cli/test_cli_realtime_zshim.py index 1fdfc5316..c8f79b150 100644 --- a/test/cli/test_cli_realtime_zshim.py +++ b/test/cli/test_cli_realtime_zshim.py @@ -58,8 +58,12 @@ def test_cli_realtime_zshim(): fname_json = os.path.join(__dir_testing__, 'realtime_zshimming_data', 'nifti', 'sub-example', 'fmap', 'sub-example_magnitude1.json') + # Path for mag anat image + fname_anat = os.path.join(__dir_testing__, 'realtime_zshimming_data', 'nifti', 'sub-example', 'anat', + 'sub-example_unshimmed_e1.nii.gz') + result = runner.invoke(realtime_zshim, ['-fmap', fname_fieldmap, '-coil', fname_coil, '-mask', fname_mask, - '-resp', fname_resp, '-json', fname_json], + '-resp', fname_resp, '-json', fname_json, '-anat', fname_anat], catch_exceptions=False) assert result.exit_code == 0 From f951c635657897e8c318964aff3a226b8d3d88d2 Mon Sep 17 00:00:00 2001 From: Alexandre D'Astous Date: Tue, 3 Nov 2020 15:25:43 -0500 Subject: [PATCH 25/51] Use nibabel instead of nilearn for resampling --- setup.py | 1 - shimmingtoolbox/cli/realtime_zshim.py | 22 +++++++++++++--------- 2 files changed, 13 insertions(+), 10 deletions(-) diff --git a/setup.py b/setup.py index aba63801e..202419ba1 100644 --- a/setup.py +++ b/setup.py @@ -39,7 +39,6 @@ "pytest~=4.6.3", "pytest-cov~=2.5.1", "sklearn~=0.0", - "nilearn~=0.6.2" ], extras_require={ 'docs': ["sphinx>=1.6", "sphinx_rtd_theme>=0.2.4"], diff --git a/shimmingtoolbox/cli/realtime_zshim.py b/shimmingtoolbox/cli/realtime_zshim.py index 4550a5707..12484133d 100644 --- a/shimmingtoolbox/cli/realtime_zshim.py +++ b/shimmingtoolbox/cli/realtime_zshim.py @@ -7,7 +7,7 @@ import nibabel as nib import json from sklearn.linear_model import LinearRegression -from nilearn.image import resample_img +from nibabel.processing import resample_from_to # TODO: remove matplotlib and dirtesting import from matplotlib.figure import Figure from shimmingtoolbox import __dir_testing__ @@ -99,7 +99,7 @@ def realtime_zshim(fname_coil, fname_fmap, fname_mask, fname_resp, fname_json, f acq_timestamps = get_acquisition_times(nii_fmap, json_data) pmu = PmuResp(fname_resp) # TODO: deal with saturation - acq_pressures = pmu.interp_resp_trace(acq_timestamps) + 2048 # [0, 4095] + acq_pressures = pmu.interp_resp_trace(acq_timestamps) # TODO: # fit PMU and fieldmap values @@ -128,14 +128,18 @@ def realtime_zshim(fname_coil, fname_fmap, fname_mask, fname_resp, fname_json, f static[i_x, i_y, i_z] = reg.intercept_ # Resample masked_fieldmaps, riro and static to target anatomical image - nii_masked_fmap = nib.Nifti1Image(masked_fieldmaps, nii_fmap.affine) + # TODO: convert to a function + masked_fmap_4d = np.zeros(anat.shape + (nt,)) + for it in range(nt): + nii_masked_fmap_3d = nib.Nifti1Image(masked_fieldmaps[..., it], nii_fmap.affine) + nii_resampled_fmap_3d = resample_from_to(nii_masked_fmap_3d, nii_anat, mode='nearest') + masked_fmap_4d[..., it] = nii_resampled_fmap_3d.get_fdata() + + nii_resampled_fmap = nib.Nifti1Image(masked_fmap_4d, nii_anat.affine) nii_riro = nib.Nifti1Image(riro, nii_fmap.affine) nii_static = nib.Nifti1Image(static, nii_fmap.affine) - - target_affine = nii_anat.affine[:-1, :-1] - nii_resampled_fmap = resample_img(nii_masked_fmap, target_affine=target_affine, interpolation='nearest') - nii_resampled_riro = resample_img(nii_riro, target_affine=target_affine, interpolation='nearest') - nii_resampled_static = resample_img(nii_static, target_affine=target_affine, interpolation='nearest') + nii_resampled_riro = resample_from_to(nii_riro, nii_anat, mode='nearest') + nii_resampled_static = resample_from_to(nii_static, nii_anat, mode='nearest') nib.save(nii_resampled_fmap, os.path.join(__dir_shimmingtoolbox__, 'resampled_fmap.nii.gz')) nib.save(nii_resampled_riro, os.path.join(__dir_shimmingtoolbox__, 'resampled_riro.nii.gz')) @@ -250,7 +254,7 @@ def realtime_zshim(fname_coil, fname_fmap, fname_mask, fname_resp, fname_json, f ax = fig.add_subplot(2, 1, 2) im = ax.imshow(nii_resampled_fmap.get_fdata()[0, :-1, :-1, 0]) fig.colorbar(im) - ax.set_title("Resampled fieldmap [0, :-1, :-1, 0]") + ax.set_title("Resampled fieldmap [:-1, :-1, 0, 0]") fname_figure = os.path.join(__dir_shimmingtoolbox__, 'reatime_zshime_anat.png') fig.savefig(fname_figure) From dbafdd774a4c031745957220fd64fc18bf77012f Mon Sep 17 00:00:00 2001 From: Alexandre D'Astous Date: Tue, 3 Nov 2020 17:23:41 -0500 Subject: [PATCH 26/51] Added gradient calculation --- shimmingtoolbox/cli/realtime_zshim.py | 29 +++++++++++++++++++++++---- 1 file changed, 25 insertions(+), 4 deletions(-) diff --git a/shimmingtoolbox/cli/realtime_zshim.py b/shimmingtoolbox/cli/realtime_zshim.py index 12484133d..732ad9abb 100644 --- a/shimmingtoolbox/cli/realtime_zshim.py +++ b/shimmingtoolbox/cli/realtime_zshim.py @@ -92,6 +92,18 @@ def realtime_zshim(fname_coil, fname_fmap, fname_mask, fname_resp, fname_json, f shimmed[..., i_t] = fieldmap[..., i_t] + np.sum(currents[:, i_t] * coil, axis=3, keepdims=False) masked_fieldmaps[..., i_t] = mask * fieldmap[..., i_t] + # Calculate gz gradient + # Image is z, y, x + # Pixdim[3] is the space between pixels in the z direction in millimeters + # TODO: This is the gradient is the voxel reference frame meaning the axes could possibly not follow the scanner's + # axes (therefore the z gradient), also investigate is axes should be 1 or 0 + g = 1000 / 42.576e6 # [mT / Hz] + gz_gradient = np.zeros_like(masked_fieldmaps) + for it in range(nt): + gz_gradient[..., 0, it] = np.gradient(g * masked_fieldmaps[:, :, 0, it], + nii_fmap.header['pixdim'][3] / 1000, + axis=1) # [mT / m] + # Fetch PMU timing # TODO: Add json to fieldmap instead of asking for another json file with open(fname_json) as json_file: @@ -248,16 +260,25 @@ def realtime_zshim(fname_coil, fname_fmap, fname_mask, fname_resp, fname_json, f # Show anatomical image fig = Figure(figsize=(10, 10)) ax = fig.add_subplot(2, 1, 1) - im = ax.imshow(anat[:-1, :-1, 10]) + im = ax.imshow(anat[:, :, 10]) fig.colorbar(im) - ax.set_title("Anatomical image [:-1, :-1, 10]") + ax.set_title("Anatomical image [:, :, 10]") ax = fig.add_subplot(2, 1, 2) - im = ax.imshow(nii_resampled_fmap.get_fdata()[0, :-1, :-1, 0]) + im = ax.imshow(nii_resampled_fmap.get_fdata()[:, :, 10, 0]) fig.colorbar(im) - ax.set_title("Resampled fieldmap [:-1, :-1, 0, 0]") + ax.set_title("Resampled fieldmap [:, :, 10, 0]") fname_figure = os.path.join(__dir_shimmingtoolbox__, 'reatime_zshime_anat.png') fig.savefig(fname_figure) + # Show Gradient + fig = Figure(figsize=(10, 10)) + ax = fig.add_subplot(1, 1, 1) + im = ax.imshow(gz_gradient[:, :, 0, 0]) + fig.colorbar(im) + ax.set_title("Gradient [:, :, 0, 0]") + fname_figure = os.path.join(__dir_shimmingtoolbox__, 'reatime_zshime_gradient.png') + fig.savefig(fname_figure) + return fname_figure # Debug From 64cae4434489a6d3a6e0958cf209684e41a7c1b4 Mon Sep 17 00:00:00 2001 From: Alexandre D'Astous Date: Tue, 3 Nov 2020 23:43:49 -0500 Subject: [PATCH 27/51] Output pmu from 0 to 4095 to be consistent with matlab code --- shimmingtoolbox/pmu.py | 4 ++-- test/test_pmu.py | 4 ++-- 2 files changed, 4 insertions(+), 4 deletions(-) diff --git a/shimmingtoolbox/pmu.py b/shimmingtoolbox/pmu.py index 522c5a261..62671313f 100644 --- a/shimmingtoolbox/pmu.py +++ b/shimmingtoolbox/pmu.py @@ -14,7 +14,7 @@ class PmuResp(object): Attributes: fname (str): Filename of the Siemens .resp file - data (numpy.ndarray): Pressure values ranging from -2048 to 2047 + data (numpy.ndarray): Pressure values ranging from 0 to 4095 start_time_mdh (int): Start time in milliseconds past midnight (mdh clock is expected to be the closest to the image header) stop_time_mdh (int): Stop time in milliseconds past midnight (mdh clock is expected to be the closest to @@ -107,7 +107,7 @@ def read_resp(self, fname_pmu): # files limiting the data to points in the range 0..4095 will give us just the first channel's data. # We are assuming that the 5000/6000 marks are inserted into the data values list rather than overwriting. # That is we assume they do *not* occupy a raster position. - data_cleaned = data[data < 4096] - 2048 + data_cleaned = data[data < 4096] attributes = { 'fname': fname_pmu, diff --git a/test/test_pmu.py b/test/test_pmu.py index fb19b6066..11c095629 100644 --- a/test/test_pmu.py +++ b/test/test_pmu.py @@ -17,8 +17,8 @@ def test_read_resp(): fname = os.path.join(__dir_testing__, 'realtime_zshimming_data', 'PMUresp_signal.resp') pmu = PmuResp(fname) - expected_first_data = np.array([1667, 1667, 1682, 1682, 1667]) - expected_last_data = np.array([-2048, -2048, -2048, -2048, -2048]) + expected_first_data = np.array([1667, 1667, 1682, 1682, 1667]) + 2048 + expected_last_data = np.array([-2048, -2048, -2048, -2048, -2048]) + 2048 expected_start_time_mdh = 44294387 expected_stop_time_mdh = 44343130 expected_start_time_mpcu = 44294295 From 671d411d3a17010bd89f6ee10e915dd5c847493c Mon Sep 17 00:00:00 2001 From: Alexandre D'Astous Date: Tue, 3 Nov 2020 23:44:23 -0500 Subject: [PATCH 28/51] Calculate mean for each slice --- shimmingtoolbox/cli/realtime_zshim.py | 8 ++++++++ 1 file changed, 8 insertions(+) diff --git a/shimmingtoolbox/cli/realtime_zshim.py b/shimmingtoolbox/cli/realtime_zshim.py index 732ad9abb..34fc9947a 100644 --- a/shimmingtoolbox/cli/realtime_zshim.py +++ b/shimmingtoolbox/cli/realtime_zshim.py @@ -157,6 +157,14 @@ def realtime_zshim(fname_coil, fname_fmap, fname_mask, fname_resp, fname_json, f nib.save(nii_resampled_riro, os.path.join(__dir_shimmingtoolbox__, 'resampled_riro.nii.gz')) nib.save(nii_resampled_static, os.path.join(__dir_shimmingtoolbox__, 'resampled_static.nii.gz')) + # Calculate the mean for riro and static for a perticular slice + n_slices = nii_resampled_fmap.get_fdata().shape[2] + static_correction = np.zeros([n_slices]) + riro_correction = np.zeros([n_slices]) + for i_slice in range(n_slices): + static_correction[i_slice] = np.mean(nii_resampled_static.get_fdata()[..., i_slice]) + riro_correction[i_slice] = np.mean(nii_resampled_riro.get_fdata()[..., i_slice]) + # ================ PLOTS ================ # Calculate masked shim for spherical harmonics plot From 9601591bab3cb0a7c183a17fb76b49984d6a4782 Mon Sep 17 00:00:00 2001 From: Alexandre D'Astous Date: Wed, 4 Nov 2020 00:03:12 -0500 Subject: [PATCH 29/51] Write correction coefficients to a file --- shimmingtoolbox/cli/realtime_zshim.py | 14 ++++++++++++-- 1 file changed, 12 insertions(+), 2 deletions(-) diff --git a/shimmingtoolbox/cli/realtime_zshim.py b/shimmingtoolbox/cli/realtime_zshim.py index 34fc9947a..9f0b448d5 100644 --- a/shimmingtoolbox/cli/realtime_zshim.py +++ b/shimmingtoolbox/cli/realtime_zshim.py @@ -95,8 +95,8 @@ def realtime_zshim(fname_coil, fname_fmap, fname_mask, fname_resp, fname_json, f # Calculate gz gradient # Image is z, y, x # Pixdim[3] is the space between pixels in the z direction in millimeters - # TODO: This is the gradient is the voxel reference frame meaning the axes could possibly not follow the scanner's - # axes (therefore the z gradient), also investigate is axes should be 1 or 0 + # TODO: This is the gradient in the voxel reference frame meaning the axes could possibly not follow the scanner's + # axes (therefore the z gradient), also investigate if axes should be 1 or 0 g = 1000 / 42.576e6 # [mT / Hz] gz_gradient = np.zeros_like(masked_fieldmaps) for it in range(nt): @@ -135,6 +135,7 @@ def realtime_zshim(fname_coil, fname_fmap, fname_mask, fname_resp, fname_json, f for i_y in range(fieldmap.shape[1]): for i_z in range(fieldmap.shape[2]): # TODO: Fit for -masked_field? + # reg = LinearRegression().fit(acq_pressures.reshape(-1, 1), -gz_gradient[i_x, i_y, i_z, :]) reg = LinearRegression().fit(acq_pressures.reshape(-1, 1), -masked_fieldmaps[i_x, i_y, i_z, :]) riro[i_x, i_y, i_z] = reg.coef_ static[i_x, i_y, i_z] = reg.intercept_ @@ -165,6 +166,15 @@ def realtime_zshim(fname_coil, fname_fmap, fname_mask, fname_resp, fname_json, f static_correction[i_slice] = np.mean(nii_resampled_static.get_fdata()[..., i_slice]) riro_correction[i_slice] = np.mean(nii_resampled_riro.get_fdata()[..., i_slice]) + # Write to a text file + fname_corrections = os.path.join(__dir_shimmingtoolbox__, 'zshim_gradients.txt') + file_gradients = open(fname_corrections, 'w') + for i_slice in range(n_slices): + file_gradients.write(f'Vector_Gz[0][{i_slice}]= {static_correction[i_slice]:.6f}\n') + file_gradients.write(f'Vector_Gz[1][{i_slice}]= {riro_correction[i_slice]:.12f}\n') + # Matlab includes the mean pressure + file_gradients.close() + # ================ PLOTS ================ # Calculate masked shim for spherical harmonics plot From 6171af39ec101fefa4c9b6398f2444c2aa9c87a4 Mon Sep 17 00:00:00 2001 From: Alexandre D'Astous Date: Wed, 4 Nov 2020 00:10:02 -0500 Subject: [PATCH 30/51] Comments --- shimmingtoolbox/cli/realtime_zshim.py | 5 +++-- 1 file changed, 3 insertions(+), 2 deletions(-) diff --git a/shimmingtoolbox/cli/realtime_zshim.py b/shimmingtoolbox/cli/realtime_zshim.py index 9f0b448d5..d7ac02b7d 100644 --- a/shimmingtoolbox/cli/realtime_zshim.py +++ b/shimmingtoolbox/cli/realtime_zshim.py @@ -167,12 +167,13 @@ def realtime_zshim(fname_coil, fname_fmap, fname_mask, fname_resp, fname_json, f riro_correction[i_slice] = np.mean(nii_resampled_riro.get_fdata()[..., i_slice]) # Write to a text file + # TODO: Add as an option to output the file to a specified location fname_corrections = os.path.join(__dir_shimmingtoolbox__, 'zshim_gradients.txt') file_gradients = open(fname_corrections, 'w') for i_slice in range(n_slices): file_gradients.write(f'Vector_Gz[0][{i_slice}]= {static_correction[i_slice]:.6f}\n') file_gradients.write(f'Vector_Gz[1][{i_slice}]= {riro_correction[i_slice]:.12f}\n') - # Matlab includes the mean pressure + # Matlab includes the mean pressure file_gradients.close() # ================ PLOTS ================ @@ -267,7 +268,7 @@ def realtime_zshim(fname_coil, fname_fmap, fname_mask, fname_resp, fname_json, f # ax.plot(pmu_times / 1000, pmu.data, label='Raw pressures') ax.plot(pmu_times_within_range / 1000, pmu_data_within_range, label='Pmu pressures') ax.legend() - ax.set_title("Pressure [-2048, 2047] vs time (s) ") + ax.set_title("Pressure [0, 4095] vs time (s) ") ax = fig.add_subplot(212) ax.plot(acq_timestamps / 1000, fieldmap_avg, label='Mean B0') ax.legend() From 5bc832f53f81b82f6a7b0bcab2c5fe98bd1c495e Mon Sep 17 00:00:00 2001 From: Alexandre D'Astous Date: Wed, 4 Nov 2020 14:41:17 -0500 Subject: [PATCH 31/51] Add mean pressure in calculation to account if bellow moves --- shimmingtoolbox/cli/realtime_zshim.py | 13 +++++++------ 1 file changed, 7 insertions(+), 6 deletions(-) diff --git a/shimmingtoolbox/cli/realtime_zshim.py b/shimmingtoolbox/cli/realtime_zshim.py index d7ac02b7d..abe33ce04 100644 --- a/shimmingtoolbox/cli/realtime_zshim.py +++ b/shimmingtoolbox/cli/realtime_zshim.py @@ -95,8 +95,7 @@ def realtime_zshim(fname_coil, fname_fmap, fname_mask, fname_resp, fname_json, f # Calculate gz gradient # Image is z, y, x # Pixdim[3] is the space between pixels in the z direction in millimeters - # TODO: This is the gradient in the voxel reference frame meaning the axes could possibly not follow the scanner's - # axes (therefore the z gradient), also investigate if axes should be 1 or 0 + # TODO: Investigate if axes should be 1 or 0 g = 1000 / 42.576e6 # [mT / Hz] gz_gradient = np.zeros_like(masked_fieldmaps) for it in range(nt): @@ -117,7 +116,7 @@ def realtime_zshim(fname_coil, fname_fmap, fname_mask, fname_resp, fname_json, f # fit PMU and fieldmap values # do regression to separate static componant and RIRO component # output coefficient with proper scaling - # field(i_vox) = a(i_vox) * acq_pressures + b(i_vox) + # field(i_vox) = a(i_vox) * (acq_pressures - mean_p) + b(i_vox) # could use: https://scikit-learn.org/stable/modules/generated/sklearn.linear_model.LinearRegression.html # Note: strong spatial autocorrelation on the a and b coefficients. Ie: two adjacent voxels are submitted to similar # static B0 field and RIRO component. --> we need to find a way to account for that @@ -129,14 +128,15 @@ def realtime_zshim(fname_coil, fname_fmap, fname_mask, fname_resp, fname_json, f # cons: (from Ryan): regularized fitting took a lot of time on Matlab # Shim using PMU + mean_p = np.mean(acq_pressures) riro = np.zeros_like(fieldmap[:, :, :, 0]) static = np.zeros_like(fieldmap[:, :, :, 0]) for i_x in range(fieldmap.shape[0]): for i_y in range(fieldmap.shape[1]): for i_z in range(fieldmap.shape[2]): # TODO: Fit for -masked_field? - # reg = LinearRegression().fit(acq_pressures.reshape(-1, 1), -gz_gradient[i_x, i_y, i_z, :]) - reg = LinearRegression().fit(acq_pressures.reshape(-1, 1), -masked_fieldmaps[i_x, i_y, i_z, :]) + # reg = LinearRegression().fit(acq_pressures.reshape(-1, 1) - mean_p, -gz_gradient[i_x, i_y, i_z, :]) + reg = LinearRegression().fit(acq_pressures.reshape(-1, 1) - mean_p, -masked_fieldmaps[i_x, i_y, i_z, :]) riro[i_x, i_y, i_z] = reg.coef_ static[i_x, i_y, i_z] = reg.intercept_ @@ -173,6 +173,7 @@ def realtime_zshim(fname_coil, fname_fmap, fname_mask, fname_resp, fname_json, f for i_slice in range(n_slices): file_gradients.write(f'Vector_Gz[0][{i_slice}]= {static_correction[i_slice]:.6f}\n') file_gradients.write(f'Vector_Gz[1][{i_slice}]= {riro_correction[i_slice]:.12f}\n') + file_gradients.write(f'Vector_Gz[2][{i_slice}]= {mean_p:.3f}\n') # Matlab includes the mean pressure file_gradients.close() @@ -228,7 +229,7 @@ def realtime_zshim(fname_coil, fname_fmap, fname_mask, fname_resp, fname_json, f fig.savefig(fname_figure) # Calculate fitted and shimmed for pressure fitted plot - fitted_fieldmap = riro * acq_pressures + static + fitted_fieldmap = riro * (acq_pressures -mean_p) + static shimmed_pressure_fitted = np.expand_dims(fitted_fieldmap, 2) + masked_fieldmaps # Plot pressure fitted fieldmap From 0e9d27f7e9ffe4488f6d0b26c90287be473d6238 Mon Sep 17 00:00:00 2001 From: Alexandre D'Astous Date: Wed, 4 Nov 2020 14:55:45 -0500 Subject: [PATCH 32/51] Add local full acquisition and refine cube mask --- test/cli/test_cli_realtime_zshim.py | 22 ++++++++++++---------- 1 file changed, 12 insertions(+), 10 deletions(-) diff --git a/test/cli/test_cli_realtime_zshim.py b/test/cli/test_cli_realtime_zshim.py index c8f79b150..4893cb412 100644 --- a/test/cli/test_cli_realtime_zshim.py +++ b/test/cli/test_cli_realtime_zshim.py @@ -20,26 +20,28 @@ def test_cli_realtime_zshim(): with tempfile.TemporaryDirectory(prefix='st_' + pathlib.Path(__file__).stem) as tmp: runner = CliRunner() + # Local fname_fieldmap = os.path.join(__dir_testing__, 'realtime_zshimming_data', 'nifti', 'sub-example', 'fmap', 'sub-example_fieldmap.nii.gz') + # fname_fieldmap = os.path.join(__dir_testing__, 'test_realtime_zshim', 'sub-example_fieldmap.nii.gz') nii_fmap = nib.load(fname_fieldmap) fmap = nii_fmap.get_fdata() affine = nii_fmap.affine # Set up mask # Cube - # nx, ny, nz, nt = fmap.shape - # mask = shapes(fmap[:, :, :, 0], 'cube', - # center_dim1=int(fmap.shape[0] / 2 - 8), - # center_dim2=int(fmap.shape[1] / 2 - 5), - # len_dim1=15, len_dim2=25, len_dim3=nz) + nx, ny, nz, nt = fmap.shape + mask = shapes(fmap[:, :, :, 0], 'cube', + center_dim1=int(fmap.shape[0] / 2 - 8), + center_dim2=int(fmap.shape[1] / 2 - 5), + len_dim1=15, len_dim2=40, len_dim3=nz) # Threshold - fname_mag = os.path.join(__dir_testing__, 'realtime_zshimming_data', 'nifti', 'sub-example', 'fmap', - 'sub-example_magnitude1.nii.gz') - mag = nib.load(fname_mag).get_fdata() - mask = threshold(mag, thr=50) + # fname_mag = os.path.join(__dir_testing__, 'realtime_zshimming_data', 'nifti', 'sub-example', 'fmap', + # 'sub-example_magnitude1.nii.gz') + # mag = nib.load(fname_mag).get_fdata()[..., 0] + # mask = threshold(mag, thr=50) - nii_mask = nib.Nifti1Image(mask.astype(int)[:, :, :, 0], affine) + nii_mask = nib.Nifti1Image(mask.astype(int), affine) fname_mask = os.path.join(tmp, 'mask.nii.gz') nib.save(nii_mask, fname_mask) From f9ccc0ec30e2d829ec56174a2e6259e949cbe57a Mon Sep 17 00:00:00 2001 From: Alexandre D'Astous Date: Wed, 4 Nov 2020 15:03:02 -0500 Subject: [PATCH 33/51] Add plot of evolution of riro and static coefficients --- shimmingtoolbox/cli/realtime_zshim.py | 11 +++++++++++ 1 file changed, 11 insertions(+) diff --git a/shimmingtoolbox/cli/realtime_zshim.py b/shimmingtoolbox/cli/realtime_zshim.py index abe33ce04..760fbe7b0 100644 --- a/shimmingtoolbox/cli/realtime_zshim.py +++ b/shimmingtoolbox/cli/realtime_zshim.py @@ -299,6 +299,17 @@ def realtime_zshim(fname_coil, fname_fmap, fname_mask, fname_resp, fname_json, f fname_figure = os.path.join(__dir_shimmingtoolbox__, 'reatime_zshime_gradient.png') fig.savefig(fname_figure) + # Show evolution of coefficients + fig = Figure(figsize=(10, 10)) + ax = fig.add_subplot(2, 1, 1) + ax.plot(range(n_slices), static_correction, label='Static correction') + ax.set_title("Static correction evolution through slices") + ax = fig.add_subplot(2, 1, 2) + ax.plot(range(n_slices), riro_correction, label='Riro correction') + ax.set_title("Riro correction evolution through slices") + fname_figure = os.path.join(__dir_shimmingtoolbox__, 'reatime_zshime_correction_slice.png') + fig.savefig(fname_figure) + return fname_figure # Debug From 8072d2a12061b9ff87902ca406d3621d643ce02d Mon Sep 17 00:00:00 2001 From: Julien Cohen-Adad Date: Wed, 4 Nov 2020 15:29:34 -0500 Subject: [PATCH 34/51] Reorganized code, added comments --- shimmingtoolbox/cli/realtime_zshim.py | 17 ++++++++++------- 1 file changed, 10 insertions(+), 7 deletions(-) diff --git a/shimmingtoolbox/cli/realtime_zshim.py b/shimmingtoolbox/cli/realtime_zshim.py index 760fbe7b0..2e6cfc4ff 100644 --- a/shimmingtoolbox/cli/realtime_zshim.py +++ b/shimmingtoolbox/cli/realtime_zshim.py @@ -131,6 +131,7 @@ def realtime_zshim(fname_coil, fname_fmap, fname_mask, fname_resp, fname_json, f mean_p = np.mean(acq_pressures) riro = np.zeros_like(fieldmap[:, :, :, 0]) static = np.zeros_like(fieldmap[:, :, :, 0]) + # TODO add tqdm progress for i_x in range(fieldmap.shape[0]): for i_y in range(fieldmap.shape[1]): for i_z in range(fieldmap.shape[2]): @@ -140,23 +141,25 @@ def realtime_zshim(fname_coil, fname_fmap, fname_mask, fname_resp, fname_json, f riro[i_x, i_y, i_z] = reg.coef_ static[i_x, i_y, i_z] = reg.intercept_ - # Resample masked_fieldmaps, riro and static to target anatomical image + # Resample masked_fieldmaps to target anatomical image # TODO: convert to a function masked_fmap_4d = np.zeros(anat.shape + (nt,)) for it in range(nt): nii_masked_fmap_3d = nib.Nifti1Image(masked_fieldmaps[..., it], nii_fmap.affine) - nii_resampled_fmap_3d = resample_from_to(nii_masked_fmap_3d, nii_anat, mode='nearest') + nii_resampled_fmap_3d = resample_from_to(nii_masked_fmap_3d, nii_anat, order=2, mode='nearest') masked_fmap_4d[..., it] = nii_resampled_fmap_3d.get_fdata() - nii_resampled_fmap = nib.Nifti1Image(masked_fmap_4d, nii_anat.affine) - nii_riro = nib.Nifti1Image(riro, nii_fmap.affine) + nib.save(nii_resampled_fmap, os.path.join(__dir_shimmingtoolbox__, 'resampled_fmap.nii.gz')) + + # Resample static to target anatomical image nii_static = nib.Nifti1Image(static, nii_fmap.affine) - nii_resampled_riro = resample_from_to(nii_riro, nii_anat, mode='nearest') nii_resampled_static = resample_from_to(nii_static, nii_anat, mode='nearest') + nib.save(nii_resampled_static, os.path.join(__dir_shimmingtoolbox__, 'resampled_static.nii.gz')) - nib.save(nii_resampled_fmap, os.path.join(__dir_shimmingtoolbox__, 'resampled_fmap.nii.gz')) + # Resample riro to target anatomical image + nii_riro = nib.Nifti1Image(riro, nii_fmap.affine) + nii_resampled_riro = resample_from_to(nii_riro, nii_anat, mode='nearest') nib.save(nii_resampled_riro, os.path.join(__dir_shimmingtoolbox__, 'resampled_riro.nii.gz')) - nib.save(nii_resampled_static, os.path.join(__dir_shimmingtoolbox__, 'resampled_static.nii.gz')) # Calculate the mean for riro and static for a perticular slice n_slices = nii_resampled_fmap.get_fdata().shape[2] From 9f60b357717cd55d7d1efa0727d09eb880eb159e Mon Sep 17 00:00:00 2001 From: Julien Cohen-Adad Date: Wed, 4 Nov 2020 15:29:53 -0500 Subject: [PATCH 35/51] Now outputting Gz instead of fieldmap --- shimmingtoolbox/cli/realtime_zshim.py | 4 ++-- 1 file changed, 2 insertions(+), 2 deletions(-) diff --git a/shimmingtoolbox/cli/realtime_zshim.py b/shimmingtoolbox/cli/realtime_zshim.py index 2e6cfc4ff..3f8f02d03 100644 --- a/shimmingtoolbox/cli/realtime_zshim.py +++ b/shimmingtoolbox/cli/realtime_zshim.py @@ -136,8 +136,8 @@ def realtime_zshim(fname_coil, fname_fmap, fname_mask, fname_resp, fname_json, f for i_y in range(fieldmap.shape[1]): for i_z in range(fieldmap.shape[2]): # TODO: Fit for -masked_field? - # reg = LinearRegression().fit(acq_pressures.reshape(-1, 1) - mean_p, -gz_gradient[i_x, i_y, i_z, :]) - reg = LinearRegression().fit(acq_pressures.reshape(-1, 1) - mean_p, -masked_fieldmaps[i_x, i_y, i_z, :]) + reg = LinearRegression().fit(acq_pressures.reshape(-1, 1) - mean_p, -gz_gradient[i_x, i_y, i_z, :]) + # reg = LinearRegression().fit(acq_pressures.reshape(-1, 1) - mean_p, -masked_fieldmaps[i_x, i_y, i_z, :]) riro[i_x, i_y, i_z] = reg.coef_ static[i_x, i_y, i_z] = reg.intercept_ From a120c2616fda5127858a1b6d25bdcf660e6e37ef Mon Sep 17 00:00:00 2001 From: Julien Cohen-Adad Date: Wed, 4 Nov 2020 15:47:42 -0500 Subject: [PATCH 36/51] Added DEBUG global var --- shimmingtoolbox/cli/realtime_zshim.py | 4 ++++ 1 file changed, 4 insertions(+) diff --git a/shimmingtoolbox/cli/realtime_zshim.py b/shimmingtoolbox/cli/realtime_zshim.py index 3f8f02d03..c3cdad4ce 100644 --- a/shimmingtoolbox/cli/realtime_zshim.py +++ b/shimmingtoolbox/cli/realtime_zshim.py @@ -17,6 +17,7 @@ from shimmingtoolbox.pmu import PmuResp from shimmingtoolbox import __dir_shimmingtoolbox__ +DEBUG = True CONTEXT_SETTINGS = dict(help_option_names=['-h', '--help']) @@ -102,6 +103,9 @@ def realtime_zshim(fname_coil, fname_fmap, fname_mask, fname_resp, fname_json, f gz_gradient[..., 0, it] = np.gradient(g * masked_fieldmaps[:, :, 0, it], nii_fmap.header['pixdim'][3] / 1000, axis=1) # [mT / m] + if DEBUG: + nii_gz_gradient = nib.Nifti1Image(gz_gradient, nii_fmap.affine) + nib.save(nii_gz_gradient, os.path.join(__dir_shimmingtoolbox__, 'tmp.gz_gradient.nii.gz')) # Fetch PMU timing # TODO: Add json to fieldmap instead of asking for another json file From 48891b42ee4fa8427e4a898ae4447fc29d65c6f4 Mon Sep 17 00:00:00 2001 From: Julien Cohen-Adad Date: Wed, 4 Nov 2020 15:50:24 -0500 Subject: [PATCH 37/51] Added fig_ prefix to all output figures --- shimmingtoolbox/cli/realtime_zshim.py | 17 +++++++++-------- 1 file changed, 9 insertions(+), 8 deletions(-) diff --git a/shimmingtoolbox/cli/realtime_zshim.py b/shimmingtoolbox/cli/realtime_zshim.py index c3cdad4ce..8fb227c77 100644 --- a/shimmingtoolbox/cli/realtime_zshim.py +++ b/shimmingtoolbox/cli/realtime_zshim.py @@ -10,6 +10,7 @@ from nibabel.processing import resample_from_to # TODO: remove matplotlib and dirtesting import from matplotlib.figure import Figure +from tqdm import tqdm from shimmingtoolbox import __dir_testing__ from shimmingtoolbox.optimizer.sequential import sequential_zslice @@ -210,7 +211,7 @@ def realtime_zshim(fname_coil, fname_fmap, fname_mask, fname_resp, fname_json, f im = ax.imshow(shimmed[:-1, :-1, 0, i_t]) fig.colorbar(im) ax.set_title("Shimmed") - fname_figure = os.path.join(__dir_shimmingtoolbox__, 'realtime_zshim_sphharm_shimmed.png') + fname_figure = os.path.join(__dir_shimmingtoolbox__, 'fig_realtime_zshim_sphharm_shimmed.png') fig.savefig(fname_figure) # Plot the coil coefs through time @@ -219,7 +220,7 @@ def realtime_zshim(fname_coil, fname_fmap, fname_mask, fname_resp, fname_json, f ax = fig.add_subplot(n_coils, 1, i_coil + 1) ax.plot(np.arange(nt), currents[i_coil, :]) ax.set_title(f"Channel {i_coil}") - fname_figure = os.path.join(__dir_shimmingtoolbox__, 'realtime_zshim_sphharm_currents.png') + fname_figure = os.path.join(__dir_shimmingtoolbox__, 'fig_realtime_zshim_sphharm_currents.png') fig.savefig(fname_figure) # Plot Static and RIRO @@ -232,7 +233,7 @@ def realtime_zshim(fname_coil, fname_fmap, fname_mask, fname_resp, fname_json, f im = ax.imshow(static[:-1, :-1, 0]) fig.colorbar(im) ax.set_title("Static") - fname_figure = os.path.join(__dir_shimmingtoolbox__, 'realtime_zshim_riro_static.png') + fname_figure = os.path.join(__dir_shimmingtoolbox__, 'fig_realtime_zshim_riro_static.png') fig.savefig(fname_figure) # Calculate fitted and shimmed for pressure fitted plot @@ -253,7 +254,7 @@ def realtime_zshim(fname_coil, fname_fmap, fname_mask, fname_resp, fname_json, f im = ax.imshow(shimmed_pressure_fitted[:-1, :-1, 0, i_t]) fig.colorbar(im) ax.set_title("Shimmed (fit + fieldmap") - fname_figure = os.path.join(__dir_shimmingtoolbox__, 'realtime_zshim_pressure_fitted.png') + fname_figure = os.path.join(__dir_shimmingtoolbox__, 'fig_realtime_zshim_pressure_fitted.png') fig.savefig(fname_figure) # Reshape pmu datapoints to fit those of the acquisition @@ -281,7 +282,7 @@ def realtime_zshim(fname_coil, fname_fmap, fname_mask, fname_resp, fname_json, f ax.plot(acq_timestamps / 1000, fieldmap_avg, label='Mean B0') ax.legend() ax.set_title("Fieldmap average over unmasked region (Hz) vs time (s)") - fname_figure = os.path.join(__dir_shimmingtoolbox__, 'realtime_zshim_pmu_vs_B0.png') + fname_figure = os.path.join(__dir_shimmingtoolbox__, 'fig_realtime_zshim_pmu_vs_B0.png') fig.savefig(fname_figure) # Show anatomical image @@ -294,7 +295,7 @@ def realtime_zshim(fname_coil, fname_fmap, fname_mask, fname_resp, fname_json, f im = ax.imshow(nii_resampled_fmap.get_fdata()[:, :, 10, 0]) fig.colorbar(im) ax.set_title("Resampled fieldmap [:, :, 10, 0]") - fname_figure = os.path.join(__dir_shimmingtoolbox__, 'reatime_zshime_anat.png') + fname_figure = os.path.join(__dir_shimmingtoolbox__, 'fig_reatime_zshim_anat.png') fig.savefig(fname_figure) # Show Gradient @@ -303,7 +304,7 @@ def realtime_zshim(fname_coil, fname_fmap, fname_mask, fname_resp, fname_json, f im = ax.imshow(gz_gradient[:, :, 0, 0]) fig.colorbar(im) ax.set_title("Gradient [:, :, 0, 0]") - fname_figure = os.path.join(__dir_shimmingtoolbox__, 'reatime_zshime_gradient.png') + fname_figure = os.path.join(__dir_shimmingtoolbox__, 'fig_realtime_zshim_gradient.png') fig.savefig(fname_figure) # Show evolution of coefficients @@ -314,7 +315,7 @@ def realtime_zshim(fname_coil, fname_fmap, fname_mask, fname_resp, fname_json, f ax = fig.add_subplot(2, 1, 2) ax.plot(range(n_slices), riro_correction, label='Riro correction') ax.set_title("Riro correction evolution through slices") - fname_figure = os.path.join(__dir_shimmingtoolbox__, 'reatime_zshime_correction_slice.png') + fname_figure = os.path.join(__dir_shimmingtoolbox__, 'fig_realtime_zshim_correction_slice.png') fig.savefig(fname_figure) return fname_figure From c135a641f6a3f84d16f6f7f9dc248f19169ae003 Mon Sep 17 00:00:00 2001 From: Julien Cohen-Adad Date: Wed, 4 Nov 2020 16:08:32 -0500 Subject: [PATCH 38/51] Added progress bar to fitting Progress bar not showing though-- will need to fix that --- shimmingtoolbox/cli/realtime_zshim.py | 6 ++++-- 1 file changed, 4 insertions(+), 2 deletions(-) diff --git a/shimmingtoolbox/cli/realtime_zshim.py b/shimmingtoolbox/cli/realtime_zshim.py index 8fb227c77..3229b4174 100644 --- a/shimmingtoolbox/cli/realtime_zshim.py +++ b/shimmingtoolbox/cli/realtime_zshim.py @@ -10,13 +10,13 @@ from nibabel.processing import resample_from_to # TODO: remove matplotlib and dirtesting import from matplotlib.figure import Figure -from tqdm import tqdm from shimmingtoolbox import __dir_testing__ from shimmingtoolbox.optimizer.sequential import sequential_zslice from shimmingtoolbox.load_nifti import get_acquisition_times from shimmingtoolbox.pmu import PmuResp from shimmingtoolbox import __dir_shimmingtoolbox__ +from shimmingtoolbox.utils import st_progress_bar DEBUG = True CONTEXT_SETTINGS = dict(help_option_names=['-h', '--help']) @@ -136,7 +136,8 @@ def realtime_zshim(fname_coil, fname_fmap, fname_mask, fname_resp, fname_json, f mean_p = np.mean(acq_pressures) riro = np.zeros_like(fieldmap[:, :, :, 0]) static = np.zeros_like(fieldmap[:, :, :, 0]) - # TODO add tqdm progress + # TODO fix progress bar not showing up + progress_bar = st_progress_bar(fieldmap[..., 0].size, desc="Fitting", ascii=False) for i_x in range(fieldmap.shape[0]): for i_y in range(fieldmap.shape[1]): for i_z in range(fieldmap.shape[2]): @@ -145,6 +146,7 @@ def realtime_zshim(fname_coil, fname_fmap, fname_mask, fname_resp, fname_json, f # reg = LinearRegression().fit(acq_pressures.reshape(-1, 1) - mean_p, -masked_fieldmaps[i_x, i_y, i_z, :]) riro[i_x, i_y, i_z] = reg.coef_ static[i_x, i_y, i_z] = reg.intercept_ + progress_bar.update(1) # Resample masked_fieldmaps to target anatomical image # TODO: convert to a function From 17f4ed8f1caa07abaa750355dba794a093bbf8d0 Mon Sep 17 00:00:00 2001 From: Alexandre D'Astous Date: Fri, 6 Nov 2020 14:15:04 -0500 Subject: [PATCH 39/51] Apply mask on anat instead of the fieldmap --- shimmingtoolbox/cli/realtime_zshim.py | 41 +++++++++++++++++---------- test/cli/test_cli_realtime_zshim.py | 33 ++++++++++----------- 2 files changed, 41 insertions(+), 33 deletions(-) diff --git a/shimmingtoolbox/cli/realtime_zshim.py b/shimmingtoolbox/cli/realtime_zshim.py index 3229b4174..b5cb7e327 100644 --- a/shimmingtoolbox/cli/realtime_zshim.py +++ b/shimmingtoolbox/cli/realtime_zshim.py @@ -32,8 +32,9 @@ "Siemens gradient/shim coils and external custom coils).") @click.option('-fmap', 'fname_fmap', required=True, type=click.Path(), help="B0 fieldmap. For realtime shimming, this should be a 4d file (4th dimension being time") -@click.option('-mask', 'fname_mask', type=click.Path(), - help="3D nifti file with voxels between 0 and 1 used to weight the spatial region to shim.") +@click.option('-mask', 'fname_mask_anat', type=click.Path(), + help="3D nifti file with voxels between 0 and 1 used to weight the spatial region to shim. " + "The coordinate system should be the same as ``anat``'s coordinate system.") @click.option('-resp', 'fname_resp', type=click.Path(), help="Siemens respiratory file containing pressure data.") @click.option('-anat', 'fname_anat', type=click.Path(), @@ -42,13 +43,13 @@ @click.option('-json', 'fname_json', type=click.Path(), help="Filename of json corresponding BIDS sidecar.") @click.option("-verbose", is_flag=True, help="Be more verbose.") -def realtime_zshim(fname_coil, fname_fmap, fname_mask, fname_resp, fname_json, fname_anat, verbose=True): +def realtime_zshim(fname_coil, fname_fmap, fname_mask_anat, fname_resp, fname_json, fname_anat, verbose=True): """ Args: fname_coil: Pointing to coil profile. 4-dimensional: x, y, z, coil. fname_fmap: - fname_mask: + fname_mask_anat: fname_resp: verbose: @@ -70,29 +71,39 @@ def realtime_zshim(fname_coil, fname_fmap, fname_mask, fname_resp, fname_json, f raise RuntimeError("fmap must be 4d (x, y, z, t)") nx, ny, nz, nt = fieldmap.shape - # Load mask - # TODO: check good practice below - if fname_mask is not None: - mask = nib.load(fname_mask).get_fdata() - else: - mask = np.ones_like(fieldmap) - # Load anat nii_anat = nib.load(fname_anat) anat = nii_anat.get_fdata() if anat.ndim != 3: raise RuntimeError("Anatomical image must be in 3d") + # Load mask + # TODO: check good practice below + if fname_mask_anat is not None: + nii_mask_anat = nib.load(fname_mask_anat) + if not np.all(nii_anat.affine == nii_mask_anat.affine) or not np.all(nii_mask_anat.shape == nii_anat.shape): + raise RuntimeError("Mask must have the same shape and affine transformation as anat") + nii_fmap_3d_temp = nib.Nifti1Image(fieldmap[..., 0], nii_fmap.affine) + nii_mask_fmap = resample_from_to(nii_mask_anat, nii_fmap_3d_temp, mode='nearest') + mask_fmap = nii_mask_fmap.get_fdata() + else: + mask_fmap = np.ones_like(fieldmap) + nii_mask_fmap = nib.Nifti1Image(mask_fmap, nii_anat.affine) + nii_mask_anat = nib.Nifti1Image(np.ones_like(anat), nii_anat.affine) + + if DEBUG: + nib.save(nii_mask_fmap, os.path.join(__dir_shimmingtoolbox__, 'tmp.mask_fmap_resample.nii.gz')) + # Shim using sequencer and optimizer n_coils = coil.shape[-1] currents = np.zeros([n_coils, nt]) shimmed = np.zeros_like(fieldmap) masked_fieldmaps = np.zeros_like(fieldmap) for i_t in range(nt): - currents[:, i_t] = sequential_zslice(fieldmap[..., i_t], coil, mask, z_slices=np.array(range(nz)), + currents[:, i_t] = sequential_zslice(fieldmap[..., i_t], coil, mask_fmap, z_slices=np.array(range(nz)), bounds=[(-np.inf, np.inf)] * n_coils) shimmed[..., i_t] = fieldmap[..., i_t] + np.sum(currents[:, i_t] * coil, axis=3, keepdims=False) - masked_fieldmaps[..., i_t] = mask * fieldmap[..., i_t] + masked_fieldmaps[..., i_t] = mask_fmap * fieldmap[..., i_t] # Calculate gz gradient # Image is z, y, x @@ -192,7 +203,7 @@ def realtime_zshim(fname_coil, fname_fmap, fname_mask, fname_resp, fname_json, f # Calculate masked shim for spherical harmonics plot masked_shimmed = np.zeros_like(shimmed) for i_t in range(nt): - masked_shimmed[..., i_t] = mask * shimmed[..., i_t] + masked_shimmed[..., i_t] = mask_fmap * shimmed[..., i_t] # Plot unshimmed vs shimmed and their mask for spherical harmonics i_t = 0 @@ -269,7 +280,7 @@ def realtime_zshim(fname_coil, fname_fmap, fname_mask, fname_resp, fname_json, f # Calc fieldmap average within mask fieldmap_avg = np.zeros([fieldmap.shape[3]]) for i_time in range(nt): - masked_array = np.ma.array(fieldmap[:, :, :, i_time], mask=mask == False) + masked_array = np.ma.array(fieldmap[:, :, :, i_time], mask=mask_fmap == False) fieldmap_avg[i_time] = np.ma.average(masked_array) # Plot pmu vs B0 in masked region diff --git a/test/cli/test_cli_realtime_zshim.py b/test/cli/test_cli_realtime_zshim.py index 4893cb412..88d95e331 100644 --- a/test/cli/test_cli_realtime_zshim.py +++ b/test/cli/test_cli_realtime_zshim.py @@ -25,31 +25,32 @@ def test_cli_realtime_zshim(): 'sub-example_fieldmap.nii.gz') # fname_fieldmap = os.path.join(__dir_testing__, 'test_realtime_zshim', 'sub-example_fieldmap.nii.gz') nii_fmap = nib.load(fname_fieldmap) - fmap = nii_fmap.get_fdata() - affine = nii_fmap.affine + + # Path for mag anat image + fname_anat = os.path.join(__dir_testing__, 'realtime_zshimming_data', 'nifti', 'sub-example', 'anat', + 'sub-example_unshimmed_e1.nii.gz') + nii_anat = nib.load(fname_anat) + anat = nii_anat.get_fdata() # Set up mask # Cube - nx, ny, nz, nt = fmap.shape - mask = shapes(fmap[:, :, :, 0], 'cube', - center_dim1=int(fmap.shape[0] / 2 - 8), - center_dim2=int(fmap.shape[1] / 2 - 5), - len_dim1=15, len_dim2=40, len_dim3=nz) + nx, ny, nz = anat.shape + mask = shapes(anat, 'cube', + center_dim1=int(nx / 2), + center_dim2=int(ny / 2), + len_dim1=30, len_dim2=30, len_dim3=nz) # Threshold - # fname_mag = os.path.join(__dir_testing__, 'realtime_zshimming_data', 'nifti', 'sub-example', 'fmap', - # 'sub-example_magnitude1.nii.gz') - # mag = nib.load(fname_mag).get_fdata()[..., 0] - # mask = threshold(mag, thr=50) + # mask = threshold(anat, thr=50) - nii_mask = nib.Nifti1Image(mask.astype(int), affine) + nii_mask = nib.Nifti1Image(mask.astype(int), nii_anat.affine) fname_mask = os.path.join(tmp, 'mask.nii.gz') nib.save(nii_mask, fname_mask) # Set up coils - coord_phys = generate_meshgrid(fmap.shape[0:3], affine) + coord_phys = generate_meshgrid(nii_fmap.get_fdata().shape[0:3], nii_fmap.affine) coil_profile = siemens_basis(coord_phys[0], coord_phys[1], coord_phys[2]) - nii_coil = nib.Nifti1Image(coil_profile, affine) + nii_coil = nib.Nifti1Image(coil_profile, nii_fmap.affine) fname_coil = os.path.join(tmp, 'coil_profile.nii.gz') nib.save(nii_coil, fname_coil) @@ -60,10 +61,6 @@ def test_cli_realtime_zshim(): fname_json = os.path.join(__dir_testing__, 'realtime_zshimming_data', 'nifti', 'sub-example', 'fmap', 'sub-example_magnitude1.json') - # Path for mag anat image - fname_anat = os.path.join(__dir_testing__, 'realtime_zshimming_data', 'nifti', 'sub-example', 'anat', - 'sub-example_unshimmed_e1.nii.gz') - result = runner.invoke(realtime_zshim, ['-fmap', fname_fieldmap, '-coil', fname_coil, '-mask', fname_mask, '-resp', fname_resp, '-json', fname_json, '-anat', fname_anat], catch_exceptions=False) From 18affd478ef80183655853c886941bee64ac154f Mon Sep 17 00:00:00 2001 From: Alexandre D'Astous Date: Fri, 6 Nov 2020 14:49:04 -0500 Subject: [PATCH 40/51] Commented out unecessary code --- shimmingtoolbox/cli/realtime_zshim.py | 98 ++++++++++++--------------- test/cli/test_cli_realtime_zshim.py | 14 ++-- 2 files changed, 50 insertions(+), 62 deletions(-) diff --git a/shimmingtoolbox/cli/realtime_zshim.py b/shimmingtoolbox/cli/realtime_zshim.py index b5cb7e327..cc88f0295 100644 --- a/shimmingtoolbox/cli/realtime_zshim.py +++ b/shimmingtoolbox/cli/realtime_zshim.py @@ -26,10 +26,10 @@ context_settings=CONTEXT_SETTINGS, help=f"Perform realtime z-shimming." ) -@click.option('-coil', 'fname_coil', required=True, type=click.Path(), - help="Coil basis to use for shimming. Enter multiple files if " - "you wish to use more than one set of shim coils (eg: " - "Siemens gradient/shim coils and external custom coils).") +# @click.option('-coil', 'fname_coil', required=True, type=click.Path(), +# help="Coil basis to use for shimming. Enter multiple files if " +# "you wish to use more than one set of shim coils (eg: " +# "Siemens gradient/shim coils and external custom coils).") @click.option('-fmap', 'fname_fmap', required=True, type=click.Path(), help="B0 fieldmap. For realtime shimming, this should be a 4d file (4th dimension being time") @click.option('-mask', 'fname_mask_anat', type=click.Path(), @@ -43,7 +43,7 @@ @click.option('-json', 'fname_json', type=click.Path(), help="Filename of json corresponding BIDS sidecar.") @click.option("-verbose", is_flag=True, help="Be more verbose.") -def realtime_zshim(fname_coil, fname_fmap, fname_mask_anat, fname_resp, fname_json, fname_anat, verbose=True): +def realtime_zshim(fname_fmap, fname_mask_anat, fname_resp, fname_json, fname_anat, verbose=True): """ Args: @@ -60,7 +60,7 @@ def realtime_zshim(fname_coil, fname_fmap, fname_mask_anat, fname_resp, fname_js # When using only z channel (corresponding to index 0) TODO:Remove # coil = np.expand_dims(nib.load(fname_coil).get_fdata()[:, :, :, 0], -1) # When using all channels TODO: Keep - coil = nib.load(fname_coil).get_fdata() + # coil = nib.load(fname_coil).get_fdata() # Load fieldmap nii_fmap = nib.load(fname_fmap) @@ -84,7 +84,7 @@ def realtime_zshim(fname_coil, fname_fmap, fname_mask_anat, fname_resp, fname_js if not np.all(nii_anat.affine == nii_mask_anat.affine) or not np.all(nii_mask_anat.shape == nii_anat.shape): raise RuntimeError("Mask must have the same shape and affine transformation as anat") nii_fmap_3d_temp = nib.Nifti1Image(fieldmap[..., 0], nii_fmap.affine) - nii_mask_fmap = resample_from_to(nii_mask_anat, nii_fmap_3d_temp, mode='nearest') + nii_mask_fmap = resample_from_to(nii_mask_anat, nii_fmap_3d_temp) mask_fmap = nii_mask_fmap.get_fdata() else: mask_fmap = np.ones_like(fieldmap) @@ -95,20 +95,19 @@ def realtime_zshim(fname_coil, fname_fmap, fname_mask_anat, fname_resp, fname_js nib.save(nii_mask_fmap, os.path.join(__dir_shimmingtoolbox__, 'tmp.mask_fmap_resample.nii.gz')) # Shim using sequencer and optimizer - n_coils = coil.shape[-1] - currents = np.zeros([n_coils, nt]) - shimmed = np.zeros_like(fieldmap) + # n_coils = coil.shape[-1] + # currents = np.zeros([n_coils, nt]) + # shimmed = np.zeros_like(fieldmap) masked_fieldmaps = np.zeros_like(fieldmap) for i_t in range(nt): - currents[:, i_t] = sequential_zslice(fieldmap[..., i_t], coil, mask_fmap, z_slices=np.array(range(nz)), - bounds=[(-np.inf, np.inf)] * n_coils) - shimmed[..., i_t] = fieldmap[..., i_t] + np.sum(currents[:, i_t] * coil, axis=3, keepdims=False) + # currents[:, i_t] = sequential_zslice(fieldmap[..., i_t], coil, mask_fmap, z_slices=np.array(range(nz)), + # bounds=[(-np.inf, np.inf)] * n_coils) + # shimmed[..., i_t] = fieldmap[..., i_t] + np.sum(currents[:, i_t] * coil, axis=3, keepdims=False) masked_fieldmaps[..., i_t] = mask_fmap * fieldmap[..., i_t] # Calculate gz gradient # Image is z, y, x # Pixdim[3] is the space between pixels in the z direction in millimeters - # TODO: Investigate if axes should be 1 or 0 g = 1000 / 42.576e6 # [mT / Hz] gz_gradient = np.zeros_like(masked_fieldmaps) for it in range(nt): @@ -201,40 +200,40 @@ def realtime_zshim(fname_coil, fname_fmap, fname_mask_anat, fname_resp, fname_js # ================ PLOTS ================ # Calculate masked shim for spherical harmonics plot - masked_shimmed = np.zeros_like(shimmed) - for i_t in range(nt): - masked_shimmed[..., i_t] = mask_fmap * shimmed[..., i_t] + # masked_shimmed = np.zeros_like(shimmed) + # for i_t in range(nt): + # masked_shimmed[..., i_t] = mask_fmap * shimmed[..., i_t] # Plot unshimmed vs shimmed and their mask for spherical harmonics - i_t = 0 - fig = Figure(figsize=(10, 10)) - ax = fig.add_subplot(2, 2, 1) - im = ax.imshow(masked_fieldmaps[:-1, :-1, 0, i_t]) - fig.colorbar(im) - ax.set_title("Masked unshimmed") - ax = fig.add_subplot(2, 2, 2) - im = ax.imshow(masked_shimmed[:-1, :-1, 0, i_t]) - fig.colorbar(im) - ax.set_title("Masked shimmed") - ax = fig.add_subplot(2, 2, 3) - im = ax.imshow(fieldmap[:-1, :-1, 0, i_t]) - fig.colorbar(im) - ax.set_title("Unshimmed") - ax = fig.add_subplot(2, 2, 4) - im = ax.imshow(shimmed[:-1, :-1, 0, i_t]) - fig.colorbar(im) - ax.set_title("Shimmed") - fname_figure = os.path.join(__dir_shimmingtoolbox__, 'fig_realtime_zshim_sphharm_shimmed.png') - fig.savefig(fname_figure) + # i_t = 0 + # fig = Figure(figsize=(10, 10)) + # ax = fig.add_subplot(2, 2, 1) + # im = ax.imshow(masked_fieldmaps[:-1, :-1, 0, i_t]) + # fig.colorbar(im) + # ax.set_title("Masked unshimmed") + # ax = fig.add_subplot(2, 2, 2) + # im = ax.imshow(masked_shimmed[:-1, :-1, 0, i_t]) + # fig.colorbar(im) + # ax.set_title("Masked shimmed") + # ax = fig.add_subplot(2, 2, 3) + # im = ax.imshow(fieldmap[:-1, :-1, 0, i_t]) + # fig.colorbar(im) + # ax.set_title("Unshimmed") + # ax = fig.add_subplot(2, 2, 4) + # im = ax.imshow(shimmed[:-1, :-1, 0, i_t]) + # fig.colorbar(im) + # ax.set_title("Shimmed") + # fname_figure = os.path.join(__dir_shimmingtoolbox__, 'fig_realtime_zshim_sphharm_shimmed.png') + # fig.savefig(fname_figure) # Plot the coil coefs through time - fig = Figure(figsize=(10, 10)) - for i_coil in range(n_coils): - ax = fig.add_subplot(n_coils, 1, i_coil + 1) - ax.plot(np.arange(nt), currents[i_coil, :]) - ax.set_title(f"Channel {i_coil}") - fname_figure = os.path.join(__dir_shimmingtoolbox__, 'fig_realtime_zshim_sphharm_currents.png') - fig.savefig(fname_figure) + # fig = Figure(figsize=(10, 10)) + # for i_coil in range(n_coils): + # ax = fig.add_subplot(n_coils, 1, i_coil + 1) + # ax.plot(np.arange(nt), currents[i_coil, :]) + # ax.set_title(f"Channel {i_coil}") + # fname_figure = os.path.join(__dir_shimmingtoolbox__, 'fig_realtime_zshim_sphharm_currents.png') + # fig.savefig(fname_figure) # Plot Static and RIRO fig = Figure(figsize=(10, 10)) @@ -332,14 +331,3 @@ def realtime_zshim(fname_coil, fname_fmap, fname_mask_anat, fname_resp, fname_js fig.savefig(fname_figure) return fname_figure - -# Debug -# fname_coil = os.path.join(__dir_testing__, 'test_realtime_zshim', 'coil_profile.nii.gz') -# fname_fmap = os.path.join(__dir_testing__, 'test_realtime_zshim', 'sub-example_fieldmap.nii.gz') -# fname_mask = os.path.join(__dir_testing__, 'test_realtime_zshim', 'mask.nii.gz') -# fname_resp = os.path.join(__dir_testing__, 'realtime_zshimming_data', 'PMUresp_signal.resp') -# fname_json = os.path.join(__dir_testing__, 'test_realtime_zshim', 'sub-example_magnitude1.json') -# # fname_coil='/Users/julien/code/shimming-toolbox/shimming-toolbox/test_realtime_zshim/coil_profile.nii.gz' -# # fname_fmap='/Users/julien/code/shimming-toolbox/shimming-toolbox/test_realtime_zshim/sub-example_fieldmap.nii.gz' -# # fname_mask='/Users/julien/code/shimming-toolbox/shimming-toolbox/test_realtime_zshim/mask.nii.gz' -# realtime_zshim(fname_coil, fname_fmap, fname_mask, fname_resp, fname_json) diff --git a/test/cli/test_cli_realtime_zshim.py b/test/cli/test_cli_realtime_zshim.py index 88d95e331..0f0d45357 100644 --- a/test/cli/test_cli_realtime_zshim.py +++ b/test/cli/test_cli_realtime_zshim.py @@ -47,12 +47,12 @@ def test_cli_realtime_zshim(): nib.save(nii_mask, fname_mask) # Set up coils - coord_phys = generate_meshgrid(nii_fmap.get_fdata().shape[0:3], nii_fmap.affine) - coil_profile = siemens_basis(coord_phys[0], coord_phys[1], coord_phys[2]) - - nii_coil = nib.Nifti1Image(coil_profile, nii_fmap.affine) - fname_coil = os.path.join(tmp, 'coil_profile.nii.gz') - nib.save(nii_coil, fname_coil) + # coord_phys = generate_meshgrid(nii_fmap.get_fdata().shape[0:3], nii_fmap.affine) + # coil_profile = siemens_basis(coord_phys[0], coord_phys[1], coord_phys[2]) + # + # nii_coil = nib.Nifti1Image(coil_profile, nii_fmap.affine) + # fname_coil = os.path.join(tmp, 'coil_profile.nii.gz') + # nib.save(nii_coil, fname_coil) # Path for resp data fname_resp = os.path.join(__dir_testing__, 'realtime_zshimming_data', 'PMUresp_signal.resp') @@ -61,7 +61,7 @@ def test_cli_realtime_zshim(): fname_json = os.path.join(__dir_testing__, 'realtime_zshimming_data', 'nifti', 'sub-example', 'fmap', 'sub-example_magnitude1.json') - result = runner.invoke(realtime_zshim, ['-fmap', fname_fieldmap, '-coil', fname_coil, '-mask', fname_mask, + result = runner.invoke(realtime_zshim, ['-fmap', fname_fieldmap, '-mask', fname_mask, '-resp', fname_resp, '-json', fname_json, '-anat', fname_anat], catch_exceptions=False) From fb75c623b4062eb987929210feb6c407e4403558 Mon Sep 17 00:00:00 2001 From: Alexandre D'Astous Date: Fri, 6 Nov 2020 16:02:09 -0500 Subject: [PATCH 41/51] Improve masking and mean calculation --- shimmingtoolbox/cli/realtime_zshim.py | 24 ++++++++++++++++-------- 1 file changed, 16 insertions(+), 8 deletions(-) diff --git a/shimmingtoolbox/cli/realtime_zshim.py b/shimmingtoolbox/cli/realtime_zshim.py index cc88f0295..bd5f98b5d 100644 --- a/shimmingtoolbox/cli/realtime_zshim.py +++ b/shimmingtoolbox/cli/realtime_zshim.py @@ -109,9 +109,9 @@ def realtime_zshim(fname_fmap, fname_mask_anat, fname_resp, fname_json, fname_an # Image is z, y, x # Pixdim[3] is the space between pixels in the z direction in millimeters g = 1000 / 42.576e6 # [mT / Hz] - gz_gradient = np.zeros_like(masked_fieldmaps) + gz_gradient = np.zeros_like(fieldmap) for it in range(nt): - gz_gradient[..., 0, it] = np.gradient(g * masked_fieldmaps[:, :, 0, it], + gz_gradient[..., 0, it] = np.gradient(g * fieldmap[:, :, 0, it], nii_fmap.header['pixdim'][3] / 1000, axis=1) # [mT / m] if DEBUG: @@ -171,20 +171,28 @@ def realtime_zshim(fname_fmap, fname_mask_anat, fname_resp, fname_json, fname_an # Resample static to target anatomical image nii_static = nib.Nifti1Image(static, nii_fmap.affine) nii_resampled_static = resample_from_to(nii_static, nii_anat, mode='nearest') - nib.save(nii_resampled_static, os.path.join(__dir_shimmingtoolbox__, 'resampled_static.nii.gz')) + nii_resampled_static_masked = nib.Nifti1Image(nii_resampled_static.get_fdata() * nii_mask_anat.get_fdata(), + nii_resampled_static.affine) + nib.save(nii_resampled_static_masked, os.path.join(__dir_shimmingtoolbox__, 'resampled_static.nii.gz')) # Resample riro to target anatomical image nii_riro = nib.Nifti1Image(riro, nii_fmap.affine) nii_resampled_riro = resample_from_to(nii_riro, nii_anat, mode='nearest') - nib.save(nii_resampled_riro, os.path.join(__dir_shimmingtoolbox__, 'resampled_riro.nii.gz')) + nii_resampled_riro_masked = nib.Nifti1Image(nii_resampled_riro.get_fdata() * nii_mask_anat.get_fdata(), + nii_resampled_riro.affine) + nib.save(nii_resampled_riro_masked, os.path.join(__dir_shimmingtoolbox__, 'resampled_riro.nii.gz')) # Calculate the mean for riro and static for a perticular slice - n_slices = nii_resampled_fmap.get_fdata().shape[2] + n_slices = nii_anat.get_fdata().shape[2] static_correction = np.zeros([n_slices]) riro_correction = np.zeros([n_slices]) for i_slice in range(n_slices): - static_correction[i_slice] = np.mean(nii_resampled_static.get_fdata()[..., i_slice]) - riro_correction[i_slice] = np.mean(nii_resampled_riro.get_fdata()[..., i_slice]) + ma_static_anat = np.ma.array(nii_resampled_static.get_fdata()[..., i_slice], + mask=nii_mask_anat.get_fdata()[..., i_slice] == False) + static_correction[i_slice] = np.ma.mean(ma_static_anat) + ma_riro_anat = np.ma.array(nii_resampled_riro.get_fdata()[..., i_slice], + mask=nii_mask_anat.get_fdata()[..., i_slice] == False) + riro_correction[i_slice] = np.ma.mean(ma_riro_anat) # Write to a text file # TODO: Add as an option to output the file to a specified location @@ -249,7 +257,7 @@ def realtime_zshim(fname_fmap, fname_mask_anat, fname_resp, fname_json, fname_an fig.savefig(fname_figure) # Calculate fitted and shimmed for pressure fitted plot - fitted_fieldmap = riro * (acq_pressures -mean_p) + static + fitted_fieldmap = riro * (acq_pressures - mean_p) + static shimmed_pressure_fitted = np.expand_dims(fitted_fieldmap, 2) + masked_fieldmaps # Plot pressure fitted fieldmap From a909449bb134eed5a4ccaec4aea694f50e752d5c Mon Sep 17 00:00:00 2001 From: Alexandre D'Astous Date: Fri, 6 Nov 2020 16:17:06 -0500 Subject: [PATCH 42/51] add outpout option to output figures and text file --- shimmingtoolbox/cli/realtime_zshim.py | 39 ++++++++++++++++----------- test/cli/test_cli_realtime_zshim.py | 13 +++++++-- 2 files changed, 34 insertions(+), 18 deletions(-) diff --git a/shimmingtoolbox/cli/realtime_zshim.py b/shimmingtoolbox/cli/realtime_zshim.py index bd5f98b5d..1342dc668 100644 --- a/shimmingtoolbox/cli/realtime_zshim.py +++ b/shimmingtoolbox/cli/realtime_zshim.py @@ -39,18 +39,22 @@ help="Siemens respiratory file containing pressure data.") @click.option('-anat', 'fname_anat', type=click.Path(), help="Filename of the anatomical image to apply the correction.") +@click.option('-output', 'fname_output', type=click.Path(), + help="Directory to output gradient text file and figures") # TODO: Remove json file as input @click.option('-json', 'fname_json', type=click.Path(), help="Filename of json corresponding BIDS sidecar.") @click.option("-verbose", is_flag=True, help="Be more verbose.") -def realtime_zshim(fname_fmap, fname_mask_anat, fname_resp, fname_json, fname_anat, verbose=True): +def realtime_zshim(fname_fmap, fname_mask_anat, fname_resp, fname_json, fname_anat, fname_output, verbose=True): """ Args: - fname_coil: Pointing to coil profile. 4-dimensional: x, y, z, coil. fname_fmap: fname_mask_anat: fname_resp: + fname_json: + fname_anat: + fname_output verbose: Returns: @@ -62,6 +66,10 @@ def realtime_zshim(fname_fmap, fname_mask_anat, fname_resp, fname_json, fname_an # When using all channels TODO: Keep # coil = nib.load(fname_coil).get_fdata() + # Look if output directory exists, if not, create it + if not os.path.exists(fname_output): + os.makedirs(fname_output) + # Load fieldmap nii_fmap = nib.load(fname_fmap) fieldmap = nii_fmap.get_fdata() @@ -92,7 +100,7 @@ def realtime_zshim(fname_fmap, fname_mask_anat, fname_resp, fname_json, fname_an nii_mask_anat = nib.Nifti1Image(np.ones_like(anat), nii_anat.affine) if DEBUG: - nib.save(nii_mask_fmap, os.path.join(__dir_shimmingtoolbox__, 'tmp.mask_fmap_resample.nii.gz')) + nib.save(nii_mask_fmap, os.path.join(fname_output, 'tmp.mask_fmap_resample.nii.gz')) # Shim using sequencer and optimizer # n_coils = coil.shape[-1] @@ -116,7 +124,7 @@ def realtime_zshim(fname_fmap, fname_mask_anat, fname_resp, fname_json, fname_an axis=1) # [mT / m] if DEBUG: nii_gz_gradient = nib.Nifti1Image(gz_gradient, nii_fmap.affine) - nib.save(nii_gz_gradient, os.path.join(__dir_shimmingtoolbox__, 'tmp.gz_gradient.nii.gz')) + nib.save(nii_gz_gradient, os.path.join(fname_output, 'tmp.gz_gradient.nii.gz')) # Fetch PMU timing # TODO: Add json to fieldmap instead of asking for another json file @@ -166,21 +174,21 @@ def realtime_zshim(fname_fmap, fname_mask_anat, fname_resp, fname_json, fname_an nii_resampled_fmap_3d = resample_from_to(nii_masked_fmap_3d, nii_anat, order=2, mode='nearest') masked_fmap_4d[..., it] = nii_resampled_fmap_3d.get_fdata() nii_resampled_fmap = nib.Nifti1Image(masked_fmap_4d, nii_anat.affine) - nib.save(nii_resampled_fmap, os.path.join(__dir_shimmingtoolbox__, 'resampled_fmap.nii.gz')) + nib.save(nii_resampled_fmap, os.path.join(fname_output, 'resampled_fmap.nii.gz')) # Resample static to target anatomical image nii_static = nib.Nifti1Image(static, nii_fmap.affine) nii_resampled_static = resample_from_to(nii_static, nii_anat, mode='nearest') nii_resampled_static_masked = nib.Nifti1Image(nii_resampled_static.get_fdata() * nii_mask_anat.get_fdata(), nii_resampled_static.affine) - nib.save(nii_resampled_static_masked, os.path.join(__dir_shimmingtoolbox__, 'resampled_static.nii.gz')) + nib.save(nii_resampled_static_masked, os.path.join(fname_output, 'resampled_static.nii.gz')) # Resample riro to target anatomical image nii_riro = nib.Nifti1Image(riro, nii_fmap.affine) nii_resampled_riro = resample_from_to(nii_riro, nii_anat, mode='nearest') nii_resampled_riro_masked = nib.Nifti1Image(nii_resampled_riro.get_fdata() * nii_mask_anat.get_fdata(), nii_resampled_riro.affine) - nib.save(nii_resampled_riro_masked, os.path.join(__dir_shimmingtoolbox__, 'resampled_riro.nii.gz')) + nib.save(nii_resampled_riro_masked, os.path.join(fname_output, 'resampled_riro.nii.gz')) # Calculate the mean for riro and static for a perticular slice n_slices = nii_anat.get_fdata().shape[2] @@ -196,13 +204,12 @@ def realtime_zshim(fname_fmap, fname_mask_anat, fname_resp, fname_json, fname_an # Write to a text file # TODO: Add as an option to output the file to a specified location - fname_corrections = os.path.join(__dir_shimmingtoolbox__, 'zshim_gradients.txt') + fname_corrections = os.path.join(fname_output, 'zshim_gradients.txt') file_gradients = open(fname_corrections, 'w') for i_slice in range(n_slices): file_gradients.write(f'Vector_Gz[0][{i_slice}]= {static_correction[i_slice]:.6f}\n') file_gradients.write(f'Vector_Gz[1][{i_slice}]= {riro_correction[i_slice]:.12f}\n') file_gradients.write(f'Vector_Gz[2][{i_slice}]= {mean_p:.3f}\n') - # Matlab includes the mean pressure file_gradients.close() # ================ PLOTS ================ @@ -253,7 +260,7 @@ def realtime_zshim(fname_fmap, fname_mask_anat, fname_resp, fname_json, fname_an im = ax.imshow(static[:-1, :-1, 0]) fig.colorbar(im) ax.set_title("Static") - fname_figure = os.path.join(__dir_shimmingtoolbox__, 'fig_realtime_zshim_riro_static.png') + fname_figure = os.path.join(fname_output, 'fig_realtime_zshim_riro_static.png') fig.savefig(fname_figure) # Calculate fitted and shimmed for pressure fitted plot @@ -274,7 +281,7 @@ def realtime_zshim(fname_fmap, fname_mask_anat, fname_resp, fname_json, fname_an im = ax.imshow(shimmed_pressure_fitted[:-1, :-1, 0, i_t]) fig.colorbar(im) ax.set_title("Shimmed (fit + fieldmap") - fname_figure = os.path.join(__dir_shimmingtoolbox__, 'fig_realtime_zshim_pressure_fitted.png') + fname_figure = os.path.join(fname_output, 'fig_realtime_zshim_pressure_fitted.png') fig.savefig(fname_figure) # Reshape pmu datapoints to fit those of the acquisition @@ -302,7 +309,7 @@ def realtime_zshim(fname_fmap, fname_mask_anat, fname_resp, fname_json, fname_an ax.plot(acq_timestamps / 1000, fieldmap_avg, label='Mean B0') ax.legend() ax.set_title("Fieldmap average over unmasked region (Hz) vs time (s)") - fname_figure = os.path.join(__dir_shimmingtoolbox__, 'fig_realtime_zshim_pmu_vs_B0.png') + fname_figure = os.path.join(fname_output, 'fig_realtime_zshim_pmu_vs_B0.png') fig.savefig(fname_figure) # Show anatomical image @@ -315,7 +322,7 @@ def realtime_zshim(fname_fmap, fname_mask_anat, fname_resp, fname_json, fname_an im = ax.imshow(nii_resampled_fmap.get_fdata()[:, :, 10, 0]) fig.colorbar(im) ax.set_title("Resampled fieldmap [:, :, 10, 0]") - fname_figure = os.path.join(__dir_shimmingtoolbox__, 'fig_reatime_zshim_anat.png') + fname_figure = os.path.join(fname_output, 'fig_reatime_zshim_anat.png') fig.savefig(fname_figure) # Show Gradient @@ -324,7 +331,7 @@ def realtime_zshim(fname_fmap, fname_mask_anat, fname_resp, fname_json, fname_an im = ax.imshow(gz_gradient[:, :, 0, 0]) fig.colorbar(im) ax.set_title("Gradient [:, :, 0, 0]") - fname_figure = os.path.join(__dir_shimmingtoolbox__, 'fig_realtime_zshim_gradient.png') + fname_figure = os.path.join(fname_output, 'fig_realtime_zshim_gradient.png') fig.savefig(fname_figure) # Show evolution of coefficients @@ -335,7 +342,7 @@ def realtime_zshim(fname_fmap, fname_mask_anat, fname_resp, fname_json, fname_an ax = fig.add_subplot(2, 1, 2) ax.plot(range(n_slices), riro_correction, label='Riro correction') ax.set_title("Riro correction evolution through slices") - fname_figure = os.path.join(__dir_shimmingtoolbox__, 'fig_realtime_zshim_correction_slice.png') + fname_figure = os.path.join(fname_output, 'fig_realtime_zshim_correction_slice.png') fig.savefig(fname_figure) - return fname_figure + return fname_corrections diff --git a/test/cli/test_cli_realtime_zshim.py b/test/cli/test_cli_realtime_zshim.py index 0f0d45357..9f1a93811 100644 --- a/test/cli/test_cli_realtime_zshim.py +++ b/test/cli/test_cli_realtime_zshim.py @@ -14,6 +14,7 @@ from shimmingtoolbox.coils.coordinates import generate_meshgrid from shimmingtoolbox.coils.siemens_basis import siemens_basis from shimmingtoolbox import __dir_testing__ +from shimmingtoolbox import __dir_shimmingtoolbox__ def test_cli_realtime_zshim(): @@ -61,8 +62,16 @@ def test_cli_realtime_zshim(): fname_json = os.path.join(__dir_testing__, 'realtime_zshimming_data', 'nifti', 'sub-example', 'fmap', 'sub-example_magnitude1.json') - result = runner.invoke(realtime_zshim, ['-fmap', fname_fieldmap, '-mask', fname_mask, - '-resp', fname_resp, '-json', fname_json, '-anat', fname_anat], + # Specify output for text file and figures + fname_output = os.path.join(__dir_shimmingtoolbox__, 'test_realtime_zshim') + + # Run the CLI + result = runner.invoke(realtime_zshim, ['-fmap', fname_fieldmap, + '-mask', fname_mask, + '-output', fname_output, + '-resp', fname_resp, + '-json', fname_json, + '-anat', fname_anat], catch_exceptions=False) assert result.exit_code == 0 From 6b5b64239016d377e3a08b9eaf872d8362a35268 Mon Sep 17 00:00:00 2001 From: Alexandre D'Astous Date: Fri, 6 Nov 2020 16:22:02 -0500 Subject: [PATCH 43/51] Clean up and added more to DEBUG --- shimmingtoolbox/cli/realtime_zshim.py | 277 +++++++++++++------------- 1 file changed, 139 insertions(+), 138 deletions(-) diff --git a/shimmingtoolbox/cli/realtime_zshim.py b/shimmingtoolbox/cli/realtime_zshim.py index 1342dc668..2b8178f0f 100644 --- a/shimmingtoolbox/cli/realtime_zshim.py +++ b/shimmingtoolbox/cli/realtime_zshim.py @@ -10,12 +10,10 @@ from nibabel.processing import resample_from_to # TODO: remove matplotlib and dirtesting import from matplotlib.figure import Figure -from shimmingtoolbox import __dir_testing__ from shimmingtoolbox.optimizer.sequential import sequential_zslice from shimmingtoolbox.load_nifti import get_acquisition_times from shimmingtoolbox.pmu import PmuResp -from shimmingtoolbox import __dir_shimmingtoolbox__ from shimmingtoolbox.utils import st_progress_bar DEBUG = True @@ -159,9 +157,7 @@ def realtime_zshim(fname_fmap, fname_mask_anat, fname_resp, fname_json, fname_an for i_x in range(fieldmap.shape[0]): for i_y in range(fieldmap.shape[1]): for i_z in range(fieldmap.shape[2]): - # TODO: Fit for -masked_field? reg = LinearRegression().fit(acq_pressures.reshape(-1, 1) - mean_p, -gz_gradient[i_x, i_y, i_z, :]) - # reg = LinearRegression().fit(acq_pressures.reshape(-1, 1) - mean_p, -masked_fieldmaps[i_x, i_y, i_z, :]) riro[i_x, i_y, i_z] = reg.coef_ static[i_x, i_y, i_z] = reg.intercept_ progress_bar.update(1) @@ -174,21 +170,25 @@ def realtime_zshim(fname_fmap, fname_mask_anat, fname_resp, fname_json, fname_an nii_resampled_fmap_3d = resample_from_to(nii_masked_fmap_3d, nii_anat, order=2, mode='nearest') masked_fmap_4d[..., it] = nii_resampled_fmap_3d.get_fdata() nii_resampled_fmap = nib.Nifti1Image(masked_fmap_4d, nii_anat.affine) - nib.save(nii_resampled_fmap, os.path.join(fname_output, 'resampled_fmap.nii.gz')) + + if DEBUG: + nib.save(nii_resampled_fmap, os.path.join(fname_output, 'resampled_fmap.nii.gz')) # Resample static to target anatomical image nii_static = nib.Nifti1Image(static, nii_fmap.affine) nii_resampled_static = resample_from_to(nii_static, nii_anat, mode='nearest') nii_resampled_static_masked = nib.Nifti1Image(nii_resampled_static.get_fdata() * nii_mask_anat.get_fdata(), nii_resampled_static.affine) - nib.save(nii_resampled_static_masked, os.path.join(fname_output, 'resampled_static.nii.gz')) + if DEBUG: + nib.save(nii_resampled_static_masked, os.path.join(fname_output, 'resampled_static.nii.gz')) # Resample riro to target anatomical image nii_riro = nib.Nifti1Image(riro, nii_fmap.affine) nii_resampled_riro = resample_from_to(nii_riro, nii_anat, mode='nearest') nii_resampled_riro_masked = nib.Nifti1Image(nii_resampled_riro.get_fdata() * nii_mask_anat.get_fdata(), nii_resampled_riro.affine) - nib.save(nii_resampled_riro_masked, os.path.join(fname_output, 'resampled_riro.nii.gz')) + if DEBUG: + nib.save(nii_resampled_riro_masked, os.path.join(fname_output, 'resampled_riro.nii.gz')) # Calculate the mean for riro and static for a perticular slice n_slices = nii_anat.get_fdata().shape[2] @@ -203,7 +203,6 @@ def realtime_zshim(fname_fmap, fname_mask_anat, fname_resp, fname_json, fname_an riro_correction[i_slice] = np.ma.mean(ma_riro_anat) # Write to a text file - # TODO: Add as an option to output the file to a specified location fname_corrections = os.path.join(fname_output, 'zshim_gradients.txt') file_gradients = open(fname_corrections, 'w') for i_slice in range(n_slices): @@ -214,135 +213,137 @@ def realtime_zshim(fname_fmap, fname_mask_anat, fname_resp, fname_json, fname_an # ================ PLOTS ================ - # Calculate masked shim for spherical harmonics plot - # masked_shimmed = np.zeros_like(shimmed) - # for i_t in range(nt): - # masked_shimmed[..., i_t] = mask_fmap * shimmed[..., i_t] - - # Plot unshimmed vs shimmed and their mask for spherical harmonics - # i_t = 0 - # fig = Figure(figsize=(10, 10)) - # ax = fig.add_subplot(2, 2, 1) - # im = ax.imshow(masked_fieldmaps[:-1, :-1, 0, i_t]) - # fig.colorbar(im) - # ax.set_title("Masked unshimmed") - # ax = fig.add_subplot(2, 2, 2) - # im = ax.imshow(masked_shimmed[:-1, :-1, 0, i_t]) - # fig.colorbar(im) - # ax.set_title("Masked shimmed") - # ax = fig.add_subplot(2, 2, 3) - # im = ax.imshow(fieldmap[:-1, :-1, 0, i_t]) - # fig.colorbar(im) - # ax.set_title("Unshimmed") - # ax = fig.add_subplot(2, 2, 4) - # im = ax.imshow(shimmed[:-1, :-1, 0, i_t]) - # fig.colorbar(im) - # ax.set_title("Shimmed") - # fname_figure = os.path.join(__dir_shimmingtoolbox__, 'fig_realtime_zshim_sphharm_shimmed.png') - # fig.savefig(fname_figure) - - # Plot the coil coefs through time - # fig = Figure(figsize=(10, 10)) - # for i_coil in range(n_coils): - # ax = fig.add_subplot(n_coils, 1, i_coil + 1) - # ax.plot(np.arange(nt), currents[i_coil, :]) - # ax.set_title(f"Channel {i_coil}") - # fname_figure = os.path.join(__dir_shimmingtoolbox__, 'fig_realtime_zshim_sphharm_currents.png') - # fig.savefig(fname_figure) - - # Plot Static and RIRO - fig = Figure(figsize=(10, 10)) - ax = fig.add_subplot(2, 1, 1) - im = ax.imshow(riro[:-1, :-1, 0]) - fig.colorbar(im) - ax.set_title("RIRO") - ax = fig.add_subplot(2, 1, 2) - im = ax.imshow(static[:-1, :-1, 0]) - fig.colorbar(im) - ax.set_title("Static") - fname_figure = os.path.join(fname_output, 'fig_realtime_zshim_riro_static.png') - fig.savefig(fname_figure) - - # Calculate fitted and shimmed for pressure fitted plot - fitted_fieldmap = riro * (acq_pressures - mean_p) + static - shimmed_pressure_fitted = np.expand_dims(fitted_fieldmap, 2) + masked_fieldmaps - - # Plot pressure fitted fieldmap - fig = Figure(figsize=(10, 10)) - ax = fig.add_subplot(3, 1, 1) - im = ax.imshow(masked_fieldmaps[:-1, :-1, 0, i_t]) - fig.colorbar(im) - ax.set_title("fieldmap") - ax = fig.add_subplot(3, 1, 2) - im = ax.imshow(fitted_fieldmap[:-1, :-1, i_t]) - fig.colorbar(im) - ax.set_title("Fit") - ax = fig.add_subplot(3, 1, 3) - im = ax.imshow(shimmed_pressure_fitted[:-1, :-1, 0, i_t]) - fig.colorbar(im) - ax.set_title("Shimmed (fit + fieldmap") - fname_figure = os.path.join(fname_output, 'fig_realtime_zshim_pressure_fitted.png') - fig.savefig(fname_figure) - - # Reshape pmu datapoints to fit those of the acquisition - pmu_times = np.linspace(pmu.start_time_mdh, pmu.stop_time_mdh, len(pmu.data)) - pmu_times_within_range = pmu_times[pmu_times > acq_timestamps[0]] - pmu_data_within_range = pmu.data[pmu_times > acq_timestamps[0]] - pmu_data_within_range = pmu_data_within_range[pmu_times_within_range < acq_timestamps[fieldmap.shape[3] - 1]] - pmu_times_within_range = pmu_times_within_range[pmu_times_within_range < acq_timestamps[fieldmap.shape[3] - 1]] - - # Calc fieldmap average within mask - fieldmap_avg = np.zeros([fieldmap.shape[3]]) - for i_time in range(nt): - masked_array = np.ma.array(fieldmap[:, :, :, i_time], mask=mask_fmap == False) - fieldmap_avg[i_time] = np.ma.average(masked_array) - - # Plot pmu vs B0 in masked region - fig = Figure(figsize=(10, 10)) - ax = fig.add_subplot(211) - ax.plot(acq_timestamps / 1000, acq_pressures, label='Interpolated pressures') - # ax.plot(pmu_times / 1000, pmu.data, label='Raw pressures') - ax.plot(pmu_times_within_range / 1000, pmu_data_within_range, label='Pmu pressures') - ax.legend() - ax.set_title("Pressure [0, 4095] vs time (s) ") - ax = fig.add_subplot(212) - ax.plot(acq_timestamps / 1000, fieldmap_avg, label='Mean B0') - ax.legend() - ax.set_title("Fieldmap average over unmasked region (Hz) vs time (s)") - fname_figure = os.path.join(fname_output, 'fig_realtime_zshim_pmu_vs_B0.png') - fig.savefig(fname_figure) - - # Show anatomical image - fig = Figure(figsize=(10, 10)) - ax = fig.add_subplot(2, 1, 1) - im = ax.imshow(anat[:, :, 10]) - fig.colorbar(im) - ax.set_title("Anatomical image [:, :, 10]") - ax = fig.add_subplot(2, 1, 2) - im = ax.imshow(nii_resampled_fmap.get_fdata()[:, :, 10, 0]) - fig.colorbar(im) - ax.set_title("Resampled fieldmap [:, :, 10, 0]") - fname_figure = os.path.join(fname_output, 'fig_reatime_zshim_anat.png') - fig.savefig(fname_figure) - - # Show Gradient - fig = Figure(figsize=(10, 10)) - ax = fig.add_subplot(1, 1, 1) - im = ax.imshow(gz_gradient[:, :, 0, 0]) - fig.colorbar(im) - ax.set_title("Gradient [:, :, 0, 0]") - fname_figure = os.path.join(fname_output, 'fig_realtime_zshim_gradient.png') - fig.savefig(fname_figure) - - # Show evolution of coefficients - fig = Figure(figsize=(10, 10)) - ax = fig.add_subplot(2, 1, 1) - ax.plot(range(n_slices), static_correction, label='Static correction') - ax.set_title("Static correction evolution through slices") - ax = fig.add_subplot(2, 1, 2) - ax.plot(range(n_slices), riro_correction, label='Riro correction') - ax.set_title("Riro correction evolution through slices") - fname_figure = os.path.join(fname_output, 'fig_realtime_zshim_correction_slice.png') - fig.savefig(fname_figure) + if DEBUG: + + # Calculate masked shim for spherical harmonics plot + # masked_shimmed = np.zeros_like(shimmed) + # for i_t in range(nt): + # masked_shimmed[..., i_t] = mask_fmap * shimmed[..., i_t] + + # Plot unshimmed vs shimmed and their mask for spherical harmonics + # i_t = 0 + # fig = Figure(figsize=(10, 10)) + # ax = fig.add_subplot(2, 2, 1) + # im = ax.imshow(masked_fieldmaps[:-1, :-1, 0, i_t]) + # fig.colorbar(im) + # ax.set_title("Masked unshimmed") + # ax = fig.add_subplot(2, 2, 2) + # im = ax.imshow(masked_shimmed[:-1, :-1, 0, i_t]) + # fig.colorbar(im) + # ax.set_title("Masked shimmed") + # ax = fig.add_subplot(2, 2, 3) + # im = ax.imshow(fieldmap[:-1, :-1, 0, i_t]) + # fig.colorbar(im) + # ax.set_title("Unshimmed") + # ax = fig.add_subplot(2, 2, 4) + # im = ax.imshow(shimmed[:-1, :-1, 0, i_t]) + # fig.colorbar(im) + # ax.set_title("Shimmed") + # fname_figure = os.path.join(__dir_shimmingtoolbox__, 'fig_realtime_zshim_sphharm_shimmed.png') + # fig.savefig(fname_figure) + + # Plot the coil coefs through time + # fig = Figure(figsize=(10, 10)) + # for i_coil in range(n_coils): + # ax = fig.add_subplot(n_coils, 1, i_coil + 1) + # ax.plot(np.arange(nt), currents[i_coil, :]) + # ax.set_title(f"Channel {i_coil}") + # fname_figure = os.path.join(__dir_shimmingtoolbox__, 'fig_realtime_zshim_sphharm_currents.png') + # fig.savefig(fname_figure) + + # Plot Static and RIRO + fig = Figure(figsize=(10, 10)) + ax = fig.add_subplot(2, 1, 1) + im = ax.imshow(riro[:-1, :-1, 0]) + fig.colorbar(im) + ax.set_title("RIRO") + ax = fig.add_subplot(2, 1, 2) + im = ax.imshow(static[:-1, :-1, 0]) + fig.colorbar(im) + ax.set_title("Static") + fname_figure = os.path.join(fname_output, 'fig_realtime_zshim_riro_static.png') + fig.savefig(fname_figure) + + # Calculate fitted and shimmed for pressure fitted plot + fitted_fieldmap = riro * (acq_pressures - mean_p) + static + shimmed_pressure_fitted = np.expand_dims(fitted_fieldmap, 2) + masked_fieldmaps + + # Plot pressure fitted fieldmap + fig = Figure(figsize=(10, 10)) + ax = fig.add_subplot(3, 1, 1) + im = ax.imshow(masked_fieldmaps[:-1, :-1, 0, i_t]) + fig.colorbar(im) + ax.set_title("fieldmap") + ax = fig.add_subplot(3, 1, 2) + im = ax.imshow(fitted_fieldmap[:-1, :-1, i_t]) + fig.colorbar(im) + ax.set_title("Fit") + ax = fig.add_subplot(3, 1, 3) + im = ax.imshow(shimmed_pressure_fitted[:-1, :-1, 0, i_t]) + fig.colorbar(im) + ax.set_title("Shimmed (fit + fieldmap") + fname_figure = os.path.join(fname_output, 'fig_realtime_zshim_pressure_fitted.png') + fig.savefig(fname_figure) + + # Reshape pmu datapoints to fit those of the acquisition + pmu_times = np.linspace(pmu.start_time_mdh, pmu.stop_time_mdh, len(pmu.data)) + pmu_times_within_range = pmu_times[pmu_times > acq_timestamps[0]] + pmu_data_within_range = pmu.data[pmu_times > acq_timestamps[0]] + pmu_data_within_range = pmu_data_within_range[pmu_times_within_range < acq_timestamps[fieldmap.shape[3] - 1]] + pmu_times_within_range = pmu_times_within_range[pmu_times_within_range < acq_timestamps[fieldmap.shape[3] - 1]] + + # Calc fieldmap average within mask + fieldmap_avg = np.zeros([fieldmap.shape[3]]) + for i_time in range(nt): + masked_array = np.ma.array(fieldmap[:, :, :, i_time], mask=mask_fmap == False) + fieldmap_avg[i_time] = np.ma.average(masked_array) + + # Plot pmu vs B0 in masked region + fig = Figure(figsize=(10, 10)) + ax = fig.add_subplot(211) + ax.plot(acq_timestamps / 1000, acq_pressures, label='Interpolated pressures') + # ax.plot(pmu_times / 1000, pmu.data, label='Raw pressures') + ax.plot(pmu_times_within_range / 1000, pmu_data_within_range, label='Pmu pressures') + ax.legend() + ax.set_title("Pressure [0, 4095] vs time (s) ") + ax = fig.add_subplot(212) + ax.plot(acq_timestamps / 1000, fieldmap_avg, label='Mean B0') + ax.legend() + ax.set_title("Fieldmap average over unmasked region (Hz) vs time (s)") + fname_figure = os.path.join(fname_output, 'fig_realtime_zshim_pmu_vs_B0.png') + fig.savefig(fname_figure) + + # Show anatomical image + fig = Figure(figsize=(10, 10)) + ax = fig.add_subplot(2, 1, 1) + im = ax.imshow(anat[:, :, 10]) + fig.colorbar(im) + ax.set_title("Anatomical image [:, :, 10]") + ax = fig.add_subplot(2, 1, 2) + im = ax.imshow(nii_resampled_fmap.get_fdata()[:, :, 10, 0]) + fig.colorbar(im) + ax.set_title("Resampled fieldmap [:, :, 10, 0]") + fname_figure = os.path.join(fname_output, 'fig_reatime_zshim_anat.png') + fig.savefig(fname_figure) + + # Show Gradient + fig = Figure(figsize=(10, 10)) + ax = fig.add_subplot(1, 1, 1) + im = ax.imshow(gz_gradient[:, :, 0, 0]) + fig.colorbar(im) + ax.set_title("Gradient [:, :, 0, 0]") + fname_figure = os.path.join(fname_output, 'fig_realtime_zshim_gradient.png') + fig.savefig(fname_figure) + + # Show evolution of coefficients + fig = Figure(figsize=(10, 10)) + ax = fig.add_subplot(2, 1, 1) + ax.plot(range(n_slices), static_correction, label='Static correction') + ax.set_title("Static correction evolution through slices") + ax = fig.add_subplot(2, 1, 2) + ax.plot(range(n_slices), riro_correction, label='Riro correction') + ax.set_title("Riro correction evolution through slices") + fname_figure = os.path.join(fname_output, 'fig_realtime_zshim_correction_slice.png') + fig.savefig(fname_figure) return fname_corrections From a3b0830fff08a17beee58c73c19853d95512feb5 Mon Sep 17 00:00:00 2001 From: Alexandre D'Astous Date: Fri, 6 Nov 2020 16:42:44 -0500 Subject: [PATCH 44/51] Removed plot that only works for fieldmap fitting --- shimmingtoolbox/cli/realtime_zshim.py | 38 +++++++++++++-------------- 1 file changed, 19 insertions(+), 19 deletions(-) diff --git a/shimmingtoolbox/cli/realtime_zshim.py b/shimmingtoolbox/cli/realtime_zshim.py index 2b8178f0f..17d9199a2 100644 --- a/shimmingtoolbox/cli/realtime_zshim.py +++ b/shimmingtoolbox/cli/realtime_zshim.py @@ -265,25 +265,25 @@ def realtime_zshim(fname_fmap, fname_mask_anat, fname_resp, fname_json, fname_an fig.savefig(fname_figure) # Calculate fitted and shimmed for pressure fitted plot - fitted_fieldmap = riro * (acq_pressures - mean_p) + static - shimmed_pressure_fitted = np.expand_dims(fitted_fieldmap, 2) + masked_fieldmaps - - # Plot pressure fitted fieldmap - fig = Figure(figsize=(10, 10)) - ax = fig.add_subplot(3, 1, 1) - im = ax.imshow(masked_fieldmaps[:-1, :-1, 0, i_t]) - fig.colorbar(im) - ax.set_title("fieldmap") - ax = fig.add_subplot(3, 1, 2) - im = ax.imshow(fitted_fieldmap[:-1, :-1, i_t]) - fig.colorbar(im) - ax.set_title("Fit") - ax = fig.add_subplot(3, 1, 3) - im = ax.imshow(shimmed_pressure_fitted[:-1, :-1, 0, i_t]) - fig.colorbar(im) - ax.set_title("Shimmed (fit + fieldmap") - fname_figure = os.path.join(fname_output, 'fig_realtime_zshim_pressure_fitted.png') - fig.savefig(fname_figure) + # fitted_fieldmap = riro * (acq_pressures - mean_p) + static + # shimmed_pressure_fitted = np.expand_dims(fitted_fieldmap, 2) + masked_fieldmaps + # + # # Plot pressure fitted fieldmap + # fig = Figure(figsize=(10, 10)) + # ax = fig.add_subplot(3, 1, 1) + # im = ax.imshow(masked_fieldmaps[:-1, :-1, 0, i_t]) + # fig.colorbar(im) + # ax.set_title("fieldmap") + # ax = fig.add_subplot(3, 1, 2) + # im = ax.imshow(fitted_fieldmap[:-1, :-1, i_t]) + # fig.colorbar(im) + # ax.set_title("Fit") + # ax = fig.add_subplot(3, 1, 3) + # im = ax.imshow(shimmed_pressure_fitted[:-1, :-1, 0, i_t]) + # fig.colorbar(im) + # ax.set_title("Shimmed (fit + fieldmap") + # fname_figure = os.path.join(fname_output, 'fig_realtime_zshim_pressure_fitted.png') + # fig.savefig(fname_figure) # Reshape pmu datapoints to fit those of the acquisition pmu_times = np.linspace(pmu.start_time_mdh, pmu.stop_time_mdh, len(pmu.data)) From 86df361893e176ef15b00a9999ecb93265dcbe9a Mon Sep 17 00:00:00 2001 From: Alexandre D'Astous Date: Fri, 6 Nov 2020 17:16:29 -0500 Subject: [PATCH 45/51] Add isclose to input masks with sligth affine variation --- shimmingtoolbox/cli/realtime_zshim.py | 3 ++- 1 file changed, 2 insertions(+), 1 deletion(-) diff --git a/shimmingtoolbox/cli/realtime_zshim.py b/shimmingtoolbox/cli/realtime_zshim.py index 17d9199a2..e2bba4eab 100644 --- a/shimmingtoolbox/cli/realtime_zshim.py +++ b/shimmingtoolbox/cli/realtime_zshim.py @@ -87,7 +87,8 @@ def realtime_zshim(fname_fmap, fname_mask_anat, fname_resp, fname_json, fname_an # TODO: check good practice below if fname_mask_anat is not None: nii_mask_anat = nib.load(fname_mask_anat) - if not np.all(nii_anat.affine == nii_mask_anat.affine) or not np.all(nii_mask_anat.shape == nii_anat.shape): + if not np.all(np.isclose(nii_anat.affine, nii_mask_anat.affine)) or\ + not np.all(nii_mask_anat.shape == nii_anat.shape): raise RuntimeError("Mask must have the same shape and affine transformation as anat") nii_fmap_3d_temp = nib.Nifti1Image(fieldmap[..., 0], nii_fmap.affine) nii_mask_fmap = resample_from_to(nii_mask_anat, nii_fmap_3d_temp) From 16e76d75d2fc8505af39e5dcf4826503933b18f4 Mon Sep 17 00:00:00 2001 From: Alexandre D'Astous Date: Mon, 9 Nov 2020 18:25:05 -0500 Subject: [PATCH 46/51] Add meshgrid and local mask for validation --- shimmingtoolbox/cli/realtime_zshim.py | 8 ++++++-- test/cli/test_cli_realtime_zshim.py | 2 ++ 2 files changed, 8 insertions(+), 2 deletions(-) diff --git a/shimmingtoolbox/cli/realtime_zshim.py b/shimmingtoolbox/cli/realtime_zshim.py index e2bba4eab..2b0201c46 100644 --- a/shimmingtoolbox/cli/realtime_zshim.py +++ b/shimmingtoolbox/cli/realtime_zshim.py @@ -15,6 +15,7 @@ from shimmingtoolbox.load_nifti import get_acquisition_times from shimmingtoolbox.pmu import PmuResp from shimmingtoolbox.utils import st_progress_bar +from shimmingtoolbox.coils.coordinates import generate_meshgrid DEBUG = True CONTEXT_SETTINGS = dict(help_option_names=['-h', '--help']) @@ -114,12 +115,15 @@ def realtime_zshim(fname_fmap, fname_mask_anat, fname_resp, fname_json, fname_an # Calculate gz gradient # Image is z, y, x - # Pixdim[3] is the space between pixels in the z direction in millimeters + # Pixdim[2] is the space between pixels in the z direction in millimeters g = 1000 / 42.576e6 # [mT / Hz] gz_gradient = np.zeros_like(fieldmap) + coord = generate_meshgrid(mask_fmap.shape, nii_fmap.affine) + + for it in range(nt): gz_gradient[..., 0, it] = np.gradient(g * fieldmap[:, :, 0, it], - nii_fmap.header['pixdim'][3] / 1000, + -nii_fmap.header['pixdim'][2] / 1000, axis=1) # [mT / m] if DEBUG: nii_gz_gradient = nib.Nifti1Image(gz_gradient, nii_fmap.affine) diff --git a/test/cli/test_cli_realtime_zshim.py b/test/cli/test_cli_realtime_zshim.py index 9f1a93811..808ae0c34 100644 --- a/test/cli/test_cli_realtime_zshim.py +++ b/test/cli/test_cli_realtime_zshim.py @@ -46,6 +46,8 @@ def test_cli_realtime_zshim(): nii_mask = nib.Nifti1Image(mask.astype(int), nii_anat.affine) fname_mask = os.path.join(tmp, 'mask.nii.gz') nib.save(nii_mask, fname_mask) + # fname_mask = os.path.join(__dir_testing__, 'test_realtime_zshim', 'gre_seg_30.nii') + # Set up coils # coord_phys = generate_meshgrid(nii_fmap.get_fdata().shape[0:3], nii_fmap.affine) From 001adf2d183a27dc484e4faec22cc7c99921f73e Mon Sep 17 00:00:00 2001 From: Alexandre D'Astous Date: Thu, 12 Nov 2020 00:37:55 -0500 Subject: [PATCH 47/51] Update gradient to use meshgrid points --- shimmingtoolbox/cli/realtime_zshim.py | 8 +++----- 1 file changed, 3 insertions(+), 5 deletions(-) diff --git a/shimmingtoolbox/cli/realtime_zshim.py b/shimmingtoolbox/cli/realtime_zshim.py index 2b0201c46..8fadebe84 100644 --- a/shimmingtoolbox/cli/realtime_zshim.py +++ b/shimmingtoolbox/cli/realtime_zshim.py @@ -114,16 +114,14 @@ def realtime_zshim(fname_fmap, fname_mask_anat, fname_resp, fname_json, fname_an masked_fieldmaps[..., i_t] = mask_fmap * fieldmap[..., i_t] # Calculate gz gradient - # Image is z, y, x - # Pixdim[2] is the space between pixels in the z direction in millimeters g = 1000 / 42.576e6 # [mT / Hz] gz_gradient = np.zeros_like(fieldmap) - coord = generate_meshgrid(mask_fmap.shape, nii_fmap.affine) - + # Get voxel coordinates. Z coordinates correspond to coord[2] + z_coord = generate_meshgrid(mask_fmap.shape, nii_fmap.affine)[2] / 1000 # [m] for it in range(nt): gz_gradient[..., 0, it] = np.gradient(g * fieldmap[:, :, 0, it], - -nii_fmap.header['pixdim'][2] / 1000, + z_coord[0, :, 0], axis=1) # [mT / m] if DEBUG: nii_gz_gradient = nib.Nifti1Image(gz_gradient, nii_fmap.affine) From 0a385405cd72589ce70f69f84648486af07e56b4 Mon Sep 17 00:00:00 2001 From: Alexandre D'Astous Date: Thu, 12 Nov 2020 01:24:35 -0500 Subject: [PATCH 48/51] Add scaling by RMS to RIRO to have comparable figures across sunjects --- shimmingtoolbox/cli/realtime_zshim.py | 13 +++++++------ 1 file changed, 7 insertions(+), 6 deletions(-) diff --git a/shimmingtoolbox/cli/realtime_zshim.py b/shimmingtoolbox/cli/realtime_zshim.py index 8fadebe84..32e8b45df 100644 --- a/shimmingtoolbox/cli/realtime_zshim.py +++ b/shimmingtoolbox/cli/realtime_zshim.py @@ -120,9 +120,7 @@ def realtime_zshim(fname_fmap, fname_mask_anat, fname_resp, fname_json, fname_an z_coord = generate_meshgrid(mask_fmap.shape, nii_fmap.affine)[2] / 1000 # [m] for it in range(nt): - gz_gradient[..., 0, it] = np.gradient(g * fieldmap[:, :, 0, it], - z_coord[0, :, 0], - axis=1) # [mT / m] + gz_gradient[..., 0, it] = np.gradient(g * fieldmap[:, :, 0, it], z_coord[0, :, 0], axis=1) # [mT / m] if DEBUG: nii_gz_gradient = nib.Nifti1Image(gz_gradient, nii_fmap.affine) nib.save(nii_gz_gradient, os.path.join(fname_output, 'tmp.gz_gradient.nii.gz')) @@ -153,6 +151,7 @@ def realtime_zshim(fname_fmap, fname_mask_anat, fname_resp, fname_json, fname_an # Shim using PMU mean_p = np.mean(acq_pressures) + pressure_rms = np.sqrt(np.mean((acq_pressures - mean_p) ** 2)) riro = np.zeros_like(fieldmap[:, :, :, 0]) static = np.zeros_like(fieldmap[:, :, :, 0]) # TODO fix progress bar not showing up @@ -161,7 +160,9 @@ def realtime_zshim(fname_fmap, fname_mask_anat, fname_resp, fname_json, fname_an for i_y in range(fieldmap.shape[1]): for i_z in range(fieldmap.shape[2]): reg = LinearRegression().fit(acq_pressures.reshape(-1, 1) - mean_p, -gz_gradient[i_x, i_y, i_z, :]) - riro[i_x, i_y, i_z] = reg.coef_ + # Multiplying by the RMS of the pressure allows to make abstraction of the tightness of the bellow + # between scans. This allows to compare results between scans. + riro[i_x, i_y, i_z] = reg.coef_ * pressure_rms static[i_x, i_y, i_z] = reg.intercept_ progress_bar.update(1) @@ -210,7 +211,7 @@ def realtime_zshim(fname_fmap, fname_mask_anat, fname_resp, fname_json, fname_an file_gradients = open(fname_corrections, 'w') for i_slice in range(n_slices): file_gradients.write(f'Vector_Gz[0][{i_slice}]= {static_correction[i_slice]:.6f}\n') - file_gradients.write(f'Vector_Gz[1][{i_slice}]= {riro_correction[i_slice]:.12f}\n') + file_gradients.write(f'Vector_Gz[1][{i_slice}]= {riro_correction[i_slice] / pressure_rms:.12f}\n') file_gradients.write(f'Vector_Gz[2][{i_slice}]= {mean_p:.3f}\n') file_gradients.close() @@ -257,7 +258,7 @@ def realtime_zshim(fname_fmap, fname_mask_anat, fname_resp, fname_json, fname_an # Plot Static and RIRO fig = Figure(figsize=(10, 10)) ax = fig.add_subplot(2, 1, 1) - im = ax.imshow(riro[:-1, :-1, 0]) + im = ax.imshow(riro[:-1, :-1, 0] / pressure_rms) fig.colorbar(im) ax.set_title("RIRO") ax = fig.add_subplot(2, 1, 2) From 70ac285557d1a1a7400a941d9fd832a6f572d7d4 Mon Sep 17 00:00:00 2001 From: Alexandre D'Astous Date: Sun, 15 Nov 2020 23:45:24 -0500 Subject: [PATCH 49/51] Update figures --- shimmingtoolbox/cli/realtime_zshim.py | 10 +++++----- 1 file changed, 5 insertions(+), 5 deletions(-) diff --git a/shimmingtoolbox/cli/realtime_zshim.py b/shimmingtoolbox/cli/realtime_zshim.py index 32e8b45df..bef3a1d9d 100644 --- a/shimmingtoolbox/cli/realtime_zshim.py +++ b/shimmingtoolbox/cli/realtime_zshim.py @@ -204,14 +204,14 @@ def realtime_zshim(fname_fmap, fname_mask_anat, fname_resp, fname_json, fname_an static_correction[i_slice] = np.ma.mean(ma_static_anat) ma_riro_anat = np.ma.array(nii_resampled_riro.get_fdata()[..., i_slice], mask=nii_mask_anat.get_fdata()[..., i_slice] == False) - riro_correction[i_slice] = np.ma.mean(ma_riro_anat) + riro_correction[i_slice] = np.ma.mean(ma_riro_anat) / pressure_rms # Write to a text file fname_corrections = os.path.join(fname_output, 'zshim_gradients.txt') file_gradients = open(fname_corrections, 'w') for i_slice in range(n_slices): file_gradients.write(f'Vector_Gz[0][{i_slice}]= {static_correction[i_slice]:.6f}\n') - file_gradients.write(f'Vector_Gz[1][{i_slice}]= {riro_correction[i_slice] / pressure_rms:.12f}\n') + file_gradients.write(f'Vector_Gz[1][{i_slice}]= {riro_correction[i_slice]:.12f}\n') file_gradients.write(f'Vector_Gz[2][{i_slice}]= {mean_p:.3f}\n') file_gradients.close() @@ -324,9 +324,9 @@ def realtime_zshim(fname_fmap, fname_mask_anat, fname_resp, fname_json, fname_an fig.colorbar(im) ax.set_title("Anatomical image [:, :, 10]") ax = fig.add_subplot(2, 1, 2) - im = ax.imshow(nii_resampled_fmap.get_fdata()[:, :, 10, 0]) + im = ax.imshow(nii_mask_anat.get_fdata()[:, :, 10]) fig.colorbar(im) - ax.set_title("Resampled fieldmap [:, :, 10, 0]") + ax.set_title("Mask [:, :, 10]") fname_figure = os.path.join(fname_output, 'fig_reatime_zshim_anat.png') fig.savefig(fname_figure) @@ -345,7 +345,7 @@ def realtime_zshim(fname_fmap, fname_mask_anat, fname_resp, fname_json, fname_an ax.plot(range(n_slices), static_correction, label='Static correction') ax.set_title("Static correction evolution through slices") ax = fig.add_subplot(2, 1, 2) - ax.plot(range(n_slices), riro_correction, label='Riro correction') + ax.plot(range(n_slices), (acq_pressures.max() - mean_p) * riro_correction, label='Riro correction') ax.set_title("Riro correction evolution through slices") fname_figure = os.path.join(fname_output, 'fig_realtime_zshim_correction_slice.png') fig.savefig(fname_figure) From d62563731838a28b55f4584230e71e161ad1a1c0 Mon Sep 17 00:00:00 2001 From: Alexandre D'Astous Date: Wed, 25 Nov 2020 13:22:19 -0500 Subject: [PATCH 50/51] Remove comments related to sequencer. Couls be relevant for #101 --- shimmingtoolbox/cli/realtime_zshim.py | 73 --------------------------- 1 file changed, 73 deletions(-) diff --git a/shimmingtoolbox/cli/realtime_zshim.py b/shimmingtoolbox/cli/realtime_zshim.py index bef3a1d9d..871e383b8 100644 --- a/shimmingtoolbox/cli/realtime_zshim.py +++ b/shimmingtoolbox/cli/realtime_zshim.py @@ -25,10 +25,6 @@ context_settings=CONTEXT_SETTINGS, help=f"Perform realtime z-shimming." ) -# @click.option('-coil', 'fname_coil', required=True, type=click.Path(), -# help="Coil basis to use for shimming. Enter multiple files if " -# "you wish to use more than one set of shim coils (eg: " -# "Siemens gradient/shim coils and external custom coils).") @click.option('-fmap', 'fname_fmap', required=True, type=click.Path(), help="B0 fieldmap. For realtime shimming, this should be a 4d file (4th dimension being time") @click.option('-mask', 'fname_mask_anat', type=click.Path(), @@ -59,11 +55,6 @@ def realtime_zshim(fname_fmap, fname_mask_anat, fname_resp, fname_json, fname_an Returns: """ - # Load coil - # When using only z channel (corresponding to index 0) TODO:Remove - # coil = np.expand_dims(nib.load(fname_coil).get_fdata()[:, :, :, 0], -1) - # When using all channels TODO: Keep - # coil = nib.load(fname_coil).get_fdata() # Look if output directory exists, if not, create it if not os.path.exists(fname_output): @@ -102,15 +93,8 @@ def realtime_zshim(fname_fmap, fname_mask_anat, fname_resp, fname_json, fname_an if DEBUG: nib.save(nii_mask_fmap, os.path.join(fname_output, 'tmp.mask_fmap_resample.nii.gz')) - # Shim using sequencer and optimizer - # n_coils = coil.shape[-1] - # currents = np.zeros([n_coils, nt]) - # shimmed = np.zeros_like(fieldmap) masked_fieldmaps = np.zeros_like(fieldmap) for i_t in range(nt): - # currents[:, i_t] = sequential_zslice(fieldmap[..., i_t], coil, mask_fmap, z_slices=np.array(range(nz)), - # bounds=[(-np.inf, np.inf)] * n_coils) - # shimmed[..., i_t] = fieldmap[..., i_t] + np.sum(currents[:, i_t] * coil, axis=3, keepdims=False) masked_fieldmaps[..., i_t] = mask_fmap * fieldmap[..., i_t] # Calculate gz gradient @@ -219,42 +203,6 @@ def realtime_zshim(fname_fmap, fname_mask_anat, fname_resp, fname_json, fname_an if DEBUG: - # Calculate masked shim for spherical harmonics plot - # masked_shimmed = np.zeros_like(shimmed) - # for i_t in range(nt): - # masked_shimmed[..., i_t] = mask_fmap * shimmed[..., i_t] - - # Plot unshimmed vs shimmed and their mask for spherical harmonics - # i_t = 0 - # fig = Figure(figsize=(10, 10)) - # ax = fig.add_subplot(2, 2, 1) - # im = ax.imshow(masked_fieldmaps[:-1, :-1, 0, i_t]) - # fig.colorbar(im) - # ax.set_title("Masked unshimmed") - # ax = fig.add_subplot(2, 2, 2) - # im = ax.imshow(masked_shimmed[:-1, :-1, 0, i_t]) - # fig.colorbar(im) - # ax.set_title("Masked shimmed") - # ax = fig.add_subplot(2, 2, 3) - # im = ax.imshow(fieldmap[:-1, :-1, 0, i_t]) - # fig.colorbar(im) - # ax.set_title("Unshimmed") - # ax = fig.add_subplot(2, 2, 4) - # im = ax.imshow(shimmed[:-1, :-1, 0, i_t]) - # fig.colorbar(im) - # ax.set_title("Shimmed") - # fname_figure = os.path.join(__dir_shimmingtoolbox__, 'fig_realtime_zshim_sphharm_shimmed.png') - # fig.savefig(fname_figure) - - # Plot the coil coefs through time - # fig = Figure(figsize=(10, 10)) - # for i_coil in range(n_coils): - # ax = fig.add_subplot(n_coils, 1, i_coil + 1) - # ax.plot(np.arange(nt), currents[i_coil, :]) - # ax.set_title(f"Channel {i_coil}") - # fname_figure = os.path.join(__dir_shimmingtoolbox__, 'fig_realtime_zshim_sphharm_currents.png') - # fig.savefig(fname_figure) - # Plot Static and RIRO fig = Figure(figsize=(10, 10)) ax = fig.add_subplot(2, 1, 1) @@ -268,27 +216,6 @@ def realtime_zshim(fname_fmap, fname_mask_anat, fname_resp, fname_json, fname_an fname_figure = os.path.join(fname_output, 'fig_realtime_zshim_riro_static.png') fig.savefig(fname_figure) - # Calculate fitted and shimmed for pressure fitted plot - # fitted_fieldmap = riro * (acq_pressures - mean_p) + static - # shimmed_pressure_fitted = np.expand_dims(fitted_fieldmap, 2) + masked_fieldmaps - # - # # Plot pressure fitted fieldmap - # fig = Figure(figsize=(10, 10)) - # ax = fig.add_subplot(3, 1, 1) - # im = ax.imshow(masked_fieldmaps[:-1, :-1, 0, i_t]) - # fig.colorbar(im) - # ax.set_title("fieldmap") - # ax = fig.add_subplot(3, 1, 2) - # im = ax.imshow(fitted_fieldmap[:-1, :-1, i_t]) - # fig.colorbar(im) - # ax.set_title("Fit") - # ax = fig.add_subplot(3, 1, 3) - # im = ax.imshow(shimmed_pressure_fitted[:-1, :-1, 0, i_t]) - # fig.colorbar(im) - # ax.set_title("Shimmed (fit + fieldmap") - # fname_figure = os.path.join(fname_output, 'fig_realtime_zshim_pressure_fitted.png') - # fig.savefig(fname_figure) - # Reshape pmu datapoints to fit those of the acquisition pmu_times = np.linspace(pmu.start_time_mdh, pmu.stop_time_mdh, len(pmu.data)) pmu_times_within_range = pmu_times[pmu_times > acq_timestamps[0]] From c51b82d56932fdfd3fcf73f176e16e830a5f1e91 Mon Sep 17 00:00:00 2001 From: Alexandre D'Astous Date: Wed, 25 Nov 2020 14:24:20 -0500 Subject: [PATCH 51/51] Rafactor cli into an API and a CLI --- shimmingtoolbox/cli/realtime_zshim.py | 212 ++--------------------- shimmingtoolbox/shim/realtime_zshim.py | 223 +++++++++++++++++++++++++ test/cli/test_cli_realtime_zshim.py | 26 +-- 3 files changed, 243 insertions(+), 218 deletions(-) create mode 100644 shimmingtoolbox/shim/realtime_zshim.py diff --git a/shimmingtoolbox/cli/realtime_zshim.py b/shimmingtoolbox/cli/realtime_zshim.py index 871e383b8..31a4c0983 100644 --- a/shimmingtoolbox/cli/realtime_zshim.py +++ b/shimmingtoolbox/cli/realtime_zshim.py @@ -6,6 +6,9 @@ import os import nibabel as nib import json + +from shimmingtoolbox.shim.realtime_zshim import realtime_zshim + from sklearn.linear_model import LinearRegression from nibabel.processing import resample_from_to # TODO: remove matplotlib and dirtesting import @@ -27,7 +30,7 @@ ) @click.option('-fmap', 'fname_fmap', required=True, type=click.Path(), help="B0 fieldmap. For realtime shimming, this should be a 4d file (4th dimension being time") -@click.option('-mask', 'fname_mask_anat', type=click.Path(), +@click.option('-mask', 'fname_mask_anat', type=click.Path(), required=False, help="3D nifti file with voxels between 0 and 1 used to weight the spatial region to shim. " "The coordinate system should be the same as ``anat``'s coordinate system.") @click.option('-resp', 'fname_resp', type=click.Path(), @@ -40,7 +43,7 @@ @click.option('-json', 'fname_json', type=click.Path(), help="Filename of json corresponding BIDS sidecar.") @click.option("-verbose", is_flag=True, help="Be more verbose.") -def realtime_zshim(fname_fmap, fname_mask_anat, fname_resp, fname_json, fname_anat, fname_output, verbose=True): +def realtime_zshim_cli(fname_fmap, fname_mask_anat, fname_resp, fname_json, fname_anat, fname_output, verbose=True): """ Args: @@ -56,225 +59,36 @@ def realtime_zshim(fname_fmap, fname_mask_anat, fname_resp, fname_json, fname_an """ - # Look if output directory exists, if not, create it - if not os.path.exists(fname_output): - os.makedirs(fname_output) - # Load fieldmap nii_fmap = nib.load(fname_fmap) - fieldmap = nii_fmap.get_fdata() - - # TODO: Error handling might move to API - if fieldmap.ndim != 4: - raise RuntimeError("fmap must be 4d (x, y, z, t)") - nx, ny, nz, nt = fieldmap.shape # Load anat nii_anat = nib.load(fname_anat) - anat = nii_anat.get_fdata() - if anat.ndim != 3: - raise RuntimeError("Anatomical image must be in 3d") - # Load mask - # TODO: check good practice below + # Load anatomical mask if fname_mask_anat is not None: nii_mask_anat = nib.load(fname_mask_anat) - if not np.all(np.isclose(nii_anat.affine, nii_mask_anat.affine)) or\ - not np.all(nii_mask_anat.shape == nii_anat.shape): - raise RuntimeError("Mask must have the same shape and affine transformation as anat") - nii_fmap_3d_temp = nib.Nifti1Image(fieldmap[..., 0], nii_fmap.affine) - nii_mask_fmap = resample_from_to(nii_mask_anat, nii_fmap_3d_temp) - mask_fmap = nii_mask_fmap.get_fdata() else: - mask_fmap = np.ones_like(fieldmap) - nii_mask_fmap = nib.Nifti1Image(mask_fmap, nii_anat.affine) - nii_mask_anat = nib.Nifti1Image(np.ones_like(anat), nii_anat.affine) + nii_mask_anat = None - if DEBUG: - nib.save(nii_mask_fmap, os.path.join(fname_output, 'tmp.mask_fmap_resample.nii.gz')) - - masked_fieldmaps = np.zeros_like(fieldmap) - for i_t in range(nt): - masked_fieldmaps[..., i_t] = mask_fmap * fieldmap[..., i_t] - - # Calculate gz gradient - g = 1000 / 42.576e6 # [mT / Hz] - gz_gradient = np.zeros_like(fieldmap) - # Get voxel coordinates. Z coordinates correspond to coord[2] - z_coord = generate_meshgrid(mask_fmap.shape, nii_fmap.affine)[2] / 1000 # [m] - - for it in range(nt): - gz_gradient[..., 0, it] = np.gradient(g * fieldmap[:, :, 0, it], z_coord[0, :, 0], axis=1) # [mT / m] - if DEBUG: - nii_gz_gradient = nib.Nifti1Image(gz_gradient, nii_fmap.affine) - nib.save(nii_gz_gradient, os.path.join(fname_output, 'tmp.gz_gradient.nii.gz')) - - # Fetch PMU timing # TODO: Add json to fieldmap instead of asking for another json file with open(fname_json) as json_file: json_data = json.load(json_file) - acq_timestamps = get_acquisition_times(nii_fmap, json_data) - pmu = PmuResp(fname_resp) - # TODO: deal with saturation - acq_pressures = pmu.interp_resp_trace(acq_timestamps) - - # TODO: - # fit PMU and fieldmap values - # do regression to separate static componant and RIRO component - # output coefficient with proper scaling - # field(i_vox) = a(i_vox) * (acq_pressures - mean_p) + b(i_vox) - # could use: https://scikit-learn.org/stable/modules/generated/sklearn.linear_model.LinearRegression.html - # Note: strong spatial autocorrelation on the a and b coefficients. Ie: two adjacent voxels are submitted to similar - # static B0 field and RIRO component. --> we need to find a way to account for that - # solution 1: post-fitting regularization. - # pros: easy to implement - # cons: fit is less robust to noise - # solution 2: accounting for regularization during fitting - # pros: fitting more robust to noise - # cons: (from Ryan): regularized fitting took a lot of time on Matlab - - # Shim using PMU - mean_p = np.mean(acq_pressures) - pressure_rms = np.sqrt(np.mean((acq_pressures - mean_p) ** 2)) - riro = np.zeros_like(fieldmap[:, :, :, 0]) - static = np.zeros_like(fieldmap[:, :, :, 0]) - # TODO fix progress bar not showing up - progress_bar = st_progress_bar(fieldmap[..., 0].size, desc="Fitting", ascii=False) - for i_x in range(fieldmap.shape[0]): - for i_y in range(fieldmap.shape[1]): - for i_z in range(fieldmap.shape[2]): - reg = LinearRegression().fit(acq_pressures.reshape(-1, 1) - mean_p, -gz_gradient[i_x, i_y, i_z, :]) - # Multiplying by the RMS of the pressure allows to make abstraction of the tightness of the bellow - # between scans. This allows to compare results between scans. - riro[i_x, i_y, i_z] = reg.coef_ * pressure_rms - static[i_x, i_y, i_z] = reg.intercept_ - progress_bar.update(1) - - # Resample masked_fieldmaps to target anatomical image - # TODO: convert to a function - masked_fmap_4d = np.zeros(anat.shape + (nt,)) - for it in range(nt): - nii_masked_fmap_3d = nib.Nifti1Image(masked_fieldmaps[..., it], nii_fmap.affine) - nii_resampled_fmap_3d = resample_from_to(nii_masked_fmap_3d, nii_anat, order=2, mode='nearest') - masked_fmap_4d[..., it] = nii_resampled_fmap_3d.get_fdata() - nii_resampled_fmap = nib.Nifti1Image(masked_fmap_4d, nii_anat.affine) - if DEBUG: - nib.save(nii_resampled_fmap, os.path.join(fname_output, 'resampled_fmap.nii.gz')) + static_correction, riro_correction, mean_p = realtime_zshim(nii_fmap, nii_anat, fname_resp, json_data, + nii_mask_anat=nii_mask_anat) - # Resample static to target anatomical image - nii_static = nib.Nifti1Image(static, nii_fmap.affine) - nii_resampled_static = resample_from_to(nii_static, nii_anat, mode='nearest') - nii_resampled_static_masked = nib.Nifti1Image(nii_resampled_static.get_fdata() * nii_mask_anat.get_fdata(), - nii_resampled_static.affine) - if DEBUG: - nib.save(nii_resampled_static_masked, os.path.join(fname_output, 'resampled_static.nii.gz')) - - # Resample riro to target anatomical image - nii_riro = nib.Nifti1Image(riro, nii_fmap.affine) - nii_resampled_riro = resample_from_to(nii_riro, nii_anat, mode='nearest') - nii_resampled_riro_masked = nib.Nifti1Image(nii_resampled_riro.get_fdata() * nii_mask_anat.get_fdata(), - nii_resampled_riro.affine) - if DEBUG: - nib.save(nii_resampled_riro_masked, os.path.join(fname_output, 'resampled_riro.nii.gz')) - - # Calculate the mean for riro and static for a perticular slice - n_slices = nii_anat.get_fdata().shape[2] - static_correction = np.zeros([n_slices]) - riro_correction = np.zeros([n_slices]) - for i_slice in range(n_slices): - ma_static_anat = np.ma.array(nii_resampled_static.get_fdata()[..., i_slice], - mask=nii_mask_anat.get_fdata()[..., i_slice] == False) - static_correction[i_slice] = np.ma.mean(ma_static_anat) - ma_riro_anat = np.ma.array(nii_resampled_riro.get_fdata()[..., i_slice], - mask=nii_mask_anat.get_fdata()[..., i_slice] == False) - riro_correction[i_slice] = np.ma.mean(ma_riro_anat) / pressure_rms + # Look if output directory exists, if not, create it + if not os.path.exists(fname_output): + os.makedirs(fname_output) # Write to a text file fname_corrections = os.path.join(fname_output, 'zshim_gradients.txt') file_gradients = open(fname_corrections, 'w') - for i_slice in range(n_slices): + for i_slice in range(static_correction.shape[-1]): file_gradients.write(f'Vector_Gz[0][{i_slice}]= {static_correction[i_slice]:.6f}\n') file_gradients.write(f'Vector_Gz[1][{i_slice}]= {riro_correction[i_slice]:.12f}\n') file_gradients.write(f'Vector_Gz[2][{i_slice}]= {mean_p:.3f}\n') file_gradients.close() - # ================ PLOTS ================ - - if DEBUG: - - # Plot Static and RIRO - fig = Figure(figsize=(10, 10)) - ax = fig.add_subplot(2, 1, 1) - im = ax.imshow(riro[:-1, :-1, 0] / pressure_rms) - fig.colorbar(im) - ax.set_title("RIRO") - ax = fig.add_subplot(2, 1, 2) - im = ax.imshow(static[:-1, :-1, 0]) - fig.colorbar(im) - ax.set_title("Static") - fname_figure = os.path.join(fname_output, 'fig_realtime_zshim_riro_static.png') - fig.savefig(fname_figure) - - # Reshape pmu datapoints to fit those of the acquisition - pmu_times = np.linspace(pmu.start_time_mdh, pmu.stop_time_mdh, len(pmu.data)) - pmu_times_within_range = pmu_times[pmu_times > acq_timestamps[0]] - pmu_data_within_range = pmu.data[pmu_times > acq_timestamps[0]] - pmu_data_within_range = pmu_data_within_range[pmu_times_within_range < acq_timestamps[fieldmap.shape[3] - 1]] - pmu_times_within_range = pmu_times_within_range[pmu_times_within_range < acq_timestamps[fieldmap.shape[3] - 1]] - - # Calc fieldmap average within mask - fieldmap_avg = np.zeros([fieldmap.shape[3]]) - for i_time in range(nt): - masked_array = np.ma.array(fieldmap[:, :, :, i_time], mask=mask_fmap == False) - fieldmap_avg[i_time] = np.ma.average(masked_array) - - # Plot pmu vs B0 in masked region - fig = Figure(figsize=(10, 10)) - ax = fig.add_subplot(211) - ax.plot(acq_timestamps / 1000, acq_pressures, label='Interpolated pressures') - # ax.plot(pmu_times / 1000, pmu.data, label='Raw pressures') - ax.plot(pmu_times_within_range / 1000, pmu_data_within_range, label='Pmu pressures') - ax.legend() - ax.set_title("Pressure [0, 4095] vs time (s) ") - ax = fig.add_subplot(212) - ax.plot(acq_timestamps / 1000, fieldmap_avg, label='Mean B0') - ax.legend() - ax.set_title("Fieldmap average over unmasked region (Hz) vs time (s)") - fname_figure = os.path.join(fname_output, 'fig_realtime_zshim_pmu_vs_B0.png') - fig.savefig(fname_figure) - - # Show anatomical image - fig = Figure(figsize=(10, 10)) - ax = fig.add_subplot(2, 1, 1) - im = ax.imshow(anat[:, :, 10]) - fig.colorbar(im) - ax.set_title("Anatomical image [:, :, 10]") - ax = fig.add_subplot(2, 1, 2) - im = ax.imshow(nii_mask_anat.get_fdata()[:, :, 10]) - fig.colorbar(im) - ax.set_title("Mask [:, :, 10]") - fname_figure = os.path.join(fname_output, 'fig_reatime_zshim_anat.png') - fig.savefig(fname_figure) - - # Show Gradient - fig = Figure(figsize=(10, 10)) - ax = fig.add_subplot(1, 1, 1) - im = ax.imshow(gz_gradient[:, :, 0, 0]) - fig.colorbar(im) - ax.set_title("Gradient [:, :, 0, 0]") - fname_figure = os.path.join(fname_output, 'fig_realtime_zshim_gradient.png') - fig.savefig(fname_figure) - - # Show evolution of coefficients - fig = Figure(figsize=(10, 10)) - ax = fig.add_subplot(2, 1, 1) - ax.plot(range(n_slices), static_correction, label='Static correction') - ax.set_title("Static correction evolution through slices") - ax = fig.add_subplot(2, 1, 2) - ax.plot(range(n_slices), (acq_pressures.max() - mean_p) * riro_correction, label='Riro correction') - ax.set_title("Riro correction evolution through slices") - fname_figure = os.path.join(fname_output, 'fig_realtime_zshim_correction_slice.png') - fig.savefig(fname_figure) - return fname_corrections diff --git a/shimmingtoolbox/shim/realtime_zshim.py b/shimmingtoolbox/shim/realtime_zshim.py new file mode 100644 index 000000000..c86f05fe9 --- /dev/null +++ b/shimmingtoolbox/shim/realtime_zshim.py @@ -0,0 +1,223 @@ +#!/usr/bin/python3 +# -*- coding: utf-8 -*- + +import numpy as np +import os +import nibabel as nib +from sklearn.linear_model import LinearRegression +from nibabel.processing import resample_from_to +# TODO: remove matplotlib and dirtesting import +from matplotlib.figure import Figure + +from shimmingtoolbox.load_nifti import get_acquisition_times +from shimmingtoolbox.pmu import PmuResp +from shimmingtoolbox.utils import st_progress_bar +from shimmingtoolbox.coils.coordinates import generate_meshgrid +from shimmingtoolbox import __dir_shimmingtoolbox__ + +DEBUG = True +fname_debug = os.path.join(__dir_shimmingtoolbox__, 'test_realtime_zshim') + +if not os.path.exists(fname_debug): + os.makedirs(fname_debug) + + +def realtime_zshim(nii_fieldmap, nii_anat, fname_resp, json_fmap, nii_mask_anat=None): + + # Make sure fieldmap has the appropriate dimensions + fieldmap = nii_fieldmap.get_fdata() + if fieldmap.ndim != 4: + raise RuntimeError("fmap must be 4d (x, y, z, t)") + nx, ny, nz, nt = nii_fieldmap.shape + + # Make sure anat has the appropriate dimensions + anat = nii_anat.get_fdata() + if anat.ndim != 3: + raise RuntimeError("Anatomical image must be in 3d") + + # Load mask + if nii_mask_anat is not None: + if not np.all(np.isclose(nii_anat.affine, nii_mask_anat.affine)) or \ + not np.all(nii_mask_anat.shape == nii_anat.shape): + raise RuntimeError("Mask must have the same shape and affine transformation as anat") + nii_fmap_3d_temp = nib.Nifti1Image(fieldmap[..., 0], nii_fieldmap.affine) + nii_mask_fmap = resample_from_to(nii_mask_anat, nii_fmap_3d_temp) + mask_fmap = nii_mask_fmap.get_fdata() + else: + mask_fmap = np.ones_like(fieldmap) + nii_mask_fmap = nib.Nifti1Image(mask_fmap, nii_anat.affine) + nii_mask_anat = nib.Nifti1Image(np.ones_like(anat), nii_anat.affine) + + if DEBUG: + nib.save(nii_mask_fmap, os.path.join(fname_debug, 'fig_mask_fmap.nii.gz')) + + masked_fieldmaps = np.zeros_like(fieldmap) + for i_t in range(nt): + masked_fieldmaps[..., i_t] = mask_fmap * fieldmap[..., i_t] + + # Calculate gz gradient + g = 1000 / 42.576e6 # [mT / Hz] + gz_gradient = np.zeros_like(fieldmap) + # Get voxel coordinates. Z coordinates correspond to coord[2] + z_coord = generate_meshgrid(mask_fmap.shape, nii_fieldmap.affine)[2] / 1000 # [m] + + for it in range(nt): + gz_gradient[..., 0, it] = np.gradient(g * fieldmap[:, :, 0, it], z_coord[0, :, 0], axis=1) # [mT / m] + if DEBUG: + nii_gz_gradient = nib.Nifti1Image(gz_gradient, nii_fieldmap.affine) + nib.save(nii_gz_gradient, os.path.join(fname_debug, 'fig_gz_gradient.nii.gz')) + + # Fetch PMU timing + acq_timestamps = get_acquisition_times(nii_fieldmap, json_fmap) + pmu = PmuResp(fname_resp) + # TODO: deal with saturation + # fit PMU and fieldmap values + acq_pressures = pmu.interp_resp_trace(acq_timestamps) + + # Shim using PMU + # field(i_vox) = a(i_vox) * (acq_pressures - mean_p) + b(i_vox) + # Note: strong spatial autocorrelation on the a and b coefficients. Ie: two adjacent voxels are submitted to similar + # static B0 field and RIRO component. --> we need to find a way to account for that + # solution 1: post-fitting regularization. + # pros: easy to implement + # cons: fit is less robust to noise + # solution 2: accounting for regularization during fitting + # pros: fitting more robust to noise + # cons: (from Ryan): regularized fitting took a lot of time on Matlab + mean_p = np.mean(acq_pressures) + pressure_rms = np.sqrt(np.mean((acq_pressures - mean_p) ** 2)) + riro = np.zeros_like(fieldmap[:, :, :, 0]) + static = np.zeros_like(fieldmap[:, :, :, 0]) + # TODO fix progress bar not showing up + progress_bar = st_progress_bar(fieldmap[..., 0].size, desc="Fitting", ascii=False) + for i_x in range(fieldmap.shape[0]): + for i_y in range(fieldmap.shape[1]): + for i_z in range(fieldmap.shape[2]): + # do regression to separate static componant and RIRO component + reg = LinearRegression().fit(acq_pressures.reshape(-1, 1) - mean_p, -gz_gradient[i_x, i_y, i_z, :]) + # Multiplying by the RMS of the pressure allows to make abstraction of the tightness of the bellow + # between scans. This allows to compare results between scans. + riro[i_x, i_y, i_z] = reg.coef_ * pressure_rms + static[i_x, i_y, i_z] = reg.intercept_ + progress_bar.update(1) + + # Resample masked_fieldmaps to target anatomical image + # TODO: convert to a function + masked_fmap_4d = np.zeros(anat.shape + (nt,)) + for it in range(nt): + nii_masked_fmap_3d = nib.Nifti1Image(masked_fieldmaps[..., it], nii_fieldmap.affine) + nii_resampled_fmap_3d = resample_from_to(nii_masked_fmap_3d, nii_anat, order=2, mode='nearest') + masked_fmap_4d[..., it] = nii_resampled_fmap_3d.get_fdata() + nii_resampled_fmap = nib.Nifti1Image(masked_fmap_4d, nii_anat.affine) + + if DEBUG: + nib.save(nii_resampled_fmap, os.path.join(fname_debug, 'fig_resampled_fmap.nii.gz')) + + # Resample static to target anatomical image + nii_static = nib.Nifti1Image(static, nii_fieldmap.affine) + nii_resampled_static = resample_from_to(nii_static, nii_anat, mode='nearest') + nii_resampled_static_masked = nib.Nifti1Image(nii_resampled_static.get_fdata() * nii_mask_anat.get_fdata(), + nii_resampled_static.affine) + if DEBUG: + nib.save(nii_resampled_static_masked, os.path.join(fname_debug, 'fig_resampled_static.nii.gz')) + + # Resample riro to target anatomical image + nii_riro = nib.Nifti1Image(riro, nii_fieldmap.affine) + nii_resampled_riro = resample_from_to(nii_riro, nii_anat, mode='nearest') + nii_resampled_riro_masked = nib.Nifti1Image(nii_resampled_riro.get_fdata() * nii_mask_anat.get_fdata(), + nii_resampled_riro.affine) + if DEBUG: + nib.save(nii_resampled_riro_masked, os.path.join(fname_debug, 'fig_resampled_riro.nii.gz')) + + # Calculate the mean for riro and static for a perticular slice + n_slices = nii_anat.get_fdata().shape[2] + static_correction = np.zeros([n_slices]) + riro_correction = np.zeros([n_slices]) + for i_slice in range(n_slices): + ma_static_anat = np.ma.array(nii_resampled_static.get_fdata()[..., i_slice], + mask=nii_mask_anat.get_fdata()[..., i_slice] == False) + static_correction[i_slice] = np.ma.mean(ma_static_anat) + ma_riro_anat = np.ma.array(nii_resampled_riro.get_fdata()[..., i_slice], + mask=nii_mask_anat.get_fdata()[..., i_slice] == False) + riro_correction[i_slice] = np.ma.mean(ma_riro_anat) / pressure_rms + + # ================ PLOTS ================ + + if DEBUG: + + # Plot Static and RIRO + fig = Figure(figsize=(10, 10)) + ax = fig.add_subplot(2, 1, 1) + im = ax.imshow(riro[:-1, :-1, 0] / pressure_rms) + fig.colorbar(im) + ax.set_title("RIRO") + ax = fig.add_subplot(2, 1, 2) + im = ax.imshow(static[:-1, :-1, 0]) + fig.colorbar(im) + ax.set_title("Static") + fname_figure = os.path.join(fname_debug, 'fig_realtime_zshim_riro_static.png') + fig.savefig(fname_figure) + + # Reshape pmu datapoints to fit those of the acquisition + pmu_times = np.linspace(pmu.start_time_mdh, pmu.stop_time_mdh, len(pmu.data)) + pmu_times_within_range = pmu_times[pmu_times > acq_timestamps[0]] + pmu_data_within_range = pmu.data[pmu_times > acq_timestamps[0]] + pmu_data_within_range = pmu_data_within_range[pmu_times_within_range < acq_timestamps[fieldmap.shape[3] - 1]] + pmu_times_within_range = pmu_times_within_range[pmu_times_within_range < acq_timestamps[fieldmap.shape[3] - 1]] + + # Calc fieldmap average within mask + fieldmap_avg = np.zeros([fieldmap.shape[3]]) + for i_time in range(nt): + masked_array = np.ma.array(fieldmap[:, :, :, i_time], mask=mask_fmap == False) + fieldmap_avg[i_time] = np.ma.average(masked_array) + + # Plot pmu vs B0 in masked region + fig = Figure(figsize=(10, 10)) + ax = fig.add_subplot(211) + ax.plot(acq_timestamps / 1000, acq_pressures, label='Interpolated pressures') + # ax.plot(pmu_times / 1000, pmu.data, label='Raw pressures') + ax.plot(pmu_times_within_range / 1000, pmu_data_within_range, label='Pmu pressures') + ax.legend() + ax.set_title("Pressure [0, 4095] vs time (s) ") + ax = fig.add_subplot(212) + ax.plot(acq_timestamps / 1000, fieldmap_avg, label='Mean B0') + ax.legend() + ax.set_title("Fieldmap average over unmasked region (Hz) vs time (s)") + fname_figure = os.path.join(fname_debug, 'fig_realtime_zshim_pmu_vs_B0.png') + fig.savefig(fname_figure) + + # Show anatomical image + fig = Figure(figsize=(10, 10)) + ax = fig.add_subplot(2, 1, 1) + im = ax.imshow(anat[:, :, 10]) + fig.colorbar(im) + ax.set_title("Anatomical image [:, :, 10]") + ax = fig.add_subplot(2, 1, 2) + im = ax.imshow(nii_mask_anat.get_fdata()[:, :, 10]) + fig.colorbar(im) + ax.set_title("Mask [:, :, 10]") + fname_figure = os.path.join(fname_debug, 'fig_reatime_zshim_anat.png') + fig.savefig(fname_figure) + + # Show Gradient + fig = Figure(figsize=(10, 10)) + ax = fig.add_subplot(1, 1, 1) + im = ax.imshow(gz_gradient[:, :, 0, 0]) + fig.colorbar(im) + ax.set_title("Gradient [:, :, 0, 0]") + fname_figure = os.path.join(fname_debug, 'fig_realtime_zshim_gradient.png') + fig.savefig(fname_figure) + + # Show evolution of coefficients + fig = Figure(figsize=(10, 10)) + ax = fig.add_subplot(2, 1, 1) + ax.plot(range(n_slices), static_correction, label='Static correction') + ax.set_title("Static correction evolution through slices") + ax = fig.add_subplot(2, 1, 2) + ax.plot(range(n_slices), (acq_pressures.max() - mean_p) * riro_correction, label='Riro correction') + ax.set_title("Riro correction evolution through slices") + fname_figure = os.path.join(fname_debug, 'fig_realtime_zshim_correction_slice.png') + fig.savefig(fname_figure) + + # TODO: output pressure_rms to scale for interperson testing + return static_correction, riro_correction, mean_p diff --git a/test/cli/test_cli_realtime_zshim.py b/test/cli/test_cli_realtime_zshim.py index 808ae0c34..8e4886acc 100644 --- a/test/cli/test_cli_realtime_zshim.py +++ b/test/cli/test_cli_realtime_zshim.py @@ -8,11 +8,8 @@ import numpy as np from click.testing import CliRunner -from shimmingtoolbox.cli.realtime_zshim import realtime_zshim +from shimmingtoolbox.cli.realtime_zshim import realtime_zshim_cli from shimmingtoolbox.masking.shapes import shapes -from shimmingtoolbox.masking.threshold import threshold -from shimmingtoolbox.coils.coordinates import generate_meshgrid -from shimmingtoolbox.coils.siemens_basis import siemens_basis from shimmingtoolbox import __dir_testing__ from shimmingtoolbox import __dir_shimmingtoolbox__ @@ -48,15 +45,6 @@ def test_cli_realtime_zshim(): nib.save(nii_mask, fname_mask) # fname_mask = os.path.join(__dir_testing__, 'test_realtime_zshim', 'gre_seg_30.nii') - - # Set up coils - # coord_phys = generate_meshgrid(nii_fmap.get_fdata().shape[0:3], nii_fmap.affine) - # coil_profile = siemens_basis(coord_phys[0], coord_phys[1], coord_phys[2]) - # - # nii_coil = nib.Nifti1Image(coil_profile, nii_fmap.affine) - # fname_coil = os.path.join(tmp, 'coil_profile.nii.gz') - # nib.save(nii_coil, fname_coil) - # Path for resp data fname_resp = os.path.join(__dir_testing__, 'realtime_zshimming_data', 'PMUresp_signal.resp') @@ -68,12 +56,12 @@ def test_cli_realtime_zshim(): fname_output = os.path.join(__dir_shimmingtoolbox__, 'test_realtime_zshim') # Run the CLI - result = runner.invoke(realtime_zshim, ['-fmap', fname_fieldmap, - '-mask', fname_mask, - '-output', fname_output, - '-resp', fname_resp, - '-json', fname_json, - '-anat', fname_anat], + result = runner.invoke(realtime_zshim_cli, ['-fmap', fname_fieldmap, + '-mask', fname_mask, + '-output', fname_output, + '-resp', fname_resp, + '-json', fname_json, + '-anat', fname_anat], catch_exceptions=False) assert result.exit_code == 0