Poisson Blending

09 9月 2011 Under: opencv2.x-samples

このサンプルは,OpenCVの機能とはあまり関係なく,差分やデータの保持にcv::Mat形式を利用している程度です.このサンプルでは,ある画像の一部を別の画像にコピーする際に,それらを滑らかにブレンディングします.
コピー後の画素値は,Dirichlet条件の下でPoisson方程式を解くことで求められます.つまり,コピー元画像の画像勾配をなるべく保ったまま,コピー境界の画素値をコピー先の画素値と合うように,コピー結果画素値を決定します.
詳しくは,SIGGRAPH2003の論文 Poisson Image Editting (PDF) を参照してください.

Souce Image, Target Image, Mask Image, Blending Result

また,Gradient Mixtureを行うことで,コピー元とコピー先の画像のうち,より強い勾配を保存することができます.
下図に,紙に書かれた文字(左)を,手のひらの合成するデモ結果を示します.
Gradient Mixtureを行わない通常のブレンディングの場合(中央),コピー元画像の背景である紙部分が上書きされ,手のシワが消えてのっぺりした画像になっています.
ここで Gradient Mixture を行うと(右),シワなどの強い勾配が保存されたブレンディング結果になります.

 

以下は,サンプル用のソースコードですが,実行には,OpenCV-2.2以降,Eigen2/3,UmfPack が必要です.Ubuntu 11.04 上でのみ動作を確認しています.
また,OpenCV自身の Blender クラスとはまったくの無関係です.
疎な連立方程式の解法部分は,UmfPackのソルバに丸投げですが,Gauss-Seidel 法,Mutigrid 法などで解いても良いでしょう.
モナリザの画像(375×505[pixel])を例に挙げると,Core i7 2.67[GHz] の VM 上 Linux で,約950[msec] ほどの処理時間が必要です.

C++

poisson_blending.hpp

#ifndef __POISSON_BLENDING_H__
#define __POISSON_BLENDING_H__

#include <map>
#include <opencv2/imgproc/imgproc.hpp>
#include <Eigen/Core>
#define EIGEN_YES_I_KNOW_SPARSE_MODULE_IS_NOT_STABLE_YET
#include <Eigen/Sparse>
#include <unsupported/Eigen/UmfPackSupport>
//#include <unsupported/Eigen/SuperLUSupport>

namespace Blend {

class PoissonBlender
{
private:
  cv::Mat _src, _target, _mask;

  cv::Rect mask_roi1;
  cv::Mat mask1;
  cv::Mat dst1;
  cv::Mat target1;
  cv::Mat drvxy;

  int ch;

  std::map<int,int> mp;

