translation, rotation, skewing (so all affine transformations)

You'll probably have more points than parameters in the affine transformation matrix, making this an overdetermined linear system. This is good: you can solve it in the least squares sense and handle a bit of the noise in the coordinates this way.

plus barrel distortion, lens distortion

This could be harder, because it's definitely not linear. Panorama stitching software like Hugin can let you optimise the radial distortion parameters at the same time as optimising the affine transform matrices; maybe you could look at its source code for some inspiration. Wikipedia lists a formula for a generalised distortion model, which you could fit the parameters for using Levenberg-Marquardt.