[Insight-users] Re: Curvature in Level Set

Fri, 05 Mar 2004 23:17:17 -0500

Hi Sah,

You shouldn't expect the curvature to
be 0.0125 everywhere. That may be the
curvature value for a circle... but..
you *don't* have a circle.

When you represent a circle as a binary
object in an image you are in the DIGITAL
WORLD. There is no such thing as a circle
anymore, only a discrete approximation.

If you want to be formal, your 'circle' only
has two distinct curvature values, they are
Zero and Infnite.

If you zoom enough in your image, you will
see that the circle is made out of squares,
if you visit the circle's boundary, you will
only find straight lines and discontinuities
at 90 degress. The lines have zero curvature
and the corners (discontinuities) have plus
and minus infinite curvature. Interestingly
enough, if you integrate those curvature
values all around the circle you will still
recover the expected 1/R value. This is because
there are infinitely more places with zero
curvature than places with infinite curvature.

The profile of the curvature function along
the boundary of a digital circle is a distribution
composed of dirac deltas located on the positions
of the corners, and zeros everywhere else.

When you look at the circle at a higher
scale, you do the equivalent of smoothing
the values of curvature over the boundary.
the higher the scale you look at, the smoother
and more continous those values will become.

When you estimate curvature using finite
differences, the size of the neighborhood
used as support for the computation imposes
a limit to how much that curvature could be
smoothed. You could further smooth the curvature
by using larger neighborhood, but then the extra
computation time will probably make the algorithm
unusable.

Note that your image itself is also a distribution,
you only have weighted dirac deltas on the nodes
of the image sampling grid, and zero values
everywhere. That's why you need to interpolate
in order to 'estimate' the image values in non-grid
positions.

---

You may want to look at the publications by
A. Rosenfeld on computation of geometrical
features from digital objects. For example:

- "Arcs and Curves in Digital Pictures"
- "Digital Straightness and Convexity"
- "Digital Geometry"

http://www.informatik.uni-trier.de/~ley/db/indices/a-tree/r/Rosenfeld:Azriel.html

You will find interesting the online paper by
  P.J.P. Pimienta, W.C. Carter, and E.J. Garboczi,
"Cellular Automaton Algorithm for Surface Mass Transport
Due to Curvature Gradients: Simulations of Sintering",
in particular the section on:
"Local Curvature Counting Algorithms"
http://ciks.cbt.nist.gov/~garbocz/paper28/node2.html
and the appendix
"Relationship between curvature and pixel counting"
http://ciks.cbt.nist.gov/~garbocz/paper28/node10.html#SECTION00080000000000000000

The paper by V. Kovalevsky
"Curvature in Digital 2D Images"
http://www.worldscinet.com/ijprai/15/1507/S0218001401001283.html

You should probably also look at the book
"Geometry of Digital Spaces"
by G.T. Herman.
http://www.birkhauser.com/detail.tpl?ISBN=0817638970

   Regards,

      Luis

-

BTW:
Dr. Rosenfeld passed away recently
http://www.cfar.umd.edu/~ar/
(February 22 2004).
Among many other things Dr. Rosenfeld was
a Pioneer on the application of Cellular
Automata to the analysis of digital images.

------------------
Sah Rayman wrote:
> I think I eliminated the possibility that my code was
> wrong. I did one more test using original ITK code. 
> 
> 1. Compile
> ITK/Examples/Segmentation/ThresholdSegmentationLevelSetImageFilter.cxx.
> 
> 
> 2. Make the input image to have 2 values, one for
> background, and the other for a R=80 binary circle.
> Calculate upper and lower threshold so that level set
> will grow inside the circle, and shrink otherwise.
> 
> 3. Set initial radius to be 60, make sure after 400
> iterations, level set converges to the true edge. 
> 
> 4. Set breakpoint in
> itkSegmentationLevelSetImageFilter.h,
> InitializeIteration() (Line 438), set it break only
> after count >=399 
> 
> 5. After 4, set break point in
> itkLevelSetFunction.txx, ComputeMeanCurvature() (Line
> 181)
> 
> 6. Watch ( curvature_term / gd->m_GradMagSqr )
> 
> Theoretical curvature is: 1/80 = 0.0125
> I observe a range from (-0.25, +0.5)
> That's true for both curvatureweight = 1e-6 and 1
> 
> Now the question is, is that exactly what the
> algorithm (Sethian Book, pp. 70, eq6.36) should behave
> like, or there is some bug in ITK implementation? 
> 
> I tried to had a look at some related ITK calculation,
> but don't have a clue yet.
> 
> The algorithm itself has some space for calculation
> errors. Especially, a numerical calculation for 2nd
> order derivative could be away from truth. But still,
> I am not really confident about this.
> 
> 
> --- Sah Rayman <sahrayman at yahoo.com> wrote:
> 
>>I found some "wierd behavior" of the curvature
>>calculated in SegmentationLevelSetFunction. 
>>
>>While I expect the curvature to be 0.0125
>>everywhere, 
>>the actual curvature has a mean of 0.0137 (which is
>>good), but a s.t.d of 0.4. In another words, when
>>the
>>tangent circle is actually of radius 80, the
>>calculation tells me the radius is 2 or 1 at some
>>edge
>>points, and even concave in some other points.
>>
>>Is it true that there are errors in curvature
>>calculation in the range of my discription? Or,
>>there
>>must be something wrong in my application?
>>
>>======================================
>>Here is what I do in the test.
>>
>>Full image is 256x256, ground truth is a circle
>>whose
>>radius is 80 pixels, initial is a 60-pixel-radius
>>circle with a same center.
>>
>>I disabled advection, curvature term, and there are
>>only speed term. Speed is +1 inside my ground truth,
>>and -1 otherwise. In fact, curvature weight is 1e-6
>>to
>>enable the curvature calculation.
>>
>>After 180 iterations, I observe the output of
>>ComputeMeanCurvature() for one iteration, and take
>>statistics.
>>
>>The output segmentation matches well with my ground
>>truth.
>>======================================
>>
>>
>>
>>
>>
>>__________________________________
>>Do you Yahoo!?
>>Yahoo! Search - Find what you’re looking for faster
>>http://search.yahoo.com
>>
> 
> 
> 
> __________________________________
> Do you Yahoo!?
> Yahoo! Search - Find what you’re looking for faster
> http://search.yahoo.com
> _______________________________________________
> Insight-users mailing list
> Insight-users at itk.org
> http://www.itk.org/mailman/listinfo/insight-users
> 