PTdecode/CImg-1.3.0/examples/image_registration.cpp

/*
 #
 #  File        : image_registration.cpp
 #                ( C++ source file )
 #
 #  Description : Compute a motion field between two images,
 #                with a multiscale and variational algorithm.
 #                This file is a part of the CImg Library project.
 #                ( http://cimg.sourceforge.net )
 #
 #  Copyright   : David Tschumperle
 #                ( http://www.greyc.ensicaen.fr/~dtschump/ )
 #
 #  License     : CeCILL v2.0
 #                ( http://www.cecill.info/licences/Licence_CeCILL_V2-en.html )
 #
 #  This software is governed by the CeCILL  license under French law and
 #  abiding by the rules of distribution of free software.  You can  use,
 #  modify and/ or redistribute the software under the terms of the CeCILL
 #  license as circulated by CEA, CNRS and INRIA at the following URL
 #  "http://www.cecill.info".
 #
 #  As a counterpart to the access to the source code and  rights to copy,
 #  modify and redistribute granted by the license, users are provided only
 #  with a limited warranty  and the software's author,  the holder of the
 #  economic rights,  and the successive licensors  have only  limited
 #  liability.
 #
 #  In this respect, the user's attention is drawn to the risks associated
 #  with loading,  using,  modifying and/or developing or reproducing the
 #  software by the user in light of its specific status of free software,
 #  that may mean  that it is complicated to manipulate,  and  that  also
 #  therefore means  that it is reserved for developers  and  experienced
 #  professionals having in-depth computer knowledge. Users are therefore
 #  encouraged to load and test the software's suitability as regards their
 #  requirements in conditions enabling the security of their systems and/or
 #  data to be ensured and,  more generally, to use and operate it in the
 #   same conditions as regards security.
 #
 #  The fact that you are presently reading this means that you have had
 #  knowledge of the CeCILL license and that you accept its terms.
 #
*/

#include "CImg.h"
using namespace cimg_library;

// The lines below are necessary when using a non-standard compiler as visualcpp6.
#ifdef cimg_use_visualcpp6
#define std
#endif
#ifdef min
#undef min
#undef max
#endif

#ifndef cimg_imagepath
#define cimg_imagepath "img/"
#endif

// animate_warp() : Create warping animation from two images and a motion field
//----------------
void animate_warp(const CImg<unsigned char>& src, const CImg<unsigned char>& dest, const CImg<>& u,
                  const bool morph, const bool imode, const char *filename,int nb, CImgDisplay& disp) {
  CImg<unsigned char> visu = CImgList<unsigned char>(src,dest,src).get_append('x'), warp(src);
  float t=0;
  for (unsigned int iter=0; !disp || (!disp.is_closed && !disp.is_keyQ); iter++) {
    if (morph) cimg_forXYV(warp,x,y,k) {
      const float dx = u(x,y,0), dy = u(x,y,1),
        I1 = (float)src.linear_atXY(x-t*dx, y-t*dy, k),
        I2 = (float)dest.linear_atXY(x+(1-t)*dx,y+(1-t)*dy,k);
      warp(x,y,k) = (unsigned char)((1-t)*I1 + t*I2);
    } else cimg_forXYV(warp,x,y,k) {
      const float dx = u(x,y,0), dy = u(x,y,1), I1 = (float)src.linear_atXY(x-t*dx, y-t*dy, 0,k);
      warp(x,y,k) = (unsigned char)I1;
    }
    if (disp) visu.draw_image(2*src.dimx(),warp).display(disp.resize().wait(30));
    if (filename && *filename && (imode || (int)iter<nb)) {
      std::fprintf(stderr,"\r  > frame %d           ",iter);
      warp.save(filename,iter);
    }
    t+=1.0f/nb;
    if (t<0) { t=0; nb=-nb; }
    if (t>1) { t=1; nb=-nb; if (filename && *filename) std::exit(0); }
  }
}

// get_warp() : Return the image src warped by the motion field u.
//------------
template<typename T> CImg<T> getwarp(const CImg<T>& src, const CImg<>& u) {
  CImg<T> warp(src);
  cimg_forXY(warp,x,y) warp(x,y) = (T)src.linear_atXY(x - u(x,y,0), y - u(x,y,1));
  return warp;
}

// optmonoflow() : Register images for one scale ( semi-implicite PDE scheme ) between I2->I1
//---------------
CImg<> optmonoflow(const CImg<>& I1, const CImg<>& I2, const CImg<>& u0,
                   const float smooth, const float precision, CImgDisplay& disp) {

  CImg<> u = u0.get_resize(I1.dimx(),I1.dimy(),1,2,3),dI(u);
  CImg_3x3(I,float);
  float dt=2,E=1e20f;

  // compute first derivatives of I2
  cimg_for3x3(I2,x,y,0,0,I) {
    dI(x,y,0) = 0.5f*(Inc-Ipc);
    dI(x,y,1) = 0.5f*(Icn-Icp);
  }

  // Main PDE iteration
  for (unsigned int iter=0; iter<100000; iter++) {
    std::fprintf(stderr,"\r- Iteration %d - E = %g",iter,E); std::fflush(stderr);
    const float Eold = E;
    E = 0;
    cimg_for3XY(u,x,y) {
      const float
        X = x + u(x,y,0),
        Y = y + u(x,y,1),
        deltaI = (float)(I2.linear_atXY(X,Y) - I1(x,y));
      float tmpf = 0;
      cimg_forV(u,k) {
        const float
          ux  = 0.5f*(u(_n1x,y,k)-u(_p1x,y,k)),
          uy  = 0.5f*(u(x,_n1y,k)-u(x,_p1y,k));
        u(x,y,k) = (float)( u(x,y,k) +
                            dt*(
                                -deltaI*dI.linear_atXY(X,Y,k) +
                                smooth* ( u(_n1x,y,k) + u(_p1x,y,k) + u(x,_n1y,k) + u(x,_p1y,k) )
                                )
                            )/(1+4*smooth*dt);
        tmpf += ux*ux + uy*uy;
      }
      E += deltaI*deltaI + smooth * tmpf;
    }
    if (cimg::abs(Eold-E)<precision) break;
    if (Eold<E) dt*=0.5;
    if (disp) disp.resize();
    if (disp && disp.is_closed) std::exit(0);
    if (disp && !(iter%300)) {
      const unsigned char white = 255;
      CImg<unsigned char> tmp = getwarp(I1,u).normalize(0,200);
      tmp.resize(disp.dimx(),disp.dimy()).draw_quiver(u,&white,0.7f,15,-14,0).display(disp);
    }
  }
  return u;
}

