Hi Andriy --

There are actually two different algorithms in that tool. One is based on inverse interpolation pretty much like we described it here:

T. Rohlfing, M. H. Rademacher, and A. Pfefferbaum, “Volume reconstruction using inverse interpolation: application to interleaved image motion correction,” in Medical Image Computing and Computer-Assisted Intervention - MICCAI 2008. 11th International Conference, New York, NY, USA, September 6-10, 2008, Proceedings, Part I, D. Metaxas, L. Axel, G. Fichtinger, and G. Székely, Eds., Berlin/Heidelberg, 2008, vol. 5241 of Lecture Notes in Computer Science, pp. 798-806, Springer-Verlag.

This algorithm is used when one of the following command line options is given: --linear, --cubic, --cosine-sinc, --hamming-sinc

The second algorithm is a more complex deblurring algorithm that also incorporates an estimate of the imaging point spread function. That's more in line with the usual superresolution literature, but it's less well tested in CMTK. I have actually been working for over a year on writing this algorithm up, but something always gets in the way ;)

Anyway, this algorithm is used when you give one of the following options: --deblurring-box (box-shaped PSF), --deblurring-gaussian (Gaussian-shaped PSF). You typically want to use these two options also to control the PSF shape and scale: --psf and --psf-scale . The former sets the "shape" of the PSF, i.e., the width in x, y, and z direction of your acquired data. The latter scales the PSF kernel with a global factor.

I am sure you'll have more questions, which I'll be happy to answer.

Best,
Torsten

> I recently came up with an idea to apply super-resolution resampling
> in one of my project.s. After some searching, I came across such tool
> in CMTK (disclaimer: have downloaded or tried to compile CMTK). If I
> am not mistaken, it's in core/apps/volume_reconstruction.cxx
>
> I am interested which algorithm is implemented in this tool? Is there
> a reference? From just a quick google search, I see there are several
> approaches proposed, so it would be great to know more information
> about yours. I don't seen any such information in the user guide.

Re 1): I think the algorithms are both fine for your application. We have used them here as well for reconstructing super-resolution data, say, from multiple repeated high-resolution acquisitions. The only difference is that you need to take care of co-registration "yourself", whereas for the interleaved imaging motion correction there is an integrated one-stop-shop tool ("film") that integrates the entire workflow from registration to reconstruction.

Re 2): when you use inverse interpolation, linear is fastest but worst, cosine-sinc is best but slowest, and cubic is in the middle. Forget Hamming sinc, that's probably broken. The jury is still out on the deblurring kernels, but I'd start with a Gaussian kernel, set "--psf" to your acquired voxel size, and "--psf-scale" to 1.0. Important detail: regularization. By default, the tool uses a nonlinear range truncation as described in the MICCAI paper. We have very good experience with this one, but you can turn it off with the "--no-truncation" switch. You can then optionally also use L2 norm Tikhonov-type regularization with the "--l-norm-weight " option. As you'd expect, finding the optimal constraint weight is non-trivial.

By the way -- you may also want to look at the initial stage of the reconstruction, which is a push-forward volume injection with a Gaussian kernel. The kernel parameters for that are specified using the "--gauss-sigma" (width) and "--radius" (truncation radius) options. Sorry for the slightly stupid option names -- this tool has evolved slowly and isn't quite mature yet.

Best,
TR

> Yes, I do have more questions!
>
> 1) I briefly looked at your MICCAI paper, and you discuss motion
> correction as the main application. My application is quite different
> -- I have three orthogonal projections of similar quality, with high
> in-slice resolution, and thick slices. I want to get one isotropic
> volume by combining these three volumes, which I will register prior
> to that. Do you think your algorithm is a suitable tool for this task?
> Which algorithm out of the two implemented you think is more suitable?
> (I have not yet tried simple averaging and such, which will be my
> first, brain-dead approach)
>
> 2) as always -- any guidelines on parameter selection strategy?

Hi Andriy -

Did you register your input volumes before running volume reconstruction?

Also, you can write the initial, volume-injected image to disk using the "--write-splatted-image" option. That may give some extra hints.

As for the anisotropic output -- the reconstruction grid is "guessed" from the input images. In your case, you probably want to either provide an existing image as the reconstruction grid (using "--recon-grid-path"), or simply generate one on the fly (using "--recon-grid").

