Smoothen the curve found by EdgeDetect

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,

(* 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}] 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.

• I had to add ImagePadding -> None to your Plot to get the same result in Mathematica 10.2 Aug 8 '15 at 16:57
• @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? Aug 9 '15 at 5:59
• 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. 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.

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..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;
Row[{
ListLinePlot[Transpose[{xValues, smoothedValues}], ImageSize -> 600,
PlotStyle -> Orange,
Epilog -> {EdgeForm[{Gray, Dashed}], Transparent,
Rectangle @@ (Transpose[{xValues,
smoothedValues}][[{zoomStartIdx, -1}]])}],
Show[
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"}]
}] 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] • 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. Aug 9 '15 at 7:10
• @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. Aug 13 '15 at 12:22
{"ft_lang_label":"__label__en","ft_lang_prob":0.8136352,"math_prob":0.97542083,"size":1568,"snap":"2022-05-2022-21","text_gpt3_token_len":449,"char_repetition_ratio":0.08695652,"word_repetition_ratio":0.0,"special_character_ratio":0.28954083,"punctuation_ratio":0.18456376,"nsfw_num_words":0,"has_unicode_error":false,"math_prob_llama3":0.9872814,"pos_list":[0,1,2,3,4,5,6],"im_url_duplicate_count":[null,7,null,7,null,7,null],"WARC_HEADER":"{\"WARC-Type\":\"response\",\"WARC-Date\":\"2022-01-22T00:00:21Z\",\"WARC-Record-ID\":\"<urn:uuid:843d117d-088a-4593-a300-71fc20fc8ee7>\",\"Content-Length\":\"148503\",\"Content-Type\":\"application/http; msgtype=response\",\"WARC-Warcinfo-ID\":\"<urn:uuid:f450b044-754e-4c7d-93e6-5058c2a3d628>\",\"WARC-Concurrent-To\":\"<urn:uuid:f120d4f7-0a2a-4c7e-96a5-1b7a1d94a560>\",\"WARC-IP-Address\":\"151.101.65.69\",\"WARC-Target-URI\":\"https://mathematica.stackexchange.com/questions/90130/smoothen-the-curve-found-by-edgedetect?noredirect=1\",\"WARC-Payload-Digest\":\"sha1:KUIH7SETVHMJX3R54MYJ3ZALWPGK3SKX\",\"WARC-Block-Digest\":\"sha1:HTU5MFPJME7RAR6MRKREQHF2LENOCU3X\",\"WARC-Identified-Payload-Type\":\"text/html\",\"warc_filename\":\"/cc_download/warc_2022/CC-MAIN-2022-05/CC-MAIN-2022-05_segments_1642320303717.35_warc_CC-MAIN-20220121222643-20220122012643-00368.warc.gz\"}"}