openCVtrilateralFilter.cpp

#include "openCVtrilateralFilter.h"
#include <math.h>
#include "opencv2/opencv.hpp"
#include "opencv/highgui.h"
#include <vector>
using std::vector;

bool OpenCVtrilateralFilter::trilateralFilter(
    cv::Mat& inputImg,
    cv::Mat& outputImg,
    const float sigmaC, 
    const float epsilon,
    const int filterType
    )
{

    // Check sizes
    if (inputImg.cols != outputImg.cols
        || inputImg.rows != outputImg.rows
        || inputImg.depth() != outputImg.depth()
        || inputImg.channels() != outputImg.channels()
        || inputImg.channels() != 1)
    {
        return false;
    }

    if (inputImg.depth() != IPL_DEPTH_32F)
        return false;

    // Adaptive Neighbourhood pixelwise image
    cv::Mat adaptiveNeighbourhood = cv::Mat(inputImg.rows, inputImg.cols, inputImg.type(),cv::Scalar(1));
    
    // x and y gradients
    cv::Mat xGradient = cv::Mat(inputImg.rows, inputImg.cols, inputImg.type(),cv::Scalar(1));
    cv::Mat yGradient = cv::Mat(inputImg.rows, inputImg.cols, inputImg.type(),cv::Scalar(1));

    // Smoothed x and y gradients
    cv::Mat xGradientSmooth = cv::Mat(inputImg.rows, inputImg.cols, inputImg.type(),cv::Scalar(1));
    cv::Mat yGradientSmooth = cv::Mat(inputImg.rows, inputImg.cols, inputImg.type(),cv::Scalar(1));
    
    // Gradient magnitude
    cv::Mat gradientMagnitude = cv::Mat(inputImg.rows, inputImg.cols, inputImg.type(),cv::Scalar(1));

    // domain variance for the two filters: sigmaC
    // range variance of the two filters: sigmaR
    // beta = emperical value
    float sigmaR; 
    const float beta = 0.15f; //Always set between 0.1 and 0.2

    //Computes X and Y gradients of the input image
    computeGradients(inputImg, xGradient, yGradient ); 

    // Computes the gradient magnitude
    computeMagnitude(xGradient, yGradient, gradientMagnitude);

    // Get min and max gradient for sigmaR
    double minGrad, maxGrad;
    minMaxLoc(gradientMagnitude, &minGrad, &maxGrad);
    sigmaR = beta * ( (float) maxGrad - (float) minGrad );

    // Level Max and Maximum LUT values
    int levelMax, maxLUT;

    // If using LUT
    if (filterType == TYPE_LUT)
    {
        // maximum level = log2(xsize) or log2(ysize)
        levelMax = log2( __min( inputImg.cols, inputImg.rows ) );
        maxLUT = pow( 2.f, levelMax-1) + 1;
    }
    else
    {
        // Using fast-trilateral
        // Find threshold for truncation
        maxLUT = int( sigmaC * sqrt( fabs( log(epsilon) ) ) ) + 1;
        // Calculate maximum level
        levelMax = log2(2 * maxLUT + 1, true);
    }

    // Calculate Adaptive Neighbourhood
    setAdaptiveNeighbourHood(gradientMagnitude, sigmaR, levelMax, adaptiveNeighbourhood);

    // Bilaterally filter the X and Y gradients of the input image
    // to produce xGradientSmooth and yGradientSmooth.

    if (filterType == TYPE_LUT)
        BilateralGradientFilterLUT(
            xGradient,yGradient, gradientMagnitude, sigmaC, sigmaR,
            xGradientSmooth, yGradientSmooth);
    else
        BilateralGradientFilter(
        xGradient, yGradient, gradientMagnitude, sigmaC, sigmaR, epsilon, 
        xGradientSmooth, yGradientSmooth );
    
    // Performs bilateral filter on the detail signal 

    if (filterType == TYPE_LUT)
        DetailBilateralFilterLUT(
            inputImg, adaptiveNeighbourhood, xGradientSmooth, yGradientSmooth, 
            sigmaC, sigmaR, maxLUT, outputImg);
    else
        DetailBilateralFilter(
            inputImg, adaptiveNeighbourhood, xGradientSmooth, yGradientSmooth, 
            sigmaC, sigmaR, maxLUT, epsilon, outputImg);

    return true;
}

