Cfield_opticalFlow_integrated.h

/**********************************************************************
 * 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 CFIELD_OPTICALFLOW_INTEGRATED
#define CFIELD_OPTICALFLOW_INTEGRATED


//-----------------------------CIMG_NETCDF-------------------------------
#define cimg_plugin  "plugins/add_fileformat.h"
#define cimg_use_netcdf

#ifdef cimg_use_netcdf
 #include "../NetCDF.Tool/struct_parameter_NetCDF.h"
// #define     cimg_plugin1 "plugins/netcdf_file_format4CImg1.h"
// #define cimglist_plugin1 "plugins/netcdf_file_format4CImgList1.h"
 #define     cimg_plugin2 "plugins/netcdf_file_format4CImg2.h"
 #define cimglist_plugin2 "plugins/netcdf_file_format4CImgList2.h"
 #include "../CImg.Tool/CImg_NetCDF.h"
#endif
//-----------------------------CIMG_NETCDF-------------------------------

#include "Cfield_opticalFlow.h"

template<typename T, typename Timg>
class Cfield_opticalFlow_integrated: public Cfield_opticalFlow_fem<T,Timg>
{
  
  public:
    
    
    Cfield_opticalFlow_integrated() : Cfield_opticalFlow_fem<T,Timg>()
    {
      this->class_name = "Field : Optical Flow Integrated";
      
    }
    
    virtual void init(CParameterNetCDF &fp)
    {
      Cfield_opticalFlow<T,Timg>::init(fp);

      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 : [%f,%f,%f]\n", this->m_ref[0], this->m_ref[1], this->m_ref[2]);}
      
    }   

    
    virtual void exec(const Cimage<Timg> &oImage, const Cmesh<T,Timg> &oMesh, CImgList<T> &modes_prev, CshapeFunction<T,Timg> &oShape_prev, const std::vector<T> &ref_prev, const CImgList<T> &nodes_prev)
    {
      this->m_mode.assign(this->m_pShape->m_dof, 1,1,1,1, 0);

      Cfield<T,Timg>::exec(oImage, oMesh, modes_prev, oShape_prev, ref_prev, nodes_prev); 
      
    }
    
    
    virtual void project(const Cmesh<T,Timg> &oMesh, const Cimage<Timg> &oImage, CImgList<T> &modes_prev, CshapeFunction<T,Timg> &oShape_prev, const std::vector<T> &ref_prev, const CImgList<T> &nodes_prev, CImgList<T> &modes_cur)
    {//project modes on new base
      
      if ((oShape_prev.m_name).compare(0,5,"OFFEM") == 0)
      {//previous step was a FEM type
        
        //project U on current modes
        Field2Mod(oImage, oMesh, nodes_prev, modes_prev, modes_cur);
        
      }
      else if((oShape_prev.m_name).compare(0,2,"IC") == 0)
      {
        
        //build U map on previous mesh
        CImgList<T> U(oMesh.num_var(), modes_prev[0].width(), modes_prev[0].height(), modes_prev[0].depth(), 1, 0);
        cimglist_for(U, mode)
        {U[mode] = modes_prev[mode+2];}
        
        //project U on current modes
        Field2Mod(oImage, oMesh, nodes_prev, U, modes_cur);
        
      }
      else if((oShape_prev.m_name).compare(0,4,"OFIB") == 0)
      {
        
        //build U map on previous mesh
        CImgList<T> U(oMesh.num_var(), modes_prev[0].width(), modes_prev[0].height(), modes_prev[0].depth(), 1, 0);
        this->OFIBU(oMesh, oImage, oShape_prev, ref_prev, nodes_prev, modes_prev, U);
        
        //project U on current modes
        Field2Mod(oImage, oMesh, nodes_prev, U, modes_cur);
        
      }
      else if((oShape_prev.m_name).compare(0,3,"OFI") == 0)
      {
      
        if (modes_cur.size() == modes_prev.size())
        {modes_cur = modes_prev;}
        else
        {printf("warning : modes have been changes for multi-scale OFI strategy!!!");}
        
      }
      else{printf("unknown name for previous correlation method");}

      
    }

    
    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++)
        {
          if (ref[dim] == -1)
          {//set the reference position at image center location
            ref[dim] = (oImage.cur_dim(dim)*oImage.m_pix[dim])/2;
          }
          else
          {//apply scaling to reference position provided by user.
            ref[dim] = ref[dim]*oImage.m_Rpxmm;
          }
        }
        
    }
    
    virtual void Field2Mod(const Cimage<Timg> &oImage, const Cmesh<T,Timg> &oMesh, const CImgList<T> &nodes_prev, const CImgList<T> &U, CImgList<T> &Ti)
    {
        std::vector<T> X;X.assign(3,0);
        CImgList<T> N(oMesh.num_var(), Ti.size());
        CImg<T> K(Ti.size(), Ti.size(),1,1,0), F(1,Ti.size(),1,1,0);
        
        cimglist_for(oMesh.m_connect, elem)
        {
        
          e_ref(oImage, oMesh, elem, this->m_ref);
          
          K.fill(0);
          F.fill(0);
    
          cimglist_for(U, dim)
          {//for every dims
            
            cimg_forXYZ(nodes_prev[dim],x,y,z)
            {//for every nodes
              
              //physical coordonates of previous nodes considering the new reference system.
              X[0] = nodes_prev(0, x,y,z) - this->m_ref[0];
              X[1] = nodes_prev(1, x,y,z) - this->m_ref[1];
              if(nodes_prev.size() == 3){X[2] = nodes_prev(2, x,y,z) - this->m_ref[2];}
              
              //shape function at such location
              this->m_pShape->exec(U[dim],X, N);
              
              K += matrixProduct(N[dim].get_transpose(),N[dim]);
              F += (N[dim].get_transpose())*U(dim, x,y,z);
              
            }
            
          }

          F.solve(K);
          cimg_forY(F,y)
          {
            Ti[y][elem] = F[y];
          }
          
        }
        
      }
    
    
};

#endif