[Insight-users] fast marching behavior

Luca Antiga luca.antiga at gmail.com
Wed Jan 13 04:48:36 EST 2010


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



More information about the Insight-users mailing list