[Insight-developers] level sets: curvature term and time step
Lydia Ng
lng at insightful . com
Wed, 16 Jul 2003 16:27:00 -0700
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
>=20
> Hi Bill,
>=20
> Yes, I have seen these. They are coming from the
> Compute3DMinimalCurvature method when LevelSetFunction is instantiated
in
> 2D.
>=20
> 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.
>=20
> 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.
>=20
> Josh.
>=20
> ______________________________
> 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
>=20
>=20
> On Wed, 16 Jul 2003, Lorensen, William E (Research) wrote:
>=20
> > 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 =3D curvature * CurvatureWeight * =
CurvatureSpeed();
> > > >
> > > >
> > > >Say for example, there is no advection term then in
> > > >ComputeGlobalTimeStep
> > > >
> > > >dt =3D 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 =3D m_DT =3D 0.25 (2D) and hence the algorithm takes to =
big a
> step.
> > > >
> > > >While if he had it the other way around: CurvatureSpeed() =3D 1 =
and
> > > >CurvatureWeight =3D 1000 then dt =3D 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
(=3D=3D
> > > >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 =3D curve;
> > > > curvature_term *=3D 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
> >