shimming-toolbox / shimming-toolbox-matlab

Code for performing real-time shimming using external MRI shim coils
GNU General Public License v3.0
16 stars 5 forks source link

Issue with MaRdI write function #56

Closed evaalonsoortiz closed 4 years ago

evaalonsoortiz commented 4 years ago

In realtime_zshim a MaRdI image object is first created for a field map time series. B0Fields = FieldEval( FM_mag_path, FM_phase_path, Params ); The field map time series is manipulated so that a Gz (dB/dz) time series is generated from it. GzFields = B0Fields; [~,GzFields.img(:,:,1,1,measNo)] = gradient(squeeze(g*B0Fields.img(:,:,1,1,measNo)),... If the Gz time series is saved as a nifti file using the MaRdI write() function: GzFields.write('nifti/Gz_tSeries/','nii',false) the saved image actually contains the field map time series instead. @rtopfer Do you have any suggestions/ideas for how I can fix this?

rtopfer commented 4 years ago

@evaalonsoortiz because doing the dicom <--> nifti transformations would really require understanding both header formats (which is beyond me), the MaRdI.write() function writes dicoms first and then defers to dcm2niix

glancing at the code, it looks like the potential for things to go wrong is that dicm2niix just takes a folder as input and looks for dicoms, and i'm not sure if either the dicom writing or the nifti conversion will overwrite existing files.

we can figure out a workaround but for now just make sure you're working in a clean directory and also nothing is in 'nifti/Gz_tSeries' before calling MaRdI.write

also: (this may have been the actual problem): the assigned filenames are based on the dicom headers. if you take the gradient of an image, i guess the dicom header isn't changed, so the image filenames will be misleading.

evaalonsoortiz commented 4 years ago

@rtopfer thanks for the feedback.

I am creating a new nifti/Gz_tSeries folder each time (the call GzFields.write('nifti/Gz_tSeries/','nii',false) is preceded by a rm -r nifti, so I don't think that is it.

I have been looking at the images themselves to verify that a Gz map actually contains voxel values that differ from a B0 map.

Looking at the data again, I see there is a mistake in my original comment, here is the revised comment:

If I create a new MaRdI structure GzFields = B0Fields; then compute the gradient [~,GzFields.img(:,:,1,1,measNo)] = gradient( ... and save that file GzFields.write('Gz_tSeries/','nii',false) the image actually contains the gradient.

If I then also save the B0 time series: B0Fields.write('B0_tSeries/','nii',false) it also contains the gradient and not a B0 time series.

For now, I've found a workaround (to having both the B0 and Gz time series saved as nifti files): I run the code without computing the gradient and save the B0 time series: the file actually contains the B0 time series. Then I run the code again, this time computing the Gz time series and saving that (the file contains the Gz time series). It seems that once I assign GzFields = B0Fields; and compute the gradient, I loose the B0 times series.

I can not save the B0 time series before computing Gz because MaRdI.write() converts data from double to uint16, which is not compatible subsequent with operations in the code.

rtopfer commented 4 years ago

ah, you can't initialize another image that way: MaRdI is a "handle class", so when you do GzFields = B0Fields you actually end up with 2 variables referring to the same data -- changes to one will be reflected in the other. Instead, useGzFields = BoFields.copy()

(the difference is akin to passing a variable by reference vs. passing it by value in C++) I'll see if there's a way to issue a warning when someone tries to do make the direct assignment.

jcohenadad commented 4 years ago

I'll see if there's a way to issue a warning when someone tries to do make the direct assignment.

yes 👍

evaalonsoortiz commented 4 years ago

Excellent, issue solved!