3
$\begingroup$

I want to undistort 4032 x 3024 image with OpenCV model using only radial distortion. This model is described here.

The main issue is the speed. I have a function distort that takes all the execution time:

distort[{x_, y_}, {k1_, k2_, k3_, k4_, k5_, k6_}] :=
    With[{r2 = x^2 + y^2},
        With[{r4 = r2^2, r6 = r2^3},
            With[{s = (1 + k1 r2 + k2 r4 + k3 r6) / (1 + k4 r2 + k5 r4 + k6 r6)}, s {x, y}]]];

Distortion of one point on my laptop takes about 0.00002 seconds. So I have about 12M points and the total time is 12M * 0.00002 s ~ 4 minutes.

Any ideas how to speed up the computations?

P. S. A piece of code to test program:

{imageWidth, imageHeight} = {4032, 3024};
dictortionCoefficientsRadial = {0.0488079, -0.140193, 0.144845, 0, 0, 0};
cameraMatrix = {{3137.35, 0, 1990.53}, {0, 3143.89, 1514.04}, {0, 0, 1}};

AbsoluteTiming[distort[{1, 1}, dictortionCoefficientsRadial]]

pointsImage = Tuples[{Range[0, imageWidth - 1], Range[0, imageHeight - 1], {1}}];
pointsZ1 = pointsImage . Transpose[Inverse[cameraMatrix][[{1, 2}]]];
pointsZ1Distorted = Curry[distort][dictortionCoefficientsRadial] /@ pointsZ1;
$\endgroup$
2
  • 1
    $\begingroup$ Have you tried using ParallelMap instead of /@? As well, using Developer`ToPackedArray may help. $\endgroup$
    – Carl Lange
    Aug 7 '19 at 11:42
  • $\begingroup$ @CarlLange Yes, I did. ParallelMap[Curry[distort][dictortionCoefficientsRadial], Developer`ToPackedArray[widePointsZ1]] worked about 9 minutes. $\endgroup$
    – klimenkov
    Aug 7 '19 at 12:40
6
$\begingroup$

Compile is your friend in this case:

cf = Compile[{{X, _Real, 1}, {k, _Real, 1}},
 Block[{r2, s, x, y},
  x = X[[1]];
  y = X[[2]];
  r2 = x^2 + y^2;
  s = Divide[
    1. + r2 (k[[1]] + r2 (k[[2]] + k[[3]] r2)), 
    1. + r2 (k[[4]] + r2 (k[[5]] + k[[6]] r2))
    ];
  {s x, s y}],
 CompilationTarget -> "C",
 RuntimeAttributes -> {Listable},
 Parallelization -> True,
 RuntimeOptions -> "Speed"
 ];

Now you use it as follows:

{imageWidth, imageHeight} = {4032, 3024};
dictortionCoefficientsRadial = Developer`ToPackedArray[N[{0.0488079, -0.140193, 0.144845, 0, 0, 0}]];
cameraMatrix = Developer`ToPackedArray[N[{{3137.35, 0, 1990.53}, {0, 3143.89, 1514.04}, {0, 0, 1}}]];

pointsImage = RandomReal[{0, 1}, {imageWidth, imageHeight, 3}];
pointsZ1 = pointsImage.Transpose[Inverse[cameraMatrix][[{1, 2}]]];
pointsZ1Distorted2 = cf[pointsZ1, dictortionCoefficientsRadial]; // 
  AbsoluteTiming // First

0.729514

Done in under a second.

Due to the option RuntimeAttributes -> {Listable}, cf threads over matrices so that we need not flatten out the image anymore. So we have

Dimensions[pointsZ1Distorted2]

{4032, 3024, 2}

Edit

I put numerator and denominator of s into HornerForm. In principle, this should save us a few floating point multiplications, but it does not really make a difference here.

$\endgroup$
3
  • $\begingroup$ Why do you have both 1 in {{X, _Real, 1}, {k, _Real, 1}}? Should it be 2, 6? $\endgroup$
    – klimenkov
    Aug 7 '19 at 12:48
  • $\begingroup$ The 1 stands for the tensor rank of a tensor. 1 means, it is a vector. Size of the tensor arguments can be requested from within the function with Dimension. $\endgroup$ Aug 7 '19 at 13:09
  • $\begingroup$ This is actually something for which GPU programming might be very useful... $\endgroup$ Aug 7 '19 at 19:13

Your Answer

By clicking “Post Your Answer”, you agree to our terms of service, privacy policy and cookie policy

Not the answer you're looking for? Browse other questions tagged or ask your own question.