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$
  • 1
    $\begingroup$ Have you tried using ParallelMap instead of /@? As well, using Developer`ToPackedArray may help. $\endgroup$ – Carl Lange Aug 7 at 11:42
  • $\begingroup$ @CarlLange Yes, I did. ParallelMap[Curry[distort][dictortionCoefficientsRadial], Developer`ToPackedArray[widePointsZ1]] worked about 9 minutes. $\endgroup$ – klimenkov Aug 7 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$
  • $\begingroup$ Why do you have both 1 in {{X, _Real, 1}, {k, _Real, 1}}? Should it be 2, 6? $\endgroup$ – klimenkov Aug 7 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$ – Henrik Schumacher Aug 7 at 13:09
  • $\begingroup$ This is actually something for which GPU programming might be very useful... $\endgroup$ – Henrik Schumacher Aug 7 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.