[Insight-users] fast marching behavior

Tue Jan 12 11:28:04 EST 2010

Reminds me of this problem Microsoft had:

http://blogs.msdn.com/oldnewthing/archive/2003/12/24/45779.aspx

I'd say the big difference is that these algorithms may be used in
mission-critical situations, ie, medical imaging where a bad
segmentation could be the difference between good and bad diagnoses.
So, rather than bury the change in an advanced compiler switch, it
might be a good idea to implement the two in parallel, and then have a
'fixed' version of the code (with a 'fixed' suffix?) and the older
version of the code remain as it is but with big warnings that it's
broken and people should change to the 'fixed' version.  That way, if
someone continues to use the wrong version, they at least know that
it's wrong.  On windows, you can throw in a #pragma warning to provide
a compile-time warning, not sure if that works for other compilers.

just my 2cts
Mark

On Tue, Jan 12, 2010 at 7:34 AM, Bill Lorensen <bill.lorensen at gmail.com> wrote:
> Luca,
>
> I think you're being even more conservative than I am.
>
> Bill
>
> On Tue, Jan 12, 2010 at 6:07 AM, Luca Antiga <luca.antiga at gmail.com> wrote:
>> Hi Bill,
>>  indeed this is a "logical" bug, a discrepancy with the way the filter
>> implements the original algorithm,
>> which leads to an output that is not what we expect by reading the
>> algorithm.
>> However, since this class is quite old and its expected use is widespread,
>> it's not unreasonable to
>> think that some of our customers might rely on its uncorrect behavior. In
>> many situations the bug
>> will not affect the solution to a stage where the results are incorrect, but
>> they will probably be
>> numerically different to some extent (for other applications, such as
>> Siqi's, things may be worse).
>> Still, such numerical differences in the solution would lead eventual tests
>> or, worse, assumptions based
>> on expected values in the solution, to fail, and in a very subtle way.
>> For this reason, I'm all for fixing the bug, but probably a CMake flag to
>> allow users to revert to the old
>> way (that might go away in a couple of releases) is in order, so to allow
>> customers a smooth transition.
>> Or I'm being too conservative here?
>>
>> Luca
>>
>>
>>
>> On Jan 12, 2010, at 2:26 AM, Bill Lorensen wrote:
>>
>>> Luca,
>>>
>>> If this is a bug it should be fixed. I don't think we need to maintain
>>> the bad behavior as long as the API remains the same. Unless I'm
>>> missing something subtle here.
>>>
>>> Perhaps others wish to comment.
>>>
>>> Bill
>>>
>>> On Fri, Jan 8, 2010 at 11:44 AM, Luca Antiga <luca.antiga at gmail.com>
>>> wrote:
>>>>
>>>> Hi Dan,
>>>>  not at all, I was just referring to addressing the latest comments by
>>>> Siqi.
>>>> I have the file on my desktop, waiting for me to handle it :-)
>>>>
>>>> Luca
>>>>
>>>>
>>>> On Jan 8, 2010, at 5:18 PM, Dan Mueller wrote:
>>>>
>>>>> Hi Luca,
>>>>>
>>>>> Is there anything wrong with the "FastMarching.patch" file attached to
>>>>> the bug entry?
>>>>>  http://public.kitware.com/Bug/view.php?id=8990
>>>>>
>>>>> 2010/1/8 Luca Antiga <luca.antiga at gmail.com>:
>>>>>>
>>>>>> 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
>>>>
>>>> _____________________________________
>>>> 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
>