Resampled mask in Python is zero on the last slice along Z axis

[tashrifbillah]

Resampled mask in Python is zero on the last slice along Z axis

Experiment performed on the following images

Suheyla MATLAB

/data/pnl/home/sq566/Suheyla_Harmonization/target_Petra/training_x/harmMATLAB_Oct4/c028dwi-xc-QC-Ed_bMapped_Resampled.nii.gz

/data/pnl/home/sq566/Suheyla_Harmonization/target_Petra/training_x/harmMATLAB_Oct4/c028_OTSUtensormask_Resampled.nii.gz

Tashrif Python

/data/pnl/home/sq566/Suheyla_Harmonization/target_Petra/training_x/c028dwi-xc-QC-Ed_resampled.nii.gz

/data/pnl/home/sq566/Suheyla_Harmonization/target_Petra/training_y/c028_OTSUtensormask_resampled.nii.gz

This zero slice makes subsequent L0 and FA features zero on that particular slicer.

@suheyla2 and @yrathi believe this error accounts for the remaining difference between MATLAB and PYTHON. @tashrifbillah argued that the anomaly is same for all cases and all sites. @yrathi refuted, the zero slice triggers mis-registration between two cases where one is affected by this anomaly and the other is not. @suheyla2 observed this anomaly in some cases where rounding during grid creation affects resampling.

In particular, scipy resize package is used to peform mask resampling in Python

https://github.com/pnlbwh/dMRIharmonization/blob/spm-bspline/lib/resampling.py#L106

highResMask= resize(lowResMask.astype('float'), (sx, sy, sz), order= 1, mode= 'constant') # order 1 for linear interpolation

On the other hand, MATLAB uses interp3 function:

https://github.com/pnlbwh/harmonization/blob/master/Utilities/harmonizationFunctions/resampleData.m#L25

highResMask = interp3(1:sz_y, 1:sz_x, 1:sz_z,double( lowResMask), y, x, z,'linear');

Possible solution:

Other than looking for equivalent resampling option in Python (according to our experiment before, there is none), we better use MATLAB executable for mask resampling as well. Just include the mask resampling code segment in https://github.com/pnlbwh/dMRIharmonization/blob/spm-bspline/lib/spm_bspline_exec/bspline.m

cc: @sbouix

[tashrifbillah]

Proposed solution 1:

use spm_bspline with order=1 (linear interpolation):

# resample the mask ---------------------------------------------------------------
inPrefix= lowResMaskPath.split('.')[0]
savemat(inPrefix+'_sp.mat', {'sp_high':sp_high,'sp_low':sp_low, 'imgDim':lowResImg.shape[:3], 'sOrder':1})
highResMask= resize_spm(lowResMask, lowResMaskPath.split('.')[0])

Most of the difference is gone, while some is observed in the skull region, again on the last slice along the Z axis.

Double checked

[tashrifbillah]

Proposed solution 2:

Use reflect mode with scipy interpolation:

highResMask= resize(lowResMask.astype('float'), (sx, sy, sz), order= 1, mode= 'reflect')

There is more difference than that of spm_bspline. Difference is observed along all axes.

[tashrifbillah]

Proposed solution 3:

Use symmetric mode with scipy interpolation:

highResMask= resize(lowResMask.astype('float'), (sx, sy, sz), order= 1, mode= 'symmetric')

Same difference as that of reflect mode.

Guess: the difference in reflect and symmetric modes are stemming from morphological operation performed in Python.

[tashrifbillah]

Decision: accept spm_bspline(order=1) as the suitable interpolation technique for mask interpolation.

[tashrifbillah]

The following commit implemented the above decision:

https://github.com/pnlbwh/dMRIharmonization/commit/7b106d4a9fce0863ed75c1aafc8ddc1d7e777eb8

[tashrifbillah]

@tashrifbillah will merge the pull request soon. Keeping it stand by for a day in case we want to make more changes:

https://github.com/pnlbwh/dMRIharmonization/pull/30

[tashrifbillah]

Merged PR #30 into spm-bspline branch.

[tashrifbillah]

Issue marked for closing