// Computes X and Y gradients of the input image
void OpenCVtrilateralFilter::computeGradients(
    cv::Mat& inputImg, cv::Mat& xGradient, cv::Mat& yGradient )
{

    // Set up convolution kernels for forward differences
    float kernel[] = { -1, 1 };
    cv::Mat xKernel = cv::Mat(1,2,CV_32FC1,kernel);
    cv::Mat yKernel = cv::Mat(2,1,CV_32FC1,kernel);
    cv::Point anchorPt = cv::Point(0,0);

    cv::filter2D(inputImg,xGradient, CV_32FC1, xKernel, anchorPt);  // x gradient
    cv::filter2D(inputImg,yGradient,CV_32FC1, yKernel, anchorPt);   // y gradient
}

// Computes the magnitude of the gradients
void OpenCVtrilateralFilter::computeMagnitude(
        cv::Mat& xGradient, cv::Mat& yGradient,
        cv::Mat& gradientMagnitude)
{
    gradientMagnitude = cv::norm(xGradient-yGradient);
    /*
    float* xGrad = (float*) xGradient.data;
    float* yGrad = (float*) yGradient.data;
    float* gradMag = (float*) gradientMagnitude.data;

    for (int x = 0; x < xGradient.cols * xGradient.rows; x++)
        gradMag[x] = cv::norm(xGrad[x], yGrad[x]);
    */
}

// Find the adaptive neighbourhood for image
void OpenCVtrilateralFilter::setAdaptiveNeighbourHood(
    cv::Mat& gradientMagnitude, const float sigmaR, const int maxLevel,
    cv::Mat& adaptiveNeighbourhood  )
{
    // Image stacks for max and min neighbourhoods
    vector<cv::Mat> minGradientStack, maxGradientStack;
    cv::Size imgSize = gradientMagnitude.size();

    // Create image stack
    for(int i = 0 ; i < maxLevel ; i++)
        minGradientStack.push_back( cv::Mat(imgSize,gradientMagnitude.depth(),1) );
    for(int i = 0 ; i < maxLevel ; i++)
        maxGradientStack.push_back( cv::Mat(imgSize,gradientMagnitude.depth(),1) );

    // Build the min-max stack
    buildMinMaxImageStack(gradientMagnitude, minGradientStack, maxGradientStack );

    // Set up image data references
    cv::Mat minImg, maxImg;
    cv::Mat magImg = cv::Mat(gradientMagnitude);
    cv::Mat adpImg = cv::Mat(adaptiveNeighbourhood);

    for(int y = 0; y < imgSize.width; y++)
    {
        for(int x = 0; x < imgSize.height; x++)
        {
            int lev;
            const float upperThreshold = magImg.ptr(y)[x] + sigmaR;
            const float lowerThreshold = magImg.ptr(y)[x] - sigmaR;

            // Compute the adaptive neighbourhood based on the similarity of
            // the neighborhood gradients
            for(lev = 0; lev < maxLevel; lev++) 
            {
                minImg = minGradientStack[lev];
                maxImg = maxGradientStack[lev];
                if ( maxImg.ptr(y)[x] > upperThreshold || minImg.ptr(y)[x] < lowerThreshold )
                    break;
            }

            // Sets the (half) size of the adaptive region
            // i.e., floor( ( pow(2.f, lev) + 1 ) / 2.f )
            adpImg.ptr(y)[x] = pow(2.f, lev-1);
        }
    }
}


