Ccorrelation_opticalFlow_integrated.h

Go to the documentation of this file.

/**********************************************************************
 * Copyright (C) 2012, The YaDICs Project Developers. 
 * See the COPYRIGHT file at the top-level directory of this distribution ./COPYRIGHT.
 * See ./COPYING file for copying and redistribution conditions.
 * 
 * This file is part of YaDICs.
 * YaDICs is free software: you can redistribute it and/or modify
 * it under the terms of the GNU General Public License as published by
 * the Free Software Foundation, either version 3 of the License, or
 * (at your option) any later version.
 * 
 * YaDICs is distributed in the hope that it will be useful,
 * but WITHOUT ANY WARRANTY; without even the implied warranty of
 * MERCHANTABILITY or FITNESS FOR A PARTICULAR PURPOSE.  See the
 * GNU General Public License for more details.
 * 
 * You should have received a copy of the GNU General Public License
 * along with YaDICs.  If not, see <http://www.gnu.org/licenses/>.
 * 
 * Information about how to use the software are provided at http://yadic.univ-lille1.fr/
 **********************************************************************/


#ifndef CCORRELATION_OPTICALFLOW_INTEGRATED
#define CCORRELATION_OPTICALFLOW_INTEGRATED


#include "Ccorrelation_opticalFlow_fem.h"



template<typename T,typename Timg>
class Ccorrelation_opticalFlow_integrated: public Ccorrelation_opticalFlow_fem<T,Timg>
{
  
  public:

    
    Ccorrelation_opticalFlow_integrated() : Ccorrelation_opticalFlow_fem<T,Timg>()
    {// constructor
      
      this->class_name = "Correlation : Optical Flow Integrated";

    }

    virtual void init(CParameterNetCDF &fp)
    {// init from parameter file
      
      Ccorrelation_opticalFlow<T,Timg>::init(fp);
       
      this->m_medianFilter = 0;

      this->m_ref.assign(3,-1);
      
      std::string attribute_name("refX");
      int error = fp.loadAttribute(attribute_name,this->m_ref[0]);
      
      attribute_name = "refY";
      error = fp.loadAttribute(attribute_name,this->m_ref[1]);
      
      attribute_name = "refZ";
      error = fp.loadAttribute(attribute_name,this->m_ref[2]);

      if (this->m_verbose){printf("\tReference system computed by user: [%f,%f,%f]\n", this->m_ref[0], this->m_ref[1], this->m_ref[2]);}

    }

    virtual void assignG(const Cmesh<T,Timg> &oMesh, CImgList<T> &K, CImgList<T> &Ke)
    {//assign approriate size for Ke and Fe tensors

        Ccorrelation_opticalFlow<T,Timg>::assignG(oMesh, K, Ke);
        
        K[0].assign(oMesh.num_mod(),oMesh.num_mod(),1,1, 0);
        K[1].assign(1,oMesh.num_mod(),1,1, 0);
        
        Ke[0].assign(oMesh.num_mod(),oMesh.num_mod(),1,this->m_threadNB, 0);
        Ke[1].assign(1,oMesh.num_mod(),1,this->m_threadNB, 0);
   
    }
    
    virtual void assignL(const Cmesh<T,Timg> &oMesh, CImgList<T> &Ve)
    {//assign approriate size for Ve = Grad[dim].N[dim][mode] product at pix

        Ve.assign(1, oMesh.num_mod(),1,1,1, 0);
   
    }
    
    virtual void reset(const int &thread, CImgList<T> &Ke)
    {//reset tensors
    
      cimglist_for(Ke, var)
      {
        Ke[var].fill(0);
      }
      
    }
    
    virtual void evalGN(const std::vector<T> &grad,const CImgList<T> &N, CImgList<T> &Ve)
    {//eval V[dim][mode] = grad[dim].N[dim][mode] product
      
      Ve[0].fill(0);
      cimglist_for(N, dim)
      {
        cimg_forX(N[dim], mode)
        {
          Ve[0][mode] += grad[dim] * N[dim][mode];
        }
      }
      
    } 

    virtual void e_ref(const Cimage<Timg> &oImage, const Cmesh<T,Timg> &oMesh, const int &elem, std::vector<T> &ref)
    {// store position of upper left corner as element position reference "0"
     
        for (int dim=0; dim<oMesh.num_var(); dim++)
        {// user reference or image centroid in physical unit
          if (this->m_ref[dim] == -1)
          {
            ref[dim] = (oImage.cur_dim(dim)*oImage.m_pix[dim])/2.;
          }
          else
          {//apply scaling to reference position provided by user.
            ref[dim] = this->m_ref[dim]*oImage.m_Rpxmm;
          }
        }
        
    }
    
