Hans:

We are currently working on a more fancy reconstruction than in the
MICCAI paper, using joint iterative deblurring, or deconvolution if you
will. It's already implemented in CMTK actually, but it has a
substantially larger computational demand than the inverse
interpolation. The benefit in terms of reconstruction accuracy is also a
bit questionable, especially given that I am having a hard time running
it on a machine with 8GB of memory, even using only single-precision
floats. Anyway, bottom line, we may have something more interesting yet
in a few months.

InvInt reconstruction is pretty darn good though if you ask me.

Let me know how things go! Ideally we should really move the further
discussion to the forums or the Wiki though.

Best,

TR

On 07/22/2009 08:21 PM, Hans Johnson wrote:
> Torsten,
>
> I when searching pub-med after I sent you the previous e-mail and found your
> recent publication " Volume Reconstruction by Inverse Interpolation:
> Application to Interleaved MR Motion Correction"
> :"http://www.pubmedcentral.nih.gov/picrender.fcgi?artid=2646840&blobtype=pdf
> "
>
> The results look very promising. We will try this out and let you know what
> we come up with.
>
> Thank you for the very quick response.
>
> Hans
>
>
> On 7/22/09 10:18 PM, "Torsten Rohlfing" wrote:
>
>
>> Hi Hans --
>>
>> On 07/22/2009 08:01 PM, Hans Johnson wrote:
>>
>>> Torsten,
>>>
>>> I am interested in testing the ability of cmtk to do interleaved image
>>> artifact correction to fix the "sawtooth" effect in brain images where motion
>>> occurred between collecting the even numbered slices and the odd number
>>> slices of an image.
>>>
>>> 1) Can CMTK address this issue?
>>>
>>>
>> Yes. Absolutely.
>>
>>> 2) Can you please give me some pointers about which tools to use, and the
>>> proper command line options? I would prefer if the wiki page were updated
>>> with this documentation, but an email response would also be appreciated.
>>>
>>>
>> The most simple tool is "film" (as in: Fix InterLeaved Motion) :) It
>> assumes that you have a 3D input stack with sawtooth effect, which it
>> will split into the passes, co-register them, and reconstruct a
>> corrected image using volume injection, followed by inverse interpolation.
>>
>> I am hoping that the command line self-description ("film --help") is
>> sufficient to get you off the ground fast, but I agree that this should
>> (and will) go on the Wiki (actually, our conversation thus far should
>> have probably gone into the CMTK forums on NITRC, but so be it...).
>>
>> I suggest you take a quick look at the film command line options, and
>> then we can clarify things that are unclear tomorrow.
>>
>> Best,
>> Torsten
>>

--
Torsten Rohlfing, PhD SRI International, Neuroscience Program
Senior Research Scientist 333 Ravenswood Ave, Menlo Park, CA 94025
Phone: ++1 (650) 859-3379 Fax: ++1 (650) 859-2743
torsten@synapse.sri.com http://www.stanford.edu/~rohlfing/

"Though this be madness, yet there is a method in't"

Torsten,

I have successfully run the program, but I don't know what file formats are supported. Apparently film will read compressed nifti files, but will not write them. How do I convert from RAW3D to nifti?


time ~/src/cmtk/releases/r101-build/bin/film -p 2 -c test.nii.gz defaults.nii.gz
f 53412.2 MSD 53412.2 MAX 1648.63 KLD 0.286495 LNORM 0
f 53412.1 MSD 53412.1 MAX 1648.63 KLD 0.286495 LNORM 0
f 120.65 MSD 120.65 MAX 263.543 KLD 0.0262974 LNORM 0
f 49.6766 MSD 49.6766 MAX 173.545 KLD 0.0263429 LNORM 0
f 19.1478 MSD 19.1478 MAX 166.245 KLD 0.0262041 LNORM 0
f 13.7889 MSD 13.7889 MAX 167.278 KLD 0.0260802 LNORM 0
f 9.99232 MSD 9.99232 MAX 169.084 KLD 0.0260529 LNORM 0
f 8.71491 MSD 8.71491 MAX 168.322 KLD 0.0259233 LNORM 0
f 7.95553 MSD 7.95553 MAX 168.392 KLD 0.025981 LNORM 0
f 7.17472 MSD 7.17472 MAX 167.882 KLD 0.0258797 LNORM 0
f 6.78608 MSD 6.78608 MAX 167.951 KLD 0.0259406 LNORM 0
f 6.40343 MSD 6.40343 MAX 167.71 KLD 0.0258565 LNORM 0
f 6.16931 MSD 6.16931 MAX 167.763 KLD 0.0259103 LNORM 0
f 5.90719 MSD 5.90719 MAX 167.642 KLD 0.0258443 LNORM 0
f 5.75607 MSD 5.75607 MAX 167.677 KLD 0.0258925 LNORM 0
f 5.55834 MSD 5.55834 MAX 167.615 KLD 0.0258336 LNORM 0
f 5.45409 MSD 5.45409 MAX 167.638 KLD 0.0258742 LNORM 0
f 5.2806 MSD 5.2806 MAX 167.609 KLD 0.0258258 LNORM 0
f 5.20751 MSD 5.20751 MAX 167.622 KLD 0.0258602 LNORM 0
f 5.02978 MSD 5.02978 MAX 167.612 KLD 0.0258282 LNORM 0
f 4.98556 MSD 4.98556 MAX 167.618 KLD 0.0258569 LNORM 0
f 4.87793 MSD 4.87793 MAX 167.615 KLD 0.0258344 LNORM 0
Fileformat not recognized; writing RAW3D instead.

