Problem with dwi2response msmt_5tt (amp2response) command for ex-vivo monkey data

1

Dear MRtrix experts,

I’m quite new to MRtrix and I’m trying to perform all the steps towards ACT for post-mortem multiple-shell monkey data (3 b=0 and b=250, 800, 1500, 4000 and 7000 s/mm²) but I’m stuck when dealing with the command dwi2response with the msmt_5tt algorithm. The 5tt file comes from a previous segmentation and I ordered it in the right way to feed the command.

It seems that the first run for the white matter went smoothly for all the shells but when it gets to the grey matter, it’s taking a lot more time than expected. When I display information with the -info or -debug option, everything is ok until the amp2response for GM (amp2response dwi.mif gm_mask.mif dirs.mif gm.txt -shell 0,250,800,1500,4000,7000 -isotropic) where I arrived to the 3 messages for the firsts shells (“shell b=0: Response function […] from average of 1564410 voxels” and so on for the two others) after 3 days of computation; and it’s still running. I thought it could be an infinite loop somewhere but it’s displaying some messages, so could it just be really long because of the averaging over more than a million voxels for each shell or could there be an error somewhere else? Will the solution might be to come back to the manual algorithm, and if so, do you any tips on how to tune the parameters?

For information, to compare I’ve also tried to launch the same command on the HCP and on one human multiple-shell (1 b=0, b=1000 and b=3000 s/mm²) datasets where it’s also taking that time at the same step…

Thanks in advance for any help you could provide.

Best regards,

Achille

2

OK, I don’t think it’s supposed to be using that many voxels for the response estimation - certainly not with the msmt_5tt algorithm. Can you post the exact command that is causing the issue, including all of its output?

3

Hi @jdtournier,

thank you for the quick response. As requested, you can find the line command and all of the outputs below :

