[Insight-developers] level sets: curvature term and time step
Joshua Cates
cates at sci . utah . edu
Wed, 16 Jul 2003 17:35:50 -0600 (MDT)
Hi Lydia,
Thanks. I'll have a look at MinMaxCurvatureFlowFunction.
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, Lydia Ng wrote:
> Hi Josh,
>
> I believe I had a similar issue with MinMaxCurvatureFlowFunction.
> Jim fixed the warnings by using Brad's Dispatch solution.
> Perhaps this might do the trick in this case too?
>
> - Lydia
>
> > -----Original Message-----
> > From: Joshua Cates [mailto:cates at sci . utah . edu]
> > Sent: Wednesday, July 16, 2003 4:13 PM
> > To: Lorensen, William E (Research)
> > Cc: Paul Yushkevich; Lydia Ng; Insight-Developers (E-mail)
> > Subject: RE: [Insight-developers] level sets: curvature term and time
> step
> >
> > 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
> > >
>
>