[Insight-users] fast marching behavior

siqi chen siqichensc at gmail.com
Fri Jan 8 11:00:28 EST 2010 

I agree with this solution and look forward to the patch. Let me know if I
can do some contribution.

Thanks
Siqi


On Fri, Jan 8, 2010 at 10:53 AM, Luca Antiga <luca.antiga at gmail.com> wrote:

> Hi Siqi,
>  I understand your points. We have to find a solution that guarantees
> backwards compatibility, we cannot switch to setting Alive points
> altogether.
> My proposal is to use the distinction between Trial points and InitialTrial
> points in Dan's patch, and to add a flag that will decide whether
> InitialTrial points should be updated or not. This will solve the problem
> maintaining backwards compatibility. With no updating allowed, InitialTrial
> points will behave just like you wish, and we won't have to write extra code
> for computing the initial layer around Alive points.
> Let me know if you agree with this.
>
> Dan, I appreciate the lack of time, and I must say we're pretty much on the
> same boat. However, this issue carries some weight. As soon as we devise an
> acceptable solution, I'm willing to put in some time to commit the patch and
> take care of the rest.
>
> Best regards
>
> Luca
>
>
>
>
> On Jan 8, 2010, at 4:06 PM, siqi chen wrote:
>
>
> The reason I suggest to let the algorithm compute the initial trial is the
> following:
>
> If we want to compute a distance map to a perfect circle. The first thing
> we do is to discretize the circle, of course, most of the points will be
> non-integer coordinates. To initialize the algorithm, we need to find the 4
> neighbor lattice of the discretized circle points and calculate the exact
> distance from these lattice to the circle. In my humble opinion, in the
> current ITK implementation, the algorithm accepts these lattice points as
> "TrialPoints". As I illustrated before, since the FMM method will update the
> values of trial points if the new value is smaller, it is very likely that
> after the FMM sweeping, the value of the these lattice (the immediate 4
> neighbors of the circle) will change. This will introduce large errors if we
> compare the zero iso curve of the result distance to the initial circle.
>
> To fix this problem, there are two options:
> 1. Distinguish between the user input trial and the algorithm generated
> trial as you guys mentioned before. If it's a user input trial, then don't
> update the value at all.
> 2. Set these immediate 4 neighbors as KnownPoints, and let the user compute
> another layer outside these knownpoints and set them as trial points. This
> is quite tedious for the user, since they need to figure out the "upwind"
> themselves. Again, according to the FMM algorithm itself, this initial trial
> computation should be done by the algorithm, not the user. That's why I
> think we should let the algorithm to compute the initial trial.
>
> Thanks
> Siqi
>
> On Thu, Jan 7, 2010 at 6:07 PM, Luca Antiga <luca.antiga at gmail.com> wrote:
>
>> Hi Siqi,
>>  this sounds like a perfect contribution for the Insight Journal.
>> I suggest to write a new class for the second-order implementation.
>>
>> As for the minheap, I agree with you, it would save some memory.
>>
>> As for setting trial points: in the current implementation, if you don't
>> set trial points you won't have anything to process in the TrialHeap.
>> I understand your point, but I am kind of in favor of the current way of
>> working. Setting alive points doesn't necessarily mean that you want
>> the solution to propagate automatically starting from their boundary
>> with far points.
>> This could indeed be an option, but it's not
>> necessarily always practically
>> convenient. Is there any particular reason why you wouldn't want to set
>> your
>> initial points as trial?
>>
>>
>> Luca
>>
>>
>> On Jan 7, 2010, at 11:47 PM, siqi chen wrote:
>>
>>
>> It definitely should be committed. I also have several suggestions about
>> fastmarching implementation in ITK.
>>
>> 1. I think it's quite misleading to say that "In order for the filter to
>> evolve, at least some trial points must be specified." I think it should be
>> "at least some alive points must be specified". The initial trial points
>> should be calculated by the algorithm itself, not by the user. For example,
>> when I start with FastMarchingImageFilter, I want to try a simple example of
>> calculating a distance map to [50,50], Mr. Kevin Hobbs told me the right way
>> is to set [50,50] as trial point. However, I think it is not right
>> description. The right way is to set [50,50] as alive points, and the
>> algorithm computes the initial trial points and start the iteration.
>>
>> 2. Another thing is about heap implementation. Currently,
>> std::priority_queue is used, and there is some issue as Doxygen described.
>> I think we can use a minheap as the basic data structure. I wrote a VTK fast
>> maching filter yesterday and I used this min heap data structure.
>>
>> 3. Possible improvements. Current implementation is based on Prof.
>> Sethian' s work, which is actually a first order approximation of distance
>> map. This implementation will introduce large errors in the diagonal
>> direction. With a little bit change of code, we can implement a higher order
>> accuracy FMM method. The detail can be found here:
>> http://www2.imm.dtu.dk/pubdb/views/publication_details.php?id=841 .
>>
>> Thanks
>> Siqi
>>
>>
>> On Thu, Jan 7, 2010 at 4:56 PM, Luca Antiga <luca.antiga at gmail.com>wrote:
>>
>>> Dan, Siqi,
>>>  the patch should be committed. I remember the original discussion, we
>>> never followed up
>>> committing the patch. It's time to do it.
>>> However, since the fast marching method is probably used by a lot of
>>> people, how about
>>> adding a cmake flag to revert to the "wrong" implementation to ensure
>>> backward compatibility?
>>>
>>>
>>> Luca
>>>
>>>
>>> On Jan 7, 2010, at 10:42 PM, siqi chen wrote:
>>>
>>>
>>> I think the patch correctly implements the FMM method described by
>>> Sethian.
>>>
>>> Thanks
>>> Siqi
>>>
>>>
>>> On Thu, Jan 7, 2010 at 1:54 PM, Dan Mueller <dan.muel at gmail.com> wrote:
>>>
>>>> Hi Siqi,
>>>>
>>>> I managed to find some time to take a look at the issue you reported.
>>>>
>>>> I think the following patch will help:
>>>>    http://public.kitware.com/Bug/view.php?id=8990
>>>>
>>>> Here is the results I get with the patch applied:
>>>>
>>>> The distance from [50,50] to [49.7,49.8] is: 0.36
>>>> The distance from [49,49] to [49.7,49.8] is: 1.06
>>>> The distance from [49,50] to [49.7,49.8] is: 0.73
>>>> The distance from [50,49] to [49.7,49.8] is: 0.85
>>>> The distance from [48,49] to [49.7,49.8] is: 2.06
>>>> The distance from [48,50] to [49.7,49.8] is: 1.73
>>>> The distance from [49,48] to [49.7,49.8] is: 2.06
>>>> The distance from [49,51] to [49.7,49.8] is: 1.73
>>>> The distance from [50,48] to [49.7,49.8] is: 1.85
>>>> The distance from [50,51] to [49.7,49.8] is: 1.36
>>>> The distance from [51,49] to [49.7,49.8] is: 1.85
>>>> The distance from [51,50] to [49.7,49.8] is: 1.36
>>>>
>>>> Is that what you expected?
>>>>
>>>> Hope this helps.
>>>>
>>>> Cheers, Dan
>>>>
>>>> 2010/1/5 siqi chen <siqichensc at gmail.com>:
>>>> >
>>>> > I read a few more papers and  thought about the question number 2
>>>> again. I
>>>> > think there is misinterpretation on my part. The 4 neighbors of
>>>> [49.7,49.8]
>>>> > along their exact distance value should be set as Alive Points instead
>>>> of
>>>> > Trial Points. Then the neighbors of the 4 neighbors (another layer
>>>> around
>>>> > the 4 neighbors) are set as trial points, their initial tentative
>>>> values are
>>>> > calculated using upwind difference method. When I check the result
>>>> distance
>>>> > map, a few things to notice:
>>>> > 1. The distance value of the 4 neighbors of [49.7, 49.8] are correct.
>>>> This
>>>> > is obvious since I set them as Known instead of Trial.
>>>> > 2. The distance value of the neighbors of the 4 neighbors (the input
>>>> trial
>>>> > points) still changed. I think this is due to the inherent FMM
>>>> accuracy and
>>>> > update method.
>>>> >
>>>> > You can try this with the following code.
>>>> > http://www.rpi.edu/~chens/download/main3.cpp<http://www.rpi.edu/%7Echens/download/main3.cpp>
>>>> >
>>>> > Thanks
>>>> > Siqi
>>>> >
>>>> > On Mon, Jan 4, 2010 at 6:03 PM, siqi chen <siqichensc at gmail.com>
>>>> wrote:
>>>> >>
>>>> >> To better illustrate my questions regarding fast marching, I put 2
>>>> example
>>>> >> code in the attachment.
>>>> >>
>>>> >> In main1.cpp, I simply compute a distance map to point [50,50]. As
>>>> you can
>>>> >> see from the output, the distances from the 4 neighbors of [50,50] to
>>>> >> [50,50] are correct, obviously the result is 1. However, the distance
>>>> from
>>>> >> [51,51] to [50,50] is 1.707 instead of 1.414, which is obviously
>>>> wrong. I
>>>> >> think this is due to the fast marching accuracy itself. If we switch
>>>> to
>>>> >> higher order FMM, the result should be improved.
>>>> >>
>>>> >> In main2.cpp, I perturb the target point a little bit. Instead, I
>>>> want to
>>>> >> compute the distance map to point [49.7,49.8]. From my point of
>>>> >> understanding, I need to initialize the 4 neighbors of [49.7, 49.8]
>>>> and put
>>>> >> them into the TrialPoints. As you can see, I compute the exact
>>>> distance from
>>>> >> these 4 neighbors to [49.7,49.8] and put them into TrialPoints.
>>>> However,
>>>> >> when I go back and check the result distance map, some thing is
>>>> different.
>>>> >> The distances from these 4 neighbors to [49.7,49.8] are changed. As
>>>> you can
>>>> >> see, the distance from [50,50] to [49.7,49.8] remains correct. This
>>>> is
>>>> >> because this value is the smallest in the TrialPoints, therefore it
>>>> is
>>>> >> pushed in to the AlivePoints heap first and the value is frozen since
>>>> then.
>>>> >> I think there is something wrong here about whether to update trial
>>>> points
>>>> >> value or not. If this trial point is user specified, then the value
>>>> should
>>>> >> not be updated. I noticed a related discussion a couple of months ago
>>>> in the
>>>> >> mailing list,
>>>> >> http://www.itk.org/pipermail/insight-users/2009-May/030282.html
>>>> >>
>>>> >> http://www.rpi.edu/~chens/download/main1.cpp<http://www.rpi.edu/%7Echens/download/main1.cpp>
>>>> >> http://www.rpi.edu/~chens/download/main2.cpp<http://www.rpi.edu/%7Echens/download/main2.cpp>
>>>> >>
>>>> >> Any input is appreciated.
>>>> >> Siqi
>>>> >>
>>>> >>
>>>> >> On Mon, Jan 4, 2010 at 3:32 PM, Dan Mueller <dan.muel at gmail.com>
>>>> wrote:
>>>> >>>
>>>> >>> Hi Siqi,
>>>> >>>
>>>> >>> Indeed I am familiar with Fast Marching. I saw your question to the
>>>> >>> mailing list, but did not respond because I have not experienced
>>>> what
>>>> >>> you describe: when I set the trial point value, that is the value in
>>>> >>> the arrival function.
>>>> >>>
>>>> >>> Perhaps you could post to the mailing list a minimal example
>>>> >>> (code+cmake+data) demonstrating your issue. That would make it
>>>> really
>>>> >>> easy for me to help you!
>>>> >>>
>>>> >>> Cheers, Dan
>>>> >>>
>>>> >>> 2010/1/4 siqi chen <siqichensc at gmail.com>:
>>>> >>> > Hi, Dan,
>>>> >>> >
>>>> >>> > Sorry to bother you. From the ITK mailing list, I noticed you
>>>> reported
>>>> >>> > a bug
>>>> >>> > about FastMarchingImageFilter couple of months ago. So I guess you
>>>> are
>>>> >>> > a
>>>> >>> > fast marching expert : )
>>>> >>> >
>>>> >>> > I am trying to use FastMarchingImageFilter to calculate a distance
>>>> map
>>>> >>> > to a
>>>> >>> > set of points which have non-integer coordinates and I want the
>>>> result
>>>> >>> > to be
>>>> >>> > as accurate as possible. Here is what I did, but the result is not
>>>> very
>>>> >>> > accurate.
>>>> >>> >
>>>> >>> > First I find the integer points which are the neighbors of the
>>>> target
>>>> >>> > points
>>>> >>> > and set these integer points as trial points. Then I use some
>>>> >>> > interpolation
>>>> >>> > method to initialize the distance from these trial points to the
>>>> target
>>>> >>> > points, which are assumed to be "exactly correct". The TrialPoints
>>>> in
>>>> >>> > the
>>>> >>> > FastMarchingImageFilter is defined as this set of trial points and
>>>> >>> > their
>>>> >>> > corresponding distances to the target points. The AlivePoints is
>>>> empty.
>>>> >>> > When
>>>> >>> > I check the result distance map, I find that the distance value of
>>>> >>> > these
>>>> >>> > trial points are changed, they are no longer what their initial
>>>> states
>>>> >>> > are.
>>>> >>> > Therefore, the iso curve deviate the original input a little bit.
>>>> I am
>>>> >>> > quite
>>>> >>> > confusing about this result.
>>>> >>> >
>>>> >>> > I noticed you mentioned on the mailing list about neighbor update,
>>>> that
>>>> >>> > is
>>>> >>> > to distinguish between user-specified trial points and
>>>> >>> > algorithm-generated
>>>> >>> > trial points.  here is the discussion,
>>>> >>> > http://www.itk.org/pipermail/insight-users/2009-May/030282.html.  I
>>>> >>> > wonder
>>>> >>> > if you have any suggestions about my problem.
>>>> >>> >
>>>> >>> > Any input is appreciated.
>>>> >>> >
>>>> >>> > Thanks
>>>> >>> > Siqi Chen
>>>> >>> >
>>>> >>
>>>> >
>>>> >
>>>> > _____________________________________
>>>> > Powered by www.kitware.com
>>>> >
>>>> > Visit other Kitware open-source projects at
>>>> > http://www.kitware.com/opensource/opensource.html
>>>> >
>>>> > Kitware offers ITK Training Courses, for more information visit:
>>>> > http://www.kitware.com/products/protraining.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
>>>> >
>>>> >
>>>>
>>>
>>> _____________________________________
>>> Powered by www.kitware.com
>>>
>>> Visit other Kitware open-source projects at
>>> http://www.kitware.com/opensource/opensource.html
>>>
>>> Kitware offers ITK Training Courses, for more information visit:
>>> http://www.kitware.com/products/protraining.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
>>>
>>>
>>>
>>
>>
>
>
-------------- next part --------------
An HTML attachment was scrubbed...
URL: <http://www.itk.org/pipermail/insight-users/attachments/20100108/ac02a81a/attachment-0001.htm>

More information about the Insight-users mailing list