[Insight-users] fast marching behavior
Luca Antiga
luca.antiga at gmail.com
Fri Jan 8 10:53:30 EST 2010
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
>>> >
>>> > 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/~chens/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/17df39a4/attachment-0001.htm>
More information about the Insight-users
mailing list