[Insight-users] invert deformable field
Luis Ibanez
luis.ibanez at kitware.com
Tue, 17 Feb 2004 02:55:12 -0500
Hi Corinne,
A) Actually I made a naive assumption when
I tell you that you could use only
landmarks with deformation larger than
a threshold. This is not correct, instead,
the right thing to do is to use a 'density'
of landmarka that is proportional to the
variation of the deformation. That is, you
need landmarks in places where the deformation
is changing abruptly.
B) Yes, that's fine. The point q is the backward
mapping of the point p, like q = B( p ).
you look at the neighbors of q, and use their
forward deformation to see if you can get closer
to p. so you look at
F( q1 )
F( q2 )
F( q3 )
F( q4 )
and for all of them you check which one is
closere to P.
C) Great, thanks for posting your graphics.
The prediction is not going to be very
precise if you use only the two previous
deformation field snapshots. You could use
three field stages in order to make a second
order estimation of value of a pixel for the
next cycle
Regards,
Luis
---------------------
Corinne Mattmann wrote:
> Hi Luis,
>
> Thanks very much for your answer!
>
> [Option A]
>
> I tried to change the DeformationFieldInverseImageFilter to just
> consider deformations larger than a certain threshold but it didn't
> work. Here is what I did:
>
> while( !ot.IsAtEnd() )
> {
> value = ot.Get();
>
> float tmp = 0;
> for(unsigned int i=0; i < ImageDimension; i++)
> {
> tmp += value[i] * value[i];
> }
> tmp = sqrt(tmp);
> if(tmp > 0.05)
> {
> sampledInput->TransformIndexToPhysicalPoint( ot.GetIndex(),
> sourcePoint );
> source->InsertElement( landmarkId, sourcePoint );
>
> for(unsigned int i=0; i < ImageDimension; i++)
> {
> targetPoint[i] = -value[i];
> }
> target->InsertElement( landmarkId, targetPoint ); // revert
> direction of displacement
>
> ++landmarkId;
> }
> ++ot;
> }
>
> There are three strange things this filter does when I input these
> lines:
> - The filter executes for quite a while and then suddenly it outputs a
> whole bunch of numbers on the screen (every second number is a zero).
> This happens just if the actual number of landmarks is smaller than the
> maximum number of landmarks (numberOfLandmarks).
> - The result I get then seems to consist only of infinite numbers even
> if the actual number of landmarks is just slightly smaller than the
> maximum number of landmarks.
> - If I change the subsampling factor to a smaller number (e.g. 4 again),
> it still crashes even though the number of landmarks is smaller than
> with a subsampling factor of 6 and taking all landmarks.
>
>
> [Option B]
>
> I think this is an interesting option. There is just one thing I don't
> understand:
> How can I get the forward mapping from exactly the points q1, q2, q3 and
> q4? If I apply F to q, I cannot choose which pixels I want to transform,
> I just get the pixel intensity from where the field F points to and
> write it to a chosen pixel in p.
> May be F has to be applied first and then B? But I guess that B then
> doesn't necessarily have to be the same as when you apply B first.
>
>
> [Option C]
>
> I had the same idea. But before implementing I wanted to get an insight
> of how the deformable field looks like: I registered two images with
> varying values just in x-direction (in y- and z-direction all values
> were the same), in other words I registered two plates with a different
> thickness. Like this I get a deformation field in just one direction
> which I can easily visualize.
> To see how the deformable field changes over time I made a second
> registration beginning with the second "plate" to a third one. In the
> attached file "images.png" you see the values in x-direction of my three
> images I registered and in the file "deffield.png" the two resulting
> deformable fields (one pixel was 1mm). If you look at those deformable
> fields I think it would be quite hard to find an appropriate prediction
> for a third one.
>
>
> Thanks very much,
> Corinne
>
>
>
> -----Original Message-----
> From: Luis Ibanez [mailto:luis.ibanez at kitware.com]
> Sent: Thursday, February 12, 2004 3:52 PM
> To: Corinne Mattmann
> Cc: insight-users at itk.org; Harald Grossauer
> Subject: Re: [Insight-users] invert deformable field
>
>
>
> Hi Corinne,
>
>
> Subsampling the input field by a small factor
> will result in a quite large number of landmarks
> that you probably don't need. It is likely that
> memory allocation is failing at that point.
>
>
> Here are some options:
>
> [Option A]
>
> Probably we should create a variant of this
> class with a smarter sampling that instead of
> using a regular grid, sets landmarks in the
> regions of large deformation. That will make
> a better use of the memory and the landmarks.
>
> A simple option is that after the regular grid
> resampling, we could visit all the generated
> landmarks and eliminate all those having deformation
> lower than a certain threshold. The KernelSpline
> has the great advantage of not requiring the landmarks
> to be in a regular grid.
>
> You could do this elimination in the method
> PrepareKernelBaseSpline() in the "while" loop
> on lines 191-207. Simply check the magnitude
> of "value" to be over the threshold before
> creating a landmark with it.
>
>
>
> [Option B]
>
> Another interesting option is what Harald suggested
> at some point which is to use this class as just
> an initializer for a more refined method that will
> iteratively increase the precision of the inverse
> field. Some sort of redistribution of error could
> do the job. That is, if F is the Forward field
> and B is the estimation of the backward field, we
> would expect for a point p in the destination
> space to satisfy
>
> p == F( B( p ) )
>
> in practice of course, there will be some error
>
> Dp = p - F( B( P ) )
>
> We could redistribute this error in B by cheking
> which one of the neighbors of B( p ) will have
> a forward mapping that is closer to p.
>
> That is, in 2D for example, you locate q = B( p )
> then check the points
>
> q1 = q + ( 1, 0)
> q2 = q + (-1, 0)
> q3 = q + ( 0,-1)
> q4 = q + ( 0, 1)
>
> and find which one has a forward mapping closer to p.
> So you test F(q1), F(q2), F(q3), F(q4), and once you
> figure out which one is a better candidate, then
> you modify the backward deformation B(p) at pixel p
> by and ammount that will bring q closer to the best
> candidate. Note that you can refine one pixel at
> a time or do an entire visit over B before comming
> back to the same pixel. You could also carry on
> corrections, that is, take the correction of the
> previous pixel as a first estimation of the correction
> needed in the current pixel.
>
>
>
> [Option C]
>
> A third option is to not follow the evolution of
> the forward field by walkig along the deformations
> but just by staying in the same pixel. In this
> approach you simply take for the same pixel all
> the previous forward deformation vector and fit
> a curve to them, then extrapolate for the next
> time step. That's probably simpler to implement
> than the two approaches above,... but it may also
> be less meaningful.
>
>
>
> -----
>
>
> I would suggest you to try first the approach
> of reducing the number or landmakrs with a
> deformation threshold (let say eliminate all
> the landmaks that deform their pixel image by
> less than 1mm). This threshold of could will
> be set by the user through a Set() method.
>
>
>
>
> Please let us know if you find this to be useful
> so we can do the equivalent modifications in the
> repository.
>
>
>
> Thanks,
>
>
>
> Luis
>
>
> -------------------------------------
>
> Corinne Mattmann wrote:
>
>
>>Hi,
>>
>>I would like to predict a future deformable field from several fields
>>from the past. It is not too difficult to calculate the future shift
>>in a pixel by going through the deformable fields, taking always the
>>closest point (as a first approach). But then I just have the vector
>>where I have to shift this pixel intensity to ("forward field").
>>However to use the warpImageFilter and to actually apply this field I
>>need a "backward field". Therefore I have to invert the field I get. I
>
>
>>tried using the DeformationFieldInverseImageFilter (is this what I
>>need?) but I had to use a subsampling factor of 6 (for a 50x50x50
>>image) to get an acceptable speed which is not at all accurate enough
>>for my images. When I chose a smaller value (4), the program crashed.
>>Do you have any ideas how else I could get an inversed deformable
>>field? Or is there even another possibility to predict a deformable
>>field from several others?
>>
>>Thanks very much,
>>Corinne
>>
>>_______________________________________________
>>Insight-users mailing list
>>Insight-users at itk.org
>>http://www.itk.org/mailman/listinfo/insight-users
>>
>
>
>
>
> ------------------------------------------------------------------------
>
>
> ------------------------------------------------------------------------
>