# Smoothing EdgeDetect - image processing

I have to make a program that analyzes the image of saw tooth and later find the function that can describe it (or just part of it). Now I have a problem with smoothing the edges detected with EdgeDetect. After trying some codes already posted on Stack I come out to this:

img = ImageRotate[Import["https://i.stack.imgur.com/2ZBFz.jpg"], π/2];
smtFactor = 0.2;
img = TotalVariationFilter[ImageClip[img], 0.2, Method -> "Gaussian"]
imgSize = ImageDimensions[img];
imgClean =
img // Blur[#, 2] & // Binarize // ColorNegate // Erosion[#, 1] &;
edges = EdgeDetect[imgClean, Method -> Metoda] ;
Koord = Map[# &, GatherBy[PixelValuePositions[edges, 1], First] ];
L1 = ListPlot[Koord, GridLines -> Automatic, ImageSize -> Large,
Frame -> True, PlotLabel -> "Edges",
PlotRange -> {{1, imgSize[[1]]}, {0, imgSize[[2]]}}]
n = Length[Koord];
vars = Array[y, n];
constraints =
Array[Koord[[#]] - smtFactor <= y[#] <= Koord[[#]] + smtFactor &, n];
smoothness = Total[Differences[vars, 2]^2];
{fit, sol} = FindMinimum[{smoothness, constraints}, vars];
smoothKoord = vars /. sol;
ListPlot[smoothKoord,
Prolog -> {Red, Opacity[0.5],
constraints /. {y0_ <= y[x_] <= y1_ -> Line[{{x, y0}, {x, y1}}]}},
GridLines -> Automatic, ImageSize -> Large, Frame -> True,
PlotLabel -> "Edges to be smoothed",
PlotRange -> {{1, imgSize[[1]]}, {0, imgSize[[2]]}}]
Show[img, L1]


Output:

In some cases, it works, but it does not describe the bottom part (Output form a bit different code) I don't know how to solve the error. Thank you in advance!

Up to now, this is only an extended comment.

The problem is that the sawtooth (at the current orientation) is not a graph of a function: some x-values have multiple y-values. So you should not use a graph of a function as ansatz but a more general curve.

When using a list of points in the plane as unknown, you should consider the following as regularization; it's a not so untypical discretization of the Euler-Bernoulli bending energy.

getEulerBernoulli = Quiet@Block[{P, PP, X, Y, rX, rY, nX, d, nY, U, V},
d = 2;
PP = Table[P[[i, j]], {i, 1, 3}, {j, 1, d}];
X = PP[[2]] - PP[[1]];
Y = PP[[3]] - PP[[2]];
rX = Sqrt[X.X];
rY = Sqrt[Y.Y];
nX = X/rX;
nY = Y/rY;
U = nX - nY;
V = nX + nY;
P \[Function] Evaluate[Together[8 U.U/V.V/(rX + rY)]]
];

vars = Transpose[{Array[x, {n}], Array[y, n]}];
smoothness = Total[getEulerBernoulli /@ Partition[vars, 3]];


The Euler-Bernoulli bending energy of a curve $\gamma \colon [a,b] \to \mathbb{R}^2$ is given by

$$\mathscr{E}(\gamma) = \frac{1}{2}\int_a^b |\kappa_\gamma(t)|^2\, |\gamma'(t)|\, \operatorname{d} t,$$

where $\kappa_\gamma \colon [a,b] \to \mathbb{R}^2$ is the curvature vector of $\gamma$.

The idea of this discretization for a polygonal line $p = (p_0,\dotsc,p_n)$ is basically to replace the integral by a Riemann sum and the curvature by turning angles $\alpha_i = \measuredangle(p_i-p_{i-1},p_{i+1}-p_i)$:

$$\mathscr{E}_n(p) = \frac{1}{2} \sum_{i=1}^n \left( \frac{\alpha_i}{\ell_i} \right)^2 \, \ell_i = \frac{1}{2} \sum_{i=1}^n \frac{\alpha_i^2}{\ell_i},$$

where $\ell_i = \frac{1}{2}( |p_i-p_{i-1}| + |p_{i+1}-p_{i}|)$ is the average of the length of two consecutive edges.

If $n$ is not too large, Mathematica's FindMinimum may find a (local) minimizer in reasonable time.

Moreover, you will have to replace the constraints...

The price for this all is that the optimization problem is changed a quadratic one (quadratic objective with linear inequality constraints) to a nonlinear one. Thus, is might be a better idea to separate the detected edges into 2 or more groups so that each of them can be approximated by a (rotated) function graph and to apply your method onto each of these groups.

• Thank you for your answer. I was somehow thinking in the same direction, to split the curve in 3 parts (like this - (imgur.com/fF7SbXZ). But I have difficulties splitting the list Koord = Map[# &, GatherBy[PixelValuePositions[edges, 1], First] ];. I think that by using Max[Koord] the code would find the tip of the tooth. But I was unable to successfully manipulate the list. Because if it is in MatrixForm it has unequal subs. Can you please help me? Also, thank you for mentioning EulerBernoulli, I will have deeper look in it later. – Matkov Apr 8 '18 at 18:17