real 69m27.034s
user 70m10.843s
sys 0m4.653s

Torsten,

I've now written out a file with ".nii" extension, but it appears that film did nothing to correct the interleaving.

I am sure that I've got something wrong on the command line:

svn co https://www.nitrc.org/svn/brainstestdata... InterleaveTestData

1) Convert from DICOM to NIFTI preserving physical image space:

InsightApplications/ConvertBetweenFileFormats InterleaveTestData input.nii.gz


2) Run a quick test to see if anything happens:
time ~/src/cmtk/releases/r101-build/bin/film --msd -L -n 10 -p 2 input.nii.gz fast30mintest.nii

The resulting output image was nearly the same as the input image.

Any help would be appreciated.

Hans:

CMTK is a bit stupid about target file extensions. If you want to write a NIFTI file, use "something.nii" as the file name. Compression will be applied automatically, unless you set the CMTK_WRITE_UNCOMPRESSED environment variable.

I admit that it would be nicer to get uncompressed when using .nii and compressed when using .nii.gz. I'll file a feature request tracker item, but I won't guarantee I'll actually have the time to implement this.

TR

Hans:

I looked at your data, and your problem is that you have a 3-pass interleave, not 2-pass. So the command line should have "-p 3", NOT "-p 2".

Hope this solves your problem.

TR

Hans:

Further to the issue of determining the number of passes in an interleaved image. I have not yet found a DICOM field that contains this information. If you find one, I'd be grateful if you'd let me know.

In the meantime, a way to confirm that you got the number right is to use CMTK's "split" tool -- simply split the input image into as many passes as you think there are. Then open the resulting pass images in a viewer (CMTK has a simple "triplanar" viewer for that purpose if you build with Qt and GUI support).

If you got the right number of passes, then all output images of "split" should be free of the sawtooth artifacts. If you still see the artifacts, you got the wrong number of passes.

TR

Torsten,

Thanks. I'll look into the DICOM header issue (we have a few DICOM header field experts around here), but I have a vague recollection that it is not handled consistently across all scanner manufacturers.

I am now doing a parameter space exploration to find the right command line options to use:
#time ~/src/cmtk/releases/r101-build/bin/film -p 3 -c test.nii.gz cdefaults.nii &
#time ~/src/cmtk/releases/r101-build/bin/film -p 3 -a test.nii.gz adefaults.nii &
#time ~/src/cmtk/releases/r101-build/bin/film -p 3 -s test.nii.gz sdefaults.nii &

time ~/src/cmtk/releases/r101-build/bin/film --nmi -S 2 -p 3 -c test.nii.gz cS2p3defaults.nii &
time ~/src/cmtk/releases/r101-build/bin/film --nmi -S 3 -p 3 -c test.nii.gz cS3p3defaults.nii &
time ~/src/cmtk/releases/r101-build/bin/film --nmi -S 4 -p 3 -c test.nii.gz cS4p3defaults.nii &

time ~/src/cmtk/releases/r101-build/bin/film --nmi -S 4 -r 7 -p 3 -c test.nii.gz cS3r7p3defaults.nii &


My initial results were not very satisfying. I'll send you pictures in a direct e-mail that show the troubles I'm having. Every 2nd and 3rd slice has a lot of black regions and complete signal dropout. It is as if there is no information at all used to estimate those voxels.

