Thread Subject: curvature and radius of curvature of a plane curve

Hello,

I want to calculate the local curvature (or radius of curvature) for
a few point on outline on a picture, but how the hell am i suppose to
do that ?

Thanks,

Kenneth

Of course i want a way to do it in matlab...

In article <ef5d840.-1@webcrossing.raydaftYaTP>, bird
<djbirdnerd@hotmail.com> wrote:

> Hello,
>
> I want to calculate the local curvature (or radius of curvature) for
> a few point on outline on a picture, but how the hell am i suppose to
> do that ?
>
> Thanks,
> Kenneth
------------------
  If your "outline" points' coordinates are given by the two column
vectors x and y, then do:

 A = [x.^2+y.^2,x,y,ones(size(x))]; % Set up least squares problem
 [U,S,V] = svd(A,0); % Use economy version sing. value decompos.
 a = V(1,4); b = V(2,4); % Choose eigenvector from V
 c = V(3,4); d = V(4,4); % with smallest eigenvalue
 xc = -b/(2*a); yc = -c/(2*a); % Find center and radius of the
 r = sqrt(xc^2+yc^2-d/a); % circle, a*(x^2+y^2)+b*x+c*y+d=0

Then the circle with center (xc,yc) and radius r will be the best-fitting
circle to the points in a certain least square sense, and 1/r will be its
curvature.

  To be specific, the expression

 a*x^2 + a*y^2 + b*x + c*y + d,

subject to the constraint

 a^2 + b^2 + c^2 + d^2 = 1,

will have a least square mean value over the points.

Roger Stafford

In article
<ellieandrogerxyzzy-1507070827520001@dialup-4.232.57.106.dial1.losangeles1.level3.net>,
ellieandrogerxyzzy@mindspring.com.invalid (Roger Stafford) wrote:

> In article <ef5d840.-1@webcrossing.raydaftYaTP>, bird
> <djbirdnerd@hotmail.com> wrote:
> > I want to calculate the local curvature (or radius of curvature) for
> > a few point on outline on a picture, but how the hell am i suppose to
> > do that ?
------------------
  You may find the following method more satisfactory than the previous
one I sent. Its results are independent of any displacement of the axes'
origin. Again x and y are your coordinate column vectors.

 mx = mean(x); my = mean(y)
 X = x - mx; Y = y - my; % Get differences from means
 dx2 = mean(X.^2); dy2 = mean(Y.^2); % Get variances
 t = [X,Y]\(X.^2-dx2+Y.^2-dy2)/2; % Solve least mean squares problem
 a0 = t(1); b0 = t(2); % t is the 2 x 1 solution array [a0;b0]
 r = sqrt(dx2+dy2+a0^2+b0^2); % Calculate the radius
 a = a0 + mx; b = b0 + my; % Locate the circle's center
 curv = 1/r; % Get the curvature

  The circle defined by center (a,b) and radius r will yield the least
mean square value for the expression (x-a)^2 + (y-b)^2 - r^2 among all
possible parameters, a, b, and r, over the points in vectors x and y. The
circle's curvature will be 1/r.

Roger Stafford