TR

On 04/06/2010 10:25 AM, Andriy Fedorov wrote:
> Hi Torsten,
>
> I compiled your tool, and ran it on my data. I actually left it
> running overnight with --window-sinc interpolator. I have three
> volumes, each with ~.25mm resolution in-slice, and 3 mm slice
> thickness. I passed these three volumes to volume_reconstruction tool.
>
> I attach the screenshot of the result. The volume it produced is not
> isotropic, and it also has some grid lines artifact.
>
> Perhaps this was because I did not specify --psf flag. I am re-running
> it now, with the --psf 0.2734,0.2734,0.2734 option.
>
> Let me know if you have any hints as of what I might be doing wrong. Thanks
>
> AF
>

On 04/06/2010 12:09 PM, Andriy Fedorov wrote:
> On Tue, Apr 6, 2010 at 13:58, Torsten Rohlfing wrote:
>> > Did you register your input volumes before running volume reconstruction?
>> >
> Yes, I did. I pass the transforms for each image (except the first
> one) in the command line.
>
>> > Also, you can write the initial, volume injuected image to disk using the
>> > "--write-splatted-image" option. That may give some extra hints.
>> >
> Yes, I tried this, splatted image has the same problem with grid lines.

Okay, then that's the root of the problem -- if the volume injection fials miserably, then deblurring cannot recover in all cases. In general, it is better to have a slightly blurry initial reconstruction, so maybe you want to try increasing the injection kernel size, e.g., "--gauss-sigma 2.0 --radius 4"

>
>> > As for the anisotropic output -- the reconstruction grid is "guessed" from
>> > the input images. In your case, you probably want to either provide an
>> > existing image as the reconstruction grid (using "--recon-grid-path"), or
>> > simply generate one on the fly (using "--recon-grid").
>> >
> Will try this next, thanks!
>

I am not using any ITK functionality for accessing anything actually. CMTK is completely independent of ITK -- doesn't even support using ITK as a backend library.

That said, I previously confirmed that volumes read into Slicer are correctly aligned in the GUI when I register them using CMTK's tools as a plugin, and the nonrigid registration in CMTK appears to correctly read Slicer's tfm files for initialization. Of course none of this is a guarantee for complete correctness.

Anyway, to solve your problem I think the best way to proceed is to use CMTK's "registration" tool for registration.

The only other thing I can think of is to confirm that you are using the ITK-generated transformations in the correct directions, and not accidentally in the inverse direction.

TR

On 04/06/2010 03:09 PM, Andriy Fedorov wrote:
> Torsten,
>
> This indeed revealed a problem. The reformatted volumes do not overlap
> with image1.
>
> I obtained the transformations by using registration modules from
> Slicer, which are based on ITK functionality. I am very confident that
> the transforms are correct, because I can visualize the aligned
> volumes in Slicer, and my pure ITK code for constructing
> super-resolution volume by simple averaging produces aligned result.
>
> The volumes I am working with are not axis-aligned though. Is CMTK
> expected to handle such volumes correctly? Is this particular
> functionality tested? I know this has been a standing issue in ITK for
> quite a while, but it's been resolved, and I am not sure how much of
> ITK you are using to access images.
>
> AF

On 04/06/2010 03:26 PM, Andriy Fedorov wrote:
> On Tue, Apr 6, 2010 at 18:14, Torsten Rohlfing wrote:
>> I am not using any ITK functionality for accessing anything actually. CMTK
>> is completely independent of ITK -- doesn't even support using ITK as a
>> backend library
>
> Oh, ok, I didn't know this. Sorry, I didn't mean to offend you.
No offense taken -- this is just important to keep in mind. Basically, whatever behaviour you see in ITK cannot strictly predict any behaviour you see in CMTK and vice versa.
>> That said, I previously confirmed that volumes read into Slicer are
>> correctly aligned in the GUI when I register them using CMTK's tools as a
>> plugin, and the nonrigid registration in CMTK appears to correctly read
>> Slicer's tfm files for initialization. Of course none of this is a guarantee
>> for complete correctness.
>>
>
> Have you done this test for any volumes that are not axis-aligned?
>
> Here's an example header:
>
> NRRD0004
> # Complete NRRD file format specification at:
> # http://teem.sourceforge.net/nrrd/format....
> type: float
> dimension: 3
> space: left-posterior-superior
> sizes: 512 512 26
> space directions: (0.273395,0.00130997,0.00105967)
> (-0.00147682,0.269232,0.0475336) (-0.00895129,-0.521637,2.95429)
> kinds: domain domain domain
> endian: little
> encoding: gzip
> space origin: (-70.309,-86.797,-52.275)
>
> Is there a better way to call such volume, other than "not axis-aligned"?
I would call it "oblique", but I see what you're saying. I just wanted to make sure you're not implying "non-orthogonal". But it seems to me from the direction vectors herein that the volume is indeed orthogonal. With non-orthogonal volumes you're have trouble in CMTK (I guess the same goes for ITK though)

