I have read the paper "Shape similarity retrieval under affine transforms", and use c++ to implement it.The step as follow;

  1. Calculate the first and second Gaussian Derivatives

enter image description here

  1. Convolving each coordinate of the curve with a Gaussian function,and calculate the first and second Derivatives of x,y

enter image description here

  1. Compute curvature of the curve

enter image description here

The paper main method as follow

enter image description here

Question:I compute the curvature of a rectangle as sigma is 0.66667. But the result what i got is strange. The curvature of four lines of rectangle are different and not zero? The curvature are 1.9, -6.3, -5.6, 1.5. Then, i roatate rectangle,the curvature of every point is different. Why? Can somebody help me ?Thank you very much

  • 1
    $\begingroup$ You derivative calculations must be off. For a straight line, the first derivative of the position should be constant, and the second derivative should be zero (assuming constant "speed"). $\endgroup$ – Niki Estner Sep 18 '15 at 11:04
  • $\begingroup$ @nikie Thank you,mikie.But I don't get your meaning.Can you give more details. $\endgroup$ – kookoo121 Sep 18 '15 at 11:41
  • $\begingroup$ All @nikie is saying is if the shape is a straight line: $y(u) = m x(u) + c$ then the first derivative will be a constant $\frac{d y}{du} = m \frac{d x}{du} = m k$ (provided $x(u) = u k$). Then the second derivative will be zero. $\endgroup$ – Peter K. Sep 18 '15 at 13:01
  • $\begingroup$ @nikie: It seems this is some kind of hyperbolic geometry underneath; not my field of expertise, but there might be different definitions of "straight"; in Euklidean spaces, the linear function is the shortest curve through two points; that's not necessarily so for every other algebraic entitity upon which one can draw a line. Example: Construct an arbitrary manifold with a map $\Rightarrow \mathbb R^2$ that is not identity. $\endgroup$ – Marcus Müller Sep 18 '15 at 13:05
  • $\begingroup$ @nikie I tried your example but my result is an inverted form of yours. Do have an idea why that is so?![X First Derivative](i.stack.imgur.com/9jyu8.png)![X Second Derivative](i.stack.imgur.com/CifWN.png)[![Curvature][3]][3] $\endgroup$ – yaro Apr 27 '17 at 14:13

My guess is you're calculating the derivatives wrong, or you're applying the whole thing to the wrong input data. I can show you the expected output for each step graphically, so you can compare it to what you get. I've implemented this in Mathematica; I'll post the code at the end, but I don't know if it'll help you much - Mathematica has built in functions for e.g. convolutions, that you would have to write yourself in C++.

First the input data for a rectangle would look like this:

enter image description here

That's 150 points along a rectangular path, using linear interpolation between the corner points.

Next, the Gaussian and it's derivatives (using $\sigma=3$ for illustration):

enter image description here

If you convolve these with the X and Y coordinates separately, you get this:

enter image description here

The first row is just for illustration - you don't actually need the convolution with the Gaussian itself, only the two derivatives.

At the sides of the rectangle, the first derivatives (second row) are constant and the second derivatives (third row) are 0. (You might want to write the results you get to a text file and look at them e.g. using Excel to check this.)

Then, using the curvature formula you quoted, you get this:

enter image description here

Along the straight sides of the rectangle, the curvature is 0 as expected. At the corners, the curvature of a rectangle is infinite; however, what we're calculating here isn't the real curvature, but the curvature of a Gaussian-smoothed rectangle, with filter size $\sigma=3$:

enter image description here

To check this result, I'll draw a circle with radius $r=\frac{1}{\kappa }$ at each point:

enter image description here

Mathematica code:

(* interpolate between corner points *)
corners = {{-2, -1}, {2, -1}, {2, 1}, {-2, 1}, {-2, -1}};

dist = Rescale@Prepend[Accumulate[Norm /@ Differences[corners]], 0];

int = Table[
   Interpolation[Transpose[{dist, corners[[All, i]]}], 
    InterpolationOrder -> 1], {i, 1, 2}];