Great useful post. Now how would you construct the solution to determine radius of a generalized 4 point parabola instead of a circle. By this I mean a parabola which is not aligned with x or y axis (i.e. one at some arbitrary angle.

Thank you.

Hi Roger,

The code is good. Can you give me any reference for the mathematics used. I could not understand the use of variances.

Regards,
Ruchir

"Ruchir Srivastava" <ruchir@nus.edu.sg> wrote in message <h418sh$p59$1@fred.mathworks.com>...
> Hi Roger,
> The code is good. Can you give me any reference for the mathematics used. I could not understand the use of variances.
> Regards,
> Ruchir
------------------------------
Hello Ruchir Srivastava,

  This two-year-old thread predates my Watch List, and I didn't discover your question until today. My apologies for the two-week delay.

  As you recall, I defined X and Y (uppercase) by X = x - mean(x) and Y = y - mean(y) where x and y (lowercase) are the original coordinate vectors. We can thus express the quantity whose mean square value is to be minimized as

 (X-a0)^2 + (Y-b0)^2 - r^2.

  For any given fixed values of a0 and b0, it can readily be seen by taking the partial derivative of its mean square with respect to r^2, the mean square will reach a minimum with respect to r variation at the value

 r^2 = mean((X-a0)^2+(Y-b0)^2).

  Substituting this for r^2 in the above gives an expression involving only the two parameters a0 and b0:

 (X-a0)^2 - mean((X-a0)^2) + (Y-b0)^2 - mean((Y-b0)^2)
 = -2*X*a0 - 2*Y*b0 + X^2 - mean(X^2) + Y^2 - mean(Y^2)

for which we must find the least mean square value (taking advantage of the fact that mean(X) = mean(Y) = 0.) This is equivalent to finding a0 and b0 for the least mean square approximation to the equations

 a0*X + b0*Y = (X^2-mean(X^2)+Y^2-mean(Y^2)) / 2

  As you can see, this is the least squares problem that I used the backslash operator to solve. The final determinations of a, b, and r then readily follow from this.

Roger Stafford

Hi Y'all,

I am working with images that have crossing clockwise and counter clockwise arcs in them.

Very few straight lines or circles and about forty open ended arcs that are well defined against a dark background. I can use old fashion onion skin paper and trace all the clockwise or counter clockwise arcs for a photograph within a few minutes.

Would like to write a matlab program that does the same thing. Find the arcs that meets the right criteria in a 1024x1024 matrix and create two new matrixes , one each for only cw or ccw arcs.

Can anyone point me in the right direction? Thinking I need a line finder function, then some logic to keep the same lines from being picked more than once, then a function that defines each line as cw or ccw.

What gets me is the human mind can do this nearly without thinking, but I’m not sure how I could code it.

Thanks,

Michael

"Roger Stafford" <ellieandrogerxyzzy@mindspring.com.invalid> wrote in message <h554pn$dmf$1@fred.mathworks.com>...
> "Ruchir Srivastava" <ruchir@nus.edu.sg> wrote in message <h418sh$p59$1@fred.mathworks.com>...
> > Hi Roger,
> > The code is good. Can you give me any reference for the mathematics used. I could not understand the use of variances.
> > Regards,
> > Ruchir
> ------------------------------
> Hello Ruchir Srivastava,
>
> This two-year-old thread predates my Watch List, and I didn't discover your question until today. My apologies for the two-week delay.
>
> As you recall, I defined X and Y (uppercase) by X = x - mean(x) and Y = y - mean(y) where x and y (lowercase) are the original coordinate vectors. We can thus express the quantity whose mean square value is to be minimized as
>
> (X-a0)^2 + (Y-b0)^2 - r^2.
>
> For any given fixed values of a0 and b0, it can readily be seen by taking the partial derivative of its mean square with respect to r^2, the mean square will reach a minimum with respect to r variation at the value
>
> r^2 = mean((X-a0)^2+(Y-b0)^2).
>
> Substituting this for r^2 in the above gives an expression involving only the two parameters a0 and b0:
>
> (X-a0)^2 - mean((X-a0)^2) + (Y-b0)^2 - mean((Y-b0)^2)
> = -2*X*a0 - 2*Y*b0 + X^2 - mean(X^2) + Y^2 - mean(Y^2)
>
> for which we must find the least mean square value (taking advantage of the fact that mean(X) = mean(Y) = 0.) This is equivalent to finding a0 and b0 for the least mean square approximation to the equations
>
> a0*X + b0*Y = (X^2-mean(X^2)+Y^2-mean(Y^2)) / 2
>
> As you can see, this is the least squares problem that I used the backslash operator to solve. The final determinations of a, b, and r then readily follow from this.
>
> Roger Stafford

"Michael Snyder" <sirzerp@mathworks.com> wrote in message <h5hun6$cb$1@fred.mathworks.com>...
> Hi Y'all,
>
> I am working with images that have crossing clockwise and counter clockwise arcs in them.
>
> Very few straight lines or circles and about forty open ended arcs that are well defined against a dark background. I can use old fashion onion skin paper and trace all the clockwise or counter clockwise arcs for a photograph within a few minutes.
>
> Would like to write a matlab program that does the same thing. Find the arcs that meets the right criteria in a 1024x1024 matrix and create two new matrixes , one each for only cw or ccw arcs.
>
> Can anyone point me in the right direction? Thinking I need a line finder function, then some logic to keep the same lines from being picked more than once, then a function that defines each line as cw or ccw.
>
> What gets me is the human mind can do this nearly without thinking, but I’m not sure how I could code it.
>
> Thanks,
>
> Michael

Hello Michael,

  What is the difference between a clockwise and a counterclockwise arc? If you move from one of its ends toward the other, that might be clockwise motion, but if you start from the opposite end, it would then be counterclockwise motion. If you rotate an arc of a certain angle through 360 degrees you have gone through all possible such arcs. At what point does it change from being clockwise to counterclockwise? What I am saying is that there is nothing in a circular arc that makes it inherently clockwise or counterclockwise; it is only motion from one end toward the other end that would render it so.

  Of course there is a difference between an arc that is convex toward the left from one convex toward the right, or between convex downwards and convex upwards, but that is a different matter.

Roger Stafford

"Michael Snyder" <sirzerp@mathworks.com> wrote in message <h5hun6$cb$1@fred.mathworks.com>...
> Hi Y'all,
>
> I am working with images that have crossing clockwise and counter clockwise arcs in them.
>
> Very few straight lines or circles and about forty open ended arcs that are well defined against a dark background. I can use old fashion onion skin paper and trace all the clockwise or counter clockwise arcs for a photograph within a few minutes.
>
> Would like to write a matlab program that does the same thing. Find the arcs that meets the right criteria in a 1024x1024 matrix and create two new matrixes , one each for only cw or ccw arcs.
>
> Can anyone point me in the right direction? Thinking I need a line finder function, then some logic to keep the same lines from being picked more than once, then a function that defines each line as cw or ccw.
>
> What gets me is the human mind can do this nearly without thinking, but I’m not sure how I could code it.
>
> Thanks,
>
> Michael

You might take a look at Hough's transform (the circle version).

Bruno

>
> Hello Michael,
>
> What is the difference between a clockwise and a counterclockwise arc? If you move from one of its ends toward the other, that might be clockwise motion, but if you start from the opposite end, it would then be counterclockwise motion. If you rotate an arc of a certain angle through 360 degrees you have gone through all possible such arcs. At what point does it change from being clockwise to counterclockwise? What I am saying is that there is nothing in a circular arc that makes it inherently clockwise or counterclockwise; it is only motion from one end toward the other end that would render it so.
>
> Of course there is a difference between an arc that is convex toward the left from one convex toward the right, or between convex downwards and convex upwards, but that is a different matter.
>
> Roger Stafford


Hi Roger,

Good points. I have worked with the images for so long that I had forgotten that there was an implicit reference point for cw and ccw. I feel silly; an arc segment would have to have beginning point for it to have a cw or ccw direction.

Assume there is a reference point directly in the center of the matrix and I want to evaluate each arc’s cw or ccw direction in regards to that center reference. In other words, the closest endpoint of an arc segment to the center, is defined as the beginning point for a segment.

Think of the physics cloud chamber photographs of electron and positron paths. It is easy just to look at the decaying spiral paths and know if it was an electron or positron. I had completely forgotten that an EM field was the implicit reference frame for such detectors.

I am not dealing with complete spirals, but with just some spiral arcs (less than 2pi of curvature) decaying into a center point. My problem being that have both possible types, cw and ccw arcs crossing each other in the same image and I want write a program to separate them into two images.

So I have to find the lines, and find out which end point of the line is closest to the center, and then characterize each line as cw or ccw, and then separate a complex image into two simpler ones. No Problem. <smile>

I think I can do most of that with a radial search, but not sure about the characterizing the cw or ccw directions. On a side note, wonder if there is already code out there from particle physics, evaluating overlaying particle spirals would be a much harder task.

Thanks,

Michael

Hi Y'all.

Think I figured it out.

All I need is two radial searches checking for inter and outer pixels referenced from a center point. If I do it at the pixel level, all I have to do is to check for offset pixel farer out and either allow or disallow that set of pixels in a new array.

I bet I can code this today.

Thanks,

Michael

"Michael Snyder" <sirzerp@mathworks.com> wrote in message <h5pcvf$903$1@fred.mathworks.com>...
> Hi Roger,
>
> Good points. I have worked with the images for so long that I had forgotten that there was an implicit reference point for cw and ccw. I feel silly; an arc segment would have to have beginning point for it to have a cw or ccw direction.
>
> Assume there is a reference point directly in the center of the matrix and I want to evaluate each arc’s cw or ccw direction in regards to that center reference. In other words, the closest endpoint of an arc segment to the center, is defined as the beginning point for a segment.
>
> Think of the physics cloud chamber photographs of electron and positron paths. It is easy just to look at the decaying spiral paths and know if it was an electron or positron. I had completely forgotten that an EM field was the implicit reference frame for such detectors.
>
> I am not dealing with complete spirals, but with just some spiral arcs (less than 2pi of curvature) decaying into a center point. My problem being that have both possible types, cw and ccw arcs crossing each other in the same image and I want write a program to separate them into two images.
>
> So I have to find the lines, and find out which end point of the line is closest to the center, and then characterize each line as cw or ccw, and then separate a complex image into two simpler ones. No Problem. <smile>
>
> I think I can do most of that with a radial search, but not sure about the characterizing the cw or ccw directions. On a side note, wonder if there is already code out there from particle physics, evaluating overlaying particle spirals would be a much harder task.
>
> Thanks,
> Michael

Hello again Michael,

  Your criterion, "the closest endpoint of an arc segment to the center, is defined as the beginning point for a segment", would do the job all right, but it seems very unnatural to me. How is such a "reference" center related to "spiral arcs ... decaying into a center point"? A spiral arc decaying toward a central point, as opposed to a strictly circular arc, does indeed have an inherent clockwise versus counterclockwise orientation, assuming the decay process is sufficiently evident in the arc's varying curvature. Why not use this direction of "decay" as your criterion rather than some artificial reference point in the "center of the matrix"?

  Of course all this begs the most fundamental question. Such determinations of clockwise versus counterclockwise orientation are vastly simpler than the truly formidable task of detecting intersecting arcs in the first place and then separating them from one another.
 
Roger Stafford

Roger,

Thanks so much for your post describing how to calculate the radius of curvature. I added a dimension for calculating in 3D and it seems to work fine. I confess that I discovered your post after trying to calculate the radius of curvature "directly" and failing. By "directly" I mean that (given three arrays representing the x,y, and z coordinates) I was attempting to calculate the curvature (dT/dS) where T is the tangent line to my points and S is the length of the curve. It's been quite a while since my calculus days & so I cracked open Thomas & Finney and although I think I get the concept... the examples all lead to nice differentiable functions which conveniently disappear (converge to one or zero) when needed. I (on the other hand) have arrays of numbers that I'm not sure what to do with.

So... although I have the "answer" (yours), I'm a little fixated as to why I can't solve this. I think that I need to take my array of numbers (representing all the incremental curvatures) and fit them to a curve... so then the "direct" approach ends up being the same as your "indirect" approach.

So... I guess I'm just looking for some guidance.
Please get back to me when you have a chance.

Thanks,

Dan

"Roger Stafford" <ellieandrogerxyzzy@mindspring.com.invalid> wrote in message <h554pn$dmf$1@fred.mathworks.com>...
> "Ruchir Srivastava" <ruchir@nus.edu.sg> wrote in message <h418sh$p59$1@fred.mathworks.com>...
> > Hi Roger,
> > The code is good. Can you give me any reference for the mathematics used. I could not understand the use of variances.
> > Regards,
> > Ruchir
> ------------------------------
> Hello Ruchir Srivastava,
>
> This two-year-old thread predates my Watch List, and I didn't discover your question until today. My apologies for the two-week delay.
>
> As you recall, I defined X and Y (uppercase) by X = x - mean(x) and Y = y - mean(y) where x and y (lowercase) are the original coordinate vectors. We can thus express the quantity whose mean square value is to be minimized as
>
> (X-a0)^2 + (Y-b0)^2 - r^2.
>
> For any given fixed values of a0 and b0, it can readily be seen by taking the partial derivative of its mean square with respect to r^2, the mean square will reach a minimum with respect to r variation at the value
>
> r^2 = mean((X-a0)^2+(Y-b0)^2).
>
> Substituting this for r^2 in the above gives an expression involving only the two parameters a0 and b0:
>
> (X-a0)^2 - mean((X-a0)^2) + (Y-b0)^2 - mean((Y-b0)^2)
> = -2*X*a0 - 2*Y*b0 + X^2 - mean(X^2) + Y^2 - mean(Y^2)
>
> for which we must find the least mean square value (taking advantage of the fact that mean(X) = mean(Y) = 0.) This is equivalent to finding a0 and b0 for the least mean square approximation to the equations
>
> a0*X + b0*Y = (X^2-mean(X^2)+Y^2-mean(Y^2)) / 2
>
> As you can see, this is the least squares problem that I used the backslash operator to solve. The final determinations of a, b, and r then readily follow from this.
>
> Roger Stafford

"Dan" <samberpe@gmail.com> wrote in message <hivv5j$pfd$1@fred.mathworks.com>...
> Roger,
>
> Thanks so much for your post describing how to calculate the radius of curvature. I added a dimension for calculating in 3D and it seems to work fine. I confess that I discovered your post after trying to calculate the radius of curvature "directly" and failing. By "directly" I mean that (given three arrays representing the x,y, and z coordinates) I was attempting to calculate the curvature (dT/dS) where T is the tangent line to my points and S is the length of the curve. It's been quite a while since my calculus days & so I cracked open Thomas & Finney and although I think I get the concept... the examples all lead to nice differentiable functions which conveniently disappear (converge to one or zero) when needed. I (on the other hand) have arrays of numbers that I'm not sure what to do with.
>
> So... although I have the "answer" (yours), I'm a little fixated as to why I can't solve this. I think that I need to take my array of numbers (representing all the incremental curvatures) and fit them to a curve... so then the "direct" approach ends up being the same as your "indirect" approach.
>
> So... I guess I'm just looking for some guidance.
> Please get back to me when you have a chance.
>
> Thanks,
>
> Dan
-------
Hello Dan.

  You have me a bit puzzled. Finding the curvature of the unique circle in three dimensions which passes through three given (non-collinear) points is not particularly difficult. However, extending the method I presented in this thread back in 15 Jul 2007 for a set of more than three points to the three-dimensional case involves some considerable changes. It requires six parameters to uniquely determine a circle in three-dimensions, and trying to do some kind of least squares best fit would involve adjusting all six of these in an optimum way. Did you have this in mind when you said, "I added a dimension for calculating in 3D and it seems to work fine"?

  On the other hand, trying to accomplish the task by finding the curvatures of circles through successive sets of just three points each, faces the difficulty that small amounts of errors or noise in point positions can introduce large amounts of errors in the successive curvatures. When you spoke of "attempting to calculate the curvature (dT/dS) ...", it seems as though something like that is what you were attempting to do. In my opinion if you are dealing with a number of points in a space curve substantially greater than three and you want the curvature in their neighborhood, you need to use some sort of "best fit" procedure to get a reasonable answer.

Roger Stafford

Hi Roger,
Thank you. I found your short sample code and explanation very useful. With your permission, I would like to use it in an open-source (GPL) analysis tool I am writing. I would acknowledge you and reference this posting. I can't seem to contact you via email directly. Please let me know how best to reach you.

Best,
Andy


"Roger Stafford" <ellieandrogerxyzzy@mindspring.com.invalid> wrote in message <h554pn$dmf$1@fred.mathworks.com>...
> "Ruchir Srivastava" <ruchir@nus.edu.sg> wrote in message <h418sh$p59$1@fred.mathworks.com>...
> > Hi Roger,
> > The code is good. Can you give me any reference for the mathematics used. I could not understand the use of variances.
> > Regards,
> > Ruchir
> ------------------------------
> Hello Ruchir Srivastava,
>
> This two-year-old thread predates my Watch List, and I didn't discover your question until today. My apologies for the two-week delay.
>
> As you recall, I defined X and Y (uppercase) by X = x - mean(x) and Y = y - mean(y) where x and y (lowercase) are the original coordinate vectors. We can thus express the quantity whose mean square value is to be minimized as
>
> (X-a0)^2 + (Y-b0)^2 - r^2.
>
> For any given fixed values of a0 and b0, it can readily be seen by taking the partial derivative of its mean square with respect to r^2, the mean square will reach a minimum with respect to r variation at the value
>
> r^2 = mean((X-a0)^2+(Y-b0)^2).
>
> Substituting this for r^2 in the above gives an expression involving only the two parameters a0 and b0:
>
> (X-a0)^2 - mean((X-a0)^2) + (Y-b0)^2 - mean((Y-b0)^2)
> = -2*X*a0 - 2*Y*b0 + X^2 - mean(X^2) + Y^2 - mean(Y^2)
>
> for which we must find the least mean square value (taking advantage of the fact that mean(X) = mean(Y) = 0.) This is equivalent to finding a0 and b0 for the least mean square approximation to the equations
>
> a0*X + b0*Y = (X^2-mean(X^2)+Y^2-mean(Y^2)) / 2
>
> As you can see, this is the least squares problem that I used the backslash operator to solve. The final determinations of a, b, and r then readily follow from this.
>
> Roger Stafford

"Andrew " <leifer@fas.harvard.edu.removethis> wrote in message <i1n85c$dkn$1@fred.mathworks.com>...
> Hi Roger,
> Thank you. I found your short sample code and explanation very useful. With your permission, I would like to use it in an open-source (GPL) analysis tool I am writing. I would acknowledge you and reference this posting. I can't seem to contact you via email directly. Please let me know how best to reach you.
>
> Best,
> Andy
- - - - - - - - -
  Andy, I am not very familiar with GPL licensing provisions, but assuming that they do not restrict the further free distribution of such computer code, I would have no objection to your plans to include my piece of matlab code from this thread in your "tool".

  It is important to realize both the strengths and weaknesses of this circle-fitting code. It does well for points that lie in an approximate arc which covers a large fraction, say a half or a third, of a complete circle even if there is considerable noise present in the points. It also does well for points that lie in an arc which constitutes only a small fraction of a circle provided there is comparatively little noise present - that is, provided the points do actually lie close to a circular arc.

  It does not do so well in cases where there is substantial noise present in the points' locations and only a small fraction of a full circle is occupied. By looking at matlab plots I have made of an ideal circle with associated noisy points on a small arc of it and comparing these with the algorithm results, it is clear that in such difficult cases, a human could very likely draw a better circle to fit the points than is generated by the code. To come up to a human's level of performance would probably require a comparatively sophisticated pattern recognition algorithm rather than the present algorithm's simple least squares approach.

  It is also important to point out to potential users that the code assumes the cartesian coordinates of points in two-dimensional space to be fitted to a circle have been previously stored in two column vectors which are referred to as x and y. Row vectors would come to grief here with matlab's backslash operator.

Roger Stafford

Hello,
I have question. I used formula to calculate the local curvature for 6 degree polynomial curve at certain x coordinate. However, When compare the calculated radius(Exmple:5 mm radius ) with experimental data(Example: 5 mm diameter). It is not radius. It always cames like diameter.

If anyone has any idea why it is not radius, I would be very thankfulll.

Yurdaer



ellieandrogerxyzzy@mindspring.com.invalid (Roger Stafford) wrote in message <ellieandrogerxyzzy-1507070827520001@dialup-4.232.57.106.dial1.losangeles1.level3.net>...
> In article <ef5d840.-1@webcrossing.raydaftYaTP>, bird
> <djbirdnerd@hotmail.com> wrote:
>
> > Hello,
> >
> > I want to calculate the local curvature (or radius of curvature) for
> > a few point on outline on a picture, but how the hell am i suppose to
> > do that ?
> >
> > Thanks,
> > Kenneth
> ------------------
> If your "outline" points' coordinates are given by the two column
> vectors x and y, then do:
>
> A = [x.^2+y.^2,x,y,ones(size(x))]; % Set up least squares problem
> [U,S,V] = svd(A,0); % Use economy version sing. value decompos.
> a = V(1,4); b = V(2,4); % Choose eigenvector from V
> c = V(3,4); d = V(4,4); % with smallest eigenvalue
> xc = -b/(2*a); yc = -c/(2*a); % Find center and radius of the
> r = sqrt(xc^2+yc^2-d/a); % circle, a*(x^2+y^2)+b*x+c*y+d=0
>
> Then the circle with center (xc,yc) and radius r will be the best-fitting
> circle to the points in a certain least square sense, and 1/r will be its
> curvature.
>
> To be specific, the expression
>
> a*x^2 + a*y^2 + b*x + c*y + d,
>
> subject to the constraint
>
> a^2 + b^2 + c^2 + d^2 = 1,
>
> will have a least square mean value over the points.
>
> Roger Stafford