  template <typename T>
  bool buildMatrix(Eigen::SparseMatrix<T> &A, Eigen::Matrix<T, Eigen::Dynamic, 1> &b, 
                   Eigen::Matrix<T, Eigen::Dynamic, 1> &u);
  template <typename T, typename Backend> 
  bool solve(const Eigen::SparseMatrix<T> &A, const Eigen::Matrix<T, Eigen::Dynamic, 1> &b,
             Eigen::Matrix<T, Eigen::Dynamic, 1> &u);
  template <typename T>
  bool copyResult(Eigen::Matrix<T, Eigen::Dynamic, 1> &u);
public:
  PoissonBlender();
  PoissonBlender(const cv::Mat &src, const cv::Mat &target, const cv::Mat &mask);
  ~PoissonBlender() {};
  bool setImages(const cv::Mat &src, const cv::Mat &target, const cv::Mat &mask);
  void copyTo(PoissonBlender &b) const;
  PoissonBlender clone() const;
  bool seamlessClone(cv::Mat &dst, int offx, int offy, bool mix);
  // bool smoothComplete(cv::Mat &dst); // implemented easily with zero boundary conditions.
};

PoissonBlender::PoissonBlender()
{
}

void PoissonBlender::copyTo(PoissonBlender &b) const
{
  b.setImages(_src, _target, _mask);
}

inline PoissonBlender PoissonBlender::clone() const
{
  PoissonBlender b;
  copyTo(b);
  return b; 
}

// constructor
PoissonBlender::PoissonBlender(const cv::Mat &src, const cv::Mat &target, const cv::Mat &mask=cv::Mat())
  : _src(src), _target(target), _mask(mask)
{
  CV_Assert(_mask.channels()==1);
  CV_Assert(_src.cols==_mask.cols && _src.rows==_mask.rows);
}

// set source, tareget and destination images
bool PoissonBlender::setImages(const cv::Mat &src, const cv::Mat &target, const cv::Mat &mask=cv::Mat())
{
  _src = src;
  _target = target;
  _mask = mask;
  
  CV_Assert(_mask.channels()==1);
  CV_Assert(_src.cols==_mask.cols && _src.rows==_mask.rows);
  
  return true;
}

// solver sparse linear system
template <typename T, typename Backend> 
bool PoissonBlender::solve(const Eigen::SparseMatrix<T> &A, const Eigen::Matrix<T, Eigen::Dynamic, 1> &b, 
                           Eigen::Matrix<T, Eigen::Dynamic, 1> &u)
{
  Eigen::SparseLU<Eigen::SparseMatrix<T>, Backend> lu_of_A(A);
  
  if(!lu_of_A.succeeded()) {
    std::cerr<< "decomposition failed" << std::endl;
    return false;
  }
  if(!lu_of_A.solve(b,&u)) {
    std::cerr<< "solving failed" << std::endl;
    return false;
  }
  
  return true;
}

// build matrix as linear system
template <typename T>
bool PoissonBlender::buildMatrix(Eigen::SparseMatrix<T> &A, Eigen::Matrix<T, Eigen::Dynamic,1> &b, 
                                 Eigen::Matrix<T, Eigen::Dynamic, 1> &u)
{
  int w = mask_roi1.width;
  int h = mask_roi1.height;
  int nz=0;
  for(int y=0; y<h-1; ++y) {
    uchar *p = mask1.ptr(y);
    for(int x=0; x<w-1; ++x, ++p) {
      if(*p==0) continue;

      int id = y*(w*ch)+(x*ch);
      mp[id] = nz++;   // r
      mp[++id] = nz++; // g
      mp[++id] = nz++; // b
    }
  }

  A = Eigen::SparseMatrix<double>(nz, nz);
  b = Eigen::VectorXd(nz);
  u = Eigen::VectorXd(nz);
  int rowA = 0;

  A.reserve(5*nz);
  for(int y=1; y<h-1; ++y) {
    uchar *p = mask1.ptr(y)+1;
    cv::Vec3d *drv = drvxy.ptr<cv::Vec3d>(y)+1;
    for(int x=1; x<w-1; ++x, ++p, ++drv) {
      if(*p==0) continue;

      int id = y*(w*ch)+(x*ch);
      int tidx=id-ch*w, lidx=id-ch, ridx=id+ch, bidx=id+ch*w;

      // to omtimize insertion
      uchar tlrb = 15; // 0b1111
      if(mask1.at<uchar>(y-1,x)==0) {
        *drv -= target1.at<cv::Vec3b>(y-1,x);
        tlrb &= 7; //0b0111
      }
      if(mask1.at<uchar>(y,x-1)==0) {
        *drv -= target1.at<cv::Vec3b>(y,x-1);
        tlrb &= 11; //0b1011
      }
      if(mask1.at<uchar>(y,x+1)==0) {
        *drv -= target1.at<cv::Vec3b>(y,x+1);
        tlrb &= 13; //0b1101
      }
      if(mask1.at<uchar>(y+1,x)==0) {
        *drv -= target1.at<cv::Vec3b>(y+1,x);
        tlrb &= 14; //0b1110
      }
      for(int k=0; k<ch; ++k) {
        A.startVec(rowA+k);
        if(tlrb&8) A.insertBack(mp[tidx+k], rowA+k) = 1.0; // top
        if(tlrb&4) A.insertBack(mp[lidx+k], rowA+k) = 1.0; // left
                   A.insertBack(mp[id  +k], rowA+k) = -4.0;// center
        if(tlrb&2) A.insertBack(mp[ridx+k], rowA+k) = 1.0; // right
        if(tlrb&1) A.insertBack(mp[bidx+k], rowA+k) = 1.0; // bottom
      }
      b(rowA+0) = cv::saturate_cast<double>((*drv)[0]);
      b(rowA+1) = cv::saturate_cast<double>((*drv)[1]);
      b(rowA+2) = cv::saturate_cast<double>((*drv)[2]);
      rowA+=ch;
    }
  }
  A.finalize();
  CV_Assert(nz==rowA);

  return true;
}

template <typename T>
bool PoissonBlender::copyResult(Eigen::Matrix<T, Eigen::Dynamic, 1> &u)
{
  int w = mask_roi1.width;
  int h = mask_roi1.height;
  for(int y=1; y<h-1; ++y) {
    uchar *pd = dst1.ptr(y);
    uchar *pm = mask1.ptr(y)+1;
    for(int x=1; x<w-1; ++x, ++pm) {
      if(*pm==0) {
        pd += 3;
      } else {
        int idx = mp[y*(w*ch)+(x*ch)];
        *pd++ = cv::saturate_cast<uchar>(u[idx+0]);
        *pd++ = cv::saturate_cast<uchar>(u[idx+1]);
        *pd++ = cv::saturate_cast<uchar>(u[idx+2]);
      }
    }
  }

  return true;
}

bool 
PoissonBlender::seamlessClone(cv::Mat &_dst, const int offx, const int offy, const bool mix=false)
{
  // 
  ch = _target.channels();
  cv::Point offset(offx, offy);
  cv::Point tl(_mask.size()), br(-1,-1);

  // calc bounding box
  for(int y=0; y<_mask.rows; ++y) {
    uchar *p = _mask.ptr(y);
    for(int x=0; x<_mask.cols; ++x,++p) {
      if(*p==0) continue;
      if(tl.x>x) tl.x=x;
      if(tl.y>y) tl.y=y;
      if(br.x<x) br.x=x;
      if(br.y<y) br.y=y;
    }
  }
  br.x += 1;
  br.y += 1;

  // add borders
  cv::Rect mask_roi(tl, br);
  cv::Rect mask_roi2(tl-cv::Point(2,2), br+cv::Point(2,2));
  cv::Mat _srcUp, _targetUp, _maskUp, _dstUp;
  cv::copyMakeBorder(_src, _srcUp, 2,2,2,2, cv::BORDER_REPLICATE);
  cv::copyMakeBorder(_target, _targetUp, 2,2,2,2, cv::BORDER_REPLICATE);
  cv::copyMakeBorder(_mask, _maskUp, 1,1,1,1, cv::BORDER_CONSTANT);

  // allocate destination image
  _dstUp = _targetUp.clone();
  _dst = cv::Mat(_dstUp, cv::Rect(2,2,_dstUp.cols-2, _dstUp.rows-2));

  mask_roi1 = cv::Rect(tl-cv::Point(1,1), br+cv::Point(1,1));  
  mask1 = cv::Mat(_mask, mask_roi1);
  target1 = cv::Mat(_targetUp, mask_roi1+offset-cv::Point(1,1));
  dst1 = cv::Mat(_dstUp, mask_roi1+offset-cv::Point(1,1));

  cv::Mat src(_srcUp, mask_roi2);
  cv::Mat target(_targetUp, mask_roi2+offset-cv::Point(2,2));
  cv::Mat mask(_mask, mask_roi2);
  cv::Mat dst(_dstUp, mask_roi2+offset-cv::Point(2,2));
  CV_Assert(src.cols==dst.cols && src.rows==dst.rows);
  
  // calc differential image
  cv::Mat src64, target64;
  int pw = mask_roi2.width-1, ph = mask_roi2.height-1;
  src.convertTo(src64, CV_64F);
  target.convertTo(target64, CV_64F);
  cv::Rect roi00(0,0,pw,ph), roi10(1,0,pw,ph), roi01(0,1,pw,ph);
  cv::Mat _src64_00(src64, roi00), _target64_00(target64, roi00);
  cv::Mat src_dx = cv::Mat(src64, roi10) - _src64_00;
  cv::Mat src_dy = cv::Mat(src64, roi01) - _src64_00;
  cv::Mat target_dx = cv::Mat(target64,roi10) - _target64_00;
  cv::Mat target_dy = cv::Mat(target64,roi01) - _target64_00;

  // gradient mixture
  cv::Mat Dx, Dy;
  if(mix) {   // with gradient mixture
    cv::Mat pdx_src[ch], pdy_src[ch], pdx_target[ch], pdy_target[ch];
    cv::split(src_dx, pdx_src);
    cv::split(src_dy, pdy_src);
    cv::split(target_dx, pdx_target);
    cv::split(target_dy, pdy_target);
    
    cv::Mat _masks_dx, _masks_dy;
    for(int i=0; i<ch; i++) {
      _masks_dx = cv::abs(pdx_src[i]) < cv::abs(pdx_target[i]);
      _masks_dy = cv::abs(pdy_src[i]) < cv::abs(pdy_target[i]);
      pdx_target[i].copyTo(pdx_src[i], _masks_dx);
      pdy_target[i].copyTo(pdy_src[i], _masks_dy);
    }
    cv::merge(pdx_src, ch, Dx);
    cv::merge(pdy_src, ch, Dy);
  } else {   // without gradient mixture
    Dx = src_dx;
    Dy = src_dy;
  }

  // lapilacian
  int w = pw-1, h = ph-1;
  drvxy = cv::Mat(Dx,cv::Rect(1,0,w,h)) - cv::Mat(Dx,cv::Rect(0,0,w,h))
    + cv::Mat(Dy,cv::Rect(0,1,w,h)) - cv::Mat(Dy,cv::Rect(0,0,w,h));

  // 
  // solve an poisson's equation
  //

  Eigen::SparseMatrix<double> A;
  Eigen::VectorXd b;
  Eigen::VectorXd u;

  // build right-hand and left-hand matrix
  buildMatrix<double>(A, b, u);

  // solve sparse linear system
  solve<double, Eigen::UmfPack>(A, b, u);

  // copy computed result to destination image
  copyResult<double>(u);

  return true;
}

}

