Image processing parallelization

Question

I'm working in a big project so I decided to take just part of it which you can copy/paste and compile in your machine. You might get weird image at the end but that's fine, that's what I want.

I'm sharing with you this part of the code for the main reason to make it run faster, but if you have any other advice or anything to say about my code feel free.

The program simply takes two images and do some processing and gives you an image at the end as a result. The program is little bit long but I believe it is easy to understand.

You can compile this program:

g++ -g -std=c++1z -Wall -Weffc++ -Ofast -march=native test5.cpp -o test5 -fopenmp `pkg-config --cflags --libs opencv`

And run the program like that:

./test5 image1.png image2.png

This is the code:

#include <opencv2/core.hpp>
#include <opencv2/highgui.hpp>
#include <opencv2/imgproc.hpp>
#include <opencv2/imgcodecs.hpp>
#include <opencv2/core/utility.hpp>

#include <iostream>
#include <chrono>

using std::chrono::high_resolution_clock;
using std::chrono::duration_cast;
using std::chrono::microseconds;


class Parallel_process : public cv::ParallelLoopBody
    {

    private:
        cv::Mat img;
        std::vector<int> A;
        int diff;

    public:
        Parallel_process(cv::Mat inputImgage, std::vector<int> AA, int diffVal)
                           : img(inputImgage), A(AA), diff(diffVal){}

        virtual void operator()(const cv::Range& range) const
        {
            for(int i = range.start; i < range.end; i++)
            {

               cv::Mat in(img, cv::Rect(0, (img.rows/diff)*i, img.cols, img.rows/diff));
               std::vector<int> AAA (A);
               in.forEach<cv::Vec3f>
                (
                  [&AAA](cv::Vec3f &pixel, const int* po) -> void
                  {
                    pixel[0]/=AAA[0];
                    pixel[1]/=AAA[1];
                    pixel[2]/=AAA[2];
                  }
                );
            }
        }
};

cv::Mat dcp(const cv::Mat&, auto, auto, const cv::Mat&, double);
auto calculateSD(const cv::Mat&,auto, auto);
void fftshift(cv::Mat&);
cv::Mat transmission(cv::Mat&, cv::Mat&);

void GammaCorrection(cv::Mat&, unsigned char*, cv::Mat&);

template <typename T, typename ... Ts>
void insert_all(std::vector<T> &vec, Ts ... ts)
{
    (vec.push_back(ts), ...);
}

typedef std::vector<std::vector<int> > Matrix;


int main(int argc, char* argv[])
{

    cv::Mat im_test = cv::imread(argv[1]);// = cv::Mat::zeros(src.rows, src.cols, CV_32FC1);

    auto rows=im_test.rows,
         cols=im_test.cols;
    cv::Mat fin_img;                
    cv::Mat src=cv::imread(argv[2]);
    cv::Mat src_temp = src.clone();

    // build look up table
    unsigned char lut[256];
    auto fGamma=0.4;
    #pragma omp for
    for (size_t i=0; i<256; i++)
        lut[i] = cv::saturate_cast<uchar>(pow((float)(i / 255.0), fGamma) * 255.0f);

    //std::cout<<cv::getBuildInformation()<<std::endl;  

        high_resolution_clock::time_point t1(high_resolution_clock::now());

    GammaCorrection(src_temp, lut, src_temp);
    std::vector<cv::Mat> rgb;
    cv::split(src_temp, rgb);
    Matrix histSum(3, std::vector<int>(256,0));   

    src_temp.forEach<cv::Vec3b>
    (
      [&histSum](cv::Vec3b &pixel, const int* po) -> void
      {
        ++histSum[0][pixel[0]];
        ++histSum[1][pixel[1]];
        ++histSum[2][pixel[2]];
      }
    );

    std::vector<int> A(3, 255);
    auto A_estim_lambda([&A, rows, cols, &histSum]{

        for (auto index=8*rows*cols/1000; index>histSum[0][A[0]]; --A[0])        
         index -= histSum[0][A[0]];
    for (auto index=8*rows*cols/1000; index>histSum[1][A[1]]; --A[1])
             index -= histSum[1][A[1]];
    for (auto index=8*rows*cols/1000; index>histSum[2][A[2]]; --A[2])
             index -= histSum[2][A[2]];
        return A;
   });
    auto AA=A_estim_lambda();

    cv::Mat srcN = src_temp.clone();
    srcN.convertTo(srcN, CV_32FC3);
    im_test.convertTo(im_test, CV_32FC3);

    cv::parallel_for_(cv::Range(0, 91), Parallel_process(srcN, AA, 91));   

    cv::Mat IllumTrans = transmission(srcN, im_test);

    std::vector<cv::Mat> rgbDCP;    
    rgbDCP.reserve(3);
    insert_all(rgbDCP, dcp(rgb[0], rows, cols, IllumTrans, A[0]), 
                       dcp(rgb[1], rows, cols, IllumTrans, A[1]), 
                       dcp(rgb[2], rows, cols, IllumTrans, A[2]));

    cv::merge(rgbDCP, fin_img);

    cv::medianBlur(fin_img, fin_img, 3); //5 c trop

    fin_img.convertTo(fin_img, CV_8UC3, 255.0);

    cv::Mat temp;
    cv::GaussianBlur(fin_img, temp, cv::Size(0, 0), 3);
    cv::addWeighted(fin_img, 1.5, temp, -0.5, 0, fin_img);

        fGamma=1.5;
        for (size_t i=0; i<256; i++)
                lut[i] = cv::saturate_cast<uchar>(pow((float)(i / 255.0), fGamma) * 255.0f);      

        GammaCorrection(fin_img, lut, fin_img);

    high_resolution_clock::time_point t2(high_resolution_clock::now());
    auto timeEnd=1.0/static_cast<double>(duration_cast<microseconds>(t2 - t1).count())*1000000;
    std::cout<<timeEnd<<std::endl;

    cv::imshow("kernel", fin_img);
    cv::waitKey();

    return 0;
}