dwi2response msmt_5tt -voxels voxels.mif -info ../01-MifConversion/dwi.mif ../01-MifConversion/5TT/5TT.mif WM_response.txt GM_response.txt CSF_response.txt
dwi2response: 
dwi2response: Note that this script makes use of commands / algorithms that have relevant articles for citation. Please consult the help page (-help option) for more information.
dwi2response: 
dwi2response: Command: '/home/ateillac/Softs/MRTRIX/mrtrix3/bin/mrinfo /scratch/ateillac/NHP/M157_5shells/Essai/01-MifConversion/dwi.mif -dwgrad' (piping data to local storage)
dwi2response: Result: (218 lines)
dwi2response: Generated temporary directory: /scratch/ateillac/NHP/M157_5shells/Essai/02-TissueResponse/dwi2response-tmp-DRG9XQ/
Command:  mrconvert /scratch/ateillac/NHP/M157_5shells/Essai/01-MifConversion/dwi.mif /scratch/ateillac/NHP/M157_5shells/Essai/02-TissueResponse/dwi2response-tmp-DRG9XQ/dwi.mif -stride 0,0,0,1
mrconvert: [100%] copying from "/scratch/a...ai/01-MifConversion/dwi.mif" to "/scratch/a...response-tmp-DRG9XQ/dwi.mif"
Command:  mrconvert /scratch/ateillac/NHP/M157_5shells/Essai/01-MifConversion/5TT/5TT.mif /scratch/ateillac/NHP/M157_5shells/Essai/02-TissueResponse/dwi2response-tmp-DRG9XQ/5tt.mif
mrconvert: [100%] copying from "/scratch/a...1-MifConversion/5TT/5TT.mif" to "/scratch/a...response-tmp-DRG9XQ/5tt.mif"
dwi2response: Changing to temporary directory (/scratch/ateillac/NHP/M157_5shells/Essai/02-TissueResponse/dwi2response-tmp-DRG9XQ/)
Command:  dwi2mask dwi.mif mask.mif
dwi2mask: [100%] finding min/max of "mean b=0 image"
dwi2mask: [done] optimising threshold
dwi2mask: [100%] thresholding
dwi2mask: [100%] finding min/max of "mean b=244 image"
dwi2mask: [done] optimising threshold
dwi2mask: [100%] thresholding
dwi2mask: [100%] finding min/max of "mean b=800 image"
dwi2mask: [done] optimising threshold
dwi2mask: [100%] thresholding
dwi2mask: [100%] finding min/max of "mean b=1536 image"
dwi2mask: [done] optimising threshold
dwi2mask: [100%] thresholding
dwi2mask: [100%] finding min/max of "mean b=4075 image"
dwi2mask: [done] optimising threshold
dwi2mask: [100%] thresholding
dwi2mask: [100%] finding min/max of "mean b=7106 image"
dwi2mask: [done] optimising threshold
dwi2mask: [100%] thresholding
dwi2mask: [done] computing dwi brain mask
dwi2mask: [done] applying mask cleaning filter
Command:  5ttcheck 5tt.mif
dwi2response: Command: '/home/ateillac/Softs/MRTRIX/mrtrix3/bin/mrinfo dwi.mif -shells' (piping data to local storage)
dwi2response: Result: 0 243.7 800.471 1536.39 4075.06 7106.2
Command:  dwi2tensor dwi.mif - -mask mask.mif | tensor2metric - -fa fa.mif -vector vector.mif
dwi2tensor: [100%] computing tensors
tensor2metric: [100%] computing metrics
Function: shutil.copy(vector.mif, dirs.mif)
Command:  mrtransform 5tt.mif 5tt_regrid.mif -template fa.mif -interp linear
mrtransform: [100%] reslicing "5tt.mif"
Command:  mrconvert 5tt_regrid.mif - -coord 3 2 -axes 0,1,2 | mrcalc - 0.95 -gt mask.mif -mult wm_mask.mif
mrconvert: [100%] copying from "5tt_regrid.mif" to "/tmp/mrtrix-tmp-2Kvuz0.mif"
mrcalc: [100%] computing: ((/tmp/mrtrix-tmp-2Kvuz0.mif > 0.95) * mask.mif)
Command:  mrconvert 5tt_regrid.mif - -coord 3 0 -axes 0,1,2 | mrcalc - 0.95 -gt fa.mif 0.2 -lt -mult mask.mif -mult gm_mask.mif
mrconvert: [100%] copying from "5tt_regrid.mif" to "/tmp/mrtrix-tmp-2Kvuz0.mif"
mrcalc: [100%] computing: (((/tmp/mrtrix-tmp-2Kvuz0.mif > 0.95) * (fa.mif < 0.2)) * mask.mif)
Command:  mrconvert 5tt_regrid.mif - -coord 3 3 -axes 0,1,2 | mrcalc - 0.95 -gt fa.mif 0.2 -lt -mult mask.mif -mult csf_mask.mif
mrconvert: [100%] copying from "5tt_regrid.mif" to "/tmp/mrtrix-tmp-2Kvuz0.mif"
mrcalc: [100%] computing: (((/tmp/mrtrix-tmp-2Kvuz0.mif > 0.95) * (fa.mif < 0.2)) * mask.mif)
dwi2response: Calling dwi2response recursively to select WM single-fibre voxels using 'tournier' algorithm
Command:  dwi2response tournier dwi.mif wm_ss_response.txt -mask wm_mask.mif -voxels wm_sf_mask.mif -tempdir /scratch/ateillac/NHP/M157_5shells/Essai/02-TissueResponse/dwi2response-tmp-DRG9XQ/
dwi2response: 
dwi2response: Note that this script makes use of commands / algorithms that have relevant articles for citation. Please consult the help page (-help option) for more information.
dwi2response: 
dwi2response: Generated temporary directory: /scratch/ateillac/NHP/M157_5shells/Essai/02-TissueResponse/dwi2response-tmp-DRG9XQ/dwi2response-tmp-BWAOF1/
Command:  mrconvert /scratch/ateillac/NHP/M157_5shells/Essai/02-TissueResponse/dwi2response-tmp-DRG9XQ/dwi.mif - -stride 0,0,0,1 | dwiextract - /scratch/ateillac/NHP/M157_5shells/Essai/02-TissueResponse/dwi2response-tmp-DRG9XQ/dwi2response-tmp-BWAOF1/dwi.mif -singleshell -no_bzero
Command:  mrconvert /scratch/ateillac/NHP/M157_5shells/Essai/02-TissueResponse/dwi2response-tmp-DRG9XQ/wm_mask.mif /scratch/ateillac/NHP/M157_5shells/Essai/02-TissueResponse/dwi2response-tmp-DRG9XQ/dwi2response-tmp-BWAOF1/mask.mif -datatype bit
dwi2response: Changing to temporary directory (/scratch/ateillac/NHP/M157_5shells/Essai/02-TissueResponse/dwi2response-tmp-DRG9XQ/dwi2response-tmp-BWAOF1/)
Command:  dwi2fod csd dwi.mif init_RF.txt iter0_FOD.mif -mask mask.mif -lmax 4
Command:  fod2fixel iter0_FOD.mif iter0_fixel -peak peaks.mif -mask mask.mif -fmls_no_thresholds
Command:  fixel2voxel iter0_fixel/peaks.mif split_data iter0_amps.mif -number 2
Command:  mrconvert iter0_amps.mif iter0_first_peaks.mif -coord 3 0 -axes 0,1,2
Command:  mrconvert iter0_amps.mif iter0_second_peaks.mif -coord 3 1 -axes 0,1,2
Command:  fixel2voxel iter0_fixel/directions.mif split_dir iter0_all_dirs.mif -number 1
Command:  mrconvert iter0_all_dirs.mif iter0_first_dir.mif -coord 3 0:2
Command:  mrcalc iter0_first_peaks.mif -sqrt 1 iter0_second_peaks.mif iter0_first_peaks.mif -div -sub 2 -pow -mult iter0_CF.mif
Command:  mrthreshold iter0_CF.mif -top 300 iter0_SF.mif
Command:  amp2response dwi.mif iter0_SF.mif iter0_first_dir.mif iter0_RF.txt -lmax 4
Command:  mrthreshold iter0_CF.mif -top 3000 - | maskfilter - dilate - -npass 1 | mrcalc mask.mif - -mult iter0_SF_dilated.mif
Command:  dwi2fod csd dwi.mif iter0_RF.txt iter1_FOD.mif -mask iter0_SF_dilated.mif
Command:  fod2fixel iter1_FOD.mif iter1_fixel -peak peaks.mif -mask iter0_SF_dilated.mif -fmls_no_thresholds
Command:  fixel2voxel iter1_fixel/peaks.mif split_data iter1_amps.mif -number 2
Command:  mrconvert iter1_amps.mif iter1_first_peaks.mif -coord 3 0 -axes 0,1,2
Command:  mrconvert iter1_amps.mif iter1_second_peaks.mif -coord 3 1 -axes 0,1,2
Command:  fixel2voxel iter1_fixel/directions.mif split_dir iter1_all_dirs.mif -number 1
Command:  mrconvert iter1_all_dirs.mif iter1_first_dir.mif -coord 3 0:2
Command:  mrcalc iter1_first_peaks.mif -sqrt 1 iter1_second_peaks.mif iter1_first_peaks.mif -div -sub 2 -pow -mult iter1_CF.mif
Command:  mrthreshold iter1_CF.mif -top 300 iter1_SF.mif
Command:  amp2response dwi.mif iter1_SF.mif iter1_first_dir.mif iter1_RF.txt
Command:  mrcalc iter1_SF.mif iter0_SF.mif -sub iter1_SF_diff.mif
Command:  mrthreshold iter1_CF.mif -top 3000 - | maskfilter - dilate - -npass 1 | mrcalc mask.mif - -mult iter1_SF_dilated.mif
Command:  dwi2fod csd dwi.mif iter1_RF.txt iter2_FOD.mif -mask iter1_SF_dilated.mif
Command:  fod2fixel iter2_FOD.mif iter2_fixel -peak peaks.mif -mask iter1_SF_dilated.mif -fmls_no_thresholds
Command:  fixel2voxel iter2_fixel/peaks.mif split_data iter2_amps.mif -number 2
Command:  mrconvert iter2_amps.mif iter2_first_peaks.mif -coord 3 0 -axes 0,1,2
Command:  mrconvert iter2_amps.mif iter2_second_peaks.mif -coord 3 1 -axes 0,1,2
Command:  fixel2voxel iter2_fixel/directions.mif split_dir iter2_all_dirs.mif -number 1
Command:  mrconvert iter2_all_dirs.mif iter2_first_dir.mif -coord 3 0:2
Command:  mrcalc iter2_first_peaks.mif -sqrt 1 iter2_second_peaks.mif iter2_first_peaks.mif -div -sub 2 -pow -mult iter2_CF.mif
Command:  mrthreshold iter2_CF.mif -top 300 iter2_SF.mif
Command:  amp2response dwi.mif iter2_SF.mif iter2_first_dir.mif iter2_RF.txt
Command:  mrcalc iter2_SF.mif iter1_SF.mif -sub iter2_SF_diff.mif
Command:  mrthreshold iter2_CF.mif -top 3000 - | maskfilter - dilate - -npass 1 | mrcalc mask.mif - -mult iter2_SF_dilated.mif
Command:  dwi2fod csd dwi.mif iter2_RF.txt iter3_FOD.mif -mask iter2_SF_dilated.mif
Command:  fod2fixel iter3_FOD.mif iter3_fixel -peak peaks.mif -mask iter2_SF_dilated.mif -fmls_no_thresholds
Command:  fixel2voxel iter3_fixel/peaks.mif split_data iter3_amps.mif -number 2
Command:  mrconvert iter3_amps.mif iter3_first_peaks.mif -coord 3 0 -axes 0,1,2
Command:  mrconvert iter3_amps.mif iter3_second_peaks.mif -coord 3 1 -axes 0,1,2
Command:  fixel2voxel iter3_fixel/directions.mif split_dir iter3_all_dirs.mif -number 1
Command:  mrconvert iter3_all_dirs.mif iter3_first_dir.mif -coord 3 0:2
Command:  mrcalc iter3_first_peaks.mif -sqrt 1 iter3_second_peaks.mif iter3_first_peaks.mif -div -sub 2 -pow -mult iter3_CF.mif
Command:  mrthreshold iter3_CF.mif -top 300 iter3_SF.mif
Command:  amp2response dwi.mif iter3_SF.mif iter3_first_dir.mif iter3_RF.txt
Command:  mrcalc iter3_SF.mif iter2_SF.mif -sub iter3_SF_diff.mif
Command:  mrthreshold iter3_CF.mif -top 3000 - | maskfilter - dilate - -npass 1 | mrcalc mask.mif - -mult iter3_SF_dilated.mif
Command:  dwi2fod csd dwi.mif iter3_RF.txt iter4_FOD.mif -mask iter3_SF_dilated.mif
Command:  fod2fixel iter4_FOD.mif iter4_fixel -peak peaks.mif -mask iter3_SF_dilated.mif -fmls_no_thresholds
Command:  fixel2voxel iter4_fixel/peaks.mif split_data iter4_amps.mif -number 2
Command:  mrconvert iter4_amps.mif iter4_first_peaks.mif -coord 3 0 -axes 0,1,2
Command:  mrconvert iter4_amps.mif iter4_second_peaks.mif -coord 3 1 -axes 0,1,2
Command:  fixel2voxel iter4_fixel/directions.mif split_dir iter4_all_dirs.mif -number 1
Command:  mrconvert iter4_all_dirs.mif iter4_first_dir.mif -coord 3 0:2
Command:  mrcalc iter4_first_peaks.mif -sqrt 1 iter4_second_peaks.mif iter4_first_peaks.mif -div -sub 2 -pow -mult iter4_CF.mif
Command:  mrthreshold iter4_CF.mif -top 300 iter4_SF.mif
Command:  amp2response dwi.mif iter4_SF.mif iter4_first_dir.mif iter4_RF.txt
Command:  mrcalc iter4_SF.mif iter3_SF.mif -sub iter4_SF_diff.mif
Command:  mrthreshold iter4_CF.mif -top 3000 - | maskfilter - dilate - -npass 1 | mrcalc mask.mif - -mult iter4_SF_dilated.mif
Command:  dwi2fod csd dwi.mif iter4_RF.txt iter5_FOD.mif -mask iter4_SF_dilated.mif
Command:  fod2fixel iter5_FOD.mif iter5_fixel -peak peaks.mif -mask iter4_SF_dilated.mif -fmls_no_thresholds
Command:  fixel2voxel iter5_fixel/peaks.mif split_data iter5_amps.mif -number 2
Command:  mrconvert iter5_amps.mif iter5_first_peaks.mif -coord 3 0 -axes 0,1,2
Command:  mrconvert iter5_amps.mif iter5_second_peaks.mif -coord 3 1 -axes 0,1,2
Command:  fixel2voxel iter5_fixel/directions.mif split_dir iter5_all_dirs.mif -number 1
Command:  mrconvert iter5_all_dirs.mif iter5_first_dir.mif -coord 3 0:2
Command:  mrcalc iter5_first_peaks.mif -sqrt 1 iter5_second_peaks.mif iter5_first_peaks.mif -div -sub 2 -pow -mult iter5_CF.mif
Command:  mrthreshold iter5_CF.mif -top 300 iter5_SF.mif
Command:  amp2response dwi.mif iter5_SF.mif iter5_first_dir.mif iter5_RF.txt
Command:  mrcalc iter5_SF.mif iter4_SF.mif -sub iter5_SF_diff.mif
Command:  mrthreshold iter5_CF.mif -top 3000 - | maskfilter - dilate - -npass 1 | mrcalc mask.mif - -mult iter5_SF_dilated.mif
Command:  dwi2fod csd dwi.mif iter5_RF.txt iter6_FOD.mif -mask iter5_SF_dilated.mif
Command:  fod2fixel iter6_FOD.mif iter6_fixel -peak peaks.mif -mask iter5_SF_dilated.mif -fmls_no_thresholds
Command:  fixel2voxel iter6_fixel/peaks.mif split_data iter6_amps.mif -number 2
Command:  mrconvert iter6_amps.mif iter6_first_peaks.mif -coord 3 0 -axes 0,1,2
Command:  mrconvert iter6_amps.mif iter6_second_peaks.mif -coord 3 1 -axes 0,1,2
Command:  fixel2voxel iter6_fixel/directions.mif split_dir iter6_all_dirs.mif -number 1
Command:  mrconvert iter6_all_dirs.mif iter6_first_dir.mif -coord 3 0:2
Command:  mrcalc iter6_first_peaks.mif -sqrt 1 iter6_second_peaks.mif iter6_first_peaks.mif -div -sub 2 -pow -mult iter6_CF.mif
Command:  mrthreshold iter6_CF.mif -top 300 iter6_SF.mif
Command:  amp2response dwi.mif iter6_SF.mif iter6_first_dir.mif iter6_RF.txt
Command:  mrcalc iter6_SF.mif iter5_SF.mif -sub iter6_SF_diff.mif
Command:  mrthreshold iter6_CF.mif -top 3000 - | maskfilter - dilate - -npass 1 | mrcalc mask.mif - -mult iter6_SF_dilated.mif
Command:  dwi2fod csd dwi.mif iter6_RF.txt iter7_FOD.mif -mask iter6_SF_dilated.mif
Command:  fod2fixel iter7_FOD.mif iter7_fixel -peak peaks.mif -mask iter6_SF_dilated.mif -fmls_no_thresholds
Command:  fixel2voxel iter7_fixel/peaks.mif split_data iter7_amps.mif -number 2
Command:  mrconvert iter7_amps.mif iter7_first_peaks.mif -coord 3 0 -axes 0,1,2
Command:  mrconvert iter7_amps.mif iter7_second_peaks.mif -coord 3 1 -axes 0,1,2
Command:  fixel2voxel iter7_fixel/directions.mif split_dir iter7_all_dirs.mif -number 1
Command:  mrconvert iter7_all_dirs.mif iter7_first_dir.mif -coord 3 0:2
Command:  mrcalc iter7_first_peaks.mif -sqrt 1 iter7_second_peaks.mif iter7_first_peaks.mif -div -sub 2 -pow -mult iter7_CF.mif
Command:  mrthreshold iter7_CF.mif -top 300 iter7_SF.mif
Command:  amp2response dwi.mif iter7_SF.mif iter7_first_dir.mif iter7_RF.txt
Command:  mrcalc iter7_SF.mif iter6_SF.mif -sub iter7_SF_diff.mif
Command:  mrthreshold iter7_CF.mif -top 3000 - | maskfilter - dilate - -npass 1 | mrcalc mask.mif - -mult iter7_SF_dilated.mif
Command:  dwi2fod csd dwi.mif iter7_RF.txt iter8_FOD.mif -mask iter7_SF_dilated.mif
Command:  fod2fixel iter8_FOD.mif iter8_fixel -peak peaks.mif -mask iter7_SF_dilated.mif -fmls_no_thresholds
Command:  fixel2voxel iter8_fixel/peaks.mif split_data iter8_amps.mif -number 2
Command:  mrconvert iter8_amps.mif iter8_first_peaks.mif -coord 3 0 -axes 0,1,2
Command:  mrconvert iter8_amps.mif iter8_second_peaks.mif -coord 3 1 -axes 0,1,2
Command:  fixel2voxel iter8_fixel/directions.mif split_dir iter8_all_dirs.mif -number 1
Command:  mrconvert iter8_all_dirs.mif iter8_first_dir.mif -coord 3 0:2
Command:  mrcalc iter8_first_peaks.mif -sqrt 1 iter8_second_peaks.mif iter8_first_peaks.mif -div -sub 2 -pow -mult iter8_CF.mif
Command:  mrthreshold iter8_CF.mif -top 300 iter8_SF.mif
Command:  amp2response dwi.mif iter8_SF.mif iter8_first_dir.mif iter8_RF.txt
Command:  mrcalc iter8_SF.mif iter7_SF.mif -sub iter8_SF_diff.mif
dwi2response: Convergence of SF voxel selection detected at iteration 8
Function: shutil.copyfile(iter8_RF.txt, response.txt)
Function: shutil.move(iter8_SF.mif, voxels.mif)
Function: shutil.copyfile(response.txt, /scratch/ateillac/NHP/M157_5shells/Essai/02-TissueResponse/dwi2response-tmp-DRG9XQ/wm_ss_response.txt)
Command:  mrconvert voxels.mif /scratch/ateillac/NHP/M157_5shells/Essai/02-TissueResponse/dwi2response-tmp-DRG9XQ/wm_sf_mask.mif
dwi2response: Changing back to original directory (/scratch/ateillac/NHP/M157_5shells/Essai/02-TissueResponse/dwi2response-tmp-DRG9XQ)
dwi2response: Deleting temporary directory /scratch/ateillac/NHP/M157_5shells/Essai/02-TissueResponse/dwi2response-tmp-DRG9XQ/dwi2response-tmp-BWAOF1/
dwi2response: Command: '/home/ateillac/Softs/MRTRIX/mrtrix3/bin/mrstats wm_sf_mask.mif -output count -mask wm_sf_mask.mif' (piping data to local storage)
dwi2response: Result: 300
dwi2response: Command: '/home/ateillac/Softs/MRTRIX/mrtrix3/bin/mrstats gm_mask.mif -output count -mask gm_mask.mif' (piping data to local storage)
dwi2response: Result: 1564410
dwi2response: Command: '/home/ateillac/Softs/MRTRIX/mrtrix3/bin/mrstats csf_mask.mif -output count -mask csf_mask.mif' (piping data to local storage)
dwi2response: Result: 134401
Command:  amp2response dwi.mif wm_sf_mask.mif dirs.mif wm.txt -shell 0,244,800,1536,4075,7106
amp2response: Shell b=0: Response function [ 20490.2344 ] from average of 300 voxels
amp2response: Shell b=244: Response function [ 19293.3223 -1157.88391 -629.721313 -75.8945007  108.068123  48.0639648 ] solved after 9 constraint iterations from 300 voxels
amp2response: Shell b=800: Response function [ 14563.3799 -2763.39209  46.4419518 -19.7071095 -32.4169006  -50.580658 ] solved after 0 constraint iterations from 300 voxels
amp2response: Shell b=1536: Response function [ 13083.3184 -3678.49072  743.963013 -121.122055  35.3785896  1.31172085 ] solved after 0 constraint iterations from 300 voxels
amp2response: Shell b=4075: Response function [ 8740.21289 -3971.81079   1678.5708 -571.056824  161.318115 -26.5825653 ] solved after 0 constraint iterations from 300 voxels
amp2response: Shell b=7106: Response function [ 6709.48633 -3382.83374  1874.79749 -918.136597  379.606537 -120.274948 ] solved after 0 constraint iterations from 300 voxels
Command:  amp2response dwi.mif gm_mask.mif dirs.mif gm.txt -shell 0,244,800,1536,4075,7106 -isotropic

