I'm working on a nice approximation for the Bilateral Filter.
I have a working code which runs pretty fast yet still I think much can be improved.
The code (C
Code, compiles with C99
) is given by (See code on Compiler Explorer):
#define _USE_MATH_DEFINES
#include <stdio.h>
#include <stdlib.h>
#include <math.h>
#include <memory.h>
#include <omp.h>
#define OFF 0
#define ON 1
#include <immintrin.h> // AVX
#define SSE_STRIDE 4
#define SSE_ALIGNMENT 16
#define AVX_STRIDE 8
#define AVX_ALIGNMENT 32
#define M_PIf (float)(M_PI)
void ImageConvolutionGaussianKernel(float* mO, float* mI, float* mTmp, int numRows, int numCols, float gaussianStd, int stdToRadiusFactor);
void InitOmegaArrays(float* mCOmega, float* mSOmega, float* mI, int numRows, int numCols, float paramOmega);
void UpdateArrays(float* mO, float* mZ, float* mC, float* mS, float* mCFiltered, float* mSFiltered, int numRows, int numCols, int iterationIdx, float paramD);
void InitArraysSC(float* mC, float* mS, float* mCOmega, float* mSOmega, int numRows, int numCols);
void UpdateArraysSC(float* mC, float* mS, float* mCOmega, float* mSOmega, int numRows, int numCols);
void UpdateOutput(float* mO, float* mZ, float* mI, int numRows, int numCols, float rangeStd, float paramL);
void BilateralFilterFastCompressive(float* mO, float* mI, int numRows, int numCols, float spatialStd, float rangeStd, int paramK)
{
int ii, paramN;
float paramL, paramTau, *vParamD, *mZ, *mT, paramOmega, *mCOmega, *mSOmega, *mC, *mS, *mCFiltered, *mSFiltered;
mZ = (float*)_mm_malloc(numRows * numCols * sizeof(float), AVX_ALIGNMENT); // Should be initialized to Zero
mT = (float*)_mm_malloc(numRows * numCols * sizeof(float), AVX_ALIGNMENT); // Buffer
mC = (float*)_mm_malloc(numRows * numCols * sizeof(float), AVX_ALIGNMENT);
mS = (float*)_mm_malloc(numRows * numCols * sizeof(float), AVX_ALIGNMENT);
mCOmega = (float*)_mm_malloc(numRows * numCols * sizeof(float), AVX_ALIGNMENT);
mSOmega = (float*)_mm_malloc(numRows * numCols * sizeof(float), AVX_ALIGNMENT);
mCFiltered = (float*)_mm_malloc(numRows * numCols * sizeof(float), AVX_ALIGNMENT);
mSFiltered = (float*)_mm_malloc(numRows * numCols * sizeof(float), AVX_ALIGNMENT);
memset(mZ, 0.0f, numRows * numCols * sizeof(float));
// Init Parameters
paramL = paramK * rangeStd;
paramTau = paramK / M_PIf;
paramN = ceilf((paramK * paramK) / M_PIf);
paramOmega = M_PIf / paramL;
vParamD = (float*)_mm_malloc(paramN * sizeof(float), AVX_ALIGNMENT);
for (ii = 1; ii <= paramN; ii++)
{
vParamD[ii - 1] = 2 * expf(-(ii * ii) / (2 * paramTau * paramTau));
}
InitOmegaArrays(mCOmega, mSOmega, mI, numRows, numCols, paramOmega);
// Iteration Number 1
ii = 1;
ImageConvolutionGaussianKernel(mCFiltered, mCOmega, mT, numRows, numCols, spatialStd, paramK);
ImageConvolutionGaussianKernel(mSFiltered, mSOmega, mT, numRows, numCols, spatialStd, paramK);
UpdateArrays(mO, mZ, mCOmega, mSOmega, mCFiltered, mSFiltered, numRows, numCols, ii, vParamD[ii - 1]);
// Iteration Number 2
ii = 2;
InitArraysSC(mC, mS, mCOmega, mSOmega, numRows, numCols);
ImageConvolutionGaussianKernel(mCFiltered, mC, mT, numRows, numCols, spatialStd, paramK);
ImageConvolutionGaussianKernel(mSFiltered, mS, mT, numRows, numCols, spatialStd, paramK);
UpdateArrays(mO, mZ, mC, mS, mCFiltered, mSFiltered, numRows, numCols, ii, vParamD[ii - 1]);
for (ii = 3; ii <= paramN; ii++)
{
UpdateArraysSC(mC, mS, mCOmega, mSOmega, numRows, numCols);
ImageConvolutionGaussianKernel(mCFiltered, mC, mT, numRows, numCols, spatialStd, paramK);
ImageConvolutionGaussianKernel(mSFiltered, mS, mT, numRows, numCols, spatialStd, paramK);
UpdateArrays(mO, mZ, mC, mS, mCFiltered, mSFiltered, numRows, numCols, ii, vParamD[ii - 1]);
}
UpdateOutput(mO, mZ, mI, numRows, numCols, rangeStd, paramL);
_mm_free(mZ);
_mm_free(mT);
_mm_free(mC);
_mm_free(mS);
_mm_free(mCOmega);
_mm_free(mSOmega);
_mm_free(mCFiltered);
_mm_free(mSFiltered);
_mm_free(vParamD);
}
// Auxiliary Functions
void InitOmegaArrays(float* mCOmega, float* mSOmega, float* mI, int numRows, int numCols, float paramOmega) {
int ii;
for (ii = 0; ii < numRows * numCols; ii++)
{
mCOmega[ii] = cosf(paramOmega * mI[ii]);
mSOmega[ii] = sinf(paramOmega * mI[ii]);
}
}
void UpdateArrays(float* mO, float* mZ, float* mC, float* mS, float* mCFiltered, float* mSFiltered, int numRows, int numCols, int iterationIdx, float paramD) {
int ii;
for (ii = 0; ii < numRows * numCols; ii++)
{
mO[ii] += (iterationIdx * paramD) * (mC[ii] * mSFiltered[ii] - mS[ii] * mCFiltered[ii]);
mZ[ii] += paramD * (mC[ii] * mCFiltered[ii] + mS[ii] * mSFiltered[ii]);
}
}
void InitArraysSC(float* mC, float* mS, float* mCOmega, float* mSOmega, int numRows, int numCols) {
int ii;
for (ii = 0; ii < numRows * numCols; ii++)
{
mS[ii] = 2.0f * mCOmega[ii] * mSOmega[ii];
mC[ii] = 2.0f * mCOmega[ii] * mCOmega[ii] - 1.0f;
}
}
void UpdateArraysSC(float* mC, float* mS, float* mCOmega, float* mSOmega, int numRows, int numCols) {
int ii;
float varTmp;
for (ii = 0; ii < numRows * numCols; ii++)
{
varTmp = mC[ii] * mSOmega[ii] + mS[ii] * mCOmega[ii];
mC[ii] = mC[ii] * mCOmega[ii] - mS[ii] * mSOmega[ii];
mS[ii] = varTmp;
}
}
void UpdateOutput(float* mO, float* mZ, float* mI, int numRows, int numCols, float rangeStd, float paramL) {
int ii;
float outFactor;
outFactor = (M_PIf * rangeStd * rangeStd) / paramL;
for (ii = 0; ii < numRows * numCols; ii++)
{
mO[ii] = mI[ii] + (outFactor * (mO[ii] / (1.0f + mZ[ii])));
}
}
Basically few iterations of Gaussian Blur and Element Wise operations.
I'm testing the code on image of size 8000 x 8000
with spatialStd = 5
, rangeStd = 10 / 255
and paramK = 5
.
This set of parameters means the number of iterations (paramN
in the code) is 8 (Gaussian Blur is done twice per iteration).
My Gaussian Blur implementation takes ~0.18 [Sec] per iteration in the settings above when tested independently.
The issues I'm having with the code:
- Per iteration, it seems the time all Element Wise operations takes more than the Gaussian Blur. It seems to me something isn't efficient with the code.
- Overhead - If I run the code skipping the Gaussian Blur it takes ~1.8 [Sec]. Thinking that 16 x 0.18 = 2.88 [Sec] (16 iterations of the Gaussian Blur) means I should expect run time of ~5 [Sec]. In practice I get ~ 7 [Sec]. It means there is huge overhead somewhere there.
What I've tried:
- Writing all small function using AVX2 SIMD intrinsics. Yet it seems I gain only 3-5% over the compiler (I use Intel Compiler).
- Using the
#pragma vector aligned
to force compiler to vectorize the code and assume no aliasing and no alignment issues. It yields results which are 3-5% slower than the hand tuned I did (See 1). - Disabling OpenMP (Didn't improve).
- Various compilers (MSVC, GCC 7.3 / 8.1, Intel Compiler 2018 [Was best]). Results above are the best achieved (Intel Compiler). I tried
-Ofast
and-O3
on GCC. UsingO3
on Intel. FP Precision is set to fast.
I'd be happy to get some feedback on that.
Remark
The code is part of a work at University and will be published on GitHub once it is ready.
UpdtaeOutput
has a typo in the function name.ImageConvolutionGaussianKernel
doesn't seem to be defined anywhere, so I couldn't look at how the whole thing compiles. \$\endgroup\$