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.

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

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 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.