void GammaCorrection(cv::Mat& src, unsigned char* lut, cv::Mat& dst)
{   
    dst.forEach<cv::Vec3b>
        (
          [&lut](cv::Vec3b &pixel, const int* po) -> void
          {
                  pixel[0] = lut[(pixel[0])];
                  pixel[1] = lut[(pixel[1])];
                  pixel[2] = lut[(pixel[2])];

          }
        );
}


auto calculateSD(const cv::Mat& src, auto rows, auto cols)
{

    double sum{0};
    double sq_sum{0};

    #pragma omp for
    for(auto j=0;j<rows;j++)
        for(auto i=0;i<cols;i++)
        {
            sum += src.at<float>(j,i);
            sq_sum += src.at<float>(j,i) * src.at<float>(j,i);
        }


    double mean = sum / (rows*cols);
    double variance = sq_sum / (rows*cols) - mean * mean;

    return sqrt(variance);
}

cv::Mat transmission(cv::Mat& srcN, cv::Mat& im_test )
{
    cv::Mat srcN_gray;
    cv::cvtColor(srcN, srcN_gray, cv::COLOR_RGB2GRAY);
    cv::cvtColor(im_test, im_test, cv::COLOR_RGB2GRAY);
    cv::Mat srcN_fft;
    dft(srcN_gray, srcN_fft, cv::DFT_COMPLEX_OUTPUT) ;
    dft(im_test, im_test, cv::DFT_COMPLEX_OUTPUT);

    cv::Mat mul_fft;    
    cv::mulSpectrums(im_test, srcN_fft, mul_fft, 0);

    cv::Mat mul_invfft;
    dft(mul_fft, mul_invfft, cv::DFT_INVERSE | cv::DFT_SCALE | cv::DFT_REAL_OUTPUT);

    fftshift(mul_invfft);

    float stddev=calculateSD(mul_invfft, im_test.rows, im_test.cols);

    return (1-(mul_invfft-stddev));

}

cv::Mat dcp(const cv::Mat& src, auto rows, auto cols, const cv::Mat& IllumTrans, double A )
{

    cv::Mat imJ=cv::Mat::zeros(rows, cols, CV_32FC1);
    #pragma omp for
    for(auto j=0;j<rows;j++)
        for(auto i=0;i<cols;i++)
            imJ.at<float>(j,i)= A+((src.at<uchar>(j,i)-A)/std::max(IllumTrans.at<float>(j,i), 0.1f));

    double minVal=0, maxVal=0;
    minMaxLoc(imJ, &minVal, &maxVal);

    return imJ/maxVal;  
}




void fftshift(cv::Mat& src)
{
    int cx = src.cols/2;
    int cy = src.rows/2;

    cv::Mat q0(src, cv::Rect(0, 0, cx, cy));   
    cv::Mat q1(src, cv::Rect(cx, 0, cx, cy));  
    cv::Mat q2(src, cv::Rect(0, cy, cx, cy));  
    cv::Mat q3(src, cv::Rect(cx, cy, cx, cy)); 

    cv::Mat tmp;                           
    q0.copyTo(tmp);
    q3.copyTo(q0);
    tmp.copyTo(q3);

    q1.copyTo(tmp);
    q2.copyTo(q1);
    tmp.copyTo(q2); 
}

Answer 1