// Building the Min-Max Stack of Image Gradients
void OpenCVtrilateralFilter::buildMinMaxImageStack(
    cv::Mat& gradientMagnitude,
    vector<cv::Mat> minStack, vector<cv::Mat> maxStack )
{
    const int imgWidth = gradientMagnitude.cols;
    const int imgHeight = gradientMagnitude.rows;

    // Set up image data references
    cv::Mat minImg1, maxImg1, minImg2, maxImg2;
    cv::Mat magImg = cv::Mat(gradientMagnitude);

    // Set up the bottom level of the pyramid
    minImg1 = minStack[0];
    maxImg1 = maxStack[0];
    
    // Loop through image setting up bottom stack
    for(int y = 0; y < imgHeight ; y++)
    {
        for (int x = 0; x < imgWidth; x++)
        {
            float outMin = 1e12;
            float outMax = -1e12;
            
            // Loop through local neighbourhood
            // To find maximum and minimum values
            for(int n = __max(y-1,0) ; n < __min(y+2, imgHeight); n++) 
            {
                for(int m=__max(x-1,0); m < __min(x+2, imgWidth) ; m++) 
                {
                    outMin = __min(magImg.ptr(n)[m], outMin);
                    outMax = __max(magImg.ptr(n)[m], outMax);
                }
            }
            minImg1.ptr(y)[x] = outMin;
            maxImg1.ptr(y)[x] = outMax;
        }
    }

    // Loop through image stack
    for (int i = 1 ; i < minStack.size(); i++)
    {
        // Lower level
        minImg1 = minStack[i-1];
        maxImg1 = maxStack[i-1];

        // Current level
        minImg2 = minStack[i];
        maxImg2 = maxStack[i];

        for(int y = 0; y < imgHeight ; y++)
        {
            for (int x = 0; x < imgWidth; x++)
            {
                float outMin = 1e12;
                float outMax = -1e12;

                // Loop through local neighbourhood
                // To find maximum and minimum values
                for(int n = __max(y-1,0) ; n < __min(y+2, imgHeight); n++) 
                {
                    for(int m =__max(x-1,0); m < __min(x+2, imgWidth) ; m++) 
                    {
                        outMin = __min(minImg1.ptr(n)[m], outMin);
                        outMax = __max(maxImg1.ptr(n)[m], outMax);
                    }
                }
                minImg2.ptr(y)[x] = outMin;
                maxImg2.ptr(y)[x] = outMax;
            }
        }
    }
}

// =======================================================================
// Specific for fast version of Trilateral Fitler
// =======================================================================


//  Bilaterally filters the X and Y gradients of the input image.
//  To produce smoothed x and y gradients

void OpenCVtrilateralFilter::BilateralGradientFilter(
        cv::Mat& xGradient, cv::Mat& yGradient,
        cv::Mat& gradientMagnitude,
        const float sigmaC, const float sigmaR, const float epsilon,
        cv::Mat& xGradientSmooth, cv::Mat& yGradientSmooth  )
{
    // Get image size
    const int imgWidth  = xGradient.cols;
    const int imgHeight = xGradient.rows;

    // Constants used for domain / range calculations
    const float domainConst = -2.f * sigmaC * sigmaC;
    const float rangeConst =  -2.f * sigmaR * sigmaR;

    // Compute the weight for the domain filter (domainWeight). 
    // The domain filter is a Gaussian lowpass filter
    const int halfSize = int(sigmaC / 2.f);
    cv::Mat domainWeightLUT =
        cv::Mat(cvSize(halfSize+1,halfSize+1),xGradient.depth(),1);
    
    // Memory reference
    cv::Mat domainWeight = cv::Mat(domainWeightLUT);

    for (int y = 0; y < domainWeightLUT.rows ; y++)
    {
        for (int x = 0; x < domainWeightLUT.cols ; x++)
        {
            // weight for the domain filter (domainWeight)
            const float diff = (float) (x*x+y*y);
            domainWeight.ptr(y)[x] = (float) exp( diff / domainConst );
        }
    }

    // Memory referencing
    cv::Mat xImg = cv::Mat(xGradient);
    cv::Mat yImg= cv::Mat(yGradient);
    cv::Mat xSmoothImg = cv::Mat(xGradientSmooth);
    cv::Mat ySmoothImg= cv::Mat(yGradientSmooth);
    cv::Mat magImg= cv::Mat(gradientMagnitude);
    

    // Loop through image
    for(int y = 0; y < imgHeight ; y++) 
    {
        for(int x = 0; x < imgWidth ; x++)  
        {       
            double normFactor = 0.f;
            double tmpX = 0.f;
            double tmpY = 0.f;

            // Calculate Middle Pixel Normalised-gradient
            const float g2 = magImg.ptr(y)[x];

            // Loop through local neighbourhood
            for (int n = -halfSize; n <= halfSize; n++) 
            { 
                for(int m = -halfSize; m <= halfSize; m++) 
                {   
                    //Compute the weight for the domain filter (domainWeight). 
                    const float dWeight = domainWeight.ptr(abs(n))[abs(m)];

                    // Only perform calculation if weight above zero
                    if( dWeight < epsilon ) continue;
                
                    // Only perform calculationg if within bounds
                    const int localX = x + m;
                    if (localX < 0) continue;
                    if (localX >= imgWidth) continue;

                    const int localY = y + n;
                    if (localY < 0) continue;
                    if (localY >= imgHeight) continue;

                    // Calculate Local Normalised Gradient
                    const float g1 = magImg.ptr(localY)[localX];

                    //Compute the gradient difference between a pixel and its neighborhood pixel 
                    const float gradDiffSq =  (float) pow(g1 - g2, 2);
                    
                    // Compute the weight for the range filter (rangeWeight). The range filter
                    // is a Gaussian filter defined by the difference in gradient magnitude.
                    const float rangeWeight = (float) exp( gradDiffSq / rangeConst );   
                    
                    // Only compute if less than epsilon
                    if (rangeWeight < epsilon) continue;

                    tmpX += xImg.ptr(localY)[localX] * dWeight * rangeWeight;
                    tmpY += yImg.ptr(localY)[localX] * dWeight * rangeWeight;

                    // Bilateral filter normalized by normFactor
                    normFactor += dWeight * rangeWeight;
                }
            }

            // Set smoothed image to normalised value
            xSmoothImg.ptr(y)[x] = tmpX / normFactor;
            ySmoothImg.ptr(y)[x] = tmpY / normFactor;
        }   
    }
}


