Solve Poisson Equation with FFT

gistfile1.cpp

// Get CImg.h from http://sourceforge.net/projects/cimg/
#include "CImg.h"
using namespace cimg_library;

int main() {
  const CImg<> img("image.jpg");

  // Estimate gradient (uses forward finite difference scheme).
  const CImgList<> gradient = img.get_gradient("xy",1);

  // Remember average value for each channel (needed for normalizing the reconstruction).
  const CImg<> average = img.get_resize(1,1,1,3,2);

  // Image estimation from gradient of 'img' starts here.
  //-----------------------------------------------------

  // Step 1.  Estimate divergence of the gradient field (use backward finite difference scheme).
  const CImg<> divergence = gradient[0].get_gradient("x",-1)[0] + gradient[1].get_gradient("y",-1)[0];

  // Step 2. Invert Laplacian operator using FFT.
  CImgList<> FFT = divergence.get_FFT();
  CImg<> factor(FFT[0].width(),FFT[0].height());
  cimg_forXY(factor,x,y) factor(x,y) = -(4-2*std::cos(2*x*cimg::PI/factor.width())-2*std::cos(2*y*cimg::PI/factor.height()));
  factor(0,0) = 1;
  FFT[0].div(factor);
  FFT[1].div(factor);
  CImg<> res0 = FFT.get_FFT(true)[0];
  res0+=average.get_resize(res0.width(),res0.height());
  res0.cut(0,255);
  res0.display("Reconstruction");

  return 0;
}