Poisson Blending II

I re-wrote the Poisson Blending code using C++ and OpenCV.

About the Algorithm, see my Previous Poisson Blending post.

This time, I just used the most stupid way, just solving the Poisson Equation. You can improve it by using advanced methods. About solving discrete Poisson Equation using Jacobi, SOR, Conjugate Gradients, and FFT, read THIS 

In this code, I used two different ways to calculate vector (actually they are the same…), The first way is calculate gradient by Convolution a Laplacian Matrix, the second way is calculate gradient directly. You can use any one you like.

CODE

// Poisson Blending
// Using stupid way (just solving Ax = b)
// Author: Eric Yuan
// Blog: http://eric-yuan.me
// You are free to use the following code for ANY purpose.
// 
// You can improve it by using advanced methods, 
// About solving discrete Poisson Equation using Jacobi, SOR, Conjugate Gradients, and FFT
// see here: http://www.cs.berkeley.edu/~demmel/cs267/lecture24/lecture24.html
// 

#include "opencv2/core/core.hpp"
#include "opencv2/imgproc/imgproc.hpp"
#include "opencv2/highgui/highgui.hpp"
#include &ltmath.h&gt
#include 
using namespace cv;
using namespace std;

#define elif else if
#define ATD at

// calculate horizontal gradient, img(i, j + 1) - img(i, j)
Mat 
getGradientXp(Mat &img){
    int height = img.rows;
    int width = img.cols;
    Mat cat = repeat(img, 1, 2);

    Rect roi = Rect(1, 0, width, height);
    Mat roimat = cat(roi);
    return roimat - img;
}

// calculate vertical gradient, img(i + 1, j) - img(i, j)
Mat 
getGradientYp(Mat &img){
    int height = img.rows;
    int width = img.cols;
    Mat cat = repeat(img, 2, 1);

    Rect roi = Rect(0, 1, width, height);
    Mat roimat = cat(roi);
    return roimat - img;
}

// calculate horizontal gradient, img(i, j - 1) - img(i, j)
Mat 
getGradientXn(Mat &img){
    int height = img.rows;
    int width = img.cols;
    Mat cat = repeat(img, 1, 2);

    Rect roi = Rect(width - 1, 0, width, height);
    Mat roimat = cat(roi);
    return roimat - img;
}

// calculate vertical gradient, img(i - 1, j) - img(i, j)
Mat 
getGradientYn(Mat &img){
    int height = img.rows;
    int width = img.cols;
    Mat cat = repeat(img, 2, 1);

    Rect roi = Rect(0, height - 1, width, height);
    Mat roimat = cat(roi);
    return roimat - img;
}

int
getLabel(int i, int j, int height, int width){
    return i * width + j;
}

// get Matrix A.
Mat
getA(int height, int width){

    Mat A = Mat::eye(height * width, height * width, CV_64FC1);
    A *= -4;
    Mat M = Mat::zeros(height, width, CV_64FC1);
    Mat temp = Mat::ones(height, width - 2, CV_64FC1);
    Rect roi = Rect(1, 0, width - 2, height);
    Mat roimat = M(roi);
    temp.copyTo(roimat);
    temp = Mat::ones(height - 2, width, CV_64FC1);
    roi = Rect(0, 1, width, height - 2);
    roimat = M(roi);
    temp.copyTo(roimat);
    temp = Mat::ones(height - 2, width - 2, CV_64FC1);
    temp *= 2;
    roi = Rect(1, 1, width - 2, height - 2);
    roimat = M(roi);
    temp.copyTo(roimat);

    for(int i=0; i<height; i++){
        for(int j=0; j<width; j++){
            int label = getLabel(i, j, height, width);
            if(M.ATD(i, j) == 0){
                if(i == 0)  A.ATD(getLabel(i + 1, j, height, width), label) = 1;
                elif(i == height - 1)   A.ATD(getLabel(i - 1, j, height, width), label) = 1;
                if(j == 0)  A.ATD(getLabel(i, j + 1, height, width), label) = 1;
                elif(j == width - 1)   A.ATD(getLabel(i, j - 1, height, width), label) = 1;
            }elif(M.ATD(i, j) == 1){
                if(i == 0){
                    A.ATD(getLabel(i + 1, j, height, width), label) = 1;
                    A.ATD(getLabel(i, j - 1, height, width), label) = 1;
                    A.ATD(getLabel(i, j + 1, height, width), label) = 1;
                }elif(i == height - 1){
                    A.ATD(getLabel(i - 1, j, height, width), label) = 1;
                    A.ATD(getLabel(i, j - 1, height, width), label) = 1;
                    A.ATD(getLabel(i, j + 1, height, width), label) = 1;
                }
                if(j == 0){
                    A.ATD(getLabel(i, j + 1, height, width), label) = 1;
                    A.ATD(getLabel(i - 1, j, height, width), label) = 1;
                    A.ATD(getLabel(i + 1, j, height, width), label) = 1;
                }elif(j == width - 1){
                    A.ATD(getLabel(i, j - 1, height, width), label) = 1;
                    A.ATD(getLabel(i - 1, j, height, width), label) = 1;
                    A.ATD(getLabel(i + 1, j, height, width), label) = 1;
                }
            }else{
                    A.ATD(getLabel(i, j - 1, height, width), label) = 1;
                    A.ATD(getLabel(i, j + 1, height, width), label) = 1;
                    A.ATD(getLabel(i - 1, j, height, width), label) = 1;
                    A.ATD(getLabel(i + 1, j, height, width), label) = 1;
            }
        }
    }
    return A;
}