// Filters the detail signal and computes the output (2nd filtering pass for trilateral filter).
void OpenCVtrilateralFilter::DetailBilateralFilter(
        cv::Mat& inputImage, cv::Mat& adaptiveRegion,
        cv::Mat& xGradientSmooth, cv::Mat& yGradientSmooth,
        const float sigmaC, const float sigmaR,
        const int maxDomainSize, const float epsilon,
        cv::Mat& outputImg)
{

    // Get image size
    const int imgWidth  = inputImage.cols;
    const int imgHeight = inputImage.rows;

    // Create constants used throughout code
    const double domainConst = -2.f * sigmaC * sigmaC;
    const double rangeConst =  -2.f * sigmaR * sigmaR;

    // Memory referencing
    cv::Mat xImg = cv::Mat(xGradientSmooth);
    cv::Mat yImg = cv::Mat(yGradientSmooth);
    cv::Mat srcImg = cv::Mat(inputImage);
    cv::Mat outImg = cv::Mat(outputImg);
    cv::Mat adpImg = cv::Mat(adaptiveRegion);

    // Define Look Up Table for weightings
    
    // Compute the weight for the domain filter (domainWeight). The domain filter
    // is a Gaussian lowpass filter
    cv::Mat domainWeightLUT =
        cv::Mat( cvSize(maxDomainSize+1,maxDomainSize+1), inputImage.depth(), 1);
        
    // Memory reference
    cv::Mat domainWeight = cv::Mat(domainWeightLUT);

    for (int y = 0; y < domainWeightLUT.rows ; y++)
    {
        for (int x = 0; x < domainWeightLUT.cols ; x++)
        {
            // weight for the domain filter (domainWeight)
            const float diff = (float) (x*x+y*y);
            domainWeight.ptr(y)[x] = (float) exp( diff / domainConst );
        }
    }

    for(int y = 0; y < imgHeight; y++) 
    {   
        for(int x = 0; x < imgWidth; x++) 
        { 
            double normFactor = 0.f;
            double tmp = 0.f;

            // filter window width is calculated from adaptive neighbourhood
            // halfsize is half of the filter window width (or maximum LUT value)
            const int halfSize = (int) __min( adpImg.ptr(y)[x], maxDomainSize );

            // Coefficients defining the centerplane is calculated
            // from the smoothed image gradients
            // Plane Equation is z = coeffA.x + coeffB.y + coeffC
            // coeffA = dI/dx, coeffB = dI/dy, coeffC = I at center pixel of the filter kernel
            const float coeffA = xImg.ptr(y)[x];
            const float coeffB = yImg.ptr(y)[x];
            const float coeffC = srcImg.ptr(y)[x];

            for (int n = -halfSize; n <= halfSize; n++) // y scan line
            {
                for(int m = -halfSize; m <= halfSize; m++) // x scan line
                {
                    // Get domain weight from LUT
                    const float dWeight = domainWeight.ptr(abs(n))[abs(m)];

                    // Only perform calculation if weight above zero
                    if( dWeight < epsilon ) continue;
                
                    // Only perform calculationg if within bounds
                    const int localX = x + m;
                    if (localX < 0) continue;
                    if (localX >= imgWidth) continue;

                    const int localY = y + n;
                    if (localY < 0) continue;
                    if (localY >= imgHeight) continue;

                    // Compute the detail signal (detail) based on the difference between a 
                    // neighborhood pixel and the centerplane passing through the center-pixel 
                    // of the filter window. 
                    const float detail = 
                        srcImg.ptr(localY)[localX] - coeffA * float(m) - coeffB * float(n) - coeffC;

                    // Compute the weight for the range filter (rangeWeight). The range filter
                    // is a Gaussian filter defined by the detail signal.
                    const float rangeWeight = exp( pow(detail,2) / rangeConst );    

                    if ( dWeight < epsilon ) continue;

                    tmp += detail * dWeight * rangeWeight;

                    //Detail Bilateral filter normalized by normFactor (eq. 9, Section 3.1)
                    normFactor += dWeight * rangeWeight;
                }
            }

            // Write result to output image
            outImg.ptr(y)[x] = tmp / normFactor + coeffC;
        }   
    }   
}

