Ccorrelation_intercor.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 CCORRELATION_INTERCOR
#define CCORRELATION_INTERCOR


#include "Ccorrelation.h"

template<typename T, typename Timg>
class Ccorrelation_intercor: public Ccorrelation<T,Timg>
{
  
  public:
  
    std::string m_peak_type;
    T m_meanVal;
    
    Ccorrelation_intercor() : Ccorrelation<T,Timg>()
    {// constructor
      
      this->class_name = "Correlation : Intercorrelation";
      
    }
 
 
    virtual void init(CParameterNetCDF &fp)
    {// init from parameter file
    
      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);}
      
      attribute_name = "peak";
      error = fp.loadAttribute(attribute_name,m_peak_type);
      if (this->m_verbose){printf("\t\t peak approximation method -> %s\n", m_peak_type.c_str());}

    }
   
    virtual void exec(const Cimage<Timg> &oImage, const Cmesh<T,Timg> &oMesh, Cfield<T,Timg> &oField)
    {// exec
    
      #ifdef cimg_use_fftw3
        pthread_mutex_t fftw_init_mutex;
      #endif
        
      this->_3D = oMesh._3D;
      
      //calcul of image mean value to fill unknown region within correlograms.
      T m_meanVal = oImage.m_curImg[0].mean();
      CImg<T> current_res(this->m_threadNB,1,1,1,0);
      CImg<int> p(this->m_threadNB,1,1,1,0);
      
      int counter(100);
      
      #pragma omp parallel
      {

        std::vector<T> offset;
        int thread = omp_get_thread_num();
        
        //initialisation of two sub_images containers and associated Fourier transform.
        // ! CImg complex = CImgList -> [0] = real part and [1] = imaginary one
        CImgList<T> correlogram(2, oMesh.m_win[0], oMesh.m_win[1],  oMesh.m_win[2], 1, 0);
        CImgList<T> F(correlogram), G(correlogram); //FFT of img1 & img2
        
        CNpeak<T> Peaks(m_peak_type, oMesh.grid_dims(3), oMesh.num_var()); //pointer vector toward m_nbPeak peak_type objects
        Peaks.m_verbose = this->m_verbose;
          
        #ifdef cimg_use_fftw3 //setup FFTw (thread safe by using mutex)
          correlogram[0].fftw_assign(F[0],&fftw_init_mutex);
        #endif

        #pragma omp for schedule(runtime)
        cimglist_for(oMesh.m_connect, elem)
        {// for each element

          p[thread]+=oMesh.m_win[0]*oMesh.m_win[1]*oMesh.m_win[2];
        
          if (cropBox(oMesh, oImage, oField, elem, correlogram, offset))
          {
            
            if(this->m_verbose and thread == 0)
            {this->count(elem,oMesh.m_connect.size(),counter);}// counter of master thread progression
            
            correl(correlogram, F, G);//intercorrelation

            if(oField.m_storeCorrelo)
            {//stores correlograms
              oField.m_correlogram.draw_image(offset[0],offset[1],offset[2],0,correlogram[0],1);
            }
            
            Peaks.exec(correlogram[0]);//n first peaks stat - correlogram[0] = real(ifft(F.G*))

            actu_field(oImage, oMesh, elem, Peaks, oField);//actualizes displacement field
            
            residual(oImage, oMesh, elem, offset, oField, current_res[thread]);//computes residues
            
          }// if elem is unmasked
          
        }// loops on elem

        #ifdef cimg_use_fftw3
          correlogram[0].fftw_clear();
        #endif //cimg_use_fftw3
      
        #pragma omp single
        {
          oField.m_meanRes.push_back(std::sqrt(current_res.sum())/p.sum()*oImage.m_curImg[0].size()/oImage.m_res*100);
          std::cout<<this->class_name.c_str()<<" -> residue : "<<oField.m_meanRes.back()<<" \%\n\n";
      
          if (this->m_medianFilter)
          {
            for(int var=0; var<oMesh.num_var(); var++)
            {// apply pixel scaling and median filter
              oField.m_mode[var+2].blur_median(this->m_medianFilter);
            }
          }
        }//barrier
      
        if (oImage.m_storeDef or oImage.m_storeRes or oImage.m_storeField)
        {
          #pragma omp for schedule(runtime)
          cimglist_for(oMesh.m_connect, elem)
          {// for each element
        
            if (cropBox(oMesh, oImage, oField, elem, correlogram, offset))
            {// if elem is unmasked
        
              if (oImage.m_storeField){field(oImage, oMesh, elem, offset, oField, current_res[thread]);}
              if (oImage.m_storeDef){deformField(oImage, oMesh, elem, offset, oField, current_res[thread]);}
              if (oImage.m_storeRes){residuField(oImage, oMesh, elem, offset, oField, current_res[thread]);}
              
            }
          }
          
        }
        
      }//parallele
      
    }

    void actu_field(const Cimage<Timg> &oImage, const Cmesh<T,Timg> &oMesh, const int &elem, const CNpeak<T> &Peaks, Cfield<T,Timg> &oField)
    {// actualizes fields adding current peak positions
      
      int x(0), y(0), z(0);
      oField.m_mode[0].contains((oField.m_mode[0])[elem], x,y,z);

      cimglist_for(oField.m_mode, L) 
      {//fill output field

        for (int c=0; c<oMesh.grid_dims(3); ++c)
        {//for desired number of peak
        
          if (L>1)
          {
            oField.m_mode(L,x,y,z,c) = (floor(oField.m_mode(L,x,y,z,c)/oImage.m_pix[L-2]) + Peaks.m_Nstat(L,0,0,0,c))*oImage.m_pix[L-2];
          }
          else
          {
            oField.m_mode(L,x,y,z,c) = Peaks.m_Nstat(L,0,0,0,c);
          }
          
        }
        
      }

    }
    
    void residual(const Cimage<Timg> &oImage, const Cmesh<T,Timg> &oMesh, const int &elem, const std::vector<T> &offset, Cfield<T,Timg> &oField, T &res)
    {//Deformed image corrected using U,V,W per block

      std::vector<T> U;U.assign(3,0);
      
      for (int dim=0; dim<oMesh.num_var(); dim++)
      {
        U[dim] = oField.m_mode[dim+2][elem]/oImage.m_pix[dim];
      }
      
      for (int x=0; x<oMesh.m_win[0]; x++)
      {
        for (int y=0; y<oMesh.m_win[1]; y++)
        {
          for (int z=0; z<oMesh.m_win[2]; z++)
          {
            
            if (oImage.m_curImg[0].contains(oImage.m_curImg(0,offset[0]+x, offset[1]+y, offset[2]+z)))
            {
              res += std::pow( (T)oImage.m_curImg(0, offset[0]+x, offset[1]+y, offset[2]+z,0) - (T)oImage.m_curImg[1].linear_atXYZ(offset[0]+x+U[0],offset[1]+y+U[1], offset[2]+z+U[2],0, oImage.m_curImg(0, offset[0]+x, offset[1]+y, offset[2]+z,0)) ,2);
            }

          }
        }
      }
      
    }
    
    void deformField(const Cimage<Timg> &oImage, const Cmesh<T,Timg> &oMesh, const int &elem, const std::vector<T> &offset, Cfield<T,Timg> &oField, T &res)
    {//compute corrected deformed image

      std::vector<T> U;U.assign(3,0);
      
      for (int dim=0; dim<oMesh.num_var(); dim++)
      {
        U[dim] = oField.m_mode[dim+2][elem];
      }
      
      for (int x=0; x<oMesh.m_win[0]; x++)
      {
        for (int y=0; y<oMesh.m_win[1]; y++)
        {
          for (int z=0; z<oMesh.m_win[2]; z++)
          {
            if (oImage.m_curImg[0].contains(oImage.m_curImg(0,offset[0]+x, offset[1]+y, offset[2]+z)))
            {
              oField.m_def(offset[0]+x, offset[1]+y, offset[2]+z) = oImage.m_curImg[1].linear_atXYZ(offset[0]+x+U[0]/oImage.m_pix[0],offset[1]+y+U[1]/oImage.m_pix[1], offset[2]+z+U[2]/oImage.m_pix[2],0,oImage.m_curImg(0,offset[0]+x, offset[1]+y, offset[2]+z));
            }
          }
        }
      }
      
      
    }
  
    void residuField(const Cimage<Timg> &oImage, const Cmesh<T,Timg> &oMesh, const int &elem, const std::vector<T> &offset, Cfield<T,Timg> &oField, T &res)
    {//compute the difference between the reference and the corrected deformed image (% of the init residue)

      std::vector<T> U;U.assign(3,0);
      
      for (int dim=0; dim<oMesh.num_var(); dim++)
      {
        U[dim] = oField.m_mode[dim+2][elem];
      }
      
      for (int x=0; x<oMesh.m_win[0]; x++)
      {
        for (int y=0; y<oMesh.m_win[1]; y++)
        {
          for (int z=0; z<oMesh.m_win[2]; z++)
          {
            if (oImage.m_curImg[0].contains(oImage.m_curImg(0,offset[0]+x, offset[1]+y, offset[2]+z)))
            {
              oField.m_res(offset[0]+x, offset[1]+y, offset[2]+z) = oImage.m_curImg(0,offset[0]+x, offset[1]+y, offset[2]+z) - oImage.m_curImg[1].linear_atXYZ(offset[0]+x+U[0]/oImage.m_pix[0],offset[1]+y+U[1]/oImage.m_pix[1], offset[2]+z+U[2]/oImage.m_pix[2],0,oImage.m_curImg(0, offset[0]+x, offset[1]+y, offset[2]+z));
            }
          }
        }
      }

    }
    
    void field(const Cimage<Timg> &oImage, const Cmesh<T,Timg> &oMesh, const int &elem, const std::vector<T> &offset, Cfield<T,Timg> &oField, T &res)
    {//Deformed image corrected using U,V,W per block

      std::vector<T> U;U.assign(3,0);
      
      for (int dim=0; dim<oMesh.num_var(); dim++)
      {
        U[dim] = oField.m_mode[dim+2][elem];
      }
      
      for (int x=0; x<oMesh.m_win[0]; x++)
      {
        for (int y=0; y<oMesh.m_win[1]; y++)
        {
          for (int z=0; z<oMesh.m_win[2]; z++)
          {
            
            if (oImage.m_curImg[0].contains(oImage.m_curImg(0,offset[0]+x, offset[1]+y, offset[2]+z)))
            {
              cimglist_for(oField.m_field, var)
              {
                oField.m_field(var,offset[0]+x, offset[1]+y, offset[2]+z) = U[var];
              }
            }
          }
        }
      }

    }
  
    int cropBox(const Cmesh<T,Timg> &oMesh, const Cimage<Timg> &oImage, const Cfield<T,Timg> &oField, const int &elem, CImgList<T> &box, std::vector<T> &coord0)
    {//crop 2 sub_image within oImage.m_img[0-1]
      
      std::vector<T> coord1;
      coord0.assign(3,0);
      coord1.assign(3,0);
      cimglist_for(box,im){box[im].fill(1);}
      int msk(0);
      
      //Upper left and lower right corners coordonates and rigid displacement of centroid
      cimglist_for(oMesh.m_node, dim)
      {
        coord0[dim] = oMesh.m_node[dim][oMesh.m_connect[elem].front()] - oMesh.m_win[dim]/2;
        coord1[dim] = coord0[dim] + oMesh.m_win[dim] - 1;
      }
      
      if (!oImage.m_curMask.is_empty())
      {box[0].draw_image(0,0,0,0,oImage.m_curMask.get_crop(coord0[0],coord0[1],coord0[2],0,coord1[0],coord1[1],coord1[2],0,0));}
      
      if (box[0].mean() > 0.5)
      {//if the box contain at least one unmasked pixel
        
        std::vector<int> U;U.assign(3,0);
        cimglist_for(oMesh.m_node, dim)
        {U[dim] = floor(oField.m_mode[dim + 2][elem]/oImage.m_pix[dim]);}
        
        cimg_forXYZ(box[0], x,y,z)
        {//if a box is analysed the entire in-bound image is taken into account
        
          if (oImage.m_curImg[0].contains(oImage.m_curImg(0,coord0[0]+x,coord0[1]+y,coord0[2]+z,0)))
          {box(0, x,y,z) = oImage.m_curImg(0, coord0[0]+x,coord0[1]+y,coord0[2]+z);}
          else
          {box(0, x,y,z) = this->m_meanVal;}
          
          if (oImage.m_curImg[1].contains(oImage.m_curImg(1,coord0[0]+U[0]+x,coord0[1]+U[1]+y,coord0[2]+U[2]+z,0)))
          {box(1, x,y,z) = oImage.m_curImg(1, coord0[0]+U[0]+x,coord0[1]+U[1]+y,coord0[2]+U[2]+z);}
          else
          {box(1, x,y,z) = box(0, x,y,z);}
          
        }
        
        msk = 1;
        
        
      }
     
      return msk;
     
    }
    

    virtual void correl(CImgList<T> &correlogram, CImgList<T> &F, CImgList<T> &G) = 0;
    
    
    void FFT(CImgList<T> &correlogram, CImgList<T> &F, CImgList<T> &G)
    {//FFT in place
    
      correlogram[0].FFT(F[0],F[1]);
      correlogram[0].FFT(G[0],G[1]);
      
    }
  
    void iFFT(CImgList<T> &correlogram)
    {//ifft
    
      correlogram[0].FFT(correlogram[0],correlogram[1],true);
      
    }
  
    void shift(CImgList<T> &correlogram)
    {//shift correlogram cadran
    
      std::vector<float> offset;
      offset.assign(4,0);
      
      offset[0] = (float)correlogram[0].width()/2;
      offset[1] = (float)correlogram[0].height()/2;
      if(this->_3D){offset[2] = (float)correlogram[0].depth()/2;}
      
      cimglist_for(correlogram,L)
      {
        correlogram[L].shift(offset[0],offset[1],offset[2],offset[3],2);//boundary condition : Fourier expansion
      }
      
    } 
    
    virtual void interCor(CImgList<T> &correlogram, CImgList<T> &F, CImgList<T> &G)
    {//inter-correlation
    
      cimg_forXYZ(correlogram[0],x,y,z)
      {
        correlogram(0,x,y,z) = F(0,x,y,z)*G(0,x,y,z) + F(1,x,y,z)*G(1,x,y,z);
        correlogram(1,x,y,z) = F(1,x,y,z)*G(0,x,y,z) - F(0,x,y,z)*G(1,x,y,z);
      }
      
    }

    void mean(CImgList<T> &correlogram)
    {//substract mean offset
      
      T mean0(correlogram[0].mean()), mean1(correlogram[1].mean());
      
      correlogram[0]-=mean0;
      correlogram[1]-=mean1;
      
    }
  
};

#endif