[MINC-users] mincresample and linear 'deformations'

Thu Aug 18 11:36:40 EDT 2016

Hi Alex,

Sure, it's possible, but I would have to be careful about the unintended
consequences. My biggest issue is how the change might affect some of the
clients of the transform functions. At the very least the change would be
the effect on some of the libminc tests - several would have to be modified
to account for the change in behavior of general_transform_point() and its
ilk.

    -bert

On Thu, Aug 18, 2016 at 11:00 AM, Alex Zijdenbos <zijdenbos at gmail.com>
wrote:

> Hi all,
>
> Thanks for all your input. I'm still trying to figure out how to do what I
> want, which is really to convert an arbitrary xfm into a deformation field
> and have everything "just work", but I have some ideas that I'll pursue.
>
> Bert, whatever the "right thing" is for mincresample, I'm pretty sure we
> don't want to change its current behaviour; however, do you think it would
> be feasible to add an option to it to control this? Would require
> implementing what you suggest, and then adding a mincresample option to
> insert the fill value as opposed to assuming '0 displacement'.
>
> -- A
>
>
>
> On Thu, Aug 18, 2016 at 10:37 AM, Robert D. Vincent <
> robert.d.vincent at mcgill.ca> wrote:
>
> > Hi all,
> >
> > The core of the issue from an implementation perspective is the function
> > evaluate_grid_volume(), which is responsible for evaluating the grid
> > transform for a point in world space. If the point lies outside the grid,
> > this function silently returns zero for all of the displacements, so the
> > point is not transformed and the original data is copied.
> >
> > It is trivial to modify this function to return an error indication if
> the
> > point is outside the grid, so that the resampling process could "do the
> > right thing", whatever that is...
> >
> >     -bert
> >
> > On Wed, Aug 17, 2016 at 8:38 PM, Andrew Janke <a.janke at gmail.com> wrote:
> >
> > > Note that xfm2def does not have to calculate the sampling of the grid
> > > from the input volume. You can define this output grid transform
> > > sampling however you choose. This will solve your current problem, but
> > > not all problems like this as the input xfm can be any xfm, which
> > > itself may have a defined sampling (eg: def grid)
> > >
> > > I'd argue that what you are seeing is correct behaviour, I use this
> > > feature of mincresample regularly when doing resampling of large (TB
> > > scale) microscopy data. In this case you perform the registration on
> > > chunks and then perform some magic to smoothly interpolate the
> > > resulting deformation grids at the edges. Move edges, repeat, iterate.
> > >
> > >
> > > a
> > >
> > > On 18 August 2016 at 08:05, Alex Zijdenbos <zijdenbos at gmail.com>
> wrote:
> > > > Thanks, Louis, Vlad!
> > > >
> > > > I was getting to the same conclusion as Vlad's explanation; but that
> > > still
> > > > leaves me with the problem that two conceptually the same operations
> > > > generate very different output:
> > > >
> > > > 1) voxel-wise transformation inversion fails (due to lack of support)
> > =>
> > > 0
> > > > deformation => insert source image value
> > > > 2) voxel-wise transformation ok/no inversion needed, points outside
> > input
> > > > image => insert 0 (or fill value).
> > > >
> > > > Not sure how to deal with this, unless by adding an option to
> > > mincresample
> > > > that would insert a 0 (fill value) when the transformation inversion
> > > fails,
> > > > as opposed to treating that as '0 transformation'. E.g.,
> > > -fill_def_missing
> > > > or some such.
> > > >
> > > > Alternatively, pad or transform the deformation grid lattice to allow
> > for
> > > > the inversion to work?
> > > >
> > > > -- A
> > > >
> > > >
> > > >
> > > > On Wed, Aug 17, 2016 at 5:24 PM, Vladimir S. FONOV <
> > > vladimir.fonov at gmail.com
> > > >> wrote:
> > > >
> > > >> xfm2def samples forward deformation in the space of input image.
> > > >> Then when you apply the transformation, mincresample will try to
> > > calculate
> > > >> the inverse (numerically) in the space of output image.
> > > >> So, if your transformation is a shift of 30mm the forward
> > transformation
> > > >> will have the same vector everywhere. However when inverse is
> > > calculated,
> > > >> mincresmple will reach the edge of the domain where vectors are
> > defined
> > > and
> > > >> will take value of 0 for the translation outside of this domain -
> > that's
> > > >> why you have a sudden jump in the result.
> > > >>
> > > >> However, if you invert the original affine transformation,  and
> sample
> > > it
> > > >> for every voxel and then invert again (applying xfminvert to a
> > > non-linear
> > > >> transformation doesn't actually recompute the transform, only adds a
> > > flag
> > > >> that it have to be inverted), then when mincresample will run, it
> will
> > > >> again apply the inverse of the transform, on top of inverted
> > transform ,
> > > >> negating the need to actually numerically recompute the
> > transformation.
> > > >>
> > > >>
> > > >>
> > > >> On 2016-08-17 04:30 PM, Alex Zijdenbos wrote:
> > > >>
> > > >>> Right - I am aware of that, but I don't think that explains any of
> > this
> > > >>> though?
> > > >>>
> > > >>> On Wed, Aug 17, 2016 at 3:21 PM, Vladimir S. FONOV <
> > > >>> vladimir.fonov at gmail.com
> > > >>>
> > > >>>> wrote:
> > > >>>>
> > > >>>
> > > >>> When you apply "simple" forward transform, mincresample have to
> > invert
> > > it
> > > >>>> internally.
> > > >>>>
> > > >>>> On Wed, Aug 17, 2016 at 1:50 PM, Alex Zijdenbos <
> > zijdenbos at gmail.com>
> > > >>>> wrote:
> > > >>>>
> > > >>>> These were simple forward resamplings. I experimented a bit with
> > > various
> > > >>>>> sequences of inversions, and found that this process:
> > > >>>>>
> > > >>>>> xfminvert( xfm2def( xfminvert( lin.xfm )))
> > > >>>>>
> > > >>>>> Yields a deformation grid that I can use in a forward resampling
> > call
> > > >>>>> and
> > > >>>>> have it produce the same result as resampling using lin.xfm
> > directly.
> > > >>>>>
> > > >>>>> Well I am puzzled.
> > > >>>>>
> > > >>>>> -- A
> > > >>>>>
> > > >>>>> For some more details, I tried these steps. Starting with
> > > >>>>>
> > > >>>>> param2xfm -translation 10 20 30 lsq3.xfm
> > > >>>>> mincresample -use_input_sampling -transform lsq3_D.xfm in.mnc
> > res.mnc
> > > >>>>>
> > > >>>>> I had done this:
> > > >>>>>
> > > >>>>> xfm2def -float [lattice from input volume] lsq3.xfm lsq3_D.xfm
> > > >>>>> mincresample -use_input_sampling -transform lsq3_D.xfm in.mnc
> > > res_D.mnc
> > > >>>>>
> > > >>>>> which "fails". Similarly, this:
> > > >>>>>
> > > >>>>> xfminvert lsq3_D.xfm lsq3_D_inv.xfm
> > > >>>>> mincresample -use_input_sampling -transform lsq3_D_inv.xfm
> -invert
> > > >>>>> in.mnc
> > > >>>>> res_D_inv.mnc
> > > >>>>>
> > > >>>>> reproduces res_D.mnc (the "bad" result). Note that the xfminvert
> > > call in
> > > >>>>> this case only sets the "Invert" flag in the xfm, it doesn't
> > actually
> > > >>>>> invert the deformation field. But now, this:
> > > >>>>>
> > > >>>>> xfminvert lsq3.xfm lsq3_inv.xfm
> > > >>>>> xfm2def -float [lattice from input volume] lsq3_inv.xfm
> > > lsq3_inv_D.xfm
> > > >>>>> mincresample -use_input_sampling -transform lsq3_inv_D.xfm
> -invert
> > > >>>>> in.mnc
> > > >>>>> res_inv_D.mnc
> > > >>>>>
> > > >>>>> Exactly reproduces res.mnc (the "good" result). So inverting the
> > > linear
> > > >>>>>
> > > >>>> xfm
> > > >>>>
> > > >>>>> first, then converting it to a deformation and resampling with
> > > >>>>> inversion,
> > > >>>>> somehow does the right thing. Then this:
> > > >>>>>
> > > >>>>> xfminvert lsq3.xfm lsq3_inv.xfm
> > > >>>>> xfm2def -float [lattice from input volume] lsq3_inv.xfm
> > > lsq3_inv_D.xfm
> > > >>>>> xfminvert lsq3_inv_D.xfm lsq3_inv_D_inv.xfm
> > > >>>>> mincresample -use_input_sampling -transform lsq3_inv_D.xfm in.mnc
> > > >>>>> res_inv_D_inv.mnc
> > > >>>>>
> > > >>>>> also does the right thing.
> > > >>>>>
> > > >>>>>
> > > >>>>> On Wed, Aug 17, 2016 at 12:28 PM, Vladimir S. FONOV <
> > > >>>>> vladimir.fonov at gmail.com> wrote:
> > > >>>>>
> > > >>>>> It looks like it's a problem of finding an inverse of the
> > > >>>>>>
> > > >>>>> transformation
> > > >>>>
> > > >>>>> outside of the domain where it's defined.
> > > >>>>>>
> > > >>>>>> When you apply the transform do you apply it directly or with
> > > >>>>>> -invert_transformation ?
> > > >>>>>>
> > > >>>>>>
> > > >>>>>> On 2016-08-17 11:54 AM, Alex Zijdenbos wrote:
> > > >>>>>>
> > > >>>>>> Hi all,
> > > >>>>>>>
> > > >>>>>>> I was always under the impression that one could turn a linear
> > xfm
> > > >>>>>>>
> > > >>>>>> into
> > > >>>>
> > > >>>>> a
> > > >>>>>
> > > >>>>>> deformation field (using xfm2def), and that the application of
> the
> > > >>>>>>> resulting deformation would be equivalent to using the linear
> > xfm.
> > > >>>>>>> However,
> > > >>>>>>> mincresample disagrees.
> > > >>>>>>>
> > > >>>>>>> In the attached image, I've created a simpe lsq3 (param2xfm
> > > >>>>>>>
> > > >>>>>> -translation
> > > >>>>
> > > >>>>> 10
> > > >>>>>>> 20 30) linear xfm, generated a 'deformation' field out of that
> > > using
> > > >>>>>>> xfm2def (same lattice as the input volume), and resampled the
> > input
> > > >>>>>>>
> > > >>>>>> volume
> > > >>>>>
> > > >>>>>> with both xfms. Clearly, in the case of a 'deformation' xfm, the
> > > >>>>>>>
> > > >>>>>> volume
> > > >>>>
> > > >>>>> is
> > > >>>>>
> > > >>>>>> filled with unmodified voxels where I didn't expect them.
> > > >>>>>>>
> > > >>>>>>> This is tied to the sampling lattice of the deformation field,
> > > which
> > > >>>>>>>
> > > >>>>>> in
> > > >>>>
> > > >>>>> this example is the same as the input volume. The resamplings are
> > > done
> > > >>>>>>> using -use_input_sampling btw. The two bottom rows show the
> > > resampling
> > > >>>>>>> using a padded input volume (no change), and using a
> > 'deformation'
> > > >>>>>>>
> > > >>>>>> grid
> > > >>>>
> > > >>>>> with a larger extent (generated on a lattice that is 10 voxels
> > larger
> > > >>>>>>>
> > > >>>>>> all
> > > >>>>>
> > > >>>>>> around).
> > > >>>>>>>
> > > >>>>>>> To pre-empt Vladimir, itk_resample behaves the same way ;-)
> > > >>>>>>>
> > > >>>>>>> Bug, or feature? I'd say the former.
> > > >>>>>>>
> > > >>>>>>>
> > > >>
> > > >> --
> > > >> Best regards,
> > > >>
> > > >>  Vladimir S. FONOV ~ vladimir.fonov <at> gmail.com
> > > >> _______________________________________________
> > > >> MINC-users at bic.mni.mcgill.ca
> > > >> http://www.bic.mni.mcgill.ca/mailman/listinfo/minc-users
> > > >>
> > > > _______________________________________________
> > > > MINC-users at bic.mni.mcgill.ca
> > > > http://www.bic.mni.mcgill.ca/mailman/listinfo/minc-users
> > > _______________________________________________
> > > MINC-users at bic.mni.mcgill.ca
> > > http://www.bic.mni.mcgill.ca/mailman/listinfo/minc-users
> > >
> > _______________________________________________
> > MINC-users at bic.mni.mcgill.ca
> > http://www.bic.mni.mcgill.ca/mailman/listinfo/minc-users
> >
> >
> _______________________________________________
> MINC-users at bic.mni.mcgill.ca
> http://www.bic.mni.mcgill.ca/mailman/listinfo/minc-users
>