Thank you for your time.
Best regards,

Achille

4

OK, two questions arise immediately:

  • is your voxels.mif input in the expected 5TT format?

  • MRtrix3 immediately breaks your data into 6 shells:

    Is this what you expected? I suspect you may be feeding in DSI data, which does not have the expected multi-shell structure, and which MRtrix3 is not designed to handle…

5

Dear @jdtournier,

thank you for your quick answer and sorry for my delay. So, yeah I double-checked and it’s definitely the expected format for the voxels.mif. Secondly, I don’t exactly know if it was DSI data (since it’s a quite old dataset and I wasn’t there for the acquisition) but I only know they did a first scan for the b~250 with 3 directions, a second scan for the b~800 with 30 directions and finally, a third scan for the b~1500, b~4000 and b~7000 with a set of 60 directions (same for each shell). I did my best to answer and hope it will clarify the situation.

One point though, I find it really strange that I can reproduce the same strange behaviour (ie. really long computational time at the same step; even if it seems to be a bit faster) on one subject of the HCP dataset.

Thanks again for your help and any advices you could provide.

Regards,
Achille

6

The main issue as far as I can tell is that it’s using 1.5 million voxels for the GM response estimation - that should be in the thousands or less. But that’s most likely why it’s slow - and why the results will most likely be poor.

