# Speed up a distortion of 4k image

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}};

pointsImage = Tuples[{Range[0, imageWidth - 1], Range[0, imageHeight - 1], {1}}];
pointsZ1 = pointsImage . Transpose[Inverse[cameraMatrix][[{1, 2}]]];

• Have you tried using ParallelMap instead of /@? As well, using DeveloperToPackedArray may help. – Carl Lange Aug 7 '19 at 11:42
• @CarlLange Yes, I did. ParallelMap[Curry[distort][dictortionCoefficientsRadial], DeveloperToPackedArray[widePointsZ1]] worked about 9 minutes. – klimenkov Aug 7 '19 at 12:40

Compile is your friend in this case:

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


Now you use it as follows:

{imageWidth, imageHeight} = {4032, 3024};
dictortionCoefficientsRadial = DeveloperToPackedArray[N[{0.0488079, -0.140193, 0.144845, 0, 0, 0}]];
cameraMatrix = DeveloperToPackedArray[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}]]];
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.

• Why do you have both 1 in {{X, _Real, 1}, {k, _Real, 1}}? Should it be 2, 6? – klimenkov Aug 7 '19 at 12:48
• 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. – Henrik Schumacher Aug 7 '19 at 13:09
• This is actually something for which GPU programming might be very useful... – Henrik Schumacher Aug 7 '19 at 19:13