-
-
Notifications
You must be signed in to change notification settings - Fork 101
New issue
Have a question about this project? Sign up for a free GitHub account to open an issue and contact its maintainers and the community.
By clicking “Sign up for GitHub”, you agree to our terms of service and privacy statement. We’ll occasionally send you account related emails.
Already on GitHub? Sign in to your account
sct_fmri_moco: Generalize motion correction for sagittal acquisitions and other improvements #2022
Conversation
There was a problem hiding this comment.
Choose a reason for hiding this comment
The reason will be displayed to describe this comment to others. Learn more.
Just a question about concurrent branches
scripts/sct_create_mask.py
Outdated
@@ -110,12 +101,11 @@ def create_mask(param): | |||
|
|||
path_tmp = sct.tmp_create(basename="create_mask", verbose=param.verbose) | |||
|
|||
sct.printv('\nCheck orientation...', param.verbose) | |||
sct.printv('\nOrientation:', param.verbose) |
There was a problem hiding this comment.
Choose a reason for hiding this comment
The reason will be displayed to describe this comment to others. Learn more.
Hmmm, which PR do you want to merge first?
There was a problem hiding this comment.
Choose a reason for hiding this comment
The reason will be displayed to describe this comment to others. Learn more.
4eefb4f
to
604c68e
Compare
This fixes issue #2017
Oherwise 4D data with singleton will become 3D data
Now, if the splitting dimension corresponds to the last dim of the image, then it reshapes to remove the last dim (e.g. 4d --> 3d). Otherwise, it keeps the singleton.
For faster computation
This is sometimes useful in case singletons need to be removed
via param squeeze_data (True by default to keep current behavior)
These modifications include the possibility to split along Z in case the input image is in sagittal orientation. The rest of the code is untouched, so that the motion correction as implemented previously is now encapsulated inside the z-splitted loop. UNFINISHED
need to be singleton
To avoid incessant warning messages.
After several unsuccessful testing using antsSliceRegularizedRegistration in 2D for moco applied to sagittal images, i realized that using antsRegistration in 2D mode, with affine transformation, gave best results. UNFINISHED: need to take care of the .mat output when using -g different from 1 (or maybe we could avoid supporting it)
Did some refactoring, notably made is_sagittal a field in param structure.
There are still a few checks to do, notably make sure the code runs when using mask with g>1 in 2D mode
...Which were introduced with the recent changes to msct_moco
Some cleaning, and miminized verbose (because it conflicts with the recent introduction of tqdm).
7efa2bd
to
272d3e2
Compare
scripts/sct_deepseg_sc.py
Outdated
|
||
# compute z_centerline in image coordinates with continuous precision | ||
voxel_coordinates = image_labels.transfo_phys2pix([[x_centerline_fit_rescorr[i], y_centerline_fit_rescorr[i], z_centerline_rescorr[i]] for i in range(len(z_centerline_rescorr))], real=False) | ||
x_centerline_voxel_cont = [coord[0] for coord in voxel_coordinates] |
There was a problem hiding this comment.
Choose a reason for hiding this comment
The reason will be displayed to describe this comment to others. Learn more.
Note that since transfo_phys2pix now returns a numpy array, you could have used something like x, y, z = voxel_coordinates.T
scripts/sct_dmri_moco.py
Outdated
@@ -369,36 +372,30 @@ def dmri_moco(param): | |||
|
|||
# DWI groups | |||
file_dwi_mean = [] | |||
tqdm_bar = tqdm(total=nb_groups, unit='iter', unit_scale=False, desc="Merge within groups", ascii=True, ncols=80) | |||
for iGroup in range(nb_groups): |
There was a problem hiding this comment.
Choose a reason for hiding this comment
The reason will be displayed to describe this comment to others. Learn more.
You could have used a tqdm iterator object instead of the range() here (avoids to update()
and close()
)
There was a problem hiding this comment.
Choose a reason for hiding this comment
The reason will be displayed to describe this comment to others. Learn more.
i was not aware of that-- i'll look into it
There was a problem hiding this comment.
Choose a reason for hiding this comment
The reason will be displayed to describe this comment to others. Learn more.
fixed in c68600d
scripts/sct_fmri_moco.py
Outdated
import time | ||
import math | ||
|
||
import shutil |
There was a problem hiding this comment.
Choose a reason for hiding this comment
The reason will be displayed to describe this comment to others. Learn more.
There is a copy()
function in sctutils
There was a problem hiding this comment.
Choose a reason for hiding this comment
The reason will be displayed to describe this comment to others. Learn more.
Fixed in 44e02e4
im_maskz.save() | ||
file_mask_splitZ.append(im_maskz.absolutepath) | ||
# initialize file list for output matrices | ||
file_mat = np.chararray([nz, nt], |
There was a problem hiding this comment.
Choose a reason for hiding this comment
The reason will be displayed to describe this comment to others. Learn more.
chararray
is deprecated and also not a good idea
https://docs.scipy.org/doc/numpy/reference/generated/numpy.chararray.html
There was a problem hiding this comment.
Choose a reason for hiding this comment
The reason will be displayed to describe this comment to others. Learn more.
after spending 10min on this i need help-- i've tried
file_mat = np.array([nz, nt], dtype="str")
but pb because string is immutable. solution?
scripts/msct_moco.py
Outdated
# Replace failed transformation with the closest good one | ||
fT = [i for i, j in enumerate(failed_transfo) if j == 1] | ||
gT = [i for i, j in enumerate(failed_transfo) if j == 0] | ||
for it in range(len(fT)): |
There was a problem hiding this comment.
Choose a reason for hiding this comment
The reason will be displayed to describe this comment to others. Learn more.
range(len(x))
-> enumerate(fT)
There was a problem hiding this comment.
Choose a reason for hiding this comment
The reason will be displayed to describe this comment to others. Learn more.
Fixed in 14f6794
scripts/msct_moco.py
Outdated
fT = [i for i, j in enumerate(failed_transfo) if j == 1] | ||
gT = [i for i, j in enumerate(failed_transfo) if j == 0] | ||
for it in range(len(fT)): | ||
abs_dist = [abs(gT[i] - fT[it]) for i in range(len(gT))] |
There was a problem hiding this comment.
Choose a reason for hiding this comment
The reason will be displayed to describe this comment to others. Learn more.
Use numpy?
There was a problem hiding this comment.
Choose a reason for hiding this comment
The reason will be displayed to describe this comment to others. Learn more.
fixed in 7b3647c
scripts/sct_fmri_moco.py
Outdated
# Copy registration matrices | ||
sct.printv('\nCopy transformations...', param.verbose) | ||
for iGroup in range(nb_groups): | ||
for data in range(len(group_indexes[iGroup])): |
There was a problem hiding this comment.
Choose a reason for hiding this comment
The reason will be displayed to describe this comment to others. Learn more.
enumerate
There was a problem hiding this comment.
Choose a reason for hiding this comment
The reason will be displayed to describe this comment to others. Learn more.
fixed in ef52dd3
There was a problem hiding this comment.
Choose a reason for hiding this comment
The reason will be displayed to describe this comment to others. Learn more.
actually we can't here-- fixed in c9aaa84
im_out.absolutepath = sct.add_suffix(im_in.absolutepath, "_{}{}".format(dim_list[dim].upper(), str(idx_img).zfill(4))) | ||
im_out_list.append(im_out) | ||
|
||
return im_out_list | ||
|
||
|
||
def concat_data(fname_in_list, dim, pixdim=None): | ||
def concat_data(fname_in_list, dim, pixdim=None, squeeze_data=False): |
There was a problem hiding this comment.
Choose a reason for hiding this comment
The reason will be displayed to describe this comment to others. Learn more.
Why do you absolutely want to have squeezing performed within this function and not elsewhere (caller or otherwise)?
There was a problem hiding this comment.
Choose a reason for hiding this comment
The reason will be displayed to describe this comment to others. Learn more.
squeeze=False, so I guess it is not performed-- or am i missing something?
There was a problem hiding this comment.
Choose a reason for hiding this comment
The reason will be displayed to describe this comment to others. Learn more.
just a question of design. why would a function named concat_data perform a squeeze.
We cannot use enumerate because group_indexes has 2 dim. So changing it back to the uglier range(len)
c9aaa84
to
d784521
Compare
Otherwise I get the following error when running sct_testing: UnicodeEncodeError: 'ascii' codec can't encode characters in position 27-28: ordinal not in range(128)
sct_fmri_moco: Generalize motion correction for sagittal acquisitions and other improvements Former-commit-id: 55756cf
Currently, sct_fmri_moco is designed such that the input image is assumed to be in axial orientation (i.e. 3rd dimension is S-I), and the algorithm estimates a set of Tx,Ty for each slice, regularized using polynom functions along the S-I direction.
However, for images acquired sagitally (or coronally), this approach does not work anymore. Beyond the simple orientation issue described in #1983, the type of motion-related corruption is not addressed by the current algorithm. More specifically, if motion occurs between sagittal slices, then the cord is non-linearly distorted in the axial plane and a simple translation within the axial plane is not sufficient. This is shown in the animation below. Without motion correction (left), a few slices are translated in the AP direction, leaving ugly distortion of the cord in the axial plane. With motion correction (right), the algorithm tries to find an average translation, but this is of course not sufficient for this type of motion.
The implemented solution is to split the data along Z (R-L dimension) and use
antsRegistration
in 2D mode for each slice independently. For now we use affine transfo, which is fast and robust, but later we could investigate the use of other transfo. What helps in particular is to use a mask to focus on a region of interest for registration (esp. because motion is non-affine). See below a comparison without (left) and with the new moco algorithm (right):This PR also provides a few other fixes, refactoring and improvements (runs faster).
Fixes #1984, Fixes #820, Fixes #2017, Fixes #1878