[Insight-users] Behavior of FastMarchingImageFilter
Amardeep Singh
amar.singh at gmx.de
Fri May 15 05:27:04 EDT 2009
Dear Dan
Thank you very much for the bugfix.
It seems to work fine. :) I will let you know if I encounter any
problems in the future.
Best regards
Amardeep
Dan Mueller wrote:
> 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
>>
>>
More information about the Insight-users
mailing list