// =======================================================================
// Specific for LUT version of Trilateral Fitler
// =======================================================================

//  Bilaterally filters the X and Y gradients of the input image.
//  To produce smoothed x and y gradients

void OpenCVtrilateralFilter::BilateralGradientFilterLUT(
        cv::Mat& xGradient, cv::Mat& yGradient,
        cv::Mat& gradientMagnitude,
        const float sigmaC, const float sigmaR,
        cv::Mat& xGradientSmooth, cv::Mat& yGradientSmooth  )
{
    // Get image size
    const int imgWidth  = xGradient.cols;
    const int imgHeight = xGradient.rows;

    // Constants used for domain / range calculations
    const float domainConst = -2.f * sigmaC * sigmaC;
    const float rangeConst =  -2.f * sigmaR * sigmaR;

    // Compute the weight for the domain filter (domainWeight). 
    // The domain filter is a Gaussian lowpass filter
    const int halfSize = int(sigmaC - 1.f / 2.f);
    cv::Mat domainWeightLUT = cv::Mat(cv::Size(halfSize+1,halfSize+1), xGradient.depth(),1);
    
    // Memory reference
    cv::Mat domainWeight = cv::Mat(domainWeightLUT);

    for (int y = 0; y < domainWeightLUT.rows ; y++)
    {
        for (int x = 0; x < domainWeightLUT.cols ; x++)
        {
            // weight for the domain filter (domainWeight)
            const float diff = (float) (x*x+y*y);
            domainWeight.ptr(y)[x] = (float) exp( diff / domainConst );
        }
    }

    // Memory referencing
    cv::Mat xImg = cv::Mat(xGradient);
    cv::Mat yImg = cv::Mat(yGradient);
    cv::Mat xSmoothImg = cv::Mat(xGradientSmooth);
    cv::Mat ySmoothImg = cv::Mat(yGradientSmooth);
    cv::Mat magImg = cv::Mat(gradientMagnitude);
    

    // Loop through image
    for(int y = 0; y < imgHeight ; y++) 
    {
        for(int x = 0; x < imgWidth ; x++)  
        {       
            double normFactor = 0.f;
            double tmpX = 0.f;
            double tmpY = 0.f;

            // Calculate Middle Pixel Normalised-gradient
            const float g2 = magImg.ptr(y)[x];

            // Loop through local neighbourhood
            for (int n = -halfSize; n <= halfSize; n++) 
            { 
                for(int m = -halfSize; m <= halfSize; m++) 
                {   
                    //Compute the weight for the domain filter (domainWeight). 
                    const float dWeight = domainWeight.ptr(abs(n))[abs(m)];
                
                    // Only perform calculationg if within bounds
                    const int localX = x + m;
                    if (localX < 0) continue;
                    if (localX >= imgWidth) continue;

                    const int localY = y + n;
                    if (localY < 0) continue;
                    if (localY >= imgHeight) continue;

                    // Calculate Local Normalised Gradient
                    const float g1 = magImg.ptr(localY)[localX];

                    // Compute the gradient difference between a pixel and its neighborhood pixel 
                    const float gradDiffSq =  (float) pow(g1 - g2, 2);
                    
                    // Compute the weight for the range filter (rangeWeight). The range filter
                    // is a Gaussian filter defined by the difference in gradient magnitude.
                    const float rangeWeight = (float) exp( gradDiffSq / rangeConst );   

                    tmpX += xImg.ptr(localY)[localX] * dWeight * rangeWeight;
                    tmpY += yImg.ptr(localY)[localX] * dWeight * rangeWeight;

                    // Bilateral filter normalized by normFactor
                    normFactor += dWeight * rangeWeight;
                }
            }

            // Set smoothed image to normalised value
            xSmoothImg.ptr(y)[x] = tmpX / normFactor;
            ySmoothImg.ptr(y)[x] = tmpY / normFactor;
        }   
    }
}


