It seems that curves found by EdgeDetect always persist $C^{0}$ continuity only. Consider the following example:

l = 10; r = Pi/2;
(*Create an image of a smooth curve *)
pic = Rasterize@Plot[ArcTan[x], {x, -l, l}, Filling -> Bottom, Axes -> None, 
                                            PlotRangePadding -> None,PlotRange -> r];

(* Recover data points from the image *)
data = ImageValuePositions[Thinning@EdgeDetect@Binarize@pic, 1];

(* Create a interpolating function with the points *)
{w, h} = ImageDimensions@pic;
trf = Last@FindGeometricTransform[{{0, 0}, {0, r}, {l, 0}}, 
                                  {{w/2, h/2}, {w/2, h}, {w, h/2}}];
func = Interpolation[DeleteDuplicates[trf /@ data, First@# == First@#2 &]];

{{lb, rb}} = func["Domain"]
(* Check the derivatives of func *)
Plot[{ArcTan'[x], func'[x]}, {x, lb, rb}, PlotRange -> {0, 1}]

enter image description here

As one can see, the recovered solution is far from the analytic one, oscillating disastrously, full of noise, in a word, bad. So the question is, with what kind of postprocessing can I get a smooth ($C^{1}$ continuity, at least) and distortionless interpolating curve? I've played with GaussianFilter and LowpassFilter for a while but the result isn't great.

  • $\begingroup$ I had to add ImagePadding -> None to your Plot to get the same result in Mathematica 10.2 $\endgroup$ Aug 8 '15 at 16:57
  • $\begingroup$ @nikie That sounds interesting, in v9 when Axes->None is set, ImagePadding->All (I think this is the default setting) and ImagePadding->None gives the same result. Maybe it's worth asking another question for this? $\endgroup$
    – xzczd
    Aug 9 '15 at 5:59
  • $\begingroup$ In v10, I see a thin white border without ImagePadding->None, and EdgeDetect finds another set of edges at the bottom and at the right. Some subtle change due to the new default themes, perhaps. $\endgroup$ Aug 9 '15 at 7:33

In a "natural" image, you'd look at each edge pixel in, use some approximation (e.g. 2nd order polynomial) of the gradients above/below that pixel and calculate the sub-pixel position of the steepest gradient.

But in your case, all EdgeDetect gets to work on is a binary image, and any the potential anti aliasing sub-pixel information is lost. So the best you can probably do is find a curve that is as smooth as possible, while still less than 0.5 pixel from the discrete pixel values EdgeDetect found. You can do that using constrained optimization.

Reusing code from this answer:

xValues = Array[# &, w, func["Domain"]];
discreteValues = func[xValues];

n = Length[discreteValues];
vars = Array[y, n];
maxDist = 0.5 Norm[trf[{0, 0}] - trf[{0, 1}]];

Here are the optimization objectives: find a list of values y[1]..y[n] s.t. the distance to the original values discreteValues[i] is below 0.5 pixels and smoothness as small as possible:

constraints = 
  Array[discreteValues[[#]] - maxDist <= y[#] <= 
     discreteValues[[#]] + maxDist &, n];
smoothness = Total[Differences[vars, 2]^2];
startValues = Array[{y[#], discreteValues[[#]]} &, n];
{fit, sol} = 
  FindMinimum[{smoothness, constraints}, startValues, 
   AccuracyGoal -> 10];
smoothedValues = (vars /. sol);

Here's a graphic visualization of the constraints and the results:

zoomStartIdx = 250;
  ListLinePlot[Transpose[{xValues, smoothedValues}], ImageSize -> 600,
    PlotStyle -> Orange, 
   Epilog -> {EdgeForm[{Gray, Dashed}], Transparent, 
     Rectangle @@ (Transpose[{xValues, 
          smoothedValues}][[{zoomStartIdx, -1}]])}],
   ListLinePlot[{Transpose[{xValues, discreteValues - maxDist}], 
      Transpose[{xValues, discreteValues + maxDist}]}[[All, 
     zoomStartIdx ;;]],
    Filling -> {1 -> {2}}, InterpolationOrder -> 0, 
    PlotStyle -> Directive[Blue, Thin], ImageSize -> 600],
   ListLinePlot[Transpose[{xValues, smoothedValues}], 
    PlotStyle -> Orange, Mesh -> All, MeshStyle -> PointSize[Medium]]],
  LineLegend[{Orange, Blue}, {"Smoothed curve", "Constraints"}]

enter image description here Using (mostly) your code to display the result

(*Create a interpolating function with the points*)
{w, h} = ImageDimensions@pic;
funcSmooth = Interpolation[Transpose[{xValues, smoothedValues}]];

{{lb, rb}} = funcSmooth["Domain"];
(*Check the derivatives of func*)
Plot[{ArcTan'[x], funcSmooth'[x]}, {x, lb, rb}, PlotRange -> All]

enter image description here

  • $\begingroup$ Seems that FindMinimum is improved (?) in v10.2, in v9 AccuracyGoal -> 10 only makes things worse, but the result without AccuracyGoal isn't bad anyway. $\endgroup$
    – xzczd
    Aug 9 '15 at 7:10
  • 1
    $\begingroup$ @nikie Could you please give an example of using of the anti-aliasing sub-pixel information for recovering data from a raster plot with maximum exactness possible? Data recovering is a very common task but usually it is made through crude binarization of the original anti-aliased image. It would be very valuable to have an example of more advanced approach. $\endgroup$ Aug 13 '15 at 12:22

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.