OK, good point. That and the previous point does make me think that maybe there is an issue with your 5TT image. One possibility is that your image might not be stored using a floating-point type, which may give odd results when regridded. What does mrinfo voxels.mif report?

It’s also worth interrupting and taking a look at the intermediate results in the temporary folder,
located here:

In particular, what do the fa.mif and 5tt_regrid.mif images look like…?

7

I thought it could be an infinite loop somewhere but it’s displaying some messages, so could it just be really long because of the averaging over more than a million voxels for each shell or could there be an error somewhere else?

For such a large number of voxels, it’ll likely be the fact that a transformation between amplitudes and spherical harmonics needs to be derived for each voxel. Another point that may slow this process down quite a lot is if your input DWI does not have data for the fourth axis (i.e. across volumes) stored contiguously in memory. But that’s only a symptom of the real problem, which is the number of voxels selected.

is your voxels.mif input in the expected 5TT format?

What does mrinfo voxels.mif report?

The -voxels option provides an output image, highlighting the voxels used for response function generation. The image that needs to conform to the 5TT format is: ../01-MifConversion/5TT/5TT.mif. However this should be guaranteed, given that the msmt_5tt algorithm calls:

Command:  5ttcheck 5tt.mif

Secondly, I don’t exactly know if it was DSI data (since it’s a quite old dataset and I wasn’t there for the acquisition) but I only know they did a first scan for the b~250 with 3 directions, a second scan for the b~800 with 30 directions and finally, a third scan for the b~1500, b~4000 and b~7000 with a set of 60 directions (same for each shell).

