Ccorrelation_opticalFlow_fem_gradient.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_FEM_GRADIENT
#define CCORRELATION_OPTICALFLOW_FEM_GRADIENT


#include "Ccorrelation_opticalFlow_fem.h"



template<typename T,typename Timg>
class Ccorrelation_opticalFlow_fem_gradient: public Ccorrelation_opticalFlow_fem<T,Timg>
{
  
  public:
    
    double m_alpha;
    int m_iter_max;

    Ccorrelation_opticalFlow_fem_gradient() : Ccorrelation_opticalFlow_fem<T,Timg>()
    {// constructor
      
      this->class_name = "Correlation : Optical Flow Finite Element Method using Gradient strategy";

    }
    
    virtual void init(CParameterNetCDF &fp)
    {// init from parameter file
      
      Ccorrelation_opticalFlow_fem<T,Timg>::init(fp);
       
      std::string attribute_name = "alpha0";
      int error = fp.loadAttribute(attribute_name,m_alpha);
      if (this->m_verbose){std::cout<<"\t\t initial coefficient -> "<<m_alpha<<"\n";}
      
      attribute_name = "iter_max";
      m_iter_max = 100;
      error = fp.loadAttribute(attribute_name,m_iter_max);
      if (this->m_verbose){std::cout<<"\t\t max number of iteration -> "<<m_iter_max<<"\n";}
      
    }
    
    virtual void assignG(const Cmesh<T,Timg> &oMesh, CImgList<T> &K, CImgList<T> &Ke)
    {//assign approriate size for Ke and Fe tensors

        K.assign(1,1,oMesh.num_node()*oMesh.num_var(),1,1, 0);
        Ke.assign(1,1,oMesh.num_mod(),oMesh.num_var(),this->m_threadNB, 0);
   
    }
   
    virtual void tensors(const int &thread, const CImgList<T> &Ve, const T &residue, CImgList<T> &Ke)
    {//eval stiffness matrix & force vector at element
      
      cimg_forYZ(Ke[0],modi,dim)
      {//for every modes or nodes
        Ke[0](0, modi,dim,thread) += Ve[dim][modi] * residue;//Fe
      }
      
    }

    virtual void assembly(const int &thread, const Cmesh<T,Timg> &oMesh, const int &elem, const CImgList<T> &Ke, CImgList<T> &K)
    {//perform the assembly of K and F
    
      cimg_forYZ(Ke[0], modi,dim)
      {
        #pragma omp critical
        {
          K[0][oMesh.m_connect[elem][modi]*Ke[0].depth()+dim] += Ke[0](0,modi,dim,thread);
        }
        
      }
      
    }

    virtual void solve(CImgList<T> &K, const Cmesh<T,Timg> &oMesh, CImgList<T> &solution)
    {// K.U = F

      double norme = std::sqrt((K[0].get_sqr()).sum())/K[0].height();
      
      for (int node=0; node<oMesh.m_node[0].size(); node++)
      {//put solution vector in oField.m_mode[dim][node] framework
        for (int dim=0; dim<oMesh.num_var(); dim++)
        {
          solution[dim][node] += m_alpha*K[0][node*oMesh.num_var()+dim]/(norme*oMesh.m_win[0]*oMesh.m_win[1]*oMesh.m_win[2]);
        }
      }
     
    } 
    
    virtual void test(const CImgList<T> &solution, const double &epsilon, Cfield<T,Timg> &oField, T &previous_res, int &iter, int &converged)
    {//convergency test
    
      if (std::abs(oField.m_meanRes.back()-previous_res) >= epsilon and iter<m_iter_max)
      {
      
        if (oField.m_meanRes.back()<=previous_res)
        {
          m_alpha*=1.15;
        }
        else
        {
          m_alpha*=-0.5;
        }
        
        iter++;
        previous_res = oField.m_meanRes.back();
        oField.m_mode = solution;
        if(this->m_verbose){std::cout<<"(iter "<<iter<<") (alpha "<<m_alpha<<"): "<<oField.m_meanRes.back()<<" \%\n\n";}
        
      }
      else
      {
        converged = 1;
        if (oField.m_meanRes.back()>previous_res){printf("Warning : last step has diverged\n");}
        else{printf("Last step has converged\n");}
        oField.m_meanRes.back() = previous_res;
        std::cout<<this->class_name.c_str()<<" -> residue ("<<iter<<" iter) : "<<previous_res<<" \%\n";
      }
      
    }
 
 
    
};
 
#endif