// Filters the detail signal and computes the output (2nd filtering pass for trilateral filter).
void OpenCVtrilateralFilter::DetailBilateralFilterLUT(
        cv::Mat& inputImage, cv::Mat& adaptiveRegion,
        cv::Mat& xGradientSmooth, cv::Mat& yGradientSmooth,
        const float sigmaC, const float sigmaR,
        const int maxDomainSize, cv::Mat& outputImg)
{

    // Get image size
    const int imgWidth  = inputImage.cols;
    const int imgHeight = inputImage.rows;

    // Create constants used throughout code
    const double domainConst = -2.f * sigmaC * sigmaC;
    const double rangeConst =  -2.f * sigmaR * sigmaR;

    // Memory referencing
    cv::Mat xImg = cv::Mat(xGradientSmooth);
    cv::Mat yImg = cv::Mat(yGradientSmooth);
    cv::Mat srcImg = cv::Mat(inputImage);
    cv::Mat outImg = cv::Mat(outputImg);
    cv::Mat adpImg = cv::Mat(adaptiveRegion);

    // Define Look Up Table for weightings
    
    // Compute the weight for the domain filter (domainWeight). The domain filter
    // is a Gaussian lowpass filter
    cv::Mat domainWeightLUT = cv::Mat( cv::Size(maxDomainSize,maxDomainSize), inputImage.depth(), 1);
        
    // Memory reference
    cv::Mat domainWeight = cv::Mat(domainWeightLUT);

    for (int y = 0; y < domainWeightLUT.rows ; y++)
    {
        for (int x = 0; x < domainWeightLUT.cols ; x++)
        {
            // weight for the domain filter (domainWeight)
            const float diff = (float) (x*x+y*y);
            domainWeight.ptr(y)[x] = (float) exp( diff / domainConst );
        }
    }


    for(int y = 0; y < imgHeight; y++) 
    {   
        for(int x = 0; x < imgWidth; x++) 
        { 
            double normFactor = 0.f;
            double tmp = 0.f;

            // filter window width is calculated from adaptive neighbourhood
            // halfsize is half of the filter window width (or maximum LUT value)
            const int halfSize = (int) adpImg.ptr(y)[x];

            // Coefficients defining the centerplane is calculated
            // from the smoothed image gradients
            // Plane Equation is z = coeffA.x + coeffB.y + coeffC
            // coeffA = dI/dx, coeffB = dI/dy, coeffC = I at center pixel of the filter kernel
            const float coeffA = xImg.ptr(y)[x];
            const float coeffB = yImg.ptr(y)[x];
            const float coeffC = srcImg.ptr(y)[x];

            for (int n = -halfSize; n <= halfSize; n++) // y scan line
            {
                for(int m = -halfSize; m <= halfSize; m++) // x scan line
                {
                    // Get domain weight from LUT
                    const float dWeight = domainWeight.ptr(abs(n))[abs(m)];
                
                    // Only perform calculationg if within bounds
                    const int localX = x + m;
                    if (localX < 0) continue;
                    if (localX >= imgWidth) continue;

                    const int localY = y + n;
                    if (localY < 0) continue;
                    if (localY >= imgHeight) continue;

                    // Compute the detail signal (detail) based on the difference between a 
                    // neighborhood pixel and the centerplane passing through the center-pixel 
                    // of the filter window. 
                    const float detail = 
                        srcImg.ptr(localY)[localX] - coeffA * float(m) - coeffB * float(n) - coeffC;

                    // Compute the weight for the range filter (rangeWeight). The range filter
                    // is a Gaussian filter defined by the detail signal.
                    const float rangeWeight = exp( pow(detail,2) / rangeConst );    

                    // Add to function
                    tmp += detail * dWeight * rangeWeight;

                    //Detail Bilateral filter normalized by normFactor 
                    normFactor += dWeight * rangeWeight;
                }
            }

            // Normalise according to weight
            outImg.ptr(y)[x] = tmp / normFactor + coeffC;
        }   
    }   


}