Why the curvature of straight line in curvature scale space is not zero?

Question

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

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

3. Compute curvature of the curve

The paper main method as follow

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

Answer 1

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:

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):

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

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:

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$:

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

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];

Row[{
  ListLinePlot[pts, AspectRatio -> Automatic, Mesh -> All, 
   ImageSize -> 400, PlotLabel -> "Points", 
   Epilog -> {Red, PointSize[Large], Point[pts[[cornerIdx]]]}],
  Spacer[50],
  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];

Grid[MapThread[
  ListLinePlot[#1, PlotRange -> All, PlotLabel -> #2, 
    GridLines -> {cornerIdx, {}}] & , 
     {derivatives,
    {
    {"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 = 
  Select[Range[Length[xs]], 
   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];