opencv vs scipy: tiff rotation performance

[decarlof]

in #117 the tiff rotation function was replaced by scipy (from opencv). The new implementation works correctly except we are receiving reports of poor performance.

Using scipy for rotation:

Namespace(exp_log='20idd_mpe_arp25.yaml', scan_id='niv1', rec_type='full', rot_axis=None, rot_step=None, mblur=None, save_sino=None, use_half=None)
parameter file successfully loaded!!
Start reading files and rotate image 0.205 degree....
filname template is /scratch/step_0p100_10f_000001.tif
Done. 1676.0738 sec
Size of projection is (3601, 1590, 1740)
Normalization....

Using Opencv"

Namespace(exp_log='20idd_mpe_arp25.yaml', scan_id='niv1', rec_type='full', rot_axis=None, rot_step=None, mblur=None, save_sino=None, use_half=None)
parameter file successfully loaded!!
Start reading files and rotate image 0.205 degree....
filname template is /scratch/step_0p100_10f_000001.tif
Done. 27.3548 sec
Size of projection is (3601, 1590, 1740)
Normalization....

@samtygier-stfc did you get similar performance?

[samtygier-stfc]

Hi, yes I can reproduce that with 3d arrays. And it is coming from that change

import time
import numpy as np
import cv2
from scipy import ndimage

def timer(f, data, angle):
    time_start = time.monotonic()
    out = f(data, angle)
    time_end = time.monotonic()
    assert out.shape == data.shape
    return time_end - time_start

def scipy_rotate(arr, angle):
    return ndimage.rotate(arr, angle, reshape=False)

def cv_rotate(arr, angle):
    scaling = 1
    (h, w) = arr.shape[:2]
    (cX, cY) = (w // 2, h // 2)
    M = cv2.getRotationMatrix2D((cX, cY), angle, scaling)
    return cv2.warpAffine(arr, M, (w, h))

angle = 0.2
N = 400
data = np.random.default_rng().standard_normal((N,N,N), dtype=np.float32)

print(f"scipy_rotate: {timer(scipy_rotate, data, angle)}")
print(f"scipy_rotate: {timer(scipy_rotate, data, angle)}")
print(f"cv_rotate: {timer(cv_rotate, data, angle)}")
print(f"cv_rotate: {timer(cv_rotate, data, angle)}")

gives

scipy_rotate: 7.875956077012233
scipy_rotate: 7.7427476240554824
cv_rotate: 0.06871002994012088
cv_rotate: 0.04576018895022571

profiling shows most of the time is in affine_transform

INFO:__main__:   ncalls  tottime  percall  cumtime  percall filename:lineno(function)
INFO:__main__:        1    0.000    0.000    7.719    7.719 /home/sam/git/dxchange/test_122_rot.py:12(scipy_rotate)
INFO:__main__:        1    0.002    0.002    7.719    7.719 /home/sam/miniforge3/envs/dxchange-dev/lib/python3.9/site-packages/scipy/ndimage/_interpolation.py:821(rotate)
INFO:__main__:      400    0.005    0.000    7.717    0.019 /home/sam/miniforge3/envs/dxchange-dev/lib/python3.9/site-packages/scipy/ndimage/_interpolation.py:461(affine_transform)
INFO:__main__:      400    5.924    0.015    5.924    0.015 {built-in method scipy.ndimage._nd_image.geometric_transform}
INFO:__main__:      400    0.002    0.000    1.778    0.004 /home/sam/miniforge3/envs/dxchange-dev/lib/python3.9/site-packages/scipy/ndimage/_interpolation.py:135(spline_filter)
INFO:__main__:      800    0.002    0.000    1.762    0.002 /home/sam/miniforge3/envs/dxchange-dev/lib/python3.9/site-packages/scipy/ndimage/_interpolation.py:47(spline_filter1d)
INFO:__main__:      800    1.756    0.002    1.756    0.002 {built-in method scipy.ndimage._nd_image.spline_filter1d}
INFO:__main__:     1601    0.002    0.000    0.016    0.000 /home/sam/miniforge3/envs/dxchange-dev/lib/python3.9/site-packages/scipy/ndimage/_ni_support.py:73(_get_output)
INFO:__main__:      401    0.013    0.000    0.013    0.000 {built-in method numpy.zeros}
INFO:__main__:     2402    0.001    0.000    0.004    0.000 /home/sam/miniforge3/envs/dxchange-dev/lib/python3.9/site-packages/numpy/core/_asarray.py:14(asarray)
INFO:__main__:     1601    0.001    0.000    0.004    0.000 <__array_function__ internals>:2(iscomplexobj)
INFO:__main__:      400    0.002    0.000    0.003    0.000 /home/sam/miniforge3/envs/dxchange-dev/lib/python3.9/site-packages/scipy/ndimage/_ni_support.py:57(_normalize_sequence)
INFO:__main__:     2403    0.003    0.000    0.003    0.000 {built-in method numpy.array}
INFO:__main__:     1601    0.001    0.000    0.002    0.000 {built-in method numpy.core._multiarray_umath.implement_array_function}

I also found scipy/scipy#19369

I suspect some of the difference is that openCV is using linear interpolation and only covers cases up to 3 dimensions. Scipy is using 3rd order spline (cubic) and can handle higher dimensions.

I get a very small improvement if I pass order=1 to scipy. I suspect there is not special case so it still calling the affine_transform.

I tried flags = cv2.INTER_CUBIC but get an error because that is not implemented for more than 4 channels (I assume its treating the 3rd dimension as image channels)

[samtygier-stfc]

Given the size of the slowdown my patch could be rolled back.

Or the code could do something like use opencv for rotation if it is available, and either fall back to scipy or give an error that the rotation requires opencv.

We could also special case rotations of multiples of 90 degrees on square arrays, which do not need interpolation. It would not help here, but is probably a common case.