It’s not a typical multi-shell acquisition, but it is indeed multi-shell as opposed to DSI.

It’s also worth interrupting and taking a look at the intermediate results in the temporary folder

:arrow_up: :+1:

8

Oops, quite right… :blush:

9

This would be helpful to show us a screenshot of. If doing so, definitely add a colour scale bar to the FA image. Depending on the spatial resolution of your data, I’m not per se surprised that a large number of voxels get selected by that particular algorithm, in particular in an ex-vivo monkey brain (but still it would have to be supported by your spatial resolution).

Unless something that we hadn’t anticipated is wrong with your 5TT image, you may get luckier with the dhollander algorithm (if I say so myself :roll_eyes:), which doesn’t need a 5TT image. For ex-vivo data you may need to lower the -fa parameter, although given your highest b-value being 7000, it may not be necessary after all. See a few posts about this over here:

But definitely screen shot those FA and 5TT images for us… who knows what we might see.

10

Dear experts,

thank you for all your answers. As @rsmith pointed out, the 5ttcheck command is done within the dwi2response command but I did it again to verify the file format: the 5tt.mif image that I give is definitely the expected one with a 32 bit float encoding. I used the 5tt2vis command (fed with the 5tt_regrid.mif) to get this image which seems as expected:

all
all.png1020×841 86.4 KB