====================
Here is my interpretation of the command line arguments:

Interleave options:

At most one of [a|s|c|x|y|z] should be used. Where [a|s|c] are with respect to the image direction coded in the file. The [x|y|z] are more appropriate for directionless file formats and correspond to "voxel adjacent direction", "row adjacent direction", and "slice adjacent direction" respectively.

-p The number of interleaved passes.

-W ***no idea what this is for, or what effect it would have on the processing.


Motion Correction/Registration:

-R This could be a T1 weighted high resolution image, and it might be a good idea to use this as a reference image to minimize the number of times that resampling the image needs to be done. It is very likely that after the interleaving step a rigid registration to the T1 image would be done, and this allows you to minimize the number of resamplings that are occuring.

[nmi|mi|cr|msd|cc] Which of these is the default? I guess the recommendation would be to use msd if the reference image is the same as the interleave image (where the first pass is the default reference), or nmi if the reference is a different modality.

[--import-xforms-path|--export-xforms-path] ?? What is this useful for?

-S I've not figured out a good value for this yet. I am assuming that a good value for this is related to how much motion occured and the between pass slice thickness. i.e. The more out of plane motion the larger the smothing is needed to force the regions to overlap. This is at the expense of blurring the entire image.

-r ????? Defines a region about the current output voxel that should be used to estimate it's value. The default value of "0" means only voxels from the transformed passes that intersect the region of the current output voxel are used.


Inverse Interpolation Options

As shown in your paper, -L is fast but not very accurate, -O is SLOW, but very accurate, -C is close to -O but much faster (I assum -H is between -C and -O).

-f ??? Just a tradeoff between speed and errors, if you are happy with second-order, then it is unlikely I'd use forth order for day to day processing. Perhaps if I had a cluster of computers and this worked great it would be worth the extra time.

-T Only used to show that you should not use it in you paper???

Hi Hans --

Thanks for writing the film manual ;)

> My initial results were not very satisfying. I'll send you pictures in a direct e-mail that show the troubles I'm having. Every 2nd and 3rd slice has a lot of black regions and complete signal dropout. It is as if there is no information at all used to estimate those voxels.

This could also be a result of too small a volume injection kernel in the initial reconstruction.

In general, my advice is to

a) use a rather larger injection kernel, with --injection-kernel-sigma about the size of the slice thickness and --injection-kernel-radius about twice that.

b) use at least cubic kernel for inverse interpolation.

Don't worry about the blurring of the image in a) for large sigma; the whole point of the inverse interpolation stage is to un-blur the initial estimate, and our experience is that inverse interpolation is better at unblurring a smooth image where needed than it is at smoothing an initially noise image.

As for your interpretation of the command line parameters, I will paste that into the Wiki and add my corrections:

https://neuro.sri.com/cmtk/wiki/index.ph...

Thanks for posting this information with the corrections!

I've started another run with your suggestions... I'll let you know the results in a few hours when it finished.

Hans

Hans:

Just to let you know -- I just checked the improved "film" source code into svn/trunk. I will probably do a minor version release later today, assuming I get everything compiled on Windows too.

Most relevant for you: film now uses anisotropic volume injection, which is more appropriate than isotropic for anisotropic input data.

However, this changed the meaning of the kernel width and truncation parameters. Kernel width is now in multiples of pixel size PER DIMENSION, i.e., for sigma=1.0 the Gaussian kernel has standard deviation 1*delta_x in x, 1*delta_y in y, and 1*delta_z in z direction. The deltas refer to your reconstructed grid, by the way, not the pass image grids, because the point spread function does not depend on the number of passes in the interleaved imaging. Sigma defaults to 0.5 by the way.

Likewise, the truncation radius is now a relative factor that is multiplied with sigma*pixelsize in each dimension, effectively generating an ellipsoid truncation surface. Default here is 2.0, i.e., truncation happens after two times sigma in each direction.

As a side note, I also made some improvements to the command line class in CMTK. All options now automatically print their default values and states in the command line --help output. Some tools may have redundant information because I used to print this information manually, but these will be cleaned up over time.

So now for your testing: I would think that the default values for sigma and truncation should pretty much be okay for most cases. You may want to increased the truncation radius a little if you have trouble, but that makes things slower of course.

Best,
Torsten

Thanks! I'm running the tests now. I'll let you know how it goes.

Hans