#endif /*  __POISSON_BLENDING_H__ */

main.cpp

#include <iostream>

#include <opencv2/core/core.hpp>
#include <opencv2/highgui/highgui.hpp>

#include "poisson_blending.hpp"

void
show_usage()
{
  // src_file: source image file.
  // target_file: target image file (to which the source image is copied).
  // mask_file: mask image file (whose size must be the same as soure image).
  // offset_x: offset x-coordinate , which indicates where the origin of source image is copied.
  // offset_y: offset y-coordinate , which indicates where the origin of source image is copied.
  // mix: flag which indicates gradient-mixture is used.
  std::cout << "./pb source_image target_image mask_image offset_x offset_y mix" << std::endl;
}

int
main(int argc, char* argv[])
{
  if(argc < 7) {
    show_usage();
    return -1;
  }
  std::string src_file = argv[1];
  std::string target_file = argv[2];
  std::string mask_file = argv[3];
  int offx = atoi(argv[4]);
  int offy = atoi(argv[5]);
  bool mix = (argv[6]=="true" || atoi(argv[6])>0) ? true : false;

  cv::Mat src_img = cv::imread(src_file, 1);
  if(!src_img.data) return -1;
  
  cv::Mat target_img = cv::imread(target_file, 1);
  if(!target_img.data) return -1;

  cv::Mat mask_img = cv::imread(mask_file, 0);
  if(mask_img.empty()) return -1;

  //cv::Ptr<Blend::PoissonBlender> pb = &Blend::PoissonBlender(src_img, target_img, mask_img);
  Blend::PoissonBlender pb = Blend::PoissonBlender(src_img, target_img, mask_img);
  
  cv::Mat dst_img;
  double f = 1000.0/cv::getTickFrequency();
  int64 time = cv::getTickCount();
  pb.seamlessClone(dst_img, offx, offy, mix);
  std::cout<<(cv::getTickCount()-time)*f<<" [ms]"<<std::endl;
  cv::imwrite("test.png", dst_img);

  return 0;
}


