[Insight-developers] RE: Demons voodoo
Lydia Ng
lng at insightful.com
Thu, 8 Apr 2004 11:00:19 -0700
This is a multi-part message in MIME format.
------_=_NextPart_001_01C41D93.5ADDCF88
Content-Type: text/plain;
charset="us-ascii"
Content-Transfer-Encoding: quoted-printable
Hi Jim,
=20
>From my understanding of your email, what you are proposing is to =
estimate
the "forward mapping" instead of the "reverse mapping".=20
Both of which do allow you to warp the moving image to the fixed image - =
the
difference as you pointed out is that one keeps tracks of where a pixel =
goes
to and the other keeps track of where a pixel comes from.
=20
Out of curiosity, could you also repeat the experiment with the circle =
and
square swapped? I would be interested at the reverse behavior as well.
=20
For your forward mapping scheme, how do you compute (f-g)? Do you =
interpolate
f?
=20
As previously discussed in the list, "reverse mapping" is essential for =
image
resampling to avoid problems of holes (where no points map to an output =
pixel
- simple example is uniform scaling) and overlap (where several points =
of
different intensity maps to the same output pixel). While is it not =
strictly
essential for registration (because you can always interpolate f instead =
of
g), I have been also doing "reverse mapping" for registration so that =
the
mapping can be used for resampling.
=20
Although I don't have Thiron's technical with me right now - I believe =
it did
mention somewhere that the "reverse mapping" might be more convenient to
estimate (I will need to find the reference.)
=20
It could well be that "reverse mapping" is my personal soapbox but since =
the
implementation have shown to be useful for me for inter-subject brain =
mapping
- I really loathe to have it yanked out. Would it be possible to somehow =
have
both forward and reverse mapping?
=20
- Lydia=20
=20
=20
-----Original Message-----
From: Miller, James V (Research) [mailto:millerjv at crd.ge.com]=20
Sent: Thursday, April 08, 2004 8:12 AM
To: Lydia Ng; 'Insight-developers (E-mail)'
Subject: Demons voodoo
=20
(Original was bounced by the message size filter on the list.)
=20
Lydia,=20
=20
We were having some problems reconciling the output of the demons
registration filter with our understanding of the demons paper. I think =
I
have things worked out but it would require changing the code/behavior. =
I
have made the modfications locally and the effect of the results are in =
the
attached image. All of this is very complicated, having to keep track =
of
what the warps fields actually do, so I'll break my discussion up into =
stages
and use the "brain index".
=20
(1 brain needed to read this section): Demon's paper
=20
According to the demon's paper (page 15 and equation 4 on page 17), the =
goal
of the demons algorithm is to "make an image g appear more like an image =
f".
Equation 4 is essentially
=20
v =3D (g -f) grad(f) / (gradmag(f)^2 + (g - f))
=20
I interpret this to mean the vector v is the displacement that is needed =
to
move a position in image g such that it will lie on a pixel in image f =
that
has the same intensity as the pixel at the original position in g (this =
last
bit is just the basic premise of the demons algorithm).
=20
In other words, the displacement field calculated will tell you for each
pixel in g, where the corresponding pixel in f is.
=20
(2 brains needed to read this section): ITK implementation=20
=20
The demons algorithm is implemented using the finite difference engine, =
where
the output is the vector field of displacements. The finite difference =
engine
iterates over the output vector image, passing the index of vector pixel =
to
the DemonsRegistrationFunction. In pseudo-code, this function computes =
the
following:
=20
Calculate:
f(pi)
grad(f(pi))
g(pi + vi)
=20
dvi =3D (f - g) grad(f) / (gradmag(f)^2 + (f - g))
=20
Mapping the code to the equations above, f is the fixed image and g is =
the
moving image. Note that this is calculating a vector at a point pi that =
is
associated with image f. This seems to be calculating a vector field =
that is
keeping track of "where a pixel in f goes" as f is warped to g.
=20
(2 brains needed to read this section): Proposed ITK implementation=20
=20
Given the equation in the demons paper, I interpret the vector field to =
map a
position in image g to where in image f the corresponding pixel lies.
Therefore the output of the filter should be a vector field where the =
pixel
indices correspond to pixel indices in the moving image g and the vector
values at these indices say where in f to sample the values. I propose,
therefore, the code should be
=20
Calculate:
g(pi)
f(pi + v)
grad(f(pi+v))
=20
dvi =3D (g - f) grad(f) / (gradmag(f)^2 + (g-f))
=20
Again, the fixed image is f and the moving image is g. The output =
vector
field is index aligned with g and each vector gives a displacement to a
corresponding pixel in f. In essence, the vector field is keeping track =
of
"where a pixel in g came from" as g is warped to f.
=20
(3 brains needed to read this section): Warping using the demons =
displacement
field
=20
ITK's WarpImageFilter takes a displacement field that is index aligned =
with
the output of the WarpImageFilter and using the vector at each pixel to
determine where in the input image to probe to produce the output.
=20
In the current implementation, we pass the moving image (g) and the =
demons
displacement filed (d) to the WarpImageFilter to produce an "estimate" =
of the
fixed image.
=20
In my proposal, the deformation field is in the other direction, so we =
would
need to pass the fixed image (f) and the demons displacement field (g) =
to the
WarpImageFilter to produce an "estimate" of the moving image.
=20
Summary/Impact
=20
My main concern with the current implementation is that the vector field
computed does not match the mapping (g->f) of the original demons paper. =
So
people familar using the paper as a reference to how the filter should =
behave
can be confused (we have already proven that). The filter is currently
calculating vectors of where a pixel in f should go and opposed to where =
a
pixel in g came from. I don't think these are quite equivalent. As an
example, see the attached png image. I took a circle and warped it to a
square using the current ITK implementation and my proposed =
modifications.
Because the sense of the vector fields is different in the two
implementations, which image is considered the fixed or moving is =
different.
However, you can accomplish the end task with both implementations. The
first row of the figure shows the current implementation. As the circle
deforms into a square it actually changes topology twice (4 small hole =
appear
in the corners of the square and are then filled in). The second row of =
the
figure shows my alterations which seems more "regularized". However, =
the new
implementation requires more iterations to reach the destination.
=20
So the new implementation needs more iterations. It also requires that =
the
grad(f(pi+vi)) be calculated. Currently, I am using the same central
difference image function as before but now I pass in a physical =
coordinate
instead of an index. This image function will simply evaluate the =
gradient
at the nearest pixel. So overall, the proposed implementation requires =
more
computations. However, it looks to me like the deformation path is more
stable and the computed deformation field seems to agree (with me =
anyway)
with the original Demons paper.=20
=20
But perhaps more disruptively, the "sense" of the deformation field =
(moving
to fixed as opposed to fixed to moving) is different.
=20
=20
=20
Jim Miller=20
_____________________________________
Visualization & Computer Vision
GE Research
Bldg. KW, Room C218B
P.O. Box 8, Schenectady NY 12301
millerjv at research.ge.com
james.miller at research.ge.com
(518) 387-4005, Dial Comm: 8*833-4005,=20
Cell: (518) 505-7065, Fax: (518) 387-6981=20
=20
------_=_NextPart_001_01C41D93.5ADDCF88
Content-Type: text/html;
charset="us-ascii"
Content-Transfer-Encoding: quoted-printable
<!DOCTYPE HTML PUBLIC "-//W3C//DTD HTML 4.0 Transitional//EN">
<html xmlns:o=3D"urn:schemas-microsoft-com:office:office" =
xmlns:w=3D"urn:schemas-microsoft-com:office:word" =
xmlns:st1=3D"urn:schemas-microsoft-com:office:smarttags" =
xmlns=3D"http://www.w3.org/TR/REC-html40">
<head>
<META HTTP-EQUIV=3D"Content-Type" CONTENT=3D"text/html; =
charset=3Dus-ascii">
<meta name=3DProgId content=3DWord.Document>
<meta name=3DGenerator content=3D"Microsoft Word 10">
<meta name=3DOriginator content=3D"Microsoft Word 10">
<link rel=3DFile-List href=3D"cid:filelist.xml at 01C41D58.AD235D00">
<o:SmartTagType =
namespaceuri=3D"urn:schemas-microsoft-com:office:smarttags"
name=3D"country-region"/>
<o:SmartTagType =
namespaceuri=3D"urn:schemas-microsoft-com:office:smarttags"
name=3D"place"/>
<o:SmartTagType =
namespaceuri=3D"urn:schemas-microsoft-com:office:smarttags"
name=3D"time"/>
<o:SmartTagType =
namespaceuri=3D"urn:schemas-microsoft-com:office:smarttags"
name=3D"date"/>
<!--[if gte mso 9]><xml>
<o:OfficeDocumentSettings>
<o:DoNotRelyOnCSS/>
</o:OfficeDocumentSettings>
</xml><![endif]--><!--[if gte mso 9]><xml>
<w:WordDocument>
<w:SpellingState>Clean</w:SpellingState>
<w:GrammarState>Clean</w:GrammarState>
<w:DocumentKind>DocumentEmail</w:DocumentKind>
<w:EnvelopeVis/>
<w:BrowserLevel>MicrosoftInternetExplorer4</w:BrowserLevel>
</w:WordDocument>
</xml><![endif]--><!--[if !mso]>
<style>
st1\:*{behavior:url(#default#ieooui) }
</style>
<![endif]-->
<style>
<!--
/* Font Definitions */
at font-face
{font-family:Tahoma;
panose-1:2 11 6 4 3 5 4 4 2 4;
mso-font-charset:0;
mso-generic-font-family:swiss;
mso-font-pitch:variable;
mso-font-signature:1627421319 -2147483648 8 0 66047 0;}
at font-face
{font-family:"Comic Sans MS";
panose-1:3 15 7 2 3 3 2 2 2 4;
mso-font-charset:0;
mso-generic-font-family:script;
mso-font-pitch:variable;
mso-font-signature:647 0 0 0 159 0;}
/* Style Definitions */
p.MsoNormal, li.MsoNormal, div.MsoNormal
{mso-style-parent:"";
margin:0in;
margin-bottom:.0001pt;
mso-pagination:widow-orphan;
font-size:12.0pt;
font-family:"Times New Roman";
mso-fareast-font-family:"Times New Roman";}
a:link, span.MsoHyperlink
{color:blue;
text-decoration:underline;
text-underline:single;}
a:visited, span.MsoHyperlinkFollowed
{color:blue;
text-decoration:underline;
text-underline:single;}
p
{mso-margin-top-alt:auto;
margin-right:0in;
mso-margin-bottom-alt:auto;
margin-left:0in;
mso-pagination:widow-orphan;
font-size:12.0pt;
font-family:"Times New Roman";
mso-fareast-font-family:"Times New Roman";}
span.EmailStyle20
{mso-style-type:personal-reply;
mso-style-noshow:yes;
mso-ansi-font-size:10.0pt;
mso-bidi-font-size:10.0pt;
font-family:Arial;
mso-ascii-font-family:Arial;
mso-hansi-font-family:Arial;
mso-bidi-font-family:Arial;
color:navy;}
span.SpellE
{mso-style-name:"";
mso-spl-e:yes;}
span.GramE
{mso-style-name:"";
mso-gram-e:yes;}
at page Section1
{size:8.5in 11.0in;
margin:1.0in 1.25in 1.0in 1.25in;
mso-header-margin:.5in;
mso-footer-margin:.5in;
mso-paper-source:0;}
div.Section1
{page:Section1;}
-->
</style>
<!--[if gte mso 10]>
<style>
/* Style Definitions */=20
table.MsoNormalTable
{mso-style-name:"Table Normal";
mso-tstyle-rowband-size:0;
mso-tstyle-colband-size:0;
mso-style-noshow:yes;
mso-style-parent:"";
mso-padding-alt:0in 5.4pt 0in 5.4pt;
mso-para-margin:0in;
mso-para-margin-bottom:.0001pt;
mso-pagination:widow-orphan;
font-size:10.0pt;
font-family:"Times New Roman";}
</style>
<![endif]-->
</head>
<body lang=3DEN-US link=3Dblue vlink=3Dblue style=3D'tab-interval:.5in'>
<div class=3DSection1>
<p class=3DMsoNormal><font size=3D2 color=3Dnavy face=3DArial><span =
style=3D'font-size:
10.0pt;font-family:Arial;color:navy'>Hi =
Jim,<o:p></o:p></span></font></p>
<p class=3DMsoNormal><font size=3D2 color=3Dnavy face=3DArial><span =
style=3D'font-size:
10.0pt;font-family:Arial;color:navy'><o:p> </o:p></span></font></p>
<p class=3DMsoNormal><font size=3D2 color=3Dnavy face=3DArial><span =
style=3D'font-size:
10.0pt;font-family:Arial;color:navy'>From my understanding of your =
email, what
you are proposing is to estimate the “forward mapping” =
instead of
the “reverse mapping”. <o:p></o:p></span></font></p>
<p class=3DMsoNormal><font size=3D2 color=3Dnavy face=3DArial><span =
style=3D'font-size:
10.0pt;font-family:Arial;color:navy'>Both of which do allow you to warp =
the
moving image to the fixed image – the difference as you pointed =
out is
that one keeps tracks of where a pixel goes to and the other keeps track =
of
where a pixel comes from.<o:p></o:p></span></font></p>
<p class=3DMsoNormal><font size=3D2 color=3Dnavy face=3DArial><span =
style=3D'font-size:
10.0pt;font-family:Arial;color:navy'><o:p> </o:p></span></font></p>
<p class=3DMsoNormal><font size=3D2 color=3Dnavy face=3DArial><span =
style=3D'font-size:
10.0pt;font-family:Arial;color:navy'>Out of curiosity, could you also =
repeat
the experiment with the circle and square swapped? I would be interested =
at the
reverse behavior as well.<o:p></o:p></span></font></p>
<p class=3DMsoNormal><font size=3D2 color=3Dnavy face=3DArial><span =
style=3D'font-size:
10.0pt;font-family:Arial;color:navy'><o:p> </o:p></span></font></p>
<p class=3DMsoNormal><font size=3D2 color=3Dnavy face=3DArial><span =
style=3D'font-size:
10.0pt;font-family:Arial;color:navy'>For your forward mapping scheme, =
how do
you compute (f-g)? Do you interpolate f?<o:p></o:p></span></font></p>
<p class=3DMsoNormal><font size=3D2 color=3Dnavy face=3DArial><span =
style=3D'font-size:
10.0pt;font-family:Arial;color:navy'><o:p> </o:p></span></font></p>
<p class=3DMsoNormal><font size=3D2 color=3Dnavy face=3DArial><span =
style=3D'font-size:
10.0pt;font-family:Arial;color:navy'>As previously discussed in the =
list, “reverse
mapping” is essential for image <span =
class=3DSpellE>resampling</span> to
avoid problems of holes (where no points map to an output pixel – =
simple example
is uniform scaling) and overlap (where several points of different =
intensity
maps to the same output pixel). While is it not strictly essential for
registration (because you can always interpolate f instead of g), I have =
been
also doing “reverse mapping” for registration so that the =
mapping
can be used for <span =
class=3DSpellE>resampling</span>.<o:p></o:p></span></font></p>
<p class=3DMsoNormal><font size=3D2 color=3Dnavy face=3DArial><span =
style=3D'font-size:
10.0pt;font-family:Arial;color:navy'><o:p> </o:p></span></font></p>
<p class=3DMsoNormal><font size=3D2 color=3Dnavy face=3DArial><span =
style=3D'font-size:
10.0pt;font-family:Arial;color:navy'>Although I don’t have <span
class=3DSpellE>Thiron’s</span> technical with me right now - I =
believe it
did mention somewhere that the “reverse mapping” might be =
more convenient
to estimate (I will need to find the =
reference.)<o:p></o:p></span></font></p>
<p class=3DMsoNormal><font size=3D2 color=3Dnavy face=3DArial><span =
style=3D'font-size:
10.0pt;font-family:Arial;color:navy'><o:p> </o:p></span></font></p>
<p class=3DMsoNormal><font size=3D2 color=3Dnavy face=3DArial><span =
style=3D'font-size:
10.0pt;font-family:Arial;color:navy'>It could well be that =
“reverse
mapping” is my personal soapbox but since the implementation have =
shown
to be useful for me for inter-subject brain mapping – I really =
loathe to
have it yanked out. Would it be possible to somehow have both forward =
and
reverse mapping?<o:p></o:p></span></font></p>
<p class=3DMsoNormal><font size=3D2 color=3Dnavy face=3DArial><span =
style=3D'font-size:
10.0pt;font-family:Arial;color:navy'><o:p> </o:p></span></font></p>
<p class=3DMsoNormal><font size=3D2 color=3Dnavy face=3DArial><span =
style=3D'font-size:
10.0pt;font-family:Arial;color:navy'>- =
</span></font><st1:country-region><st1:place><font
size=3D2 color=3Dnavy face=3DArial><span =
style=3D'font-size:10.0pt;font-family:Arial;
color:navy'>Lydia</span></font></st1:place></st1:country-region><font =
size=3D2
color=3Dnavy face=3DArial><span =
style=3D'font-size:10.0pt;font-family:Arial;
color:navy'> <o:p></o:p></span></font></p>
<p class=3DMsoNormal><font size=3D2 color=3Dnavy face=3DArial><span =
style=3D'font-size:
10.0pt;font-family:Arial;color:navy'><o:p> </o:p></span></font></p>
<p class=3DMsoNormal><font size=3D2 color=3Dnavy face=3DArial><span =
style=3D'font-size:
10.0pt;font-family:Arial;color:navy'><o:p> </o:p></span></font></p>
<div style=3D'border:none;border-left:solid blue 1.5pt;padding:0in 0in =
0in 4.0pt'>
<p class=3DMsoNormal><font size=3D2 face=3DTahoma><span =
style=3D'font-size:10.0pt;
font-family:Tahoma'>-----Original Message-----<br>
<b><span style=3D'font-weight:bold'>From:</span></b> Miller, James V =
(Research)
[mailto:millerjv at crd.ge.com<span class=3DGramE>] <br>
<b><span style=3D'font-weight:bold'>Sent</span></b></span><b><span
style=3D'font-weight:bold'>:</span></b> </span></font><st1:date =
Month=3D"4" Day=3D"8"
Year=3D"2004"><font size=3D2 face=3DTahoma><span =
style=3D'font-size:10.0pt;font-family:
Tahoma'>Thursday, April 08, 2004</span></font></st1:date><font size=3D2
face=3DTahoma><span style=3D'font-size:10.0pt;font-family:Tahoma'> =
</span></font><st1:time
Hour=3D"8" Minute=3D"12"><font size=3D2 face=3DTahoma><span =
style=3D'font-size:10.0pt;
font-family:Tahoma'>8:12 AM</span></font></st1:time><font size=3D2 =
face=3DTahoma><span
style=3D'font-size:10.0pt;font-family:Tahoma'><br>
<b><span style=3D'font-weight:bold'>To:</span></b> Lydia Ng; =
'Insight-developers
(E-mail)'<br>
<b><span style=3D'font-weight:bold'>Subject:</span></b> Demons =
voodoo</span></font></p>
<p class=3DMsoNormal><font size=3D3 face=3D"Times New Roman"><span =
style=3D'font-size:
12.0pt'><o:p> </o:p></span></font></p>
<div>
<p class=3DMsoNormal><font size=3D2 face=3D"Times New Roman"><span =
style=3D'font-size:
10.0pt'>(Original was bounced by the message size filter on the =
list.)</span></font><o:p></o:p></p>
</div>
<div>
<p class=3DMsoNormal><font size=3D3 face=3D"Times New Roman"><span =
style=3D'font-size:
12.0pt'> <o:p></o:p></span></font></p>
</div>
<div>
<p class=3DMsoNormal><font size=3D2 face=3D"Times New Roman"><span =
style=3D'font-size:
10.0pt'>Lydia, </span></font><o:p></o:p></p>
</div>
<div>
<p class=3DMsoNormal><font size=3D3 face=3D"Times New Roman"><span =
style=3D'font-size:
12.0pt'> <o:p></o:p></span></font></p>
</div>
<div>
<p class=3DMsoNormal><font size=3D2 face=3D"Times New Roman"><span =
style=3D'font-size:
10.0pt'>We were having some problems reconciling the output of the =
demons
registration filter with our understanding of the demons paper. I think =
I have
things worked out but it would require changing the code/behavior. =
I have
made the modfications locally and the effect of the results are in the =
attached
image. All of this is very complicated, having to keep track of =
what the
warps fields actually do, so I'll break my discussion up into stages and =
use
the "brain index".</span></font><o:p></o:p></p>
</div>
<div>
<p class=3DMsoNormal><font size=3D3 face=3D"Times New Roman"><span =
style=3D'font-size:
12.0pt'> <o:p></o:p></span></font></p>
</div>
<div>
<p class=3DMsoNormal><strong><b><font size=3D2 face=3D"Times New =
Roman"><span
style=3D'font-size:10.0pt'>(1 brain needed to read this section): =
Demon's
paper</span></font></b></strong><o:p></o:p></p>
</div>
<div>
<p class=3DMsoNormal><font size=3D3 face=3D"Times New Roman"><span =
style=3D'font-size:
12.0pt'> <o:p></o:p></span></font></p>
</div>
<div>
<p class=3DMsoNormal><font size=3D2 face=3D"Times New Roman"><span =
style=3D'font-size:
10.0pt'>According to the demon's paper (page 15 and equation 4 on page =
17), the
goal of the demons algorithm is to "make an image g appear more =
like an
image f". Equation 4 is =
essentially</span></font><o:p></o:p></p>
</div>
<div>
<p class=3DMsoNormal><font size=3D3 face=3D"Times New Roman"><span =
style=3D'font-size:
12.0pt'> <o:p></o:p></span></font></p>
</div>
<div>
<p class=3DMsoNormal><font size=3D3 face=3D"Times New Roman"><span =
style=3D'font-size:
12.0pt'> </span></font><font size=3D2><span =
style=3D'font-size:
10.0pt'>v =3D (g -f) grad(f) / (gradmag(f)^2 + (g - =
f))</span></font><o:p></o:p></p>
</div>
<div>
<p class=3DMsoNormal><font size=3D3 face=3D"Times New Roman"><span =
style=3D'font-size:
12.0pt'> <o:p></o:p></span></font></p>
</div>
<div>
<p class=3DMsoNormal><font size=3D2 face=3D"Times New Roman"><span =
style=3D'font-size:
10.0pt'>I interpret this to mean the vector v is the =
displacement that is
needed to move a position in image g such that it will lie on a =
pixel in
image f that has the same intensity as the pixel at the =
original
position in g (this last bit is just the basic premise of the demons
algorithm).</span></font><o:p></o:p></p>
</div>
<div>
<p class=3DMsoNormal><font size=3D3 face=3D"Times New Roman"><span =
style=3D'font-size:
12.0pt'> <o:p></o:p></span></font></p>
</div>
<div>
<p class=3DMsoNormal><font size=3D2 face=3D"Times New Roman"><span =
style=3D'font-size:
10.0pt'>In other words, the displacement field calculated will tell you =
for
each pixel in g, where the corresponding pixel in f =
is.</span></font><o:p></o:p></p>
</div>
<div>
<p class=3DMsoNormal><font size=3D3 face=3D"Times New Roman"><span =
style=3D'font-size:
12.0pt'> <o:p></o:p></span></font></p>
</div>
<div>
<p class=3DMsoNormal><strong><b><font size=3D2 face=3D"Times New =
Roman"><span
style=3D'font-size:10.0pt'>(2 brains needed to read this section): ITK =
implementation</span></font> </b></strong><o:p></o:p></p>
</div>
<div>
<p class=3DMsoNormal><font size=3D3 face=3D"Times New Roman"><span =
style=3D'font-size:
12.0pt'> <o:p></o:p></span></font></p>
</div>
<div>
<p class=3DMsoNormal><font size=3D2 face=3D"Times New Roman"><span =
style=3D'font-size:
10.0pt'>The demons algorithm is implemented using the finite difference =
engine,
where the output is the vector field of displacements. The finite =
difference
engine iterates over the output vector image, passing the index of =
vector pixel
to the DemonsRegistrationFunction. In pseudo-code, this function =
computes
the following:</span></font><o:p></o:p></p>
</div>
<div>
<p class=3DMsoNormal><font size=3D3 face=3D"Times New Roman"><span =
style=3D'font-size:
12.0pt'> <o:p></o:p></span></font></p>
</div>
<div>
<p class=3DMsoNormal><font size=3D2 face=3D"Times New Roman"><span =
style=3D'font-size:
10.0pt'> Calculate:</span></font><o:p></o:p></p>
</div>
<div>
<p class=3DMsoNormal><font size=3D2 face=3D"Times New Roman"><span =
style=3D'font-size:
10.0pt'> f(pi)</span></font><o:p></o:p></p>
</div>
<div>
<p class=3DMsoNormal><font size=3D2 face=3D"Times New Roman"><span =
style=3D'font-size:
10.0pt'> grad(f(pi))</span></font><o:p></o:p></p>
</div>
<div>
<p class=3DMsoNormal><font size=3D3 face=3D"Times New Roman"><span =
style=3D'font-size:
12.0pt'> </span></font><font size=3D2><span =
style=3D'font-size:
10.0pt'>g(pi + vi)</span></font><o:p></o:p></p>
</div>
<div>
<p class=3DMsoNormal><font size=3D3 face=3D"Times New Roman"><span =
style=3D'font-size:
12.0pt'> <o:p></o:p></span></font></p>
</div>
<div>
<p class=3DMsoNormal><font size=3D3 face=3D"Times New Roman"><span =
style=3D'font-size:
12.0pt'> </span></font><font size=3D2><span =
style=3D'font-size:
10.0pt'>dvi =3D</span></font> <font size=3D2><span =
style=3D'font-size:10.0pt'> (f
- g) grad(f) / (gradmag(f)^2 + (f - g))</span></font><o:p></o:p></p>
</div>
<div>
<p class=3DMsoNormal><font size=3D3 face=3D"Times New Roman"><span =
style=3D'font-size:
12.0pt'> <o:p></o:p></span></font></p>
</div>
<div>
<p class=3DMsoNormal><font size=3D2 face=3D"Times New Roman"><span =
style=3D'font-size:
10.0pt'>Mapping the code to the equations above, f is the fixed image =
and g is
the moving image. Note that this is calculating a vector at a point pi =
that is
associated with image f. This seems to be calculating a vector =
field that
is keeping track of "where a pixel in f goes" as f is warped =
to g.</span></font><o:p></o:p></p>
</div>
<div>
<p class=3DMsoNormal><font size=3D3 face=3D"Times New Roman"><span =
style=3D'font-size:
12.0pt'> <o:p></o:p></span></font></p>
</div>
<div>
<div>
<p class=3DMsoNormal><strong><b><font size=3D2 face=3D"Times New =
Roman"><span
style=3D'font-size:10.0pt'>(2 brains needed to read this section): =
Proposed ITK
implementation </span></font></b></strong><font size=3D2><span
style=3D'font-size:10.0pt'><o:p></o:p></span></font></p>
</div>
<div>
<p class=3DMsoNormal><font size=3D2 face=3D"Times New Roman"><span =
style=3D'font-size:
10.0pt'> <o:p></o:p></span></font></p>
</div>
<div>
<p class=3DMsoNormal><font size=3D2 face=3D"Times New Roman"><span =
style=3D'font-size:
10.0pt'>Given the equation in the demons paper, I interpret the vector =
field to
map a position in image g to where in image f the corresponding pixel
lies. Therefore the output of the filter should be a vector field =
where
the pixel indices correspond to pixel indices in the moving image g and =
the
vector values at these indices say where in f to sample the values. I =
propose,
therefore, the code should be<o:p></o:p></span></font></p>
</div>
<div>
<p class=3DMsoNormal><font size=3D2 face=3D"Times New Roman"><span =
style=3D'font-size:
10.0pt'> <o:p></o:p></span></font></p>
</div>
<div>
<p class=3DMsoNormal><font size=3D2 face=3D"Times New Roman"><span =
style=3D'font-size:
10.0pt'> Calculate:<o:p></o:p></span></font></p>
</div>
<div>
<p class=3DMsoNormal><font size=3D2 face=3D"Times New Roman"><span =
style=3D'font-size:
10.0pt'> g(pi)<o:p></o:p></span></font></p>
</div>
<div>
<p class=3DMsoNormal><font size=3D2 face=3D"Times New Roman"><span =
style=3D'font-size:
10.0pt'> f(pi + v)<o:p></o:p></span></font></p>
</div>
<div>
<p class=3DMsoNormal><font size=3D2 face=3D"Times New Roman"><span =
style=3D'font-size:
10.0pt'> grad(f(pi+v))<o:p></o:p></span></font></p>
</div>
<div>
<p class=3DMsoNormal><font size=3D2 face=3D"Times New Roman"><span =
style=3D'font-size:
10.0pt'> <o:p></o:p></span></font></p>
</div>
<div>
<p class=3DMsoNormal><font size=3D2 face=3D"Times New Roman"><span =
style=3D'font-size:
10.0pt'> dvi =3D (g - f) grad(f) / (gradmag(f)^2 + =
(g-f))<o:p></o:p></span></font></p>
</div>
<div>
<p class=3DMsoNormal><font size=3D2 face=3D"Times New Roman"><span =
style=3D'font-size:
10.0pt'> <o:p></o:p></span></font></p>
</div>
<div>
<p class=3DMsoNormal><font size=3D2 face=3D"Times New Roman"><span =
style=3D'font-size:
10.0pt'>Again, the fixed image is f and the moving image is g. The =
output
vector field is index aligned with g and each vector gives a =
displacement to a
corresponding pixel in f. In essence, the vector field is keeping track =
of
"where a pixel in g came from" as g is warped to =
f.<o:p></o:p></span></font></p>
</div>
<div>
<p class=3DMsoNormal><font size=3D2 face=3D"Times New Roman"><span =
style=3D'font-size:
10.0pt'> <o:p></o:p></span></font></p>
</div>
<div>
<div>
<p class=3DMsoNormal><strong><b><font size=3D2 face=3D"Times New =
Roman"><span
style=3D'font-size:10.0pt'>(3 brains needed to read this section): =
Warping using
the demons displacement field</span></font></b></strong><font =
size=3D2><span
style=3D'font-size:10.0pt'><o:p></o:p></span></font></p>
</div>
<div>
<p class=3DMsoNormal><font size=3D2 face=3D"Times New Roman"><span =
style=3D'font-size:
10.0pt'> <o:p></o:p></span></font></p>
</div>
<div>
<p class=3DMsoNormal><font size=3D2 face=3D"Times New Roman"><span =
style=3D'font-size:
10.0pt'>ITK's WarpImageFilter takes a displacement field that is index =
aligned
with the output of the WarpImageFilter and using the vector at each =
pixel to
determine where in the input image to probe to produce the =
output.<o:p></o:p></span></font></p>
</div>
<div>
<p class=3DMsoNormal><font size=3D2 face=3D"Times New Roman"><span =
style=3D'font-size:
10.0pt'> <o:p></o:p></span></font></p>
</div>
<div>
<p class=3DMsoNormal><font size=3D2 face=3D"Times New Roman"><span =
style=3D'font-size:
10.0pt'>In the current implementation, we pass the moving image (g) and =
the
demons displacement filed (d) to the WarpImageFilter to produce an
"estimate" of the fixed image.<o:p></o:p></span></font></p>
</div>
<div>
<p class=3DMsoNormal><font size=3D2 face=3D"Times New Roman"><span =
style=3D'font-size:
10.0pt'> <o:p></o:p></span></font></p>
</div>
<div>
<p class=3DMsoNormal><font size=3D2 face=3D"Times New Roman"><span =
style=3D'font-size:
10.0pt'>In my proposal, the deformation field is in the other direction, =
so we
would need to pass the fixed image (f) and the demons displacement field =
(g) to
the WarpImageFilter to produce an "estimate" of the moving =
image.<o:p></o:p></span></font></p>
</div>
<div>
<p class=3DMsoNormal><font size=3D2 face=3D"Times New Roman"><span =
style=3D'font-size:
10.0pt'> <o:p></o:p></span></font></p>
</div>
<div>
<p class=3DMsoNormal><strong><b><font size=3D2 face=3D"Times New =
Roman"><span
style=3D'font-size:10.0pt'>Summary/Impact</span></font></b></strong><font=
size=3D2><span
style=3D'font-size:10.0pt'><o:p></o:p></span></font></p>
</div>
<div>
<p class=3DMsoNormal><font size=3D2 face=3D"Times New Roman"><span =
style=3D'font-size:
10.0pt'> <o:p></o:p></span></font></p>
</div>
<div>
<p class=3DMsoNormal><font size=3D2 face=3D"Times New Roman"><span =
style=3D'font-size:
10.0pt'>My main concern with the current implementation is =
that the
vector field computed does not match the mapping (g->f) of the =
original
demons paper. So people familar using the paper as a reference to =
how the
filter should behave can be confused (we have already proven that). The =
filter
is currently calculating vectors of where a pixel in f should go and =
opposed to
where a pixel in g came from. I don't think these are quite =
equivalent.
As an example, see the attached png image. I took a circle and =
warped it
to a square using the current ITK implementation and my proposed =
modifications.
Because the sense of the vector fields is different in the two =
implementations,
which image is considered the fixed or moving is different. =
However, you
can accomplish the end task with both implementations. The first =
row of
the figure shows the current implementation. As the circle deforms =
into a
square it actually changes topology twice (4 small hole appear in the =
corners
of the square and are then filled in). The second row of the =
figure shows
my alterations which seems more "regularized". However, =
the new
implementation requires more iterations to reach the =
destination.<o:p></o:p></span></font></p>
</div>
<div>
<p class=3DMsoNormal><font size=3D2 face=3D"Times New Roman"><span =
style=3D'font-size:
10.0pt'> <o:p></o:p></span></font></p>
</div>
<div>
<p class=3DMsoNormal><font size=3D2 face=3D"Times New Roman"><span =
style=3D'font-size:
10.0pt'>So the new implementation needs more iterations. It also =
requires that
the grad(f(pi+vi)) be calculated. Currently, I am using the same =
central
difference image function as before but now I pass in a physical =
coordinate
instead of an index. This image function will simply evaluate the
gradient at the nearest pixel. So overall, the proposed implementation =
requires
more computations. However, it looks to me like the deformation =
path is
more stable and the computed deformation field seems to agree (with me =
anyway)
with the original Demons paper. <o:p></o:p></span></font></p>
</div>
<div>
<p class=3DMsoNormal><font size=3D2 face=3D"Times New Roman"><span =
style=3D'font-size:
10.0pt'> <o:p></o:p></span></font></p>
</div>
<div>
<p class=3DMsoNormal><font size=3D2 face=3D"Times New Roman"><span =
style=3D'font-size:
10.0pt'>But perhaps more disruptively, the "sense" of the =
deformation
field (moving to fixed as opposed to fixed to moving) is =
different.<o:p></o:p></span></font></p>
</div>
<div>
<p class=3DMsoNormal><font size=3D3 face=3D"Times New Roman"><span =
style=3D'font-size:
12.0pt'> <o:p></o:p></span></font></p>
</div>
</div>
</div>
<div>
<p class=3DMsoNormal><font size=3D3 face=3D"Times New Roman"><span =
style=3D'font-size:
12.0pt'> <o:p></o:p></span></font></p>
</div>
<div>
<p class=3DMsoNormal><font size=3D3 face=3D"Times New Roman"><span =
style=3D'font-size:
12.0pt'> <o:p></o:p></span></font></p>
</div>
<div>
<p style=3D'margin:0in;margin-bottom:.0001pt'><b><font size=3D3 =
color=3Dnavy
face=3D"Comic Sans MS"><span =
style=3D'font-size:12.0pt;font-family:"Comic Sans MS";
color:navy;font-weight:bold'>Jim Miller</span></font></b> <br>
<b><i><font size=3D2 color=3Dred face=3DArial><span =
style=3D'font-size:10.0pt;
font-family:Arial;color:red;font-weight:bold;font-style:italic'>_________=
____________________________</span></font></i></b><br>
<em><i><font size=3D1 color=3Dblack face=3DArial><span =
style=3D'font-size:7.5pt;
font-family:Arial;color:black'>Visualization & Computer =
Vision</span></font></i></em><i><font
size=3D1 color=3Dblack face=3DArial><span =
style=3D'font-size:7.5pt;font-family:Arial;
color:black;font-style:italic'><br>
<em><i><font face=3DArial><span style=3D'font-family:Arial'>GE =
Research</span></font></i></em><br>
<em><i><font face=3DArial><span style=3D'font-family:Arial'>Bldg. KW, =
Room C218B</span></font></i></em><br>
<em><i><font face=3DArial><span style=3D'font-family:Arial'>P.O. Box 8, =
Schenectady
NY 12301</span></font></i></em><br>
<br>
</span></font></i><em><i><u><font size=3D1 color=3Dblue face=3D"Times =
New Roman"><span
style=3D'font-size:7.5pt;color:blue'><a =
href=3D"mailto:millerjv at research.ge.com">millerjv at research.ge.com</a></sp=
an></font></u></i></em><o:p></o:p></p>
<p style=3D'margin:0in;margin-bottom:.0001pt'><em><i><u><font size=3D1 =
color=3Dblue
face=3D"Times New Roman"><span =
style=3D'font-size:7.5pt;color:blue'>james.miller at research.ge.com</span><=
/font></u></i></em><br>
<i><font size=3D1 color=3Dblack face=3DArial><span =
style=3D'font-size:7.5pt;font-family:
Arial;color:black;font-style:italic'>(518) 387-4005, Dial Comm: =
8*833-4005, </span></font></i><br>
<i><font size=3D1 color=3Dblack face=3DArial><span =
style=3D'font-size:7.5pt;font-family:
Arial;color:black;font-style:italic'>Cell: (518) 505-7065, Fax: (518) =
387-6981</span></font></i>
<o:p></o:p></p>
</div>
<div>
<p class=3DMsoNormal><font size=3D3 face=3D"Times New Roman"><span =
style=3D'font-size:
12.0pt'> <o:p></o:p></span></font></p>
</div>
</div>
</div>
</body>
</html>
------_=_NextPart_001_01C41D93.5ADDCF88--