With the mrinfo command, I obtained for the 5tt_regrid.mif dimensions of 262 x 256 x 222 x 5 and a voxel size of 0.21 x 0.258 x 0.207 x ? (could this be linked to the storing memory you mentionned @rsmith?). However, the 5ttcheck command fails when I feed the 5tt_regrid.mif file as input, and I get:

5ttcheck: “[WARNING] Image 5tt_regrid.mif contains 141756 brain voxels with non-unity sum of partial volume fractions”
5ttcheck: “[ERROR] Input image does not conform to 5TT format”.

I suppose the 5tt_regrid.mif file should be conform with the 5TT format so, it seems that the first impression was the good one and that there is something wrong with my 5tt.mif file (which is strange as the 5ttcheck command was ok with it) or that the transformation from the 5tt.mif file to the 5tt_regrid.mif one is failing somehow…

On the other hand, for the fa.mif file, I got the same dimensions and resolutions (except in the 4th one where I didn’t have the question mark) as the 5tt_regrid.mif. There is a screenshot of it below (@ThijsDhollander what did you mean by the colour scale bar?):

fa
fa.jpg1134×1145 187 KB

I hope you’ll get all the information you needed. @ThijsDhollander, I’ll look into the -fa parameter and will definitely give the dhollander algorithm a try :wink:

Thank you again all for your help,

Achille

11

I should precise that the pathological tissue in my case is the cerebellum as you probably have guessed from the 5tt image.

Achille

12

Just the (gray scale in your case) bar that’s at the bottom right corner of your FA screenshot, so I could get an idea of the magnitude of the FA intensities I’m looking at. Thanks for that! From that screenshot, I can already assure you that you don’t need the -fa option to the dhollander algorithm; the default value (of 0.2) will work just fine out of the box!

The only thing I forgot to ask to show was an ADC map, and/or some information about the environment the brain sits in: you mentioned it was post-mortem, but does the brain still sit in the animal, e.g. suspended in CSF? Or does it sit in some other fluid? It seems to still take on a shape and even position that hints at it at least being suspended in a fluid; but it would be good to know nonetheless (and if not CSF: what fluid?). But again, and ADC map would give me the information I’m after just as well.

So far, all conditions are in place for the dhollander algorithm to work out of the box, without any special options. If you could run it, and specify the -voxels option to get the voxel output, and potentially screenshot that in the way described here: Multi-tissue CSD - #2 by ThijsDhollander , that would immediately let me confirm (or deny) that it worked well, and provide you with advice on how usable the output is, and in what way.

But just to be sure, and ADC map would definitely be useful as well to confirm a thing or two.

1 Like
13

I used the 5tt2vis command … to get this image which seems as expected:

:heart_eyes: Will be nice to get human in vivo 5TT segmentations looking like that…

could this be linked to the storing memory you mentioned @rsmith?

For the contiguity of memory storage, you want to be looking specifically at the strides of the image data as reported by mrinfo.

