Ccorrelation_opticalFlow_integrated_block_A.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_A
#define CCORRELATION_OPTICALFLOW_INTEGRATED_BLOCK_A



#include "Ccorrelation_opticalFlow_integrated_block_A.h"

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

  public:
    
    Ccorrelation_opticalFlow_integrated_block_A() : Ccorrelation_opticalFlow_integrated_block<T,Timg>()
    {// constructor
      
      this->class_name = "Correlation : Optical Flow Integrated per Block version AA";
      
    }
    
    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 exec(const Cimage<Timg> &oImage, const Cmesh<T,Timg> &oMesh, Cfield<T,Timg> &oField)
    {// exec optical flow correlation

      //-----------------------------INIT--------------------------------------
      this->_3D = oMesh._3D;

      int iter(0), converged(1), thread(0);
      T previous_res(1e10);
      
      CImg<T> box;
      CImgList<T> Ke,Fe,Ve;//elementary tensors
      CImgList<T> solution(oField.m_mode);
      assign(oMesh, this->m_threadNB, Ke, Fe, Ve);
      
      CImgList<T> N(oMesh.num_var(), oMesh.num_mod(),1,1,1, 0);//Shape function at pix

      CImg<T> current_res(this->m_threadNB,1,1,1,0);
      
      CImg<T> tmp1(Fe[0].width(), Fe[0].height(), Fe[0].depth(),1,0), tmp2(Ke[0].width(), Ke[0].height(), Ke[0].depth(),1,0);
      
      T residue(0), deformed(0);
      
      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);
      //-----------------------------------------------------------------------

      do
      {//while the convergency velocity is higher than a certain level

      
//      #pragma omp parallel for if (oMesh.m_connect.size()>=oMesh.m_win[0]*oMesh.m_win[1]*oMesh.m_win[2]) firstprivate(ref, N, Ve, tmp1, tmp2, grad, U, residue, deformed, lCoor, gCoor, offset, box, thread) schedule(runtime)
        #pragma omp parallel for firstprivate(ref, N, Ve, tmp1, tmp2, grad, U, residue, deformed, lCoor, gCoor, offset, box, thread) schedule(runtime)
        cimglist_for(oMesh.m_connect, elem)
        {//for every elements
          
          //extract cropped box
          if (this->cropBox(oMesh, oImage, elem, box, offset))
          {//if the box is not entirely masked
          
            thread = omp_get_thread_num();
          
            this->metric(thread, oImage,oMesh,elem, oField, N,Ve,Ke,Fe, current_res,ref,grad,U, residue, deformed, lCoor, gCoor, box, offset);
          
            tmp1 = Fe[0].get_channel(thread);
            tmp2 = Ke[0].get_channel(thread);

            solve_CImg(tmp2,tmp1,oMesh,solution);

            cimg_forY(tmp1,y){solution[y][elem] += tmp1[y];}
            
          }

        }//implied barrier

        //sum of square residues over threads and normalization
        oField.m_meanRes.push_back(std::sqrt(current_res.sum())/oMesh.m_numPix);
        current_res.fill(0);

        if (this->m_medianFilter){this->regularize(oMesh, oImage, solution);}

        //actualized residue
        this->actu(thread, oImage, oMesh, oField, solution, converged, N, current_res, ref,U, residue, deformed, lCoor, gCoor, box, offset);

        //condition test
        this->test(solution, this->m_current_res, this->m_epsilon, oField, previous_res, iter, converged);

      }
      while(converged);
      
      //actualize field from converged solution vector
      if (oImage.m_storeActiv)
      {actu(thread, oImage, oMesh, oField, oField.m_mode, converged, N, current_res, ref,U, residue, deformed, lCoor, gCoor, box, offset);}
      
    }
 
};

#endif