[Insight-users] fast marching behavior

Wed Jan 13 10:15:23 EST 2010

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