<br>The reason I suggest to let the algorithm compute the initial trial is the following:<br><br>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.<br>
<br>To fix this problem, there are two options:<br>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.<br>
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.<br>
<br>Thanks<br>Siqi<br><br><div class="gmail_quote">On Thu, Jan 7, 2010 at 6:07 PM, Luca Antiga <span dir="ltr"><<a href="mailto:luca.antiga@gmail.com">luca.antiga@gmail.com</a>></span> wrote:<br><blockquote class="gmail_quote" style="border-left: 1px solid rgb(204, 204, 204); margin: 0pt 0pt 0pt 0.8ex; padding-left: 1ex;">
<div style="word-wrap: break-word;">Hi Siqi,<div><div> this sounds like a perfect contribution for the Insight Journal. </div><div>I suggest to write a new class for the second-order implementation.</div><div><br></div><div>
As for the minheap, I agree with you, it would save some memory.</div><div><br></div><div>As for setting trial points: in the current implementation, if you don't </div><div>set trial points you won't have anything to process in the TrialHeap.</div>
<div>I understand your point, but I am kind of in favor of the current way of </div><div>working. Setting alive points doesn't necessarily mean that you want </div><div>the solution to propagate automatically starting from their boundary </div>
<div>with far points. </div><div>This could indeed be an option, but it's not necessarily always practically </div><div>convenient. Is there any particular reason why you wouldn't want to set your</div><div>initial points as trial?</div>
<div><br></div><font color="#888888"><div><br></div><div>Luca</div></font><div><div></div><div class="h5"><div><br></div><div><br><div><div>On Jan 7, 2010, at 11:47 PM, siqi chen wrote:</div><br><blockquote type="cite"><br>
It definitely should be committed. I also have several suggestions about fastmarching implementation in ITK. <br><br>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.<br>
<br>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.<br>
<br>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:<br>
<a href="http://www2.imm.dtu.dk/pubdb/views/publication_details.php?id=841" target="_blank">http://www2.imm.dtu.dk/pubdb/views/publication_details.php?id=841</a> . <br><br>Thanks<br>Siqi<br><br><br><div class="gmail_quote">
On Thu, Jan 7, 2010 at 4:56 PM, Luca Antiga <span dir="ltr"><<a href="mailto:luca.antiga@gmail.com" target="_blank">luca.antiga@gmail.com</a>></span> wrote:<br> <blockquote class="gmail_quote" style="border-left: 1px solid rgb(204, 204, 204); margin: 0pt 0pt 0pt 0.8ex; padding-left: 1ex;">
<div style="word-wrap: break-word;"><div>Dan, Siqi,</div><div> the patch should be committed. I remember the original discussion, we never followed up</div> <div>committing the patch. It's time to do it.</div><div>However, since the fast marching method is probably used by a lot of people, how about</div>
<div>adding a cmake flag to revert to the "wrong" implementation to ensure backward compatibility?</div> <div><br></div><font color="#888888"><div><br></div><div>Luca</div></font><div><div></div><div><div><br></div>
<br><div><div>On Jan 7, 2010, at 10:42 PM, siqi chen wrote:</div><br><blockquote type="cite"><br>I think the patch correctly implements the FMM method described by Sethian. <br> <br>Thanks<br>Siqi<br><br><br><div class="gmail_quote">
On Thu, Jan 7, 2010 at 1:54 PM, Dan Mueller <span dir="ltr"><<a href="mailto:dan.muel@gmail.com" target="_blank">dan.muel@gmail.com</a>></span> wrote:<br> <blockquote class="gmail_quote" style="border-left: 1px solid rgb(204, 204, 204); margin: 0pt 0pt 0pt 0.8ex; padding-left: 1ex;">
Hi Siqi,<br> <br> I managed to find some time to take a look at the issue you reported.<br> <br> I think the following patch will help:<br> <a href="http://public.kitware.com/Bug/view.php?id=8990" target="_blank">http://public.kitware.com/Bug/view.php?id=8990</a><br>
<br> Here is the results I get with the patch applied:<br> <br> The distance from [50,50] to [49.7,49.8] is: 0.36<br> The distance from [49,49] to [49.7,49.8] is: 1.06<br> The distance from [49,50] to [49.7,49.8] is: 0.73<br>
The distance from [50,49] to [49.7,49.8] is: 0.85<br> The distance from [48,49] to [49.7,49.8] is: 2.06<br> The distance from [48,50] to [49.7,49.8] is: 1.73<br> The distance from [49,48] to [49.7,49.8] is: 2.06<br> The distance from [49,51] to [49.7,49.8] is: 1.73<br>
The distance from [50,48] to [49.7,49.8] is: 1.85<br> The distance from [50,51] to [49.7,49.8] is: 1.36<br> The distance from [51,49] to [49.7,49.8] is: 1.85<br> The distance from [51,50] to [49.7,49.8] is: 1.36<br> <br>
Is that what you expected?<br> <br> Hope this helps.<br> <br> Cheers, Dan<br> <br> 2010/1/5 siqi chen <<a href="mailto:siqichensc@gmail.com" target="_blank">siqichensc@gmail.com</a>>:<br> <div><div></div><div>><br>
> I read a few more papers and thought about the question number 2 again. I<br> > think there is misinterpretation on my part. The 4 neighbors of [49.7,49.8]<br> > along their exact distance value should be set as Alive Points instead of<br>
> Trial Points. Then the neighbors of the 4 neighbors (another layer around<br> > the 4 neighbors) are set as trial points, their initial tentative values are<br> > calculated using upwind difference method. When I check the result distance<br>
> map, a few things to notice:<br> > 1. The distance value of the 4 neighbors of [49.7, 49.8] are correct. This<br> > is obvious since I set them as Known instead of Trial.<br> > 2. The distance value of the neighbors of the 4 neighbors (the input trial<br>
> points) still changed. I think this is due to the inherent FMM accuracy and<br> > update method.<br> ><br> > You can try this with the following code.<br> > <a href="http://www.rpi.edu/%7Echens/download/main3.cpp" target="_blank">http://www.rpi.edu/~chens/download/main3.cpp</a><br>
><br> > Thanks<br> > Siqi<br> ><br> > On Mon, Jan 4, 2010 at 6:03 PM, siqi chen <<a href="mailto:siqichensc@gmail.com" target="_blank">siqichensc@gmail.com</a>> wrote:<br> >><br> >> To better illustrate my questions regarding fast marching, I put 2 example<br>
>> code in the attachment.<br> >><br> >> In main1.cpp, I simply compute a distance map to point [50,50]. As you can<br> >> see from the output, the distances from the 4 neighbors of [50,50] to<br>
>> [50,50] are correct, obviously the result is 1. However, the distance from<br> >> [51,51] to [50,50] is 1.707 instead of 1.414, which is obviously wrong. I<br> >> think this is due to the fast marching accuracy itself. If we switch to<br>
>> higher order FMM, the result should be improved.<br> >><br> >> In main2.cpp, I perturb the target point a little bit. Instead, I want to<br> >> compute the distance map to point [49.7,49.8]. From my point of<br>
>> understanding, I need to initialize the 4 neighbors of [49.7, 49.8] and put<br> >> them into the TrialPoints. As you can see, I compute the exact distance from<br> >> these 4 neighbors to [49.7,49.8] and put them into TrialPoints. However,<br>
>> when I go back and check the result distance map, some thing is different.<br> >> The distances from these 4 neighbors to [49.7,49.8] are changed. As you can<br> >> see, the distance from [50,50] to [49.7,49.8] remains correct. This is<br>
>> because this value is the smallest in the TrialPoints, therefore it is<br> >> pushed in to the AlivePoints heap first and the value is frozen since then.<br> >> I think there is something wrong here about whether to update trial points<br>
>> value or not. If this trial point is user specified, then the value should<br> >> not be updated. I noticed a related discussion a couple of months ago in the<br> >> mailing list,<br> >> <a href="http://www.itk.org/pipermail/insight-users/2009-May/030282.html" target="_blank">http://www.itk.org/pipermail/insight-users/2009-May/030282.html</a><br>
>><br> >> <a href="http://www.rpi.edu/%7Echens/download/main1.cpp" target="_blank">http://www.rpi.edu/~chens/download/main1.cpp</a><br> >> <a href="http://www.rpi.edu/%7Echens/download/main2.cpp" target="_blank">http://www.rpi.edu/~chens/download/main2.cpp</a><br>
>><br> >> Any input is appreciated.<br> >> Siqi<br> >><br> >><br> >> On Mon, Jan 4, 2010 at 3:32 PM, Dan Mueller <<a href="mailto:dan.muel@gmail.com" target="_blank">dan.muel@gmail.com</a>> wrote:<br>
>>><br> >>> Hi Siqi,<br> >>><br> >>> Indeed I am familiar with Fast Marching. I saw your question to the<br> >>> mailing list, but did not respond because I have not experienced what<br>
>>> you describe: when I set the trial point value, that is the value in<br> >>> the arrival function.<br> >>><br> >>> Perhaps you could post to the mailing list a minimal example<br> >>> (code+cmake+data) demonstrating your issue. That would make it really<br>
>>> easy for me to help you!<br> >>><br> >>> Cheers, Dan<br> >>><br> >>> 2010/1/4 siqi chen <<a href="mailto:siqichensc@gmail.com" target="_blank">siqichensc@gmail.com</a>>:<br>
>>> > Hi, Dan,<br> >>> ><br> >>> > Sorry to bother you. From the ITK mailing list, I noticed you reported<br> >>> > a bug<br> >>> > about FastMarchingImageFilter couple of months ago. So I guess you are<br>
>>> > a<br> >>> > fast marching expert : )<br> >>> ><br> >>> > I am trying to use FastMarchingImageFilter to calculate a distance map<br> >>> > to a<br> >>> > set of points which have non-integer coordinates and I want the result<br>
>>> > to be<br> >>> > as accurate as possible. Here is what I did, but the result is not very<br> >>> > accurate.<br> >>> ><br> >>> > First I find the integer points which are the neighbors of the target<br>
>>> > points<br> >>> > and set these integer points as trial points. Then I use some<br> >>> > interpolation<br> >>> > method to initialize the distance from these trial points to the target<br>
>>> > points, which are assumed to be "exactly correct". The TrialPoints in<br> >>> > the<br> >>> > FastMarchingImageFilter is defined as this set of trial points and<br> >>> > their<br>
>>> > corresponding distances to the target points. The AlivePoints is empty.<br> >>> > When<br> >>> > I check the result distance map, I find that the distance value of<br> >>> > these<br>
>>> > trial points are changed, they are no longer what their initial states<br> >>> > are.<br> >>> > Therefore, the iso curve deviate the original input a little bit. I am<br> >>> > quite<br>
>>> > confusing about this result.<br> >>> ><br> >>> > I noticed you mentioned on the mailing list about neighbor update, that<br> >>> > is<br> >>> > to distinguish between user-specified trial points and<br>
>>> > algorithm-generated<br> >>> > trial points. here is the discussion,<br> >>> > <a href="http://www.itk.org/pipermail/insight-users/2009-May/030282.html" target="_blank">http://www.itk.org/pipermail/insight-users/2009-May/030282.html</a> . I<br>
>>> > wonder<br> >>> > if you have any suggestions about my problem.<br> >>> ><br> >>> > Any input is appreciated.<br> >>> ><br> >>> > Thanks<br> >>> > Siqi Chen<br>
>>> ><br> >><br> ><br> ><br> </div></div>> _____________________________________<br> > Powered by <a href="http://www.kitware.com" target="_blank">www.kitware.com</a><br> ><br> > Visit other Kitware open-source projects at<br>
> <a href="http://www.kitware.com/opensource/opensource.html" target="_blank">http://www.kitware.com/opensource/opensource.html</a><br> ><br> > Kitware offers ITK Training Courses, for more information visit:<br>
> <a href="http://www.kitware.com/products/protraining.html" target="_blank">http://www.kitware.com/products/protraining.html</a><br> ><br> > Please keep messages on-topic and check the ITK FAQ at:<br> > <a href="http://www.itk.org/Wiki/ITK_FAQ" target="_blank">http://www.itk.org/Wiki/ITK_FAQ</a><br>
><br> > Follow this link to subscribe/unsubscribe:<br> > <a href="http://www.itk.org/mailman/listinfo/insight-users" target="_blank">http://www.itk.org/mailman/listinfo/insight-users</a><br> ><br> ><br> </blockquote>
</div><br> _____________________________________<br>Powered by <a href="http://www.kitware.com" target="_blank">www.kitware.com</a><br><br>Visit other Kitware open-source projects at<br><a href="http://www.kitware.com/opensource/opensource.html" target="_blank">http://www.kitware.com/opensource/opensource.html</a><br>
<br>Kitware offers ITK Training Courses, for more information visit:<br><a href="http://www.kitware.com/products/protraining.html" target="_blank">http://www.kitware.com/products/protraining.html</a><br><br>Please keep messages on-topic and check the ITK FAQ at:<br>
<a href="http://www.itk.org/Wiki/ITK_FAQ" target="_blank">http://www.itk.org/Wiki/ITK_FAQ</a><br><br>Follow this link to subscribe/unsubscribe:<br><a href="http://www.itk.org/mailman/listinfo/insight-users" target="_blank">http://www.itk.org/mailman/listinfo/insight-users</a><br>
</blockquote></div><br></div></div></div></blockquote></div><br></blockquote></div><br></div></div></div></div></div></blockquote></div><br>