Now for the question of whether CMTK handles these ... indeed, it should. Here's the problem though: internally, for historic reasons, CMTK reorients every image into RAS space, and all transformations are computed and interpreted in this space.

That said, when a transformation is read or written in ITK format, then CMTK automatically converts between the native image space (as defined in the NRRD header for example) and the internal space. I did indeed test that, because Slicer will use non-trivial coordinate systems for images it reads, and all my tests with such images worked.

Nonetheless, I think your best bet remains using CMTK itself for registration. That's as simple as running

registration --auto-multi-levels 4 -o image12.xform image1.nrrd image2.nrrd

In your case it may make sense to first initialize the transformations using the "make_initial_affine" tool based on the NRRD coordinate system info:

make_initial_affine --direction-vectors image1.nrrd image2.nrrd image12_init.xform
registration --auto-multi-levels 4 -o image12.xform --initial image12_init.xform image1.nrrd image2.nrrd

By the way -- you can already use the aforementioned reformatx command line to see whether the result of the make_initial_affine command is reasonable (i.e., whether the coordinate information in the NRRDs actually gives valid alignments).

TR

Andriy:

So after some thinking, I believe the answer to your problem is this: ITK's tfm transformations are not converted into CMTK's coordinate space in the volume_reconstruction (or, for that matter, reformatx) tool. Should be relatively straight forward for the v_r tool, though, so I'll give it a shot tomorrow.

TR

Hi Andriy:

I just implemented what could be a fix for the ITK tfm import problem in volume_reconstruction. It's in the latest SVN/trunk. Unfortunately I don't have time to test it right now, and more unfortunately the reformatx tool will take a little longer to fix, but if you want to give this a spin...

Another question: it seems almost like you have a tagging grid on your image data. If that's the case, the grid will be reconstructed also by the deblurring algorithm. If you wanted to reconstruct non-tagged data, we'd have to do some more work.

Best,
TR

Okay, I think I am getting a better idea of what data you're working with.

My first observation is that you seem to be reconstructing at a substantially higher resolution than what you acquired. That's fine per se, but since you also seem to have only three acquisitions to work with, the whole thing breaks down when you reconstruct away from the acquired data locations. That would be my guess for why there's the observed 3D "grid."

Also, it seems to me your data were acquired in the presence of nonrigid motion, which leads to misalignment and thus inconsistencies between the acquisitions, which then affects the reconstruction.

Finally, is there any reason why there would be intensity changes during the acquisition, i.e., is there contrast enhancement or diffusion weighting going on?

TR

On 04/08/2010 01:32 PM, Andriy Fedorov wrote:
> Yes, I did get rid of the holes!
>
> I put the result in my dropbox, maybe if you take a look, you can get
> a better idea of what's going on:
>

>
> On Thu, Apr 8, 2010 at 16:22, Torsten Rohlfing wrote:
>
>> Hi Andriy.
>>
>> Did you manage to get rid of the "holes" in the volume injection stage?
>>
>> TR

On 04/23/2010 08:22 AM, Andriy Fedorov wrote:
> On Thu, Apr 8, 2010 at 16:40, Torsten Rohlfing wrote:
>> My first observation is that you seem to be reconstructing at a
>> substantially higher resolution than what you acquired. That's fine per se,
>> but since you also seem to have only three acquisitions to work with, the
>> whole thing breaks down when you reconstruct away from the acquired data
>> locations. That would be my guess for why there's the observed 3D "grid."
>>
>
> Torsten, I agree. So in my case, each image is 0.2 mm in-plane, and 2
> mm out of plane.
>
> What resolution in your opinion I can safely use for the output image
> to make sure there is enough data to reconstruct isotropically?
> Is it 2 mm?
>
Hi Andriy -

