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. 

![enter image description here][1]

On this image you apply `ImageLines` and you get exactly two lines. But even if you don't 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]

![enter image description here][2]

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]

![enter image description here][3]

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

![enter image description here][4]

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 Andy's 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.


  [1]: https://i.sstatic.net/KEgBt.png
  [2]: https://i.sstatic.net/qGotk.png
  [3]: https://i.sstatic.net/zSW2f.png
  [4]: https://i.sstatic.net/YzAlS.png