    virtual void evalU(const CImgList<T> &solution, const Cmesh<T,Timg> &oMesh, const int &elem, const CImgList<T> &N, std::vector<T> &U)
    {//eval U[dim] at pix
    
      U[2] = 0;//init required if dimension is lower than 3
      cimglist_for(N, dim)
      {
        
        U[dim] = 0;
        
        cimg_forX(N[0], mode)
        {
          //for OFI (elem = 0) and OFIB (elem = centroid)
          U[dim] += solution[mode][elem] * N[dim][mode];
        }

      }
      
      
    }
   
    virtual void gradient(const Cimage<Timg> &oImage, const std::vector<T> &x, std::vector<T> &grad)
    {// image gradient calcul at "x, y, z" image location
  
      if (x[0] != 0 and x[0] != oImage.cur_dim(0)-1)
      {
        grad[0] = ((T)oImage.m_curImg[0](x[0]+1,x[1],x[2], 0) - (T)oImage.m_curImg[0](x[0]-1,x[1],x[2], 0))/(2.*oImage.m_pix[0]);
      }
      else{grad[0] = 0;}
      
      if (x[1] != 0 and x[1] != oImage.cur_dim(1)-1)
      {     
        grad[1] = ((T)oImage.m_curImg[0](x[0],x[1]+1,x[2], 0) - (T)oImage.m_curImg[0](x[0],x[1]-1,x[2], 0))/(2.*oImage.m_pix[1]);
      }
      else{grad[1] = 0;}
      
      if (this->_3D)
      {
        if (x[2] != 0 and x[2] != oImage.cur_dim(2)-1)
        {
          grad[2] = ((T)oImage.m_curImg[0](x[0],x[1],x[2]+1, 0) - (T)oImage.m_curImg[0](x[0],x[1],x[2]-1, 0))/(2.*oImage.m_pix[2]);
        }
        else{grad[2] = 0;}
      }

    }

    virtual void evalShape(const CImg<Timg> &list, const Cimage<Timg> &oImage, const Cfield<T,Timg> &oField, const std::vector<T> &ref, const std::vector<T> &lCoor, const std::vector<T> &gCoor, CImgList<T> &N)
    { 
      
      std::vector<T> coor;coor.assign(3,0);
      for (int i=0; i<lCoor.size(); i++)
      {
        coor[i] = gCoor[i]*oImage.m_pix[i] - ref[i];
      }
      
      
//       std::cout<<"evalShape : "<<list.width()<<"; "<<list.height()<<", "<<list.depth()<<"\n";
      (oField.m_pShape)->exec(list, coor, N);
      
    }
   
    virtual void assembly(const int &thread, const Cmesh<T,Timg> &oMesh, const int &elem, const CImgList<T> &Ke, CImgList<T> &K)
    {//perform assembly of K and F

      //reduction
      cimglist_for(Ke,var)
      {
        cimg_forXYC(Ke[var],x,y,c)
        {
          K[var](x,y) += Ke[var](x,y,0,c);
        }
        
      }
      
    }
   
    virtual void metric(const Cimage<Timg> &oImage, const Cmesh<T,Timg> &oMesh, CImgList<T> &solution, Cfield<T,Timg> &oField, CImgList<T> &N, CImgList<T> &Ve, CImgList<T> &Ke, std::vector<T> &ref, std::vector<T> &offset, std::vector<T> &lCoor, std::vector<T> &gCoor, std::vector<T> &U, std::vector<T> &grad,  T &residue, T &deformed, CImg<Timg> &list, CImg<int> &p, CImg<T> &current_res, CImgList<T> &K)
    {//calcul K
    
      int thread(0), enter(1);
    
      #pragma omp single
      {//init
      
        current_res.fill(0);
        cimglist_for(K,var){K[var].fill(0);}

        //reset tensors
        reset(thread,Ke);
        
      }
      
      //look for the reference system of the element
      e_ref(oImage, oMesh, 0, ref);

      #pragma omp for schedule(runtime)
      for (int pix=0; pix<oImage.m_curImg[0].size(); pix++)
      {//eval Ke over pix

        //define global coordonates : find back "x,y,z" triplet
        oImage.m_curImg[0].contains(oImage.m_curImg[0][pix],lCoor[0],lCoor[1],lCoor[2]);
        for (int i=0; i<gCoor.size(); i++){gCoor[i]=lCoor[i];}
        
        enter = 1;
        //if pix is unmasked
        if (!oImage.m_mask.is_empty())
        {enter = oImage.m_curMask(gCoor[0], gCoor[1], gCoor[2]);}
        
        if (enter!=0)
        {//if pixels are unmasked and in bounds
            
          //check the thread number
          thread = omp_get_thread_num();
            
          // Evaluates grad[dim] at (x,y,z) in physical unit
          gradient(oImage, gCoor, grad);

          // Evaluates N[dim][mode] at (x,y,z) in physical unit
//        std::cout<<"metric : "<<oImage.m_curImg[0].width()<<"; "<<oImage.m_curImg[0].height()<<", "<<oImage.m_curImg[0].depth()<<"\n";
          evalShape(oImage.m_curImg[0], oImage, oField, ref, lCoor, gCoor, N);
          
          // Evaluates Ve[dim][mode] at (x,y,z) in physical unit
          evalGN(grad,N, Ve);

          // Evaluates U[dim] in physical unit
          evalU(solution, oMesh, 0, N, U);

          // Evaluates residue between ref and corrected deformed image at (x,y,z)
          this->evalRes_loop(oImage, U, gCoor, residue, deformed);

          //sum of square residues over pixels and elements
          current_res[thread] += std::pow(residue,2);

          //evaluate Ke
          this->tensors(thread, Ve, residue, Ke);
          
        }//endif unmasked or out-of-bound
        
      }//for implied barrier

      #pragma omp single
      {//sum of square residues over threads and normalization and store residue of previous increment
        assembly(thread, oMesh, 0, Ke, K);
        oField.m_meanRes.push_back(std::sqrt(current_res.sum())*100/oImage.m_res);
      }

    }