// Get the following Laplacian matrix
// 0  1  0
// 1 -4  1
// 0  1  0
Mat
getLaplacian(){
    Mat laplacian = Mat::zeros(3, 3, CV_64FC1);
    laplacian.ATD(0, 1) = 1.0;
    laplacian.ATD(1, 0) = 1.0;
    laplacian.ATD(1, 2) = 1.0;
    laplacian.ATD(2, 1) = 1.0;
    laplacian.ATD(1, 1) = -4.0; 
    return laplacian;
}

// Calculate b
// using convolution.
Mat
getB1(Mat &img1, Mat &img2, int posX, int posY, Rect ROI){
    Mat Lap;
    filter2D(img1, Lap, -1, getLaplacian());
    int roiheight = ROI.height;
    int roiwidth = ROI.width;
    Mat B = Mat::zeros(roiheight * roiwidth, 1, CV_64FC1);
    for(int i=0; i<roiheight; i++){
        for(int j=0; j<roiwidth; j++){
            double temp = 0.0;
            temp += Lap.ATD(i + ROI.y, j + ROI.x);
            if(i == 0)              temp -= img2.ATD(i - 1 + posY, j + posX);
            if(i == roiheight - 1)  temp -= img2.ATD(i + 1 + posY, j + posX);
            if(j == 0)              temp -= img2.ATD(i + posY, j - 1 + posX);
            if(j == roiwidth - 1)   temp -= img2.ATD(i + posY, j + 1 + posX);
            B.ATD(getLabel(i, j, roiheight, roiwidth), 0) = temp;
        }
    }
    return B;
}

// Calculate b
// using getGradient functions.
Mat
getB2(Mat &img1, Mat &img2, int posX, int posY, Rect ROI){
    Mat grad = getGradientXp(img1) + getGradientYp(img1) + getGradientXn(img1) + getGradientYn(img1);
    int roiheight = ROI.height;
    int roiwidth = ROI.width;
    Mat B = Mat::zeros(roiheight * roiwidth, 1, CV_64FC1);
    for(int i=0; i<roiheight; i++){
        for(int j=0; j<roiwidth; j++){
            double temp = 0.0;
            temp += grad.ATD(i + ROI.y, j + ROI.x);
            if(i == 0)              temp -= img2.ATD(i - 1 + posY, j + posX);
            if(i == roiheight - 1)  temp -= img2.ATD(i + 1 + posY, j + posX);
            if(j == 0)              temp -= img2.ATD(i + posY, j - 1 + posX);
            if(j == roiwidth - 1)   temp -= img2.ATD(i + posY, j + 1 + posX);
            B.ATD(getLabel(i, j, roiheight, roiwidth), 0) = temp;
        }
    }
    return B;
}

// Solve equation and reshape it back to the right height and width.
Mat
getResult(Mat &A, Mat &B, Rect &ROI){
    Mat result;
    solve(A, B, result);
    result = result.reshape(0, ROI.height);
    return  result;
}

// img1: 3-channel image, we wanna move something in it into img2.
// img2: 3-channel image, dst image.
// ROI: the position and size of the block we want to move in img1.
// posX, posY: where we want to move the block to in img2
Mat
poisson_blending(Mat &img1, Mat &img2, Rect ROI, int posX, int posY){

    int roiheight = ROI.height;
    int roiwidth = ROI.width;
    Mat A = getA(roiheight, roiwidth);

    // we must do the poisson blending to each channel.
    vector rgb1;
    split(img1, rgb1);
    vector rgb2;
    split(img2, rgb2);

    vector result;
    Mat merged, res, Br, Bg, Bb;
    // For calculating B, you can use either getB1() or getB2()
    Br = getB1(rgb1[0], rgb2[0], posX, posY, ROI);
    //Br = getB2(rgb1[0], rgb2[0], posX, posY, ROI);
    res = getResult(A, Br, ROI);
    result.push_back(res);
    cout<<"R channel finished..."<<endl;
    Bg = getB1(rgb1[1], rgb2[1], posX, posY, ROI);
    //Bg = getB2(rgb1[1], rgb2[1], posX, posY, ROI);
    res = getResult(A, Bg, ROI);
    result.push_back(res);
    cout<<"G channel finished..."<<endl;
    Bb = getB1(rgb1[2], rgb2[2], posX, posY, ROI);
    //Bb = getB2(rgb1[2], rgb2[2], posX, posY, ROI);
    res = getResult(A, Bb, ROI);
    result.push_back(res);
    cout<<"B channel finished..."<<endl;

    // merge the 3 gray images into a 3-channel image 
    merge(result,merged);
    return merged; 
}

