Adding median3d and dezinger3d filters into tomopy.misc.corr

[dkazanc]

in this commit:

• Median3d filter implemented in C language, accessible with from tomopy.misc.corr import median_filter3d for float32 and uint16 data types
• Dezinger3d filter implemented in C language, accessible with from tomopy.misc.corr import remove_outlier3d for float32 and uint16 data types
• unit tests for median and dezinger3d filters

The C functions are implemented using OpenMP API and a special care taken with TomoPy ncore multithreading option. C modules are also accepting ncore parameter and schedule computing resources accordingly. By default all cores are enabled.

[carterbox]

Thanks for your contribution! I look forward to seeing these functions merged.

[carterbox]

Wrote a small benchmark comparing this PR with scipy. Shows good improvement if you use enough cores. Note, this is for the 3D case only, I did not compare this PR with the existing 2D filtering in tomopy. Though, we should probably investgate if we can replace scipy with this PR in all cases.

import numpy as np
import timeit
import matplotlib.pyplot as plt

tomopy8_times = []
tomopy3_times = []
tomopy1_times = []
scipy_times = []
widths = np.unique(np.logspace(start=0, stop=6, base=2).astype(int))

for width in widths:

    times = timeit.repeat(
        number=1,
        stmt='tomopy.median_filter3d(x, size=7, ncore=1)',
        setup=f'''
import numpy as np
x = np.random.rand({width}, {width}, {width})
import tomopy
''',
    )
    tomopy1_times.append((np.mean(times), np.std(times)))

    times = timeit.repeat(
        number=1,
        stmt='tomopy.median_filter3d(x, size=7, ncore=3)',
        setup=f'''
import numpy as np
x = np.random.rand({width}, {width}, {width})
import tomopy
''',
    )
    tomopy3_times.append((np.mean(times), np.std(times)))

    times = timeit.repeat(
        number=1,
        stmt='tomopy.median_filter3d(x, size=7, ncore=8)',
        setup=f'''
import numpy as np
x = np.random.rand({width}, {width}, {width})
import tomopy
''',
    )
    tomopy8_times.append((np.mean(times), np.std(times)))

    times = timeit.repeat(
        number=1,
        stmt='scipy.ndimage.median_filter(x, size=7)',
        setup=f'''
import numpy as np
x = np.random.rand({width}, {width}, {width})
import scipy
    ''',
    )
    scipy_times.append((np.mean(times), np.std(times)))

plt.figure()
y, yerr = zip(*tomopy1_times)
plt.errorbar(
    x=widths,
    y=y,
    yerr=yerr,
)
y, yerr = zip(*tomopy3_times)
plt.errorbar(
    x=widths,
    y=y,
    yerr=yerr,
)
y, yerr = zip(*tomopy8_times)
plt.errorbar(
    x=widths,
    y=y,
    yerr=yerr,
)
y, yerr = zip(*scipy_times)
plt.errorbar(
    x=widths,
    y=y,
    yerr=yerr,
)
plt.legend(
    ['tomopy (ncore=1)', 'tomopy (ncore=3)', 'tomopy (ncore=8)', 'scipy'])
plt.title('Wall time of median_filter for cube of random float32')
plt.xlabel('Cube edge length [pixels]')
plt.ylabel('Wall time [s]')
plt.savefig('compare-median.svg')

[carterbox]

Please do one last check of your code. Note that I did remove some headers from the new C implementations because they seemed to be unused.

[carterbox]

Ope, I broke Windows.

[dkazanc]

Many thanks @carterbox for your help with this PR. It's been a good tutorial for us and will try to keep the suggested formatting and style in our future PRs. We've got few things that we would like to add to TomoPy.

Regarding the PR, all seems to behave for me and the multithreading performs as expected. I'm getting a good speed-up compare to scipy.ndimage.median_filter which does seems to run serially. Happy for it to be merged. Thanks!