Indeed, I think 2mm is probably the limit in your case. In general, you can't expect the super-resolution recon to make good data out of nothing. So in order to get true super-resolution, you have to increase the amount of available actual data by repeated acquisition. Basically, super-resolution helps you deal with irregular sampling locations due to, say, motion of the subject etc.

In your case, if you wanted to go better than 2mm, I would suggest dropping the orthogonal acquisitions, which really don't help much, and instead acquire more planes between the 2mm planes of the first volume. That way, you still keep the .2mm in-plane and increase the effective through-plane resolution after SR recon to a degree that depends on how many additional slices you acquire in between.

TR

On 04/23/2010 01:51 PM, Andriy Fedorov wrote:
> On Fri, Apr 23, 2010 at 16:18, Torsten Rohlfing wrote:
>> Indeed, I think 2mm is probably the limit in your case. In general, you
>> can't expect the super-resolution recon to make good data out of nothing. So
>> in order to get true super-resolution, you have to increase the amount of
>> available actual data by repeated acquisition. Basically, super-resolution
>> helps you deal with irregular sampling locations due to, say, motion of the
>> subject etc.
>>
>
> I see. Actually, I am getting more or less reasonable result by doing
> averaging over the values interpolated from the input images. I am
> also trying to use bspline interpolation of the input data without
> averaging, but this is not quite working yet.
Yes, that is an option of course. The reason why we went SR for our original application (interleaved image motion correction) is that by averaging interpolated images your resolution is limited by the acquired resolution (here: mostly slice spacing) even if the actual sampling after multiple passes is actually much more dense. Interpolation and averaging also does not allow you to take advantage of the proximity between acquired and reconstructed samples, i.e., if you acquired a pixel at exactly a location of a particular recon pixel, you still get a blurred average there, too.

As a substitute for interpolation and averaging, you may want to liik into using the volume injection, which is the initial step of our recon procedure (also, there's a separate tool "volume_injection" for that in CMTK). Basically, you get the locality and dense acquisition benefits where they apply without the numerical issues of iterative recon when data isn't dense enough for the deblurring. The only downside then is that you get locally different resolutions in your reconstructed images, but the alternative to that would be to have them blurry everywhere, even in places where it's not strictly required.

TR

On 04/23/2010 03:50 PM, Andriy Fedorov wrote:
> On Fri, Apr 23, 2010 at 16:58, Torsten Rohlfing wrote:
>> Yes, that is an option of course. The reason why we went SR for our original
>> application (interleaved image motion correction) is that by averaging
>> interpolated images your resolution is limited by the acquired resolution
>> (here: mostly slice spacing) even if the actual sampling after multiple
>> passes is actually much more dense. Interpolation and averaging also does
>> not allow you to take advantage of the proximity between acquired and
>> reconstructed samples, i.e., if you acquired a pixel at exactly a location
>> of a particular recon pixel, you still get a blurred average there, too.
>>
>
> Right, but if you don't average, but interpolate directly from the
> input points, this will not be a problem, right?

True, but then you're stuck using only a single image stack, no?

>
>> As a substitute for interpolation and averaging, you may want to liik into
>> using the volume injection, which is the initial step of our recon procedure
>> (also, there's a separate tool "volume_injection" for that in CMTK).
>> Basically, you get the locality and dense acquisition benefits where they
>> apply without the numerical issues of iterative recon when data isn't dense
>> enough for the deblurring. The only downside then is that you get locally
>> different resolutions in your reconstructed images, but the alternative to
>> that would be to have them blurry everywhere, even in places where it's not
>> strictly required.
>>
>
> I am not sure this would benefit my application. I am interested in
> using this volume for visualization (volume rendering) and
> registration. I am not quite sure I would gain much from having
> non-uniform resolution. I need to try, I guess. Thanks for the
> suggestion!

I think my suggestion is probably more appropriate for registration -- when you render, the effect of "ugliness" is important. The only problem I can see with registration is that the grid-like structure of high- vs. low-res regions might get the registration caught. I am not necessarily convinced that's going to happen though, unless you use image gradient differences (not image gradients) to drive the registration.