Correcting barrel lens distortion in an x,y series

Question

Edit: Ideally, I would like to use an identical algorithm to ImageMagick's -distort function, but implemented in R so as to transform raw x,y coordinates, rather than images. Using ImageMagick would also allow me to gain a visual proof of the effectiveness of the undistortion.

The question

I have an x,y series, produced from analysis of a video by a tracking software. The original images in the video have some barrel lens distortion that I would like to correct in order to improve the accuracy of the tracking coordinates.

I can perform what looks like a suitable correction using the distort operator in ImageMagick, e.g.:

convert input.jpg -distort barrel '0 -0.022 0' output.jpg

Now I know I could apply this correction to each frame in the video (before tracking), but this doesn't seem like the best option as I have dozens of videos, each of which consists of >7e4 frames. It seems like it would be much easier to apply the correction to the x,y coordinates themselves after tracking.

From the ImageMagick documentation, the barrel distortion equation is:

Rsrc = r * ( Ar^3 + Br^2 + C*r + D )

Where "r" is the destination radius and "Rsrc" is the source pixel to get the pixel color from. the radii are normalized so that radius = '1.0' for the half minimum width or height of the input image."

But I'm unaware of how to implement this in R in order to transform an x,y series. Could anyone offer any help? Thanks!

What I've tried so far

I've experimented with my own function, modifying a simple algorithm I found here, but this seems to introduce a greater barrel distortion if anything (and the polarity of which does not seem to be able to be reversed):

undistortFun <- function(X, Y, strength) {

    imWidth <- 640
    imHeight <- 480

    radius_u <- sqrt(imWidth^2 + imHeight^2) / strength

    normX <- X - (imWidth / 2)
    normY <- Y - (imHeight / 2)

    distance <- sqrt(normX^2 + normY^2)
    r <- distance / radius_u

    theta <- ifelse(r == 0, 1, atan(r) / r)

    newX <- (imWidth / 2) + theta * normX
    newY <- (imHeight / 2) + theta * normY

    return(data.frame(X = newX, Y = newY))

}

Similar questions

I found two similar questions, here and here, but these are concerned with undistorting images rather than raw x,y coordinates, and are implemented in Java and C++, which I'm not familiar with.

Answer 1

I've got a somewhat satisfactory solution by modifying the function provided in the question, using the ImageMagick distortion algorithm:

Rsrc = r * ( Ar^3 + Br^2 + C*r + D )

Like so:

undistortFun <- function(x, y, a, b, c, d = 1, imWidth = 640, imHeight = 480) {

    normX <- X - (imWidth / 2)
    normY <- Y - (imHeight / 2)

    radius_u <- sqrt(imWidth^2 + imHeight^2)

    r <- sqrt(normX^2 + normY^2) /  radius_u

    Rsrc <- r * (a*r^3 + b*r^2 + c*r + d)

    theta <- ifelse(Rsrc == 0, 1, 1 / atan(Rsrc) * Rsrc)

    newX <- (imWidth / 2) + theta * normX
    newY <- (imHeight / 2) + theta * normY

    return(data.frame(X = newX, Y = newY))
}

Although I'm uncertain about the calculation of theta, and specifying zero distortion (a=0, b=0, c=0) still leads to a transformation. Nevertheless, it appears to do what I wanted in these samples:

Plot of original x,y data:

Plot of adjusted x,y data (where: a =0, b= -0.02, c = 0):