A quick, partial review (just one function) [maybe I'll have more time later to figure out what your code does].

auto calculateSD(const cv::Mat& src, auto rows, auto cols)
{
    double sum{0};
    double sq_sum{0};

    #pragma omp for
    for(auto j=0;j<rows;j++)
        for(auto i=0;i<cols;i++)
        {
            sum += src.at<float>(j,i);
            sq_sum += src.at<float>(j,i) * src.at<float>(j,i);
        }

    double mean = sum / (rows*cols);
    double variance = sq_sum / (rows*cols) - mean * mean;

    return sqrt(variance);
}

This code will probably give you wrong results. You have multiple threads all updating the same two variables. sum += x is the same as sum = sum + x. Each thread reads sum, updates the value, and writes it back. Before it's written back, the value of sum has likely changed, those changes will be gone.

OpenMP has a reduction clause just for this type of loop:

#pragma omp parallel for reduction(+:sum,sq_sum)

With this clause, each thread gets its own copy of the variables sum and sq_sum, and at the end of the loop these local variables are added together into the function's copy of these variables.

The algorithm you are using here to compute the variance is fast, but it's unstable. If the mean is large w.r.t. the variance, you'll get nonsensical results. See this recent post on Code Review for some discussions regarding this issue. There are links there to resources about Welford's Algorithm, which is a stable one-pass algorithm to compute mean and variance.

Lastly, I honestly don't see the advantage of declaring a function auto. Writing double there only requires 2 more keystrokes. In exchange, your function's interface is explicit, which makes using it so much easier. The same is true for the auto in the input arguments.* I don't know what integer type OpenCV uses for row and column sizes, but if it's not an int, then your code here

for(auto j=0;j<rows;j++)

probably doesn't do what you want, because j will be declared int, not the type of rows (0 is an int).

EDIT:

*: Reading more about auto parameters in a function declaration, I learned that it is a Concepts TS thing, and is a shortcut for a template specification. That means that it is not standard C++ yet. It also means that this function is a function template. I don't like the idea of a template that is not obviously so, and don't like the idea of a template that will only ever be used with one type. It prevents the C++ type system to do its thing, and it could mean you end up compiling more versions of your function than necessary. Type is an important part of C++, I am not in favor of hiding it away from the reader of the code, and especially not from the user of your functions. I would prefer explicit types in the function declaration.

Answer 2

This sounds like a really interesting project! Here are a few thoughts on your code.

Naming

I have to admit that your naming is a mix of good and bad. The constructor for Parallel_process uses inputImgage, which, misspelling aside, tells me what it is. But then you have a vector named AA. When I see AA in image processing code, I usually assume it means either "axis-aligned" (as in AABB = "Axis-aligned bounding box"), or "antialiasing". It doesn't appear to mean either in this case. What does it mean? And you assign AA to a member variable named A. What is A?

The name diffVal and diff are almost as bad. At least I know it's the difference between two things. (Or is it a differential?) Whenever you find yourself using the word "value" (or some abbreviation of it) in a name, you probably need to rethink the name. What is it the difference of? You're dividing by it, so it seems like maybe it's actually a range - like the difference between the minimum and maximum values of… something. It would be nice to let a reader of your code know what that something is.

Then you have a value inside operator() named AAA which appears to be a copy of A, despite the fact that you don't modify it at all. Why are you not just using A wherever you have AAA? It would speed things up as you wouldn't need to copy A.

It's not immediately clear what dcp stands for.

I assume calculateSD() is calculating the standard deviation? I would clarify SD to be standardDeviation. You can shorten it if you like, but SD is overloaded, too. (Standard Definition, super density, standard deviation, Jacobi's function.)

Functions

Your main() function could be made much simpler and easier to read if you broke it into functions with descriptive names. This:

// build look up table
unsigned char lut[256];
auto fGamma=0.4;
#pragma omp for
for (size_t i=0; i<256; i++)
    lut[i] = cv::saturate_cast<uchar>(pow((float)(i / 255.0), fGamma) * 255.0f);

//std::cout<<cv::getBuildInformation()<<std::endl;  

    high_resolution_clock::time_point t1(high_resolution_clock::now());

GammaCorrection(src_temp, lut, src_temp);

could all be put into a function named convertToLinearRGB(). That said, are you sure you want to use a gamma of 0.4? If you're dealing with most normal image (sRGB) or video (Rec. 709) formats, 1.0 / 2.2 = .4545… would be a better choice. The linear offset near 0 makes 2.4 or 2.5 suboptimal choices for a conversion, despite the fact that they are used in the actual calculation.

Next, you should put the histogram calculation into a function named histogram().

For A_estim_lambda, why do you define the lambda and then immediately define a variable that is simply assigned to the lambda? Why not just do it in 1 line?

Reduce Complexity

I see at least 2 different parallel computation systems in use - OpenMP and OpenCV's parallel structures. If you get some material advantage out of them, then maybe it's worth it, but using so many different systems makes it more complex for maintenance and understanding.

Also, the first loop in main() is unlikely to be helped by being parallelized. For 256 values, the overhead of creating multiple threads is very likely to be more than the time it takes to just do the calculations.

Performance

Unfortunately, I can't run your code and profile it because I don't have OpenCV installed on my machine. But you should profile your code to see which specific lines in which specific functions are costing the most time. Without doing that, it doesn't make sense to start optimizing it because you don't know which parts are the slowest. You appear to be compiling with GCC, so you can probably use gprof for profiling. Depending on what OS you're on, there may be additional tools for profiling. I recommend you check them out.