But having looked again at the amp2response command, there’s a couple of non-optimised aspects of the code that have not really caused issues in the past due to the relatively low overall computational expense of this command, but are being highly problematic in your specific use case (very high resolution & very clean delineation of the cortex). I’ll fix those issues at some point; but using instead the dwi2response dhollander algorithm may bypass this issue for you entirely.

5ttcheck: "[WARNING] Image 5tt_regrid.mif contains 141756 brain voxels with non-unity sum of partial volume fractions"
5ttcheck: “[ERROR] Input image does not conform to 5TT format”.

That’s not concerning. The issue here is that the re-gridding operation is a basic linear interpolation, which should preserve a unity sum within the mask but will introduce voxels with non-unity sum at the outer edges of the brain. If you re-run 5ttcheck on this image and provide the -masks option, I’m pretty sure the output mask image will highlight the outer edge of the brain mask. This image is never explicitly interpreted as a 5TT image - it’s just a shortcut within the script to re-sample all tissue partial volume images into diffusion space in a single step - so its non-conformance isn’t actually consequential.

15

Dear @ThijsDhollander and @rsmith,

thanks for your time.

That’s what I thought by looking at the fa map but thanks for the confirmation! You can find below a screenshot of the adc map (while I’ll ask how the brain was fixed and in which fluid):

adc
adc930×982 138 KB

I launched the dwi2response with the dhollander algorithm and the default -fa parameter and it seems to work just fine as you can see from the voxels.mif screenshot below:

voxels_RF
voxels_RF1033×1138 173 KB

I’m really impressed by how robust this algorithm is!

@rsmith thanks a lot for the clarification of what mrtransform did to get the 5tt_regrid.mif with the linear interpolation based on the fa map. I looked at what you mentionned with the use of the option -masks in the 5ttcheck command and it did exactly what you expected (mask of the outer edge of the brain mask)! I’ll bookmark the future fix conversation and look into the strides for the data image.

Thanks a lot to all you for your reactivity, advices and explanations!

Regards,
Achille

16

Nice, that looks all right. I was (before you showed the ADC map) potentially concerned there wouldn’t be any fluid with a decent SNR (or signal altogether) in there, as is sometimes the case when tissue is fixed like this. But it looks all good, and the high resolution allows for more than enough “pure” voxels in each tissue class. If you do happen to figure out what the fluid is, I’d certainly be interested to hear. :slightly_smiling_face:

Beautiful! This result is as perfect as we could’ve hoped for. :star_struck:
Robustness is what I always aim for indeed, together with an easy to use command and interface. Glad to see it pays off. :wink:

As a side note: also very interesting to see which single-fibre white matter voxels are able to be selected at this exceptional spatial resolution! A very interesting, and quite nicely symmetric (left-right) selection.

So all of this looking perfectly fine, you can safely proceed to use these response functions for multi-tissue CSD. If you’re keen to share, it’d also be interesting to see (the contents of the text files of) your WM, GM and CSF response functions.

1 Like
17

Dear @ThijsDhollander,

for the fluid, it seems it was fixed within agarose but they don’t recall exactly what treatment they used.

For the response functions, you can find below a screenshot of the WM response functions for the 6 shells done with shview:

responseFunctionsWM
responseFunctionsWM.png988×661 34.8 KB

The CSD and GM responses have both a sphere behaviour with different sizes. You can find all the responses here:

WM:
17244.31383624893 0 0 0 0 0

16083.85586735213 -707.1519892633875 -176.4925364634225 12.48616424129028 -46.17463567012031 -30.11985191385659

12029.73724642379 -2021.178439326104 245.4735857455719 -116.7632838165835 59.55977725413623 -94.43227920279337

10626.1636013819 -2779.075666971507 548.9684350598563 -78.34954648527987 4.119423785010004 4.834080059488596

7120.870916854469 -2952.120614158137 1225.406541217849 -415.5491774251993 113.9783382770016 -20.93654460734571

5519.608575399732 -2539.422338785555 1371.690297947774 -649.3539488662808 254.2827193554359 -80.07607048340495

GM:
24566.45298356919
22914.16788546077
16457.49962354114
13731.52010004403
6826.701371636201
4015.373006406009

CSF:
70694.78019725633
45944.70532795467
10561.5608648784
2694.780816132201
1715.681299083774
1607.920714526232

Thanks again all the team for your help!

Regards,

Achille

18

No problems there, shouldn’t cause any issues for response function estimation as well as multi-tissue CSD.

Again, looks great! It’s good that you went all the way up to that highest b-value of 7000 for an ex-vivo acquisition like this. That should really come in handy for a good quality CSD result.

Looks all good as well! I did a quick plot on my end to have a look at their relative qualities and contrast: again, looks really great! :+1:

1 Like