[Insight-users] Behavior of FastMarchingImageFilter
Dan Mueller
dan.muel at gmail.com
Thu May 14 14:17:13 EDT 2009
Hi Amardeep,
Find attached a patch which should solve the issue you reported.
Please let us know if this works for you. If all is good, I will try
to get the fix (and some tests :P) into the archive this weekend.
Regards, Dan
Index: itkFastMarchingImageFilter.h
===================================================================
RCS file: /cvsroot/Insight/Insight/Code/Algorithms/itkFastMarchingImageFilter.h,v
retrieving revision 1.40
diff -u -r1.40 itkFastMarchingImageFilter.h
--- itkFastMarchingImageFilter.h 23 Apr 2009 03:53:35 -0000 1.40
+++ itkFastMarchingImageFilter.h 13 May 2009 17:29:12 -0000
@@ -162,7 +162,7 @@
* away points; TrialPoints represent points within a narrowband of the
* propagating front; and AlivePoints represent points which have already
* been processed. */
- enum LabelType { FarPoint, AlivePoint, TrialPoint };
+ enum LabelType { FarPoint, AlivePoint, TrialPoint, InitialTrialPoint };
/** LabelImage typedef support. */
typedef Image<unsigned char, itkGetStaticConstMacro(SetDimension)>
LabelImageType;
Index: itkFastMarchingImageFilter.txx
===================================================================
RCS file: /cvsroot/Insight/Insight/Code/Algorithms/itkFastMarchingImageFilter.txx,v
retrieving revision 1.52
diff -u -r1.52 itkFastMarchingImageFilter.txx
--- itkFastMarchingImageFilter.txx 21 Dec 2008 19:13:11 -0000 1.52
+++ itkFastMarchingImageFilter.txx 13 May 2009 17:28:42 -0000
@@ -239,7 +239,7 @@
}
// make this a trial point
- m_LabelImage->SetPixel( node.GetIndex(), TrialPoint );
+ m_LabelImage->SetPixel( node.GetIndex(), InitialTrialPoint );
outputPixel = node.GetValue();
output->SetPixel( node.GetIndex(), outputPixel );
@@ -289,7 +289,8 @@
}
// is this node already alive ?
- if ( m_LabelImage->GetPixel( node.GetIndex() ) != TrialPoint )
+ double label = m_LabelImage->GetPixel( node.GetIndex() );
+ if ( label == AlivePoint )
{
continue;
}
@@ -349,7 +350,8 @@
{
neighIndex[j] = index[j] - 1;
}
- if ( m_LabelImage->GetPixel( neighIndex ) != AlivePoint )
+ double label = m_LabelImage->GetPixel( neighIndex );
+ if ( label != AlivePoint && label != InitialTrialPoint )
{
this->UpdateValue( neighIndex, speedImage, output );
}
@@ -359,7 +361,8 @@
{
neighIndex[j] = index[j] + 1;
}
- if ( m_LabelImage->GetPixel( neighIndex ) != AlivePoint )
+ label = m_LabelImage->GetPixel( neighIndex );
+ if ( label != AlivePoint && label != InitialTrialPoint )
{
this->UpdateValue( neighIndex, speedImage, output );
}
2009/5/12 Luca Antiga <luca.antiga at gmail.com>:
> Hi Dan,
> I'll be happy to interact with you on the fixes next week, so they
> can be committed after the cvs goes back to unfrozen.
> Thanks in the meantime and keep in touch
>
> Luca
>
> 2009/5/12, Dan Mueller <dan.muel at gmail.com>:
>> Hi Amardeep,
>>
>> The initial patch I proposed missed a use case, so: no, the results
>> you experienced are not normal. The solution Luca proposed (adding an
>> InitialTrial type) should correctly solve your issue.
>>
>> I'm snowed under at work at the moment, so I won't get around to
>> submitting the fix until the weekend...
>>
>> Regards, Dan
>>
>> 2009/5/12 Amardeep Singh <amar.singh at gmx.de>:
>>> Hello
>>>
>>> I would just like to let you know what I have experienced with the bugfix
>>> as
>>> it is
>>> now.
>>> The "jaggednes" of the result has definitely disappeared. However, I am
>>> confused
>>> by the fact that the isolines of the fast marching output a little further
>>> away from the object look like edges of a square
>>> rotated by 45 degrees (see the snapshot at the isoline 12 and the attached
>>> output image; I have applied the filter to the
>>> binary circle again).
>>> Is this normal? I understand that the current bugfix is a preliminary one.
>>>
>>> I am looking forward to the final bugfix and thanks a lot for correcting
>>> this bug.
>>>
>>> Best regards
>>> Amar
>>>
>>>
>>> Dan Mueller wrote:
>>>>
>>>> Hi Luca,
>>>>
>>>> Sure -- adding the target stuff to the base class at the same time sounds
>>>> great.
>>>>
>>>> Do you want me to do it? I probably won't make the Friday 15th May
>>>> deadline for 3.14 though...
>>>>
>>>> Cheers, Dan
>>>>
>>>> 2009/5/11 Luca Antiga <luca.antiga at gmail.com>:
>>>>
>>>>>
>>>>> Hi Dan,
>>>>> yep, that would work.
>>>>> The user-supplied Trial points should be tagged as InitialTrial in the
>>>>> LabelImage,
>>>>> put in the TrialHeap but not be updated as neighbors. And taking care of
>>>>> condition
>>>>> if ( m_LabelImage->GetPixel( node.GetIndex() ) != TrialPoint )
>>>>> at line 292 in itkFastMarchingImageFilter.txx.
>>>>> Actually, since we are at it, it may also be good to move the Target
>>>>> business from
>>>>> FastMarchingUpwindGradientImageFilter to the base class (see your Apr 14
>>>>> message "Trying to monitor progression of fast march operation").
>>>>> Actually, we may also improve it by not looping around the target point
>>>>> list
>>>>> in
>>>>> the UpdateNeighbors method, but defining an extra label in the enum
>>>>> (Target)
>>>>> and storing the Target label in the LabelImage.
>>>>> I didn't do it back then because I just wanted to write an extra class
>>>>> without
>>>>> touching the original filter.
>>>>> What do you think?
>>>>>
>>>>> Luca
>>>>>
>>>>>
>>>>> On May 9, 2009, at 8:19 PM, Dan Mueller wrote:
>>>>>
>>>>> shouldn't you check if the update value computed for the trial point
>>>>>
>>>>> is smaller than the already existing value for that trial point and if
>>>>>
>>>>> so update it regardless?
>>>>>
>>>>> Ah yes, I see the issue with the fix I proposed -- once a trial point
>>>>> is given a solution value, it will never be improved upon.
>>>>>
>>>>> But what if the trial point is user specified? Shouldn't the value
>>>>> given by the user be used, even it is larger than the calculated one?
>>>>> Otherwise, why does the user even specify the node value?
>>>>>
>>>>> So I guess we need some way to distinguish between user specified and
>>>>> computed trial points... If the neighbour is a user specified trial
>>>>> point or alive, then don't update it...?
>>>>>
>>>>> Regards, Dan
>>>>>
>>>>> 2009/5/9 Luca Antiga <luca.antiga at gmail.com>:
>>>>>
>>>>> Hi Dan,
>>>>>
>>>>> shouldn't you check if the update value computed for the trial point
>>>>>
>>>>> is smaller than the already existing value for that trial point and if
>>>>>
>>>>> so update it regardless?
>>>>>
>>>>> Luca
>>>>>
>>>>> 2009/5/9, Dan Mueller <dan.muel at gmail.com>:
>>>>>
>>>>> Hi Amar,
>>>>>
>>>>> Congratulations! It seems you managed to find a bug in the
>>>>>
>>>>> FastMarchingImageFilter :P
>>>>>
>>>>> In itkFastMarchingImageFilter.txx the UpdateNeighbors method checks if
>>>>>
>>>>> the candidate neighbour is an alive point. If so, then it is skipped
>>>>>
>>>>> from further processing. If not, then it is processed (ie. output
>>>>>
>>>>> value updated). This method should also test if the candidate
>>>>>
>>>>> neighbour is a *trial* point. If it is a trial point, then the output
>>>>>
>>>>> value has already been set, and does not need to be updated. Doing so
>>>>>
>>>>> will actually produce incorrect values, as you have encountered with
>>>>>
>>>>> your "jagged" edges.
>>>>>
>>>>> I have submitted a bug for this issue:
>>>>>
>>>>> http://public.kitware.com/Bug/view.php?id=8990
>>>>>
>>>>> The bug report includes a patch, which I have also attached to this
>>>>> mail.
>>>>>
>>>>> Surprisingly this issue has been around since version 1.16 of the file
>>>>>
>>>>> (nearly 8 years)...
>>>>>
>>>>> Please let us know if this solves your issue.
>>>>>
>>>>> Regards, Dan
>>>>>
>>>>>
>>>>> 2009/5/8 Dan Mueller <dan.muel at gmail.com>:
>>>>>
>>>>> Hi Amar,
>>>>>
>>>>> I can reproduce your results. It is very strange; I can't see what is
>>>>>
>>>>> wrong with the code. I will look into it further on the weekend.
>>>>>
>>>>> Cheers, Dan
>>>>>
>>>>> 2009/5/8 Amardeep Singh <amar.singh at gmx.de>:
>>>>>
>>>>> Hello everybody:
>>>>>
>>>>> Dear Dan:
>>>>>
>>>>> Thank you very much for your help!
>>>>>
>>>>> I have attached some sample code and a binary image of a circle for
>>>>>
>>>>> testing.
>>>>>
>>>>> The program can be called like the following:
>>>>>
>>>>> fast_marching_itk_mailing_list ./circle.nii.gz
>>>>>
>>>>> ./fast_marching_output.nii.gz
>>>>>
>>>>> 200
>>>>>
>>>>> The last number is the stopping value which I usually set to large
>>>>>
>>>>> values. I
>>>>>
>>>>> use
>>>>>
>>>>> the fast marching filter in this example with a constant speed of 1 so
>>>>>
>>>>> that
>>>>>
>>>>> the result should be some kind of a distance map.
>>>>>
>>>>> Certainly, one important point is my calculation of the alive and trial
>>>>>
>>>>> points. All the points of the circle that have a background voxel within
>>>>>
>>>>> their
>>>>>
>>>>> 6-neighborhood are considered trial points (so the last layer around the
>>>>>
>>>>> circle becomes a layer of trial points). All other points of the circle
>>>>>
>>>>> are
>>>>>
>>>>> alive points.
>>>>>
>>>>> I have attached a screenshot of the 0 isoline of the result (the
>>>>> complete
>>>>>
>>>>> output image was too large for the mailing list, unfortunately). When
>>>>> you
>>>>>
>>>>> look at it, you can see the jagged edge again. But still: Shouldn't it
>>>>>
>>>>> propagate evenly from the trial points?
>>>>>
>>>>> The intialization (0 isoline of the result) looks jagged already.
>>>>>
>>>>> Concerning my previous snapshot:
>>>>>
>>>>> The output image was of type float. Just for the snapshopt, I mapped all
>>>>>
>>>>> values > 0 to 1 in order to visualize the jagged edge.
>>>>>
>>>>> Best regards
>>>>>
>>>>> Amar
>>>>>
>>>>>
>>>>> Dan Mueller wrote:
>>>>>
>>>>> Hi Amar,
>>>>>
>>>>> Perhaps you can post a minimal example (code, input image, output
>>>>>
>>>>> image) demonstrating the issue.
>>>>>
>>>>> I have not experienced this issue myself, but the
>>>>>
>>>>> FastMarchingImageFilter does propagate the front using
>>>>>
>>>>> 4-/6-connectivity, so after a small timestep it may be possible to
>>>>>
>>>>> experience the "jagged" edge you are seeing (though not the checked
>>>>>
>>>>> pattern inside the hole). How many iterations are you letting the
>>>>>
>>>>> filter run for? (ie. what is the stopping value?)
>>>>>
>>>>> Also, the FastMarchingImageFilter requires a real (float, double)
>>>>>
>>>>> output image. What is the pixel type of your output image?
>>>>>
>>>>> Cheers, Dan
>>>>>
>>>>> 2009/5/7 Amardeep Singh <amar.singh at gmx.de>:
>>>>>
>>>>>
>>>>> Hello everyone
>>>>>
>>>>> I am wondering about the behavior of the FastMarchingImageFilter.
>>>>>
>>>>> I want to let the algorithm march from a certain surface which is
>>>>>
>>>>> represented by a binary image. I
>>>>>
>>>>> pass the points of the object as alive points to the marcher. The layer
>>>>>
>>>>> outside of the object (in my case: a brain)
>>>>>
>>>>> is passed as trial points. When I look at the output of the fast
>>>>>
>>>>> marching filter, there is something strange -> see attached screenshot.
>>>>>
>>>>> The area within the yellow line is my initial object. So they are all
>>>>>
>>>>> alive points whose value is 0 in the output. But when you look at the
>>>>>
>>>>> layer
>>>>>
>>>>> outside of the inner object, you see that it is jagged, black (value 0)
>>>>>
>>>>> and white (value != 0) voxels take turns. Why is this the case?
>>>>>
>>>>> What do I need to do, so that the marcher marches equally away from the
>>>>>
>>>>> alive points. The speed image values on the layer around the alive
>>>>>
>>>>> points do have the same value, so this is not the reason.
>>>>>
>>>>> Thanks a lot in advance!
>>>>>
>>>>> Best Regards
>>>>>
>>>>> Amar
>>>>>
>>>>
>>>>
>>>
>> _____________________________________
>> Powered by www.kitware.com
>>
>> Visit other Kitware open-source projects at
>> http://www.kitware.com/opensource/opensource.html
>>
>> Please keep messages on-topic and check the ITK FAQ at:
>> http://www.itk.org/Wiki/ITK_FAQ
>>
>> Follow this link to subscribe/unsubscribe:
>> http://www.itk.org/mailman/listinfo/insight-users
>>
>
> --
> Inviato dal mio dispositivo mobile
>
> Luca Antiga, PhD
> Biomedical Technologies Laboratory
> Biomedical Engineering Department,
> Mario Negri Institute
> mail: Villa Camozzi, 24020, Ranica (BG), Italy
> phone: +39 035 4535-381
> email: antiga at marionegri.it
> web: http://villacamozzi.marionegri.it/~luca
>
-------------- next part --------------
A non-text attachment was scrubbed...
Name: FastMarching.patch
Type: application/octet-stream
Size: 2515 bytes
Desc: not available
URL: <http://www.itk.org/pipermail/insight-users/attachments/20090514/b76759df/attachment-0001.obj>
More information about the Insight-users
mailing list