    virtual void solve(CImgList<T> &K, const Cmesh<T,Timg> &oMesh, CImgList<T> &solution)
    {// K.U = F
    
      K[1].solve(K[0]);
      
      cimg_forY(K[1],y){solution[y][0] += K[1][y];}
      
      if(this->m_verbose){for(int i=0; i<solution.size(); i++){std::cout<<solution[i][0]<<" ";}std::cout<<"\n";}
      
    }

    virtual void test(const CImgList<T> &solution, const double &epsilon, Cfield<T,Timg> &oField, T &previous_res, int &iter, int &converged)
    {
      
        Ccorrelation_opticalFlow<T,Timg>::test(solution, epsilon, oField, previous_res, iter, converged);
        if(converged)
        {
          std::cout<<"Solution vector : [ ";
          for(int i=0; i<oField.m_mode.size(); i++){std::cout<<oField.m_mode[i][0]<<" ";}std::cout<<"]\n\n";
        }
      
    }
   
    virtual void export_res(const Cimage<Timg> &oImage, const Cmesh<T,Timg> &oMesh, Cfield<T,Timg> &oField, CImgList<T> &N, std::vector<T> &ref, std::vector<T> &offset, std::vector<T> &lCoor, std::vector<T> &gCoor, std::vector<T> &U,  T &residue, T &deformed, CImg<Timg> &list)
    {//actualize residue

      int thread(0),enter(1);

      //look for the reference system of the element
      e_ref(oImage, oMesh, 0, ref);
          
      #pragma omp for schedule(runtime)
      for (int pix=0; pix<oImage.m_curImg[0].size(); pix++)
      {//eval Ke over pix
    
        //define global coordonates : find back "x,y,z" triplet
        oImage.m_curImg[0].contains(oImage.m_curImg[0][pix],lCoor[0],lCoor[1],lCoor[2]);
        for (int i=0; i<gCoor.size(); i++){gCoor[i]=lCoor[i];}

        enter = 1;
        //if pix is unmasked
        if (!oImage.m_mask.is_empty()){enter = oImage.m_curMask(gCoor[0], gCoor[1], gCoor[2]);}
        
        if (enter!=0)
        {//if pixels are unmasked and in bounds
      
          //check the thread number
          thread = omp_get_thread_num();
      
          // Evaluates N[dim][mode] at (x,y,z) in physical unit
          evalShape(oImage.m_curImg[0], oImage, oField, ref, lCoor, gCoor, N);

          // Evaluates U[dim] in physical unit
          evalU(oField.m_mode, oMesh, 0, N, U);
      
          // Evaluates residue between ref and corrected deformed image at (x,y,z)
          this->evalRes_final(oImage, U, gCoor, residue, deformed);
      
          //export data
          store(oField,oImage,gCoor,U,deformed,residue);
          
        }//if
      
      }//for
    
    }
   
    virtual void store(Cfield<T,Timg> &oField, const Cimage<Timg> &oImage, const std::vector<T> gCoor, const std::vector<T> &U,const T &deformed,const T &residue)
    {//store output data
      
      if (oImage.m_storeField)
      {//store field a iter - 1 to avoid storing the last diverging increment
        cimglist_for(oField.m_field, dim)
        {
          oField.m_field(dim, gCoor[0],gCoor[1],gCoor[2]) = U[dim];
        }
      }

      if (oImage.m_storeDef)
      {//store imDef a iter - 1
        oField.m_def(gCoor[0],gCoor[1],gCoor[2]) = deformed;
      }

      if (oImage.m_storeRes)
      {//store imRes a iter - 1
        oField.m_res(gCoor[0],gCoor[1],gCoor[2]) = residue;
      }
      
    }
  
  
  
};
   
 
#endif