[Insight-developers] itkIsoContourDistanceImageFilter, Win32, Fixed, but additional th reading issues

Karl Krissian karl at bwh . harvard . edu
Wed, 20 Aug 2003 18:26:26 -0400


--------------070109060505020600090504
Content-Type: text/plain; charset=us-ascii; format=flowed
Content-Transfer-Encoding: 7bit


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
> (518) 387-4005, Dial Comm: 8*833-4005,
> Cell: (518) 505-7065, Fax: (518) 387-6981
>
>  



--------------070109060505020600090504
Content-Type: text/html; charset=us-ascii
Content-Transfer-Encoding: 7bit

<!DOCTYPE html PUBLIC "-//W3C//DTD HTML 4.01 Transitional//EN">
<html>
<head>
  <title></title>
</head>
<body>
<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>
&nbsp; 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 type="cite"
 cite="midFBE90DFC240BA541B38A43F39913A16D04461347 at xmb02crdge . crd . ge . com"> 
 
  <meta http-equiv="Content-Type" content="text/html; ">
   
  <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.&nbsp; 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.&nbsp;  In the narrow band case, not all
the output pixels were visited so there the  comparisons were using uninitialized
memory.&nbsp; This caused the test to  fail.&nbsp; 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>&nbsp;</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.&nbsp; The filter divides the narrow band (list of indices)
into a set  of groups, specifying one group for each thread.&nbsp; This is nice.&nbsp;
 However, the algorithm may set pixels in a 3x3x... neighborhood of a narrow
band  pixel.&nbsp; When running multiple threads this may be a problem.&nbsp; 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.&nbsp;  This can
cause a problem with memory integrity.&nbsp; 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>&nbsp;</div>
 
  <div><span class="865111014-20082003"><font size="2">&nbsp;&nbsp;&nbsp;&nbsp;&nbsp;&nbsp;&nbsp;&nbsp;&nbsp;&nbsp;&nbsp;  if(fabs(val1_new)
&lt;  fabs(outNeigIt.GetNext(n,1)))<br>
&nbsp;&nbsp;&nbsp;&nbsp;&nbsp;&nbsp;&nbsp;&nbsp;&nbsp;&nbsp;&nbsp;&nbsp;&nbsp;  outNeigIt.SetNext(n,1,static_cast&lt;PixelType&gt;(val1_new)
 );<br>
  </font></span></div>
 
  <div><span class="865111014-20082003"><font size="2">Consider&nbsp;a thread
that is  interrupted between the call to GetNext() and SetNext().&nbsp; 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.&nbsp; 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.&nbsp; Then
 when the first thread awakes, it has already satisfied its conditional so
it  writes it version of val1_new.&nbsp; 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>&nbsp;</div>
 
  <div><span class="865111014-20082003"><font size="2">&nbsp;</font></span></div>
 
  <div><span class="865111014-20082003"></span>&nbsp;</div>
 
  <div>&nbsp;</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
&amp;  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>&nbsp;</div>
</blockquote>
<br>
</body>
</html>

--------------070109060505020600090504--