“Poisson Blending” への15件のコメント

  1. panovr より:

    uchar tlrb = 15; // 0b1111
    if(mask1.at(y-1,x)==0) {
    *drv -= target1.at(y-1,x);
    tlrb &= 7; //0b0111
    }

    May you interpret these magic number (15, 7, etc.) more?
    Thanks!

  2. miros より:

    d:\eigen\unsupported\eigen\umfpacksupport(9): fatal error C1083: 无法打开包括文件:“umfpack.h”: No such file or directory

    Could u please help me about that ? thank!

    Environment:win7+MSVC2010

  3. student より:

    プログラムの提供ありがとうございます。

    実装してみたところ二つのエラーが発生しました。
    一つは他の人も言っているように、“UmfPackSupport”がないこと。(Eigen3.1.3を使用)
    もう一つは、265行目で”ch”に値が格納されていないので、宣言時にエラーが発生してしまいます。

    どうすれば解決するか教えてください。
    返信お持ちしてます。

    • teacher より:

      説明にもある通り、疎な連立方程式を解くにあたってUMFPACKを使う必要は必ずしもありません。試してみたところ、Eigenでスパース行列をコレスキー分解して普通に解いても全然大丈夫でしたよ。
      もう1つの問題ですが、265行目で”ch”に値が格納されていないためにエラーを吐いているというわけではなく、たぶん、配列の要素数を変数で宣言していることが問題なのでは? 試しに、動的に確保してみたら大丈夫でした。

  4. Nadin より:

    Dear sir,

    I have been trying to use the package in my android project, I use Native Development Kit (NDK) with OpneCV 2.4.1 library under windows.

    My question is how to use the packages: Eigen & UMFPack in my project, given that U said it’s build only under Ubunto.

    Thanks in advance,

    • Nadin より:

      Dear Sir,

      Kindly reply to my message below as it’s an urgent project for me.

      Thanks in advance,
      ——————————
      “Nadin より:
      2013/01/04 00:18

      Dear sir,

      I have been trying to use the package in my android project, I use Native Development Kit (NDK) with OpneCV 2.4.1 library under windows.

      My question is how to use the packages: Eigen & UMFPack in my project, given that U said it’s build only under Ubunto.

      Thanks in advance,”

  5. dai より:

    I have met with the same problem,then I found the solution.Maybe it’s your eigen pack is too old,download the newest editon will help solve the problem.Since I look in to the website below thttp://reference.mrpt.org/svn/_umf_pack_support_source.html#l00026,found out the statement”struct UmfPack {};”was what missing in my umfpacksupport.

  6. Chris より:

    Hi, thanks for the code. Works great! I propose two small improvements:
    1. the output image is 2 pixel too wide and two pixel to tall
    2. you aren’t be able to use a mask that goes up until the border of the input image, it has to be at least 4pixels away from the image border.

    To account for this, just change lines 224 to 242 to this:

    cv::Rect mask_roi2(tl, br+cv::Point(4,4));
    cv::Mat _srcUp, _targetUp, _maskUp, _dstUp;
    cv::copyMakeBorder(_src, _srcUp, 2,2,2,2, cv::BORDER_REPLICATE);
    cv::copyMakeBorder(_target, _targetUp, 2,2,2,2, cv::BORDER_REPLICATE);
    cv::copyMakeBorder(_mask, _maskUp, 1,1,1,1, cv::BORDER_CONSTANT);

    // allocate destination image
    _dstUp = _targetUp.clone();
    _dst = cv::Mat(_dstUp, cv::Rect(2,2,_dstUp.cols-4, _dstUp.rows-4));

    mask_roi1 = cv::Rect(tl, br+cv::Point(2,2));
    mask1 = cv::Mat(_maskUp, mask_roi1);
    target1 = cv::Mat(_targetUp, mask_roi2+offset);
    dst1 = cv::Mat(_dstUp, mask_roi2+offset);

    cv::Mat src(_srcUp, mask_roi2);
    cv::Mat target(_targetUp, mask_roi2+offset);
    cv::Mat dst(_dstUp, mask_roi2+offset);

    How the output image is the right size and the mask can be anywhere in the image.

    Cheers,
    Chris

    • dai より:

      Hi,can you help me solve the complie problem?All error occurs in the statement:
      cv::Mat pdx_src[ch], pdy_src[ch], pdx_target[ch], pdy_target[ch];
      I found out that the variable ch has not been assigned with a value.

      How you deal with this?please replay.thanks in advance.

      • dai より:

        I’ve ask a really stupid question.It is all because my complier does not support using a variable as the array dimentions.
        But I also want to know the mask image is generated by what method.

  7. mercadee より:

    You will need to composite two (or more) images to make a panorama. It is your choice on whether to incrementally warp one image to another or warp all the images to one image. It is also your choice on how to composite. You can simply do over, averaging, feathering, GraphCut seams, Poisson image blending or any other technique you feel like. You will not be penalized for doing simple composite techniques, but your results may not look as good.

  8. Gil Megidish より:

    Awesome work! What is the license for this code?

  9. lu より:

    oh to add on to the previous comment, in poisson_blending.hpp i changed a line
    in the beginning from

    #include
    to
    #include

    Because in the new version of Eigen package, there is change of directory where the file “UmfPackSupport” is located.

    Please kindly help.

    Thanks again.

  10. lu より:

    in the header file, this line

    // solve sparse linear system
    solve(A, b, u);

    caused a compile problem
    poisson_blending.hpp(297): error C2039: ‘UmfPack’ : is not a member of ‘Eigen’

    Could you help please. thanks!

コメントをどうぞ