# Recovering data points from an image

Unfortunately, some data can only be obtained in picture form (Japanese publications anyone?). Since this cannot be worked with, it has to be converted to a dataset that can; I was wondering whether this could be done in Mathematica. Consider the following example image:

The task consists of two steps:

1. Extract a point set from the graph
2. Finding a function that takes out the distortion and applying it to the data set

Is there some streamlined way of doing this in Mathematica? The result should be a point set that is as accurate as possible, think of something whose Interpolation could be easily and accurately ploted. Bonus points for no Get Coordinates to trace the graph.

Here's the code used to generate the image above:

img = ImagePerspectiveTransformation[
Rasterize[
Plot[x^((x - 2)^2 E^-x) + E^-x, {x, 0, 10}, PlotStyle -> Thick],
ImageSize -> 400
],
{{1, 0.1, 0}, {0.1, 1, 0}, {0, 0.1, 1}},
]

-
Wait, you want to sample points from a printed curve? – J. M. Feb 8 '12 at 23:43
Exactly. I would like to be able to re-generate the plot in Mathematica with reasonable precision, for example to compare my data with the results pictured in a paper. – David Feb 8 '12 at 23:46
As a start: you could use FindGeometricTransform[] to try to find the transformation matrix you can use as the second argument of ImagePerspectiveTransformation[]; this is assuming that the distorted version has axes that you know to be actually orthogonal. – J. M. Feb 8 '12 at 23:47
I believe this tool does what you want: getdata-graph-digitizer.com/index.php ; nevertheless it's not free, and since you already have a Mathematica license, and this is a Mathematica site... – P. Fonseca Feb 9 '12 at 0:17
@P.Fonseca It's Windows only. Bad idea for a research tool ;( – David Feb 9 '12 at 0:38

I started with the image you provide and called it img. This solution isn't perfect but it might serve as a starting point.

Get some known points:

I right clicked the image and selected "Get Coordinates". I then clicked as closely as possible to the origin, and the points {0,1.3} and {10.,.82}. On Windows hold Ctrl+C to copy those points. And then Ctrl+V to paste them into the notebook...

{o, y, x} = {{36.5173, 206.72}, {17.5824, 17.3711}, {391.209, 54.9028}};


Find a transformation that will return the proper points:

Here I use FindGeometricTransform and feed it the known values for the selected points along with their image coordinates. This produces a TransformationFunction to use later.

trans = FindGeometricTransform[
{{0, .82}, {0, 1.3}, {10, .82}},
{o,      y,        x}
][[2]];


Obtain and process the image data:

Here I round the RGB color values in the ImageData so that the blue curve is coded as {0,0,1}. This will allow me to extract the curve.

data = Round[ImageData[img], 1];

col = DeleteDuplicates[Flatten[Round[ImageData[img], 1], 1]];

