[Insight-users] fast marching behavior

[siqi chen]

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/f56aa33e/attachment-0001.htm>