[Insight-developers] level sets: curvature term and time step
Joshua Cates
cates at sci . utah . edu
Wed, 16 Jul 2003 17:13:15 -0600 (MDT)
Hi Bill,
Yes, I have seen these. They are coming from the
Compute3DMinimalCurvature method when LevelSetFunction is instantiated in
2D.
The minimal curvature calculation involves finding eigenvalues and it was
necessary for us to specialize for the 3D case (vnl ND eigensolver is not
thread safe and is not optimized for our 3D case). Since we are using a
switch statement to spawn the 3D case (template specialization is still
not an option, right?), the compiler is instantiating the 3D function in
the 2D case even though it will never be used and therefore we get the
subscript out of range message.
I'm not sure what the best solution to this is. Assuming the switch
statement is functioning properly, there should be no danger of
Compute3DMinimalCurvature ever being invoked in 2D. But a better
workaround for template specialization would be nice in this case.
Josh.
______________________________
Josh Cates
School of Computer Science
University of Utah
Email: cates at sci . utah . edu
Phone: (801) 587-7697
URL: http://www . sci . utah . edu/~cates
On Wed, 16 Jul 2003, Lorensen, William E (Research) wrote:
> Josh,
> Look at the viggen.crd build and you will see warnings about subscripts out
> of range in itkLevelSetFunction.txx. This may be why the test is crashing on
> Borland. It also shows the importance of getting the warnings down to zero.
>
> Bill
>
>
> -----Original Message-----
> From: Joshua Cates [mailto:cates at sci . utah . edu]
> Sent: Thursday, July 10, 2003 6:03 PM
> To: Paul Yushkevich
> Cc: Lydia Ng; Insight-Developers (E-mail)
> Subject: Re: [Insight-developers] level sets: curvature term and time
> step
>
>
> Hi Paul and Lydia,
>
> Yes I agree that the CurvatureSpeed() should be taken account in the
> time step calculation. I'll have a look at this and get it resolved.
>
> Josh.
>
> ______________________________
> Josh Cates
> School of Computer Science
> University of Utah
> Email: cates at sci . utah . edu
> Phone: (801) 587-7697
> URL: http://www . sci . utah . edu/~cates
>
>
> On Thu, 10 Jul 2003, Paul Yushkevich wrote:
>
> > I agree, in SNAP we use the curvature speed, so it would be better if it
> > were taken into consideration.
> >
> > Paul
> >
> > Lydia Ng wrote:
> >
> > >Hi Josh,
> > >
> > >I think it is possible (but may not be likely) that all the level set
> > >filters (ie. using itk::LevelSetFunction) are affected not just the
> > >spatially varying curvature ones.
> > >
> > >We have,
> > >Curvature_term = curvature * CurvatureWeight * CurvatureSpeed();
> > >
> > >
> > >Say for example, there is no advection term then in
> > >ComputeGlobalTimeStep
> > >
> > >dt = m_DT / vnl_math_abs( m_CurvatureWeight );
> > >
> > >
> > >By default, CurvatureSpeed() returns 1 - so things might go wrong if
> > >curvature is much larger than 1 / m_DT ( > 4 for 2D ).
> > >
> > >But I don't know how likely this scenario is.
> > >
> > >________________
> > >
> > >The problem Nils ( email attached ) is having is because his
> > >CurvatureSpeed() returns 1000 and he has a CurvatureWeight of 1.
> > >
> > >Say if curvature is approx 1 then curvature_term is 1000
> > >But dt = m_DT = 0.25 (2D) and hence the algorithm takes to big a step.
> > >
> > >While if he had it the other way around: CurvatureSpeed() = 1 and
> > >CurvatureWeight = 1000 then dt = 0.25/1000 and the algorithm remains
> > >stable.
> > >
> > >So I think the time step calculation to be general it needs to take into
> > >account the whole curvature_term rather than just the CurvatureWeight.
> > >
> > >- Lydia
> > >
> > >
> > >
> > >>-----Original Message-----
> > >>From: Joshua Cates [mailto:cates at sci . utah . edu]
> > >>Sent: Thursday, July 10, 2003 10:32 AM
> > >>To: Lydia Ng
> > >>Cc: Insight-Developers (E-mail)
> > >>Subject: Re: [Insight-developers] level sets: curvature term and time
> > >>
> > >>
> > >step
> > >
> > >
> > >>Hi Lydia,
> > >>
> > >>Just to clarify: you are talking about the spatially varying curvature
> > >>term, and not the curvature calculation itself? So only filters
> > >>which use this term should be affected, right?
> > >>
> > >>Josh.
> > >>
> > >>______________________________
> > >> Josh Cates
> > >> School of Computer Science
> > >> University of Utah
> > >> Email: cates at sci . utah . edu
> > >> Phone: (801) 587-7697
> > >> URL: http://www . sci . utah . edu/~cates
> > >>
> > >>
> > >>On Thu, 10 Jul 2003, Lydia Ng wrote:
> > >>
> > >>
> > >>
> > >>>Question for those interested in level sets:
> > >>>
> > >>>Currently the itk::LevelSetFunction class does not take the
> > >>>
> > >>>
> > >curvature
> > >
> > >
> > >>>term into account when computing the time step to meet the CFL
> > >>>
> > >>>
> > >criterion
> > >
> > >
> > >>>to ensure numerical stability.
> > >>>
> > >>>This is causing some grief with users playing around with the
> > >>>
> > >>>
> > >scaling
> > >
> > >
> > >>>parameter in each of the term in particular when they set the
> > >>>
> > >>>
> > >curvature
> > >
> > >
> > >>>scaling to be high relative to the propagation and advection terms.
> > >>>
> > >>>Can someone point me to a paper/book that talks about have to
> > >>>
> > >>>
> > >compute
> > >
> > >
> > >>>the CFL time step for the curvature term?
> > >>>
> > >>>Alternatively, (hacking versus mathematically sound) could we just
> > >>>restrict the time step such that the change in level set value due
> > >>>
> > >>>
> > >to
> > >
> > >
> > >>>the curvature term is no larger that half a pixel? Would this have
> > >>>
> > >>>
> > >the
> > >
> > >
> > >>>desired affect?
> > >>>
> > >>>The second alternative is for me to change my filters
> > >>>ShapeDetectionLevelSetImageFilter and
> > >>>GeodesicActiveContourLevelSetImageFilter just to use a constant time
> > >>>step - making it the user's problem to ensure they choose the time
> > >>>
> > >>>
> > >step
> > >
> > >
> > >>>small enough. This is how it was before, when I got complaints from
> > >>>other users about how it was difficult to set this parameter...
> > >>>
> > >>>Cheers,
> > >>>Lydia
> > >>>_______________________________________________
> > >>>Insight-developers mailing list
> > >>>Insight-developers at itk . org
> > >>>http://www . itk . org/mailman/listinfo/insight-developers
> > >>>
> > >>>
> > >>>
> > >
> > >
> > >
> > >
> > > ------------------------------------------------------------------------
> > >
> > > Subject:
> > > RE: [Insight-users] Crash in ShapeDetectionLevelSetImageFilter
> > > From:
> > > "Nils Hanssen" <hanssen at caesar . de>
> > > Date:
> > > Thu, 10 Jul 2003 00:16:22 -0700
> > > To:
> > > "Luis Ibanez" <luis . ibanez at kitware . com>
> > >
> > >
> > >Hi Luis,
> > >
> > >thank you for the detailed answer.
> > >
> > >
> > >My initial level set conforms to the "negative inside" convention.
> > >Now, I am using a curvature weighting of 1000 (just to see what happens
> in
> > >the extreme case) and a propagation weighting of 0.001. When running the
> > >ShapeDetectionLSIF, the contour is not collapsing at all. Maybe that's
> > >because of the feature image that dictates the curvature speed (==
> > >propagation speed in the ShapeDetectionLSIF)?
> > >
> > >Since the curvature speed in the ShapeDetectionLSF equals the propagation
> > >speed, I derived my own classes called "TubularLevelSetFunction" and
> > >"TubularLevelSetImageFilter". These are basically clones of the
> > >ShapeDetection* classes, but with the difference that I set the curvature
> > >speed to (constant) one in the Level-set function.
> > >Now, the contour is collapsing according to the mean curvature but all
> this
> > >happens _very_ slowly (the curvature weighting is 1000!).
> > >
> > >In itk::LevelSetFunction::ComputeUpdate(...) I see the following:
> > >---
> > > curvature_term = curve;
> > > curvature_term *= m_CurvatureWeight * this->CurvatureSpeed(it, offset);
> > >---
> > >so it should make no difference when I set the curvature weight to 1 and
> the
> > >curvature speed to (constant) 1000 (and not vice versa), right? But when
> I
> > >do this, I get numerical instabilities during the evolution. What is
> wrong
> > >with my assumption?
> > >
> > >Luis, what was the curvature weighting and curvature speed when you
> tested
> > >the collapsing contour with the syntethic sphere image? Was the curvature
> > >speed constant or dependent on the image?
> > >
> > >For now, I do all calculations in 2D. Could that be a problem? I have
> > >installed a observer that shows me the progress of the contour after each
> > >update event.
> > >
> > >
> > >Regards,
> > >Nils
> > >
> > >
> > >
> > >
> > >>-----Original Message-----
> > >>From: Luis Ibanez [mailto:luis . ibanez at kitware . com]
> > >>Sent: Wednesday, July 09, 2003 6:14 PM
> > >>To: Nils Hanssen
> > >>Cc: Insight-users at public . kitware . com
> > >>Subject: Re: [Insight-users] Crash in
> > >>ShapeDetectionLevelSetImageFilter
> > >>
> > >>
> > >>
> > >>Hi Nils,
> > >>
> > >>You are right in your interpretation of the
> > >>behavior for the ShapeDetectionLevelSetFilter.
> > >>
> > >>Increasing the curvature scaling while decreasing
> > >>the propagation scaling should make the contour
> > >>collapse (contract).
> > >>
> > >>I just verified this behavior using a syntethic
> > >>image of a sphere.
> > >>
> > >>Note that in version 1.2 this filter is expecting
> > >>the input level set to conform to the convention
> > >>of "negative inside" which means that the level
> > >>set has negative values 'inside' the contour and
> > >>positive values outside the contour.
> > >>
> > >>You may have to push the curvature scaling to
> > >>a value higher than 1.0. (e.g. 5.0 or so).
> > >>
> > >>I would avoid to use 0.0 for the propagation
> > >>scaling ( in part just for superstition against
> > >>zeros values as parameters). You may want to
> > >>try something like 0.1 for the propagation
> > >>scaling.
> > >>
> > >>Please verify the convention used by your initial
> > >>level set. That may be the cause for the contour
> > >>not behaving as expected.
> > >>
> > >>
> > >>Regards,
> > >>
> > >>
> > >> Luis
> > >>
> > >>
> > >>
> > >>---------------------
> > >>Nils Hanssen wrote:
> > >>
> > >>
> > >>>Hi Luis,
> > >>>
> > >>>I am using the 1.2.0 web release.
> > >>>
> > >>>Regards,
> > >>>Nils
> > >>>
> > >>>
> > >>>
> > >>>
> > >>>>-----Original Message-----
> > >>>>From: insight-users-admin at itk . org
> > >>>>[mailto:insight-users-admin at itk . org]On
> > >>>>Behalf Of Luis Ibanez
> > >>>>Sent: Tuesday, July 08, 2003 3:55 PM
> > >>>>To: Nils Hanssen
> > >>>>Cc: Insight-users at public . kitware . com
> > >>>>Subject: Re: [Insight-users] Crash in
> > >>>>ShapeDetectionLevelSetImageFilter
> > >>>>
> > >>>>
> > >>>>
> > >>>>Hi Nils,
> > >>>>
> > >>>>Are you using a very recent checkout of ITK ?
> > >>>>
> > >>>>The level set filters have been reworked in
> > >>>>recent days, so it will help us to know if you
> > >>>>are experiencing this behavior with the current
> > >>>>version or with the version as it was two weeks
> > >>>>ago (or before that).
> > >>>>
> > >>>>Please let us know,
> > >>>>
> > >>>> Thanks
> > >>>>
> > >>>>
> > >>>> Luis
> > >>>>
> > >>>>
> > >>>>------------------------
> > >>>>Nils Hanssen wrote:
> > >>>>
> > >>>>
> > >>>>
> > >>>>>Hi,
> > >>>>>
> > >>>>>I am trying to understand the behaviour of the ShapeDetectionLSIF.
> > >>>>>Therefore, I set the propagation-weighting to zero. By setting the
> > >>>>>curvature-weighting to a value of one I would expect that
> > >>>>>
> > >>>>>
> > >>>>the inital
> > >>>>
> > >>>>
> > >>>>
> > >>>>>surface is shrinking to a point (I set the MaxRMSError to
> > >>>>>
> > >>>>>
> > >>>>zero) and the
> > >>>>
> > >>>>
> > >>>>
> > >>>>>number of iterations very high.
> > >>>>>
> > >>>>>The filter is crashing in
> > >>>>>SegmentationLevelSetFunction<TImageType, TFeatureImageType>
> > >>>>>::PropagationSpeed(const NeighborhoodType &neighborhood, const
> > >>>>>FloatOffsetType &offset) const
> > >>>>>[...]
> > >>>>>--> else return (
> > >>>>>
> > >>>>>
> > >>>>static_cast<ScalarValueType>(m_SpeedImage->GetPixel(idx)) );
> > >>>>// crashing here
> > >>>>
> > >>>>
> > >>>>
> > >>>>>[...]
> > >>>>>
> > >>>>>When I set the propagation-weighting to - for instance -
> > >>>>>
> > >>>>>
> > >>0.0001 the
> > >>
> > >>
> > >>>>>filter is not crashing, but the contour is not shrinking
> > >>>>>
> > >>>>>
> > >>>>according to
> > >>>>
> > >>>>
> > >>>>
> > >>>>>the mean curvature.
> > >>>>>
> > >>>>>Is that the correct behaviour of the filter?
> > >>>>>
> > >>>>>I would appreciate any suggestions. Thanks!
> > >>>>>
> > >>>>>
> > >>>>>Regards,
> > >>>>>Nils
> > >>>>>
> > >>>>>
> > >>>>>-------------------------
> > >>>>>Nils Hanssen
> > >>>>>Surgical Systems Laboratory
> > >>>>>research center c ae sa r
> > >>>>>Ludwig-Erhard-Allee 2
> > >>>>>53175 Bonn
> > >>>>>fon: +49-228-9656-197
> > >>>>>fax: +49-228-9656-117
> > >>>>>___http://www . caesar . de/ssl_
> > >>>>>
> > >>>>>
> > >>>>>
> > >>>>>
> > >>>>>
> > >>>>
> > >>>>_______________________________________________
> > >>>>Insight-users mailing list
> > >>>>Insight-users at itk . org
> > >>>>http://www . itk . org/mailman/listinfo/insight-users
> > >>>>
> > >>>>
> > >>>>
> > >>>_______________________________________________
> > >>>Insight-users mailing list
> > >>>Insight-users at itk . org
> > >>>http://www . itk . org/mailman/listinfo/insight-users
> > >>>
> > >>>
> > >>>
> > >>
> > >>
> > >>
> > >>
> > >
> > >_______________________________________________
> > >Insight-users mailing list
> > >Insight-users at itk . org
> > >http://www . itk . org/mailman/listinfo/insight-users
> > >
> > >
> > >
> > >
> >
> >
> > --
> > --------------------------------
> > Paul A. Yushkevich, Ph.D.
> > President, Cognitica Corporation
> >
> > 17 Flemington Rd
> > Chapel Hill, NC 27517
> > Tel: 1-919-929-7652
> > --------------------------------
> >
> >
> >
>
> _______________________________________________
> Insight-developers mailing list
> Insight-developers at itk . org
> http://www . itk . org/mailman/listinfo/insight-developers
>