[Insight-users] fast marching behavior

[Dan Mueller]

Hi Siqi,

Responses inline below.

2010/1/7 siqi chen <siqichensc at gmail.com>:
>
> 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.

I believe with the patch applied, you can again just add trial points;
no need to add alive points for your situation.

> 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 .

An interesting paper. We look forward to your Insight Journal
submission with the higher-order implementation :P

> 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