Cfield_opticalFlow_fem.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_FEM
#define CFIELD_OPTICALFLOW_FEM


//-----------------------------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_fem: public Cfield_opticalFlow<T,Timg>
{
  
  public:

    
    Cfield_opticalFlow_fem() : Cfield_opticalFlow<T,Timg>()
    {
      this->class_name = "Field : Optical Flow Finite Element Method";
      
    }

    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(oMesh.num_var(), oMesh.grid_dims(0),oMesh.grid_dims(1),oMesh.grid_dims(2), 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)
      {//basic linear interpolation from previous the current grid of modes [U, V, W]

        Field2Field(modes_prev, modes_cur);
        
      }
      else if((oShape_prev.m_name).compare(0,2,"IC") == 0)
      {//basic linear interpolation from previous the current grid of modes [U, V, W]
      
        //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];}
        
        Field2Field(U, modes_cur);
     
      }
      else if((oShape_prev.m_name).compare(0,4,"OFIB") == 0)
      {//basic linear interpolation from previous the current grid and building of U, V and W at previous element centroid
        
        CImgList<T> U(oMesh.num_var(), modes_prev[0].width(), modes_prev[0].height(), modes_prev[0].depth(), 1, 0);
        OFIBU(oMesh, oImage, oShape_prev, ref_prev, nodes_prev, modes_prev, U);
        
        Field2Field(U, modes_cur);
        
      }
      else if((oShape_prev.m_name).compare(0,3,"OFI") == 0)
      {

        OFIU(oMesh, oImage, oShape_prev, ref_prev, nodes_prev, modes_prev, modes_cur);
        
      }
      else{printf("unknown name for previous correlation method");}
      
    }
    
    void OFIU(const Cmesh<T,Timg> &oMesh, const Cimage<Timg> &oImage, CshapeFunction<T,Timg> &oShape_prev, std::vector<T> ref_prev, const CImgList<T> &nodes_prev, CImgList<T> &modes_prev, CImgList<T> &modes_cur)
    {//build U from OFI
     
      for (int dim=0; dim<oMesh.num_var(); dim++)
      {
        if (ref_prev[dim] == -1)
        {
          ref_prev[dim] = (oImage.cur_dim(dim)*oImage.m_pix[dim])/2;
        }
      }
      
      std::vector<T> X;X.assign(3,0);
      CImgList<T> N(oMesh.num_var(), (oShape_prev).m_dof,1,1,1, 0);
      
      cimglist_for(modes_cur, dim)
      {//linear interpolation
        cimglist_for(modes_prev, mode)
        {//linear interpolation
          cimg_forXYZ(oMesh.m_node[dim],x,y,z)
          {
            
            X[0] = oMesh.m_node(0, x,y,z)*oImage.m_pix[0] - ref_prev[0];
            X[1] = oMesh.m_node(1, x,y,z)*oImage.m_pix[1] - ref_prev[1];
            if (this->_3D){X[2] = oMesh.m_node(2, x,y,z)*oImage.m_pix[2] - ref_prev[2];}
            
            oShape_prev.exec(modes_cur[dim], X, N);
            modes_cur(dim, x,y,z) += modes_prev[mode][0]*N[dim][mode];
          }
        }
      }

    }
    
    void OFIBU(const Cmesh<T,Timg> &oMesh, const Cimage<Timg> &oImage, CshapeFunction<T,Timg> &oShape_prev, std::vector<T> ref_prev, const CImgList<T> &nodes_prev, const CImgList<T> &modes_prev, CImgList<T> &U)
    {// build U from OFIB
      
      std::vector<T> X; X.assign(3,0);
      CImgList<T> N(oMesh.num_var(), (oShape_prev).m_dof,1,1,1, 0);
    
      cimglist_for(U, dim)
      {//linear interpolation
        cimglist_for(modes_prev, mode)
        {//linear interpolation
          cimg_forXYZ(nodes_prev[0],x,y,z)
          {
            
            if (ref_prev[0] != -1){X[0] = nodes_prev(0,x,y,z) - ref_prev[0];}
            if (ref_prev[1] != -1){X[1] = nodes_prev(1,x,y,z) - ref_prev[1];}
            if (ref_prev[2] != -1 and nodes_prev.size()==3){X[2] = nodes_prev(2,x,y,z) - ref_prev[2];}
            
            oShape_prev.exec(U[dim], X, N);
            U(dim, x,y,z) += modes_prev(mode, x,y,z)*N[dim][mode];
          }
        }
      }

    }

    void Field2Field(const CImgList<T> &U, CImgList<T> &modes_cur)
    {
     
      cimglist_for(modes_cur, mode)
      {
        modes_cur[mode] = U[mode].get_resize(this->grid_dims(0),this->grid_dims(1),this->grid_dims(2),this->grid_dims(3),3);
      }
      
    }
    
    
    
};

#endif