[Insight-developers] itkIsoContourDistanceImageFilter, Win32,
Fixed, but additional th reading issues
Miller, James V (Research)
millerjv at crd . ge . com
Thu, 21 Aug 2003 09:20:21 -0400
This message is in MIME format. Since your mail reader does not understand
this format, some or all of this message may not be legible.
------_=_NextPart_001_01C367E6.F8BDA7EE
Content-Type: text/plain
Karl,
I think the method of looking out 2n neighbors and only updating itself is
the best solution right now. This is a difficult issue for the threading
models to address. The safest solution is always to set up the threads so
there is no possibility that multiple threads will try to write to the same
output pixel.
I suppose you could break of the narrow band into several groups of pixels.
For instance, groups A1, A2, A3, A4, B.
The A* groups would be pixels in the narrow band that can safely be
processed in separated threads with no possibility of separate threads
writing to the same output pixel. The B group will be all those pixels in
the "border" regions between say A1 and A2. The A* pixels would be sent off
to multiple threads and once they completed, the B group would processed
serially. The question is whether you can define a simple algorithm for
assigning the pixels in the narrow band to the appropriate groups. Some
sort of K-means clustering might give a first pass for the grouping, then
you could erode away the interface between groups and put in the eroded
pixels in the B group.
We have been tracking down a number of these problems recently. I've been
thinking about whether it would be possible to construct an image class (or
adaptor) or maybe a set of iterators that would naively safe in a
multithreading environment. As you mentioned, you could have a mutex per
pixel but that would be prohibitive. Another option is for an NxM image
have N+M mutexes, one for each row and one for each column. To write to a
pixel, you would have to lock a row and a column. I'm not sure whether this
would result in deadlocks if a thread could lock a row but not a column. So
you would have fewer mutex but have to lock more than one to write to a
pixel. Another option is to break the image into a series of blocks (say
with the number of blocks equal the number of threads) and have a mutex per
block.
Jim
-----Original Message-----
From: Karl Krissian [mailto:karl at bwh . harvard . edu]
Sent: Wednesday, August 20, 2003 6:26 PM
To: Miller, James V (Research)
Cc: Insight-developers (E-mail); 'Raul San Jose Estepar'
Subject: Re: [Insight-developers] itkIsoContourDistanceImageFilter, Win32,
Fixed, but additional th reading issues
Hi James,
Thank you for looking and fixing these problems.
The code for setting the pixel as either +FarValue, -FarValue or 0 is fine.
It will create a processing overhead that I was avoiding in my initial
level set implementation in VTK by using the same image for input and output
(the input is initialized on the whole image when starting the level set).
One solution to keep the filter fast could be to inherit from an
itkInPlaceImageFilter.
I tried to do it for the itkFastChamferDistanceImageFilter but ran into some
problem with the neighborhood iterators.
Concerning the thread problem, you are completely right to notice
that two threads may update the same voxel.
The level set evolution would probably not suffer a lot from this
problem because there is an inherent difficulty in approximating
the distance from the interpolated surface to the grid points.
When a point has several local estimates of the distance
to the surface, I decided to keep the minimal distance,
but whichever is chosen will displace the interpolated
surface from its original position.
Alternatives could be:
- to put semaphore per each voxel but I am afraid it would require too much
memory
for a little gain.
- to change the code so that each voxel looks at 2n neighbors instead of n
and only updates itself (it would be twice slower).
Karl Krissian, PhD
Laboratory of Mathematics in Imaging (LMI)
Harvard Medical School, Brigham and Women's Hospital
Tel:617-278-0639, Fax:617-264-6887
Miller, James V (Research) wrote:
I modified the itkIsoContourDistanceImageFilter so the test will pass on
Windows. The test for this class checks to make sure the pixels outside the
levelset and inside the level before and after the reinitialization have the
same sign. In the narrow band case, not all the output pixels were visited
so there the comparisons were using uninitialized memory. This caused the
test to fail. I put code in the BeforeThreadedGenerateData() method to set
each pixel in the output levelset as either +FarValue, -FarValue or 0.
There is another problem with this algorithm that could appear when running
on multiple threads in the narrow band case. The filter divides the narrow
band (list of indices) into a set of groups, specifying one group for each
thread. This is nice. However, the algorithm may set pixels in a 3x3x...
neighborhood of a narrow band pixel. When running multiple threads this may
be a problem. If threads of examining adjacent pixels (or pixels with a
radius of 2), the two threads may end up trying to set the same output pixel
at the same time. This can cause a problem with memory integrity. The
situation here could actually be more subtle because before each write
access to an output pixel, a comparison is done with a current calculated
value.
if(fabs(val1_new) < fabs(outNeigIt.GetNext(n,1)))
outNeigIt.SetNext(n,1,static_cast<PixelType>(val1_new) );
Consider a thread that is interrupted between the call to GetNext() and
SetNext(). If the conditional is evaluated as true,
then the code wants to overwrite the pixel. But if the thread is
interrupted after the conditional is satisfied but before the call to
SetNext(), another thread may come in and evaluate the same output pixel,
pass its conditional with a value of val1_new that is less than the first
thread, write the pixel out. Then when the first thread awakes, it has
already satisfied its conditional so it writes it version of val1_new. The
the first thread had a value of val1_new that is greater than the second
thread's, the output will be incorrect.
Jim Miller
_____________________________________
Visualization & Computer Vision
GE Research
Bldg. KW, Room C218B
P.O. Box 8, Schenectady NY 12301
millerjv at research . ge . com <mailto:millerjv at research . ge . com>
james . miller at research . ge . com <mailto:james . miller at research . ge . com>
(518) 387-4005, Dial Comm: 8*833-4005,
Cell: (518) 505-7065, Fax: (518) 387-6981
------_=_NextPart_001_01C367E6.F8BDA7EE
Content-Type: text/html
<!DOCTYPE HTML PUBLIC "-//W3C//DTD HTML 4.0 Transitional//EN">
<HTML><HEAD>
<META HTTP-EQUIV="Content-Type" CONTENT="text/html; charset=US-ASCII">
<TITLE></TITLE>
<META content="MSHTML 6.00.2800.1141" name=GENERATOR></HEAD>
<BODY>
<DIV><SPAN class=139115912-21082003><FONT color=#0000ff
size=2>Karl,</FONT></SPAN></DIV>
<DIV><SPAN class=139115912-21082003><FONT color=#0000ff
size=2></FONT></SPAN> </DIV>
<DIV><SPAN class=139115912-21082003><FONT color=#0000ff size=2>I think the
method of looking out 2n neighbors and only updating itself is the best solution
right now. This is a difficult issue for the threading models to address. The
safest solution is always to set up the threads so there is no possibility that
multiple threads will try to write to the same output pixel.
</FONT></SPAN></DIV>
<DIV><SPAN class=139115912-21082003><FONT color=#0000ff
size=2></FONT></SPAN> </DIV>
<DIV><SPAN class=139115912-21082003><FONT color=#0000ff size=2>I suppose you
could break of the narrow band into several groups of pixels. For
instance, groups A1, A2, A3, A4, B.</FONT></SPAN></DIV>
<DIV><SPAN class=139115912-21082003><FONT color=#0000ff size=2>The A* groups
would be pixels in the narrow band that can safely be processed in separated
threads with no possibility of separate threads writing to the same output
pixel. The B group will be all those pixels in the "border" regions
between say A1 and A2. The A* pixels would be sent off to multiple threads
and once they completed, the B group would processed serially. The
question is whether you can define a simple algorithm for assigning the pixels
in the narrow band to the appropriate groups. Some sort of K-means
clustering might give a first pass for the grouping, then you could erode away
the interface between groups and put in the eroded pixels in the B
group.</FONT></SPAN></DIV>
<DIV><SPAN class=139115912-21082003><FONT color=#0000ff
size=2></FONT></SPAN> </DIV>
<DIV><SPAN class=139115912-21082003><FONT color=#0000ff size=2>We have been
tracking down a number of these problems recently. I've been thinking
about whether it would be possible to construct an image class (or adaptor) or
maybe a set of iterators that would naively safe in a multithreading
environment. As you mentioned, you could have a mutex per pixel but that
would be prohibitive. Another option is for an NxM image have N+M mutexes,
one for each row and one for each column. To write to a pixel, you would
have to lock a row and a column. I'm not sure whether this would result in
deadlocks if a thread could lock a row but not a column. So you would have
fewer mutex but have to lock more than one to write to a pixel. Another
option is to break the image into a series of blocks (say with the number of
blocks equal the number of threads) and have a mutex per block.
</FONT></SPAN></DIV>
<DIV><SPAN class=139115912-21082003><FONT color=#0000ff
size=2></FONT></SPAN> </DIV>
<DIV><SPAN class=139115912-21082003><FONT color=#0000ff
size=2>Jim</FONT></SPAN></DIV>
<DIV><SPAN class=139115912-21082003><FONT color=#0000ff
size=2></FONT></SPAN> </DIV>
<DIV><SPAN class=139115912-21082003><FONT color=#0000ff
size=2></FONT></SPAN> </DIV>
<BLOCKQUOTE
style="PADDING-LEFT: 5px; MARGIN-LEFT: 5px; BORDER-LEFT: #0000ff 2px solid">
<DIV class=OutlookMessageHeader dir=ltr align=left><FONT face=Tahoma
size=2>-----Original Message-----<BR><B>From:</B> Karl Krissian
[mailto:karl at bwh . harvard . edu]<BR><B>Sent:</B> Wednesday, August 20, 2003 6:26
PM<BR><B>To:</B> Miller, James V (Research)<BR><B>Cc:</B> Insight-developers
(E-mail); 'Raul San Jose Estepar'<BR><B>Subject:</B> Re: [Insight-developers]
itkIsoContourDistanceImageFilter, Win32, Fixed, but additional th reading
issues<BR><BR></FONT></DIV><BR>Hi James,<BR><BR>Thank you for looking and
fixing these problems.<BR><BR>The code for setting the pixel as either
+FarValue, -FarValue or 0 is fine.<BR><BR>It will create a processing overhead
that I was avoiding in my initial<BR>level set implementation in VTK by using
the same image for input and output<BR>(the input is initialized on the whole
image when starting the level set).<BR><BR>One solution to keep the filter
fast could be to inherit from an itkInPlaceImageFilter.<BR>I tried to do it
for the itkFastChamferDistanceImageFilter but ran into some<BR>problem with
the neighborhood iterators.<BR><BR>Concerning the thread problem, you are
completely right to notice<BR>that two threads may update the same
voxel.<BR><BR>The level set evolution would probably not suffer a lot from
this<BR>problem because there is an inherent difficulty in
approximating<BR>the distance from the interpolated surface to the grid
points.<BR>When a point has several local estimates of the distance <BR>to the
surface, I decided to keep the minimal distance,<BR>but whichever is chosen
will displace the interpolated<BR>surface from its original
position.<BR><BR>Alternatives could be: <BR>- to put semaphore per each voxel
but I am afraid it would require too much memory<BR>for a little gain.<BR>- to
change the code so that each voxel looks at 2n neighbors instead of
n<BR> and only updates itself (it would be twice
slower).<BR><BR><BR>Karl Krissian, PhD<BR>Laboratory of Mathematics in Imaging
(LMI)<BR>Harvard Medical School, Brigham and Women's
Hospital<BR>Tel:617-278-0639, Fax:617-264-6887<BR><BR><BR><BR>Miller, James V
(Research) wrote:<BR>
<BLOCKQUOTE
cite=midFBE90DFC240BA541B38A43F39913A16D04461347 at xmb02crdge . crd . ge . com
type="cite">
<META content="MSHTML 6.00.2800.1141" name=GENERATOR>
<DIV><SPAN class=865111014-20082003><FONT size=2>I modified the
itkIsoContourDistanceImageFilter so the test will pass on Windows. The
test for this class checks to make sure the pixels outside the levelset and
inside the level before and after the reinitialization have the same
sign. In the narrow band case, not all the output pixels were visited
so there the comparisons were using uninitialized memory. This caused
the test to fail. I put code in the BeforeThreadedGenerateData()
method to set each pixel in the output levelset as either +FarValue,
-FarValue or 0.</FONT></SPAN></DIV>
<DIV><SPAN class=865111014-20082003></SPAN> </DIV>
<DIV><SPAN class=865111014-20082003><FONT size=2>There is another problem
with this algorithm that could appear when running on multiple threads in
the narrow band case. The filter divides the narrow band (list of
indices) into a set of groups, specifying one group for each thread.
This is nice. However, the algorithm may set pixels in a 3x3x...
neighborhood of a narrow band pixel. When running multiple threads
this may be a problem. If threads of examining adjacent pixels (or
pixels with a radius of 2), the two threads may end up trying to set the
same output pixel at the same time. This can cause a problem with
memory integrity. The situation here could actually be more subtle
because before each write access to an output pixel, a comparison is done
with a current calculated value.</FONT></SPAN></DIV>
<DIV><SPAN class=865111014-20082003></SPAN> </DIV>
<DIV><SPAN class=865111014-20082003><FONT
size=2>
if(fabs(val1_new) <
fabs(outNeigIt.GetNext(n,1)))<BR>
outNeigIt.SetNext(n,1,static_cast<PixelType>(val1_new)
);<BR></FONT></SPAN></DIV>
<DIV><SPAN class=865111014-20082003><FONT size=2>Consider a thread that
is interrupted between the call to GetNext() and SetNext(). If the
conditional is evaluated as true, </FONT></SPAN></DIV>
<DIV><SPAN class=865111014-20082003><FONT size=2>then the code wants to
overwrite the pixel. But if the thread is interrupted after the
conditional is satisfied but before the call to SetNext(), another thread
may come in and evaluate the same output pixel, pass its conditional with a
value of val1_new that is less than the first thread, write the pixel
out. Then when the first thread awakes, it has already satisfied its
conditional so it writes it version of val1_new. The the first thread
had a value of val1_new that is greater than the second thread's, the output
will be incorrect.</FONT></SPAN></DIV>
<DIV><SPAN class=865111014-20082003></SPAN> </DIV>
<DIV><SPAN class=865111014-20082003><FONT size=2></FONT></SPAN> </DIV>
<DIV><SPAN class=865111014-20082003></SPAN> </DIV>
<DIV> </DIV>
<DIV>
<P style="MARGIN: 0in 0in 0pt"><B><SPAN
style="COLOR: navy; FONT-FAMILY: 'Comic Sans MS'">Jim Miller</SPAN></B>
<BR><B><I><SPAN
style="FONT-SIZE: 10pt; COLOR: red; FONT-FAMILY: Arial">_____________________________________</SPAN></I></B><BR><EM><SPAN
style="FONT-SIZE: 7.5pt; COLOR: black; FONT-FAMILY: Arial">Visualization
& Computer Vision</SPAN></EM><I><SPAN
style="FONT-SIZE: 7.5pt; COLOR: black; FONT-FAMILY: Arial"><BR><EM>GE
Research</EM><BR><EM>Bldg. KW, Room C218B</EM><BR><EM>P.O. Box 8,
Schenectady NY 12301</EM><BR><BR></SPAN></I><EM><U><SPAN
style="FONT-SIZE: 7.5pt; COLOR: blue"><A
href="mailto:millerjv at research . ge . com">millerjv at research . ge . com</A></SPAN></U></EM></P>
<P style="MARGIN: 0in 0in 0pt"><EM><U><SPAN
style="FONT-SIZE: 7.5pt; COLOR: blue"><A class=moz-txt-link-abbreviated
href="mailto:james . miller at research . ge . com">james . miller at research . ge . com</A></SPAN></U></EM><BR><I><SPAN
style="FONT-SIZE: 7.5pt; COLOR: black; FONT-FAMILY: Arial">(518) 387-4005,
Dial Comm: 8*833-4005, </SPAN></I><BR><I><SPAN
style="FONT-SIZE: 7.5pt; COLOR: black; FONT-FAMILY: Arial">Cell: (518)
505-7065, Fax: (518) 387-6981</SPAN></I> </P></DIV>
<DIV> </DIV></BLOCKQUOTE><BR></BLOCKQUOTE></BODY></HTML>
------_=_NextPart_001_01C367E6.F8BDA7EE--