pts = Outer[#1[#2] &, int, Range[0., 1 - 10^-9, 1/150]]\[Transpose];

cornerIdx = Most[dist*Length[pts] + 1];

  ListLinePlot[pts, AspectRatio -> Automatic, Mesh -> All, 
   ImageSize -> 400, PlotLabel -> "Points", 
   Epilog -> {Red, PointSize[Large], Point[pts[[cornerIdx]]]}],
  ListLinePlot[Transpose[pts], PlotRange -> All, 
   GridLines -> {cornerIdx, {}}, PlotLabel -> "X/Y Coordinates", 
   PlotLegends -> {"X", "Y"}, ImageSize -> 400]

(* Gaussian & derivatives *)
\[Sigma] = 2;

g0 = Table[
   1/(Sqrt[2 \[Pi]] \[Sigma])*Exp[-x^2/(2 \[Sigma]^2)], {x, -10., 10}];
g1 = Table[-x/(Sqrt[2 \[Pi]] \[Sigma]^3)*
    Exp[-x^2/(2 \[Sigma]^2)], {x, -10., 10}];
g2 = Table[(-\[Sigma]^2 + x^2)/(Sqrt[2 \[Pi]] \[Sigma]^5)*
    Exp[-x^2/(2 \[Sigma]^2)], {x, -10., 10}];

ListLinePlot[{g0, g1, g2}, PlotRange -> All, 
 PlotLegends -> {"Gaussian", "1st derivative of Gaussian", 
   "2st derivative of Gaussian"}]

(* Convolve x and y with gaussian&derivatives *)
derivatives = 
  Outer[ListConvolve[#1, #2, Round[Length[g1]/2] + 1] &, {g0, g1, g2},
    Transpose[pts], 1];

  ListLinePlot[#1, PlotRange -> All, PlotLabel -> #2, 
    GridLines -> {cornerIdx, {}}] & , 
    {"X (smoothed)", "Y (smoothed)"},
    {Overscript[X, \[Bullet]], 
     Overscript[Y, \[Bullet]]}, {Overscript[X, \[Bullet]\[Bullet]], 
     Overscript[Y, \[Bullet]\[Bullet]]}}}, 2]]

{{xs, ys}, {dx1, dy1}, {dx2, dy2}} = derivatives;

(* calculate curvature *)
\[Kappa] = dx1*dy2 - dx2*dy1/((dx1^2 + dy1^2)^(3/2));
ListLinePlot[\[Kappa], PlotRange -> All, 
 PlotLabel -> "Curvature \[Kappa]", GridLines -> {cornerIdx, {}}]

ListLinePlot[Transpose[{xs, ys}], Mesh -> All]

(* osculating circle animation *)
animationIdx = 
   Abs[\[Kappa][[#]]] > 10^-2 || (Mod[#, 5] == 0) &];

Monitor[frames = Table[Rasterize@Row[
       ListLinePlot[Transpose[{xs, ys}], Mesh -> All, 
        ImageSize -> 400, 
        PlotLabel -> "Osculating circle r=1/\[Kappa]",
        Epilog -> Module[{r, pt, n},
           pt = {xs[[i]], ys[[i]]};
           r = Clip[1/\[Kappa][[i]], {-1, 1}*10^4];
           n = Normalize[{dy1[[i]], -dx1[[i]]}];
           {Red, PointSize[Large], Point[pt], Circle[pt - n*r, r]}
       ListLinePlot[\[Kappa], PlotRange -> All, 
        PlotLabel -> "Curvature \[Kappa]", 
        GridLines -> {cornerIdx, {}}, ImageSize -> 400, 
        Epilog -> {Red, Line[{{i, 0}, {i, 10}}]}]
       }], {i, animationIdx}], i];
  • $\begingroup$ Thanks for your comment.I show my code here.and the full code is here.code $\endgroup$ – kookoo121 Sep 19 '15 at 8:36
  • $\begingroup$ I don't kown why the x'' fig is stange $\endgroup$ – kookoo121 Sep 19 '15 at 10:01

Your Answer

By clicking “Post Your Answer”, you agree to our terms of service, privacy policy and cookie policy

Not the answer you're looking for? Browse other questions tagged or ask your own question.