[Insight-users] Re: Gradient of vector fields -> how to compute eigensystem efficiently?

[Luis Ibanez]

Hi Lisa,


Thanks for your very well structured question.


Here are some suggestions on how you could implement this
processing using ITK components.


K = number of components per pixel
N = number of dimensions in the image
V = number of pixels


1) As you pointed out from the previous posting the mailing list, you
    can use ImageAdaptors or the VectorIndexSelectionCastImageFilter
    to extract the components K of the input vector image, and feed them
    into the GradientRecursiveGaussianImageFilter. This will give you
    the K rows of the "D" matrix, each one with N elements.

2) Using Iterators to visit this set of K images, you can, pixel-wise,
    compute the itkMatrix [D^T.D] of dimensions  NxN. The eigen analysis
    of this matrix can be solved using the SymmetricEigenAnalysis class:
http://www.itk.org/Insight/Doxygen/html/classitk_1_1SymmetricEigenAnalysis.html

    Note that if (N==2) then there is a closed form solution to the
    eigen analysis problem, and therefore you don't need to use the
    SymmetricEigenAnalysis class.


3) If you need to run this processing iteratively, you may find
    more convenient to modify the DerivativeImageFilter for taking
    an image of Vectors as input and producing an image of Matrix
    pixels as output. Then, in the same main loop of the GenerateData
    method, you could insert the computation of the eigen analysis.

    For an example on how to compute derivatives using finite differences
    you could look at the code of the itkLevelSetFunction and
    itkLevelSetImageFilter.




   Please let us know if you have questions about the details of
   how to use the classes suggested above.



      Regards,


         Luis



--------------------
lisat at sfu.ca wrote:
> Hello Luis,
> 
> I would like to calculate the gradient of a vector-valued image to be used for calculating the similarity between 2 vector images.
> 
> According to [2-3], we need to solve for the eigensystem of D*D' where D is the derivative of the image. Given k components per pixel, we could calculate the gradient per component per image dimension as suggested in your previous post [1]. Given the output from join-filter, what should be done next? 
> 
> Would an "ad-hoc" [2] approach be simply taking the average of the join-filter's output? 
> 
> Taking the correct approach, this requires solving for the eigenvector of D*D' per pixel. My questions then are:
> 
> 1. Could we make use of ITK filters to calculate/approximate D (i.e. DerivativeImageFilter, RecursiveGradient, etc.)? 
> 
> 2. How could I solve this eigensystem efficiently? (suppose I need to calculate this gradient image iteratively, this would require solving the system per iteration). For an k component, n dimension image of size v, this would require solving for [k x n]*[k x n] for v times, correct? Are there fast ways to approximate these calculations?
> 
> Many, many thanks in advance for your help.
> 
> Best,
> Lisa
> 
> References
> 
> [1] http://public.kitware.com/pipermail/insight-users/2003-October/005197.html
> 
> [2] http://csdl2.computer.org/persagen/DLAbsToc.jsp?resourcePath=/dl/proceedings/&toc=comp/proceedings/cvpr/1996/7258/00/7258toc.xml&DOI=10.1109/CVPR.1996.517146
> 
> [3] http://www.cs.sfu.ca/~hamarneh/ecopy/spiemi2005d.pdf
>