// optflow() : multiscale version of the image registration algorithm
//-----------
CImg<> optflow(const CImg<>& xsrc, const CImg<>& xdest,
               const float smooth, const float precision, const unsigned int pnb_scale, CImgDisplay& disp) {
  const CImg<>
    src  = xsrc.get_pointwise_norm(1).resize(xdest.dimx(),xdest.dimy(),1,1,3).normalize(0,1),
    dest = xdest.get_pointwise_norm(1).resize(xdest.dimx(),xdest.dimy(),1,1,3).normalize(0,1);
  CImg<> u = CImg<>(src.dimx(),src.dimy(),1,2).fill(0);

  const unsigned int nb_scale = pnb_scale>0?pnb_scale:(unsigned int)(2*std::log((double)(cimg::max(src.dimx(),src.dimy()))));
  for (int scale=nb_scale-1; scale>=0; scale--) {
    const CImg<> I1 = src.get_resize((int)(src.dimx()/std::pow(1.5,scale)), (int)(src.dimy()/std::pow(1.5,scale)) ,1,1,3);
    const CImg<> I2 = dest.get_resize((int)(src.dimx()/std::pow(1.5,scale)), (int)(src.dimy()/std::pow(1.5,scale)) ,1,1,3);
    std::fprintf(stderr," * Scale %d\n",scale);
    u*=1.5;
    u = optmonoflow(I1,I2,u,smooth,(float)(precision/std::pow(2.25,1+scale)),disp);
    std::fprintf(stderr,"\n");
  }
  return u;
}

/*------------------------

  Main function

  ------------------------*/

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

  // Read command line parameters
  cimg_usage("Compute an optical flow between two 2D images, and create a warped animation");
  const char
    *name_i1   = cimg_option("-i",cimg_imagepath "sh0r.pgm","Input Image 1 (Destination)"),
    *name_i2   = cimg_option("-i2",cimg_imagepath "sh1r.pgm","Input Image 2 (Source)"),
    *name_o    = cimg_option("-o",(const char*)NULL,"Output 2D flow (inrimage)"),
    *name_seq  = cimg_option("-o2",(const char*)NULL,"Output Warping Sequence");
  const float
    smooth    = cimg_option("-s",0.1f,"Flow Smoothness"),
    precision = cimg_option("-p",0.9f,"Convergence precision");
  const unsigned int
    nb        = cimg_option("-n",40,"Number of warped frames"),
    nbscale   = cimg_option("-scale",0,"Number of scales (0=auto)");
  const bool
    normalize = cimg_option("-equalize",true,"Histogram normalization of the images"),
    morph     = cimg_option("-m",true,"Morphing mode"),
    imode     = cimg_option("-c",true,"Complete interpolation (or last frame is missing)"),
    dispflag = !cimg_option("-novisu",false,"Visualization");

  // Init images and display
  std::fprintf(stderr," - Init images.\n");
  const CImg<>
    src(name_i1),
    dest(CImg<>(name_i2).resize(src,3)),
    src_blur  = normalize?src.get_blur(0.5f).equalize(256):src.get_blur(0.5f),
    dest_blur = normalize?dest.get_blur(0.5f).equalize(256):dest.get_blur(0.5f);

  CImgDisplay disp;
  if (dispflag) {
    unsigned int w = src.dimx(), h = src.dimy();
    const unsigned int dmin = cimg::min(w,h), minsiz = 512;
    if (dmin<minsiz) { w=w*minsiz/dmin; h=h*minsiz/dmin; }
    const unsigned int dmax = cimg::max(w,h), maxsiz = 1024;
    if (dmax>maxsiz) { w=w*maxsiz/dmax; h=h*maxsiz/dmax; }
    disp.assign(w,h,"Estimated Motion",0);
  }

  // Run Motion estimation algorithm
  std::fprintf(stderr," - Compute optical flow.\n");
  const CImg<> u = optflow(src_blur,dest_blur,smooth,precision,nbscale,disp);
  if (name_o) u.save(name_o);
  u.print("Computed flow");

  // Do morphing animation
  std::fprintf(stderr," - Create warped animation.\n");
  CImgDisplay disp2;
  if (dispflag) {
    unsigned int w = src.dimx(), h = src.dimy();
    const unsigned int dmin = cimg::min(w,h), minsiz = 100;
    if (dmin<minsiz) { w=w*minsiz/dmin; h=h*minsiz/dmin; }
    const unsigned int dmax = cimg::max(w,h), maxsiz = 1024/3;
    if (dmax>maxsiz) { w=w*maxsiz/dmax; h=h*maxsiz/dmax; }
    disp2.assign(3*w,h,"Source/Destination images and Motion animation",0);
  }

  animate_warp(src.get_normalize(0,255),dest.get_normalize(0,255),u,morph,imode,name_seq,nb,disp2);

  std::exit(0);
  return 0;
}