1

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.

Community
  • 1
  • 1
jogall
  • 651
  • 6
  • 21
  • 1
    There's little to no published `R` code for optical correction algorithms AFAIK. But it would be easier to comment if you provided the source of the barrel-correction algorithm you refer to as "elsewhere." (Note that I'm failing to resist suggesting that, if your first attempt increases the apparent distortion, you should **reverse the polarity** :-) ) – Carl Witthoft Aug 06 '14 at 12:34
  • 2
    This question appears to be off-topic because it is about image processing methods developments and it's not really a specific programming question (since your code runs without error, it just doesn't produce the results you want.) I might check if there is [a different stackexchange site](http://stackexchange.com/sites) that might be a better fit. Or at least be clear on the details of the method you wish to implement. – MrFlick Aug 06 '14 at 13:30
  • @MrFlick I've edited the question to hopefully make it clearer. – jogall Aug 06 '14 at 19:53
  • This really doesn't help unless you can tell us exactly what algorithm ImageMagick's `-distort` function uses. – MrFlick Aug 06 '14 at 19:56
  • @CarlWitthoft I've added the source -- reversing the polarity was my first thought (believe it or not :-P ), but the author of the function indicates that the transformation cannot be inversed, and I'm afraid I don't understand what the function actually does well enough to work around it! – jogall Aug 06 '14 at 19:57
  • @MrFlick added algorithm to question: `Rsrc = r * ( A*r^3 + B*r^2 + C*r + D )` – jogall Aug 06 '14 at 20:08
  • Jogal, I meant just change the sign of the correction terms :-( – Carl Witthoft Aug 07 '14 at 12:24
  • @jogall Are X and Y that are in: normX <- X - (imWidth / 2) normY <- Y - (imHeight / 2) supposed to be vectors of the image width (1:640) and length (1:480) or a single number of the image width and length each? – SqueakyBeak Mar 14 '19 at 19:39

1 Answers1

2

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:

original

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

jogall
  • 651
  • 6
  • 21
  • Wonderful solution. In this case, are X and Y that are in: normX <- X - (imWidth / 2) normY <- Y - (imHeight / 2) supposed to be vectors of the image width (1:640) and length (1:480) or a single number of the image width and length each? – SqueakyBeak Mar 14 '19 at 15:42