[Insight-users] fast marching behavior

Wed Jan 13 11:27:00 EST 2010

Hi Bill,

The reason the existing test misses the bug is that it does not add
two seed points side-by-side...

HTH

Cheers, Dan

2010/1/13 Bill Lorensen <bill.lorensen at gmail.com>:
> Maybe our testing of fast marching is not complete. I'll snoop around.
>
> Bill
>
> On Wed, Jan 13, 2010 at 10:12 AM, Luca Antiga <luca.antiga at gmail.com> wrote:
>> Sadly enough, yes :-)
>> Should we open a bug?
>>
>> Luca
>>
>>
>> On Jan 13, 2010, at 4:10 PM, Bill Lorensen wrote:
>>
>>> Luca,
>>>
>>> Are you saying that all tests passed with modifying any regression images?
>>>
>>> Bill
>>>
>>> On Wed, Jan 13, 2010 at 9:35 AM, Luca Antiga <luca.antiga at gmail.com>
>>> wrote:
>>>>
>>>> Ok, the fix is in. In the end I did stick to the original idea of adding
>>>> an
>>>> advanced CMake
>>>> variable (ITK_USE_DEPRECATED_FAST_MARCHING), OFF by default. If this is
>>>> not
>>>> ok I'll remove it, just let me know.
>>>> All tests are passing on my machine, I'll monitor the dashboard.
>>>> Please let me know about any issue.
>>>> Best regards
>>>>
>>>> Luca
>>>>
>>>>
>>>> On Jan 13, 2010, at 10:48 AM, Luca Antiga wrote:
>>>>
>>>>> Hey Bill,
>>>>>>>>
>>>>>>>> Luca,
>>>>>>>>
>>>>>>>> I think you're being even more conservative than I am.
>>>>>>>>
>>>>>>>> Bill
>>>>>>>
>>>>>
>>>>>
>>>>> :-)
>>>>>
>>>>> Mark, thanks a lot for your comment and for the link. However, I see
>>>>> Bill's point
>>>>> on maintenance if this strategy becomes a policy. Tests and baseline
>>>>> images
>>>>> would also have be duplicated, in the long run I see this turning into
>>>>> an
>>>>> intricate
>>>>> web of possibilities.
>>>>>
>>>>> I'll start working on tests and stuff related to the commit. I'll add
>>>>> #ifdef's to keep
>>>>> both versions available for the time being, if anybody has any vote on
>>>>> the
>>>>> various
>>>>> ways of handling this tricky matter, please drop a note on the thread.
>>>>>
>>>>>
>>>>> Luca
>>>>>
>>>>>
>>>>> On Jan 12, 2010, at 8:38 PM, Bill Lorensen wrote:
>>>>>
>>>>>> This can turn into a maintenance nightmare. This looks like it is
>>>>>> clearly a bug. We can't create new classes for each bug we fix.
>>>>>>
>>>>>> Bill
>>>>>>
>>>>>> On Tue, Jan 12, 2010 at 11:28 AM, Mark Roden <mmroden at gmail.com> wrote:
>>>>>>>
>>>>>>> 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
>>>>>>>>
>>>>>>>
>>>>>
>>>>
>>>>
>>
>>
> _____________________________________
> 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
>