Why is un-distorting so mathematically complex? Isn't it just a curve? Simple code to un-distort?

[Benjamin Hill]

My goal is to take a series of night-sky photographs and stack them.  To do this, I need to un-distort my images.  I've located the stars in the image, so I need to run each star through

stars.forEach { star ->

  val (realLocationX, realLocationY) = magicUnDistortFunction( star.x, star.y, magicDistortionConstant1,...  )

}

Which means I'm searching for magicUnDistortFunction(), and then I can try a few million values for the magicDistortionConstant until things match, and then I will have my camera lens "calibrated" and re-use those constants next time.  (Hopefully I can avoid needing magicDistortionConstant2, 3, 4, etc.)

I've found some code snippets to un-distort photographs, but all are quite complex, and I'm not sure why.

1. http://lensfun.sourceforge.net/calibration-tutorial/lens-distortion.html
   
2. http://www.tannerhelland.com/4743/simple-algorithm-correcting-lens-distortion/
   
3. https://github.com/Mihaylov93/UnBarrel
   
4. http://hugin.sourceforge.net/docs/manual/Lens_correction_model.html#Determine_lens_correction
   

But... these all seem quite complex, and require knowing about how to factor a third-order polynomial to go from a regular captured photo to an undistorted model.  I'll admit - I don't get why a simplified model doesn't work "close enough".  Or maybe I haven't found the right code snippet yet.

For most distortion, the center of the image is "ok", then gets warped as you move towards the edge. 

As you move from the center of the image to the right along the x axis:

x= your pixel's location

y= the actual star location

then

blue line is a perfect "your pixels = reality 1:1 correspondence, no distortion to be corrected"

red line is "your pixels = farther from reality as you move towards the edge"

And because it isn't a perfectly square camera, maybe a third line for moving along the Y axis, but it would be nice if it was a similar line, and even nicer if you could figure out the offset for any location using the same formula.

So... I think I need that formula.  

Option 1: From https://github.com/Mihaylov93/UnBarrel

rSrc = (a * rDest^3 + b * rDest^2 + c * rDest + d) * rDest

which is rSrc and rDest?  rSrc is the source image... the reality version?  Great, so I don't need to invert the function.  But that has 4 constants.

Option 2: Rd = k1 * Ru^3 + (1 - k1) * Ru

I really hope Rd == rSrc, or else I have to invert the function...  

And then I found that this only works if you "normalize" the coordinates, where radius = the biggest possible circle inside the photo along the shortest dimension.  ack.

[Dustin Lang]

In Astrometry.net, we compute "SIP" distortion coefficients -- and their inverses -- and save them in the wcs.fits file.

In SIP, the "distort" and "undistort" functions are both polynomials.  Why polynomials?  They're simple and general.

Why is it hard?  In general, it's hard because you can't always determine which stars match -- so you don't necessarily know that a dot in your image corresponds to a star with known RA,Dec.  (And you have to fit the center, scale, and rotation of the image as well.)

The polynomials, if they're of order N (3 probably works fine), have about N*(N-1) terms for each of X and Y -- the magicUndistortFunction is

Xu = AP00 + AP01*X + AP10*Y + AP02*X**2 + AP11*X*Y + AP20*Y**2 + AP03*X**3 + AP12*X**2*Y + AP21*X*Y**2 + AP30*Y**3

Yu = BP.... samesies

If you wanted something simpler, you could try a radial function -- instead of fitting X and Y separately, fit for a factor that scales both X and Y away from the image center by a factor that varies with radius.

Either way, inverting that function is a pain -- in Astrometry.net we just instantiate the distorted values on a grid and then fit for the undistorted coefficients by a least-squares fit.

cheers,

--dustin

[Eric Sibert]

> My goal is to take a series of night-sky photographs and stack them.

I took a sequence of wide field highly distorted night-sky photographs
and wanted to stack them using Iris
(http://www.astrosurf.com/buil/iris-software.html).

Indeed, due to high distortion and poor tracking, Iris didn't wanted
to align the sequence. So, after pretreatment, I exported pictures to
Hugin. I aligned pictures in Hugin and reprojected them (incl.
distortion correction). Then, I imported back reprojected pictures in
Iris and stacked them.

Result:

https://www.webastro.net/forums/topic/160980-la-voie-lact%C3%A9e/?do=findComment&comment=2504869

May this solve uour problem?

Eric

[Benjamin Hill]

First: thank you for the reply!  Very clear comment.

Yes, I don't know which stars match.  I'm doing a VERY slow optimization algo that calculates "Across all stars in this image compared to all stars in another image, for each star in this image, get the closest distance to a star in the other image, then take the median of all those distances" to check how "correct" my rotation and distortion is.  It isn't great, but it is the best I have currently.

I think I've got the radial version working, if I haven't gone and inverted everything.

(Not shown: scaleAroundAnchor and rotateAroundAnchor, which I'm more confident in.)

fun updateStarRealLocations() {
    stars.forEach { star ->
        star.vLoc.set(star.loc) // Stars start out with their pixel location in the image
        
        // First undo distortion
        undistort(star.vLoc)

        // Second undo rotation
        star.vLoc.rotateAroundAnchor(rotationAnchorNorthStar, rotationRadians)
    }
}

fun undistort(loc:Point2D) {
    val centerX = width / 2.0
    val centerY = height / 2.0
    val radiusNormalizationScale = 1.0 / Math.min(centerX, centerY) // I think necessary to have normalized coords to -1..1 for the formula to work.

    val distToCenterX = loc.x - centerX
    val distToCenterY = loc.y - centerY

    val radiusDistortedNormalized = radiusNormalizationScale *
            Math.sqrt(distToCenterX * distToCenterX + distToCenterY * distToCenterY)

    // poly3 model - The part I'm most concerned about.
    val radiusOriginalNormalized = distortionK1 * Math.pow(radiusDistortedNormalized, 3.0) +
            (1-distortionK1) * radiusDistortedNormalized

    loc.scaleAroundAnchor(Point2D(centerX, centerY), radiusOriginalNormalized/radiusDistortedNormalized)
}

[stub mandrel]

If you have something like Paint Shop Pro you can use the lens correction function to correct the images enough that normal stacking software will stack them.