int 
main(int argc, char** argv)
{
    long start, end;
    start = clock();

    Mat img1, img2;
    Mat in1 = imread("face.png");
    Mat in2 = imread("lisa.png");
    imshow("src", in1);
    imshow("dst", in2);
    in1.convertTo(img1, CV_64FC3);
    in2.convertTo(img2, CV_64FC3);

    Rect rc = Rect(25, 40, 80, 20);
    Mat result = poisson_blending(img1, img2, rc, 15, 50);
    result.convertTo(result, CV_8UC1);
    Rect rc2 = Rect(15, 50, 80, 20);
    Mat roimat = in2(rc2);
    result.copyTo(roimat);

    end = clock();
    cout<<"used time: "<<((double)(end - start)) / CLOCKS_PER_SEC<<" second"<<endl;
    imshow("roi", result);
    imshow("result", in2);

    waitKey(0);
    return 0;
}

Have fun with it 🙂

Here is Dobiasd’s version which fixed several bug, thanks Dobiasd.



												
					
This entry was posted in Algorithm, Maths, OpenCV and tagged , , , , . Bookmark the permalink. Post a comment or leave a trackback: Trackback URL.

15 Comments

  1. Marco
    Posted June 17, 2014 at 7:08 pm | Permalink

    Hello, thanks for sharing, your code looks great because it doesn’t require additional libraries then opencv’s ones, but on line #15 and #16 there are 2 incompleted include, fixing this error gives another error on ATD, on the compiler i got (No Matching Member Function for call to ‘at’), looks like you have forgotten the , can you provide the fixed code? Thanks a lot and have a nice day!

    • Eric
      Posted June 17, 2014 at 9:46 pm | Permalink

      Hi Marco,

      Thanks for asking, I found that angle brackets (and code in angle brackets) in code do not show correctly and that caused the problem. The two missing includes are math.h and iostream; and the ATD part is #define ATD at(angle bracket)double(closing angle bracket). Hope this comment helps.

      • Marco
        Posted June 22, 2014 at 11:08 am | Permalink

        There was another error on the vector starting from on line 216, they should be vector of Mat object, thanks again and have a nice day!

        • Eric
          Posted June 23, 2014 at 10:20 am | Permalink

          Oh yeah, another angle brackets problem. U r right, vector of Mat. Cheers!

  2. Marco
    Posted June 17, 2014 at 7:16 pm | Permalink

    #define ATD at fixed some problem but there are others like “vector rgb1;” vector does not exist…

  3. Marco
    Posted June 17, 2014 at 7:18 pm | Permalink

    sorry the type looks like is filtered by your html tag, the define should look like this #define ATD at

  4. Marco
    Posted June 17, 2014 at 7:18 pm | Permalink

    sorry it’s filtered even in the code tags, i hope you have understood the problem

  5. Posted January 7, 2015 at 7:34 am | Permalink

    Hi,
    that is a very cool algorithm. Thanks for sharing.

    Just some remarks to your code:
    – OpenCV loads color images in BGR, not in RGB. So the channel variables are labeled incorrectly. But is has no consqeuences in your case.
    – ctime has to be included to use clock()

    The fixed code (also with corrected brackets problem) can be found here: http://codepad.org/PU8ezzdl

    Perhaps you want to include a link to the correct raw code in your post, so not every of your readers has to fix this independently.

    Dobiasd

    • Eric
      Posted January 10, 2015 at 1:25 am | Permalink

      Thanks Dobiasd

  6. Josep Bosch
    Posted April 9, 2015 at 4:53 am | Permalink

    Dobiasd’s link looks dead…
    I posted the working code in a new one:

    http://codepad.org/ANqtikKR

    • Eric
      Posted April 13, 2015 at 4:18 pm | Permalink

      Thanks!

    • Lee
      Posted May 24, 2015 at 10:19 pm | Permalink

      thank you for posting the code

      • mrhloom
        Posted July 6, 2015 at 4:39 am | Permalink

        We we use only rectangular areas and no binary mask is specified?

    • Muo
      Posted June 26, 2015 at 5:43 am | Permalink

      Why I can’t process the roi bigger than 100 pixels each edge.
      This line will occur error:
      “Mat A = Mat::eye(height * width, height * width, CV_64FC1);”
      The opencv show the error message “The total matrix size does not fit to “size_t” type”
      Is there anybody know this problem.

  7. mrhloom
    Posted July 6, 2015 at 10:25 am | Permalink

    Also what if PosX=0 or PosY=0( in GetB )program breaks?

One Trackback

Post a Comment

Your email is never published nor shared. Required fields are marked *

You may use these HTML tags and attributes <a href="" title=""> <abbr title=""> <acronym title=""> <b> <blockquote cite=""> <cite> <code> <del datetime=""> <em> <i> <q cite=""> <s> <strike> <strong>

*
*