Graphics[{RGBColor[#], Disk[]}, ImageSize -> Tiny] & /@ col


The nice blue color I'm wanting to extract is the third color in the list. Now I binarize the image. I convert non-blue pixels to black and the blue to white.

binImage = Image@Replace[data, {col[[3]] -> 1, _ :> 0}, {2}]


But this has some spurious points I'd like to remove so I only have the curve remaining. I'll use a GaussianFilter to create a binary mask that will allow me to filter those points out. This should give me the curve I want.

curve = ImageApply[{0, 0, 0} &, binImage,
Masking -> ColorNegate[Binarize[GaussianFilter[binImage, 5]]]]


That's much cleaner! Now to extract the locations of the white pixels while maintaining the proper orientation.

curvLoc = (Reverse /@
Position[ImageData[curve, DataReversed -> True], {1., 1., 1.}]);


Apply the transformation before to the curve points and show it with the original plot before distortion. I called this plot...

Show[ListPlot[trans@curvLoc, PlotRange -> All], plot]


Its not perfect, but it should be a start.

EDIT: I realized that the coordinates of the origin were actually {0,.82} rather than {0,.8}. With this realization we get an even better approximation. Note that I've also employed an interpolating function. Using various smoothing techniques on the function values prior to interpolating should further improve things.

pts = Sort[trans@curvLoc];

g = Interpolation[pts, InterpolationOrder -> 1]

Show[Plot[g[x], {x, .05, 10}, PlotStyle->Red], plot]


-
Methinks this is the right approach to take for David's problem. Note that scans from papers will sometimes have a few smudges in addition to the distortion, so one might have to experiment with various filters for best results. – J. M. Feb 9 '12 at 2:34
Very nice! I guess you could also use DeleteSmallComponents instead of the Gaussian filter – R. M. Feb 9 '12 at 4:20
"Using various smoothing techniques on the function values prior to interpolating should further improve things." - for instance, With[{m = 3}, ListLinePlot[Transpose[MapAt[Join[Take[#, m], MovingAverage[#, 2 m + 1], Take[#, -m]] &, Transpose[SortBy[curvLoc, First]], {2}]], Frame -> True]] – J. M. Feb 9 '12 at 4:42
@JM Could you explain what that does? It's not very obvious to me. Is it some hack to apply MovingAverage to the data so the line gets thinner? If so, in what way, for example does it take the direction of the graph into account or is it simply averaging over y values? – David Feb 9 '12 at 17:32
@David: moving averages are a conventional, simple-minded technique for smoothing noisy ordinates. As implied by the name, it just takes averages of consecutive ordinates, and replaces the ordinate in the middle of the "window" of size 2m+1 with this average value. A lot of times, they're fine (as long as the window parameter m is chosen properly), but some situations might demand the use of a more elaborate smoothing filter, like Savitzky-Golay. – J. M. Feb 10 '12 at 23:56

Let me emphasize what IMO are the key-points in the image-processing here. First of all, if your images are not so bad there is no requirement to manually find the inverse transformation. What you should try is (as @kguler already mentioned) a Hough-transform which detects lines. An equivalent filter in Mathematica is given by ImageLines. So what you do is, you invert the colors of your image and binarize it with a high threshold.

On this image you apply ImageLines and you get exactly two lines. But even if you dont get only two lines, it should be possible to make an educated guess which are the right ones automatically.

lines = ImageLines[Binarize[ColorNegate[img], 0.8]]


These two lines can now be used to calculate the backward transformation because, lucky enough they represent your transformed system. So taking them, calculating the inverse and scaling it with your image-dimensions should do what you want

m = (Subtract @@ Reverse[#]) & /@ lines;
minv = DiagonalMatrix[ImageDimensions[img]*{1, -1}].Inverse[Transpose[m]]
orig=ImagePerspectiveTransformation[img, minv, Padding -> White]


But you don't want to transform your disturbed image back before you use your lines to remove the original axes. This happens simply by creating a mask and using ImageMultiply. The mask is created the same way you would draw the axis-lines you already extracted:

mask = Graphics[{Thickness[0.04], Black, Line /@ lines},
Background -> White,
PlotRange -> Transpose[{{1, 1}, ImageDimensions[img]}]];
axesFree = ImageMultiply[ColorNegate[img], mask]


What you see now is, that you have small objects (the rests of the labels) and the large curve. So why not using ImageComponents and it's buddies to select the curve. Basically it's one call to ImageComponents and then you select the image mask of the largest component:

axesFreeOrig =
ImagePerspectiveTransformation[axesFree, minv, Padding -> Black]
comp = MorphologicalComponents[axesFreeOrig];
curve = Thinning[Image[SelectComponents[comp, "Count", -1], "Bit"]]


Now having this image it is easy to extract all points with Position. While the output of this is often enough, it is never guarantied that the points are in the right order. For this you could use FindCurvePath

points = #[[First@FindCurvePath[#]]] &@
Position[Transpose@ImageData[curve, "Bit", DataReversed -> True],
1];


Since I only wanted to add something to the image processing, I'm done here. What is left open is the transformation into your data-range. Doing this automatically is not easy and therefore, I would suggest to follow Andys approach. Or you combine the best and use MorphologicalComponents for the curve extraction and FindCurvePath for the order and the rest you take from Andy.

-
you beat me by several days - or, rather, saved me as many days:) Short of going into detecting some numbers on the axes to get data-range info, this is probably the limit of what can be done in an unmanned-fashion.+1 – kglr Feb 11 '12 at 2:59

Not an answer but a comment too long for comment box on some ideas as starting points:

For a semi-manual approach, barChartDigitizer from Will DeBeest at MathGroup archive may be good starting point:

barChartDigitizer[g_Image] :=
DynamicModule[{min = {0, 0}, max = {0, 0}, xmin = -1., xmax = 1.,
pt = {0, 0}, data = {}, img = ImageDimensions[g], output},
Deploy@Column[{Row[{Column[{Button["Y Axis Min", min = pt],
InputField[Dynamic[xmin], Number, ImageSize -> 70]}],
Column[{Button["Y Axis Max", max = pt],
InputField[Dynamic[xmax], Number, ImageSize -> 70]}],
Column[{Button["Add point", AppendTo[data, pt]],
Button["Remove Last", data = Quiet@Check[Most@data, {}]]}],
Column[{Button["Start Over", data = {}],
Button["Print Output",
Print@Column[{BarChart[
output =
Rescale[#, {Last@min, Last@max}, {xmin, xmax}] & /@
data[[All, 2]],
PlotRange -> {Automatic, {xmin, xmax}},
ImageSize -> 400], output}],
Enabled -> Dynamic[data =!= {}]]}]}],
Row[{Graphics[{Inset[
Image[g, ImageSize -> img]], {Tooltip[Locator[Dynamic[pt]],
Dynamic[pt]]}}, ImageSize -> img, PlotRange -> 1,
AspectRatio -> img[[2]]/img[[1]]]}]}]]


It works nicely withbarcharts and it can be adapted to work with arbitrary graph images. When applied to the image

the palette

It can be component of a solution combined with pre-processing on the source image, like:

With

plt = Plot[x^((x - 2)^2 E^-x) + E^-x, {x, 0, 10}, PlotStyle -> Thick];
imgr = ImagePerspectiveTransformation[
Rasterize[plt,  ImageSize -> 400], {{1, .1, 0}, {.1, 1, 0}, {0, .1, 1}},

GraphicsRow@{EdgeDetect[imgr],
`