[Insight-developers] Re: Kernel transforms. K symmetric?

Luis Ibanez luis.ibanez@kitware.com
Thu, 12 Dec 2002 15:03:49 -0500


Hi Jim,

Sorry I missed the point about K.
Thanks for clarifying it.

You are right,
ComputeG() is invoked for the computation
of every coefficient in K. (in ComputeK())

ComputeG() actually receives as parameter
the vector that results from the difference
between the two points (landmarks). So the
notation

             G(p1,p2)

actually becomes

             G(p1 - p2)

Looking again, here is what we
can probably say about K:


ElasticBodySpline           symmetric
ElasticBodyReciprocalSpline symmetric
VolumeSpline                symmetric
ThinPlateSpline             symmetric
ThinPlateR2LogSpline        symmetric


Your additions sound quite interesting.
In particular if the transforms are used for
deformable registration since in this way
the landmark will only have to be approximations
to the final position.

How about adding a family of Approximating
splines as you proposed them and simply overload
the ComputeK() method on them.

It could be for example:

TPS aproximating Spline, deriving from the
current TPS.


   Luis



========================================

Miller, James V (Research) wrote:
> Actually, I was thinking of the optimization more in constructing 
> the K matrix which calls ComputeG() for every pairing of landmarks.
> If all the kernels of symmetric, then for G for p1->p2 should be the
> same as G for p2->p1 and we can save "some" computation.
> 
> This is a side issue for what I am going to do. The code currently
> calculates K as 
> 
> G(p1, p1)  G(p1, p2), G(p1, p3), ....
> G(p2, p1)  G(p2, p2), G(p2, p3), ....
> 
> For a thin plate spline and it looks like for the volume spline, 
> the diagonal on the K matrix is zero.
> 
> If we replace the diagonal with some other constant, lambda, then
> the splines are no longer forced to interpolate the landmarks but
> instead are allowed to approximate the landmarks.  The lambda term
> is standard lagrange multiplier to trade off fidelity to the data 
> and a penalty for smoothness.
> 
> This can be generalized further so that instead of the same constant
> along the diagonal, each dxd block around the diagonal can be set to
> reflect a weight of confidence (or covariance) in the particular landmark.
> 
> The result is a weighted approximating thin plate (or other) spline.  I'd need
> to look at the EBS closer to determine if the same can done for it.
> 
> So I am thinking of adding a method that computes G(pi, pi) that by default
> calls G(pi, pj) and can be overridden in the subclass.
> 
> In ITKs implementation, the K matrix is (Nxd)x(Nxd).  I have some other 
> thin plate spline code where the K matrix is only (NxN). While this later code
> allows me to implement approximating splines, I can only weight each landmark
> by a single scalar, in other words, I cannot put a covariance matrix of a landmark
> along the diagonal.  So while ITKs implementation uses more memory, it is more
> general.  I was thinking of whether it would be worth it for the subclasses of the 
> KernelTransform to specify the number of "degrees of freedom" in the kernel.  The
> elastic body spline and would have d degrees of freedom in the kernel and the 
> thin plate and volume spline would only have 1 degree of freedom.  So the former
> would need the (Nxd)x(Nxd) space for the K matrix while the latter two splines
> would (usually) only need NxN space for the K matrix.
> 
> Jim
> 
> 
> 
> 
>>-----Original Message-----
>>From: Luis Ibanez [mailto:luis.ibanez@kitware.com]
>>Sent: Thursday, December 12, 2002 2:18 PM
>>To: Miller, James V (Research)
>>Cc: Insight-developers (E-mail)
>>Subject: Re: Kernel transforms. K symmetric?
>>
>>
>>Hi Jim,
>>
>>The Kernel matrix is noted as "G" in the code
>>
>>Here is what I can see:
>>
>>- ElasticBodySpline            is symmetric
>>- ElasticBodyReciprocalSpline  is symmetric
>>- VolumeSpline                 is diagonal
>>- ThinPlateSpline              is diagonal
>>- ThinPlateR2LogSpline         is diagonal
>>
>>
>>The ComputG( vector ) method is virtual
>>in KernelTranform. It is overloaded in
>>every particular KernelSplineTransform.
>>
>>The method is only invoked in
>>
>>    ComputeDeformationContribution()
>>
>>Here seems to be the place where we could
>>take advantage of the matrix being symmetric
>>or diagonal in order to reduce the number
>>of operations.
>>
>>ComputeDeformationContribution() is however,
>>already overloaded in :
>>
>>- ThinPlateR2LogSpline
>>- VolumeSpline
>>- ThinPlateSpline
>>
>>To take advantage of the diagonal property.
>>
>>It may be a matter then, of overloading
>>this method in the ElasticBody and
>>ElasticBodySpline in order to take advantage
>>of the symmetry.
>>
>>
>>Thanks for improving this code.
>>
>>
>>    Luis
>>
>>
>>----------------------------
>>
>>Miller, James V (Research) wrote:
>>
>>>Luis,
>>>
>>> 
>>>
>>>Is the K matrix of the various kernel transforms (thin 
>>>
>>plate, elastic 
>>
>>>body, volume spline) always symmetric?
>>>
>>> 
>>>
>>>I think it is for the TPS.  If it is always symmetric, I am 
>>>
>>going to 
>>
>>>rewrite it so that it only evaluates the upper triagular 
>>>
>>part of K and 
>>
>>>copies values into the lower triangle.
>>>
>>> 
>>>
>>>Jim   
>>>
>>> 
>>>
>>>
>>
>>
>