Divide every channel of image by weight image in opencv

Question

I have an 3 channel image that is basically a composition of weighted patches. After adding all contributions to my sum image, I want to divide it by the weights. Right now, I use the following code:

Mat weighted_sum, weights;
// for (Patch p: patches)
//   weighted_sum[loc_p] += w * p
//   weights[loc_p] += w;

vector<Mat> channels(3);
split(weighted_sum, channels);
for (Mat chan: channels) {
    divide(chan, weights, chan);
}
merge(channels, reconstructed);

Is there any more efficient solution for dividing all 3 channels of an image by the same 1-channel image?

Answer 1

I tested your method against:

1. split - divide - merge, @panmari
2. Linearize matrices, and apply element wise division @Miki
3. Linearize matrices, and apply element wise division (with pointers) @Miki and multiplication instead of division @Micka
4. Make weights a 3 channel matrix, and apply divide on whole matrices @Micka.

Results (in ms):

Size            Method1         Metdhod2        Method3         Method4
[2 x 2]         0.0359212       0.00183271      0.000733086     1.77333
[10 x 10]       0.0117294       0.00293234      0.00109963      0.0051316
[100 x 100]     0.422624        0.241918        0.0751413       0.319625
[1000 x 1000]   20.757          20.3673         7.28284         18.4389
[2000 x 2000]   83.6238         82.942          28.4353         74.2132

NOTES

• My method works faster, but the speed-up is related to the size of the matrices.
• Method 3 (with pointers) is the fastest
• Method 2 and Method 3 modify the original matrix. If you need not to change original matrix, you need to do a deep copy (clone(), see commented lines in code below). With deep copy of the original matrix both methods are slower, but Method 3 is still the fastest.
• Method 4 won't work with double matrices, since cvtColor doesn't accepts double).
• Preallocating matrices for Method 4 is just a small improvement.
• Tested with compiler msvc12 (Visual Studio 2013). With gcc 4.8 Method 4 seems to be faster.

Here the code I used. I'm testing on float matrices, but it's easy to port to other types.

    #include <opencv2/opencv.hpp>
#include <iostream>
#include <vector>
using namespace std;
using namespace cv;


int main()
{
    vector<Size> sizes{Size(2,2), Size(10,10), Size(100,100), Size(1000,1000), Size(2000,2000)};

    cout << "Size \t\tMethod1 \tMetdhod2 \tMethod3 \tMethod4" << endl;

    for (int is = 0; is < sizes.size(); ++is)
    {

        Size sz = sizes[is];
        Mat3f weighted_sum(sz);
        randu(weighted_sum, 0, 200);

        Mat1f weights(sz);
        randu(weights, 0, 10);

        Mat3f ws1 = weighted_sum.clone();
        Mat3f ws2 = weighted_sum.clone();
        Mat3f ws3 = weighted_sum.clone();
        Mat3f ws4 = weighted_sum.clone();



        // Method 1 @parmari
        double tic1 = double(getTickCount());
        Mat3f rec1;

        vector<Mat> channels(3);
        split(ws1, channels);
        for (Mat chan : channels) {
            divide(chan, weights, chan);
        }
        merge(channels, rec1);

        double toc1 = (double(getTickCount() - tic1)) * 1000. / getTickFrequency();

        // Method 2 @Miki
        double tic2 = double(getTickCount());
        Mat3f rec2 = ws2.reshape(3, 1);
        //Mat3f rec2 = ws2.reshape(3, 1).clone(); // To not override original image
        Mat1f ww2 = weights.reshape(1, 1);
        for (int i = 0; i < rec2.cols; ++i)
        {
            double w = ww2(0, i);
            Vec3f& v = rec2(0, i);
            v[0] /= w;
            v[1] /= w;
            v[2] /= w;
        }
        rec2 = rec2.reshape(3, ws2.rows);

        double toc2 = (double(getTickCount() - tic2)) * 1000. / getTickFrequency();


        // Method 3 @Miki (+ @Micka)
        double tic3 = double(getTickCount());
        Mat3f rec3 = ws3.reshape(3, 1);
        //Mat3f rec3 = ws3.reshape(3, 1).clone(); // To not override original image
        Mat1f ww3 = weights.reshape(1, 1);

        Vec3f* prec3 = rec3.ptr<Vec3f>(0);
        float* pww = ww3.ptr<float>(0);

        for (int i = 0; i < rec3.cols; ++i)
        {
            float scale = 1. / (*pww);
            (*prec3)[0] *= scale;
            (*prec3)[1] *= scale;
            (*prec3)[2] *= scale;

            ++prec3; ++pww;
        }
        rec3 = rec3.reshape(3, ws3.rows);

        double toc3 = (double(getTickCount() - tic3)) * 1000. / getTickFrequency();


        // Method 4 @Micka
        double tic4 = double(getTickCount());
        Mat3f rec4;
        Mat3f w3ch4;
        cvtColor(weights, w3ch4, COLOR_GRAY2BGR);
        divide(ws4, w3ch4, rec4);

        double toc4 = (double(getTickCount() - tic4)) * 1000. / getTickFrequency();

        cout << sz << " \t" << toc1 << " \t" << toc2 << " \t" << toc3 << " \t" << toc4 <<  endl;

    }

    getchar();

    return 0;
}