Ccorrelation_opticalFlow_integrated_block.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_BLOCK
#define CCORRELATION_OPTICALFLOW_INTEGRATED_BLOCK



#include "Ccorrelation_opticalFlow_integrated.h"

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

  public:
    
    Ccorrelation_opticalFlow_integrated_block() : Ccorrelation_opticalFlow_integrated<T,Timg>()
    {// constructor
      
      this->class_name = "Correlation : Optical Flow Integrated per Block";
      
    }
    
    virtual void init(CParameterNetCDF &fp)
    {// init from parameter file
      
      Ccorrelation_opticalFlow_integrated<T,Timg>::init(fp);
      
      std::string attribute_name = "median_filter";
      int error = fp.loadAttribute(attribute_name,this->m_medianFilter);
      if (this->m_verbose){printf("\t\t median filter over following domain -> %i\n",this->m_medianFilter);}

    }

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

        Ke.assign(2);
        Ke[0].assign(oMesh.num_mod(),oMesh.num_mod(),1,1, 0);
        Ke[1].assign(1,oMesh.num_mod(),1,1, 0);
   
    }
    
    virtual void reset(const int &thread, CImgList<T> &Ke)
    {//reset tensors
    
      cimglist_for(Ke, var)
      {
        Ke[var].fill(0);
      }
      
    }
    
    virtual void e_ref(const Cimage<Timg> &oImage, const Cmesh<T,Timg> &oMesh, const int &elem, std::vector<T> &ref)
    {// omp_get_num_threads()store position of the centroid 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)
        {
          //centroid coordonates in physical unit
          ref[dim] = oMesh.m_node[dim][oMesh.m_connect[elem].front()]*oImage.m_pix[dim] ;
        }
        else
        {//apply scaling to reference position provided by user.
          ref[dim] = this->m_ref[dim]*oImage.m_Rpxmm;
        }
        
      }
        
    }
  
    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)
    {//evaluate the shape function at pix (x, y, z)
      
      std::vector<T> coor;coor.assign(3,0);
      for (int i=0; i<gCoor.size(); i++)
      {
        coor[i] = gCoor[i]*oImage.m_pix[i] - ref[i];
      }
      
      (oField.m_pShape)->exec(list, coor, N);
      
    }
    
    virtual void tensors(const int &thread, const CImgList<T> &Ve, const T &residue, CImgList<T> &Ke)
    {//eval stiffness matrix & force vector at element
      
      cimg_forZ(Ke[0], dim)
      {//in the case of OFI dim = 0
        cimg_forX(Ke[0],modi)
        {//for every modes or nodes
          cimg_forY(Ke[0],modj)
          {
            Ke[0](modi,modj,dim) += Ve[dim][modi] * Ve[dim][modj];//Ke
          }
          Ke[1](0, modi,dim) += Ve[dim][modi] * residue;//Fe
        }
      }
      
    }
 
    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)
    {//calcul K

      #pragma omp single
      {//init
        current_res.fill(0);
        p.fill(0);
      }

      int thread(0), enter(1);

      #pragma omp for schedule(runtime)
      cimglist_for(oMesh.m_connect, elem)
      {//for every elements
        
        if (oMesh.cropList(oImage, elem, list, offset))
        {//if the list is not entirely masked

          //check the thread number
          thread = omp_get_thread_num();
      
          //look for the reference system of the element
          e_ref(oImage, oMesh, elem, ref);

          //reset tensors
          reset(thread,Ke);

          for (int pix=0; pix<list.size(); pix++)
          {//eval Ke over pix

            //define global coordonates : find back "x,y,z" triplet
            list.contains(list[pix],lCoor[0],lCoor[1],lCoor[2]);
            for (int i=0; i<gCoor.size(); i++){gCoor[i]=lCoor[i]+offset[i];}
            
            //if pix is in img bounds
            enter = oImage.m_img[0].contains(oImage.m_img[0](gCoor[0],gCoor[1],gCoor[2]));
            //if pix is unmasked
            if (!oImage.m_mask.is_empty() and enter){enter = oImage.m_curMask(gCoor[0], gCoor[1], gCoor[2]);}
            
            p[thread]++;

            if (enter!=0)
            {//if pixels are unmasked and in bounds
              
              // Evaluates grad[dim] at (x,y,z) in physical unit
              Ccorrelation_opticalFlow<T,Timg>::gradient(oImage, gCoor, grad);

              // Evaluates N[dim][mode] at (x,y,z) in physical unit
              evalShape(list, oImage, oField, ref, lCoor, gCoor, N);

              // Evaluates Ve[dim][mode] at (x,y,z) in physical unit
              this->evalGN(grad,N, Ve);

              // Evaluates U[dim] in physical unit
              this->evalU(solution, oMesh, elem, 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
              tensors(thread, Ve, residue, Ke);

             }//endif unmasked or out-of-bound
        
          }//for implied barrier
         
          Ke[1].solve(Ke[0]);
          cimg_forY(Ke[1],y){solution[y][elem] += Ke[1][y];}
          
        }
        
      }
     
      #pragma omp single
      {//sum of square residues over threads and normalization and store residue of previous increment
        oField.m_meanRes.push_back(std::sqrt(current_res.sum())/p.sum()*oImage.m_curImg[0].size()*100/oImage.m_res);
      }

    }
 
    virtual void exec(const Cimage<Timg> &oImage, const Cmesh<T,Timg> &oMesh, Cfield<T,Timg> &oField)
    {// exec optical flow correlation
    
      //-----------------------------------------------------------------------
      //monoprocessing asignments
      this->_3D = oMesh._3D;

      int iter(0); //number of iterations
      int converged(0); //convergence parameter
      T previous_res(1e10); //residue at previous step
      
      CImg<int> p(this->m_threadNB,1,1,1,0);
      CImgList<T> solution(oField.m_mode);//intermediate solution vector
      CImg<T> current_res(this->m_threadNB,1,1,1,0);//global mean residue
      //-----------------------------------------------------------------------

      #pragma omp parallel
      {//beginning of parallel section
  
        //-----------------------------------------------------------------------
        //multithread initializations
        //for elements
        CImgList<T> Ke;//global and local tensors for linear system K[0].U=K[1] with K[0] = sum_(V.Vt) and K[1] = sum_((Img1(X)-Img2(X-U)).V)
        assignG(oMesh, Ke);//global assignments
        CImg<Timg> list; //elementary pixel list
        
        //for pixels
        CImgList<T> Ve; //shape function to image gradient product at pix
        CImgList<T> N(oMesh.num_var(), oMesh.num_mod(),1,1,1, 0);//Shape function at pix
        T residue(0), deformed(0); //residue and corrected GL at pix
        std::vector<T> ref, offset, lCoor, gCoor, grad, U;
        ref.assign(3,0);offset.assign(3,0);lCoor.assign(3,0); gCoor.assign(3,0); U.assign(3,0); grad.assign(oMesh.num_var(),0);
        
        this->assignL(oMesh, Ve);
      //-----------------------------------------------------------------------

        do
        {//while the convergency velocity is higher than a certain level
        
          metric(oImage,oMesh,solution,oField,N,Ve,Ke,ref,offset,lCoor,gCoor,U,grad,residue,deformed,list, p,current_res);

          #pragma omp single
          {
            if (this->m_medianFilter){this->regularize(oMesh, oImage, solution);}
            test(solution, this->m_epsilon, oField, previous_res, iter, converged);
          }
          
        }
        while(!converged);
      
        if (oImage.m_storeActiv)
        {Ccorrelation_opticalFlow<T,Timg>::export_res(oImage,oMesh,oField, N,ref,offset,lCoor,gCoor,U,residue,deformed,list);}
        
      }// pragma omp parallel
      
    }
  
    
    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)
        {
          #pragma omp critical
          {
            oField.m_field(dim, gCoor[0],gCoor[1],gCoor[2]) = U[dim];
          }
        }
      }

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

      if (oImage.m_storeRes)
      {//store imRes a iter - 1
          #pragma omp critical
          {
            oField.m_res(gCoor[0],gCoor[1],gCoor[2]) = residue;
          }
      }
      
    }
   
    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);
      
    } 

};

#endif