complexdata.h

/***************************************************************************
                          complexdata.h  -  description
                             -------------------
    begin                : Fri Apr 6 2001
    copyright            : (C) 2000-2021 by Thies Jochimsen & Michael von Mengershausen
    email                : thies@jochimsen.de  mengers@cns.mpg.de
 ***************************************************************************/
 
/***************************************************************************
 *                                                                         *
 *   This program 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 2 of the License, or     *
 *   (at your option) any later version.                                   *
 *                                                                         *
 ***************************************************************************/
 
#ifndef COMPLEXDATA_H
#define COMPLEXDATA_H
 
 
#include <odindata/data.h>
#include <odindata/utils.h>
 
 
 
 
template <int N_rank>
class ComplexData : public Data<STD_complex,N_rank> {
 
 public:
 
  ComplexData() {}
 
 
  ComplexData(const TinyVector<int,N_rank>& dimvec) : Data<STD_complex,N_rank>(dimvec) {(*this)=STD_complex(0.0);}
 
 
  // resolve ambiguities due to other constructors
  ComplexData(int extent1) : Data<STD_complex,N_rank>(extent1)  {}
  ComplexData(int extent1, int extent2) : Data<STD_complex,N_rank>(extent1,extent2)  {}
  ComplexData(int extent1, int extent2,int extent3) : Data<STD_complex,N_rank>(extent1,extent2,extent3)  {}
  ComplexData(int extent1, int extent2,int extent3,int extent4) : Data<STD_complex,N_rank>(extent1,extent2,extent3,extent4)  {}
  ComplexData(int extent1, int extent2,int extent3,int extent4,int extent5) : Data<STD_complex,N_rank>(extent1,extent2,extent3,extent4,extent5)  {}
//  ComplexData(int extent1, int extent2,int extent3,int extent4,int extent5,int extent6) : Data<STD_complex,N_rank>(extent1,extent2,extent3,extent4,extent5,extent6)  {}
//  ComplexData(int extent1, int extent2,int extent3,int extent4,int extent5,int extent6,int extent7) : Data<STD_complex,N_rank>(extent1,extent2,extent3,extent4,extent5,extent6,extent7)  {}
//  ComplexData(int extent1, int extent2,int extent3,int extent4,int extent5,int extent6,int extent7,int extent8) : Data<STD_complex,N_rank>(extent1,extent2,extent3,extent4,extent5,extent6,extent7,extent8)  {}
//  ComplexData(int extent1, int extent2,int extent3,int extent4,int extent5,int extent6,int extent7,int extent8,int extent9) : Data<STD_complex,N_rank>(extent1,extent2,extent3,extent4,extent5,extent6,extent7,extent8,extent9)  {}
 
 
  ComplexData(const ComplexData<N_rank>& cd) : Data<STD_complex,N_rank>(cd) {}
 
 
  ComplexData(const Data<STD_complex,N_rank>& a) : Data<STD_complex,N_rank>(a) {}
 
 
#ifdef BLITZ_REPLACEMENT
  template<int M_rank>
  ComplexData(const Array<STD_complex,M_rank>& a) : Data<STD_complex,N_rank>(a) {}
#else
  template<class T_expr>
  ComplexData(BZ_ETPARM(_bz_ArrayExpr<T_expr>) expr) : Data<STD_complex,N_rank>(expr) {}
#endif
 
  ComplexData<N_rank>& operator = (const ComplexData<N_rank>& d) {Data<STD_complex,N_rank>::operator=(d); return *this;}
 
 
  ComplexData<N_rank>& operator = (const Array<STD_complex,N_rank>& a) {Array<STD_complex,N_rank>::operator=(a); return *this;}
 
 
  ComplexData<N_rank>& operator = (const STD_complex& val) {Data<STD_complex,N_rank>::operator=(val); return *this;}
 
 
#ifdef BLITZ_REPLACEMENT
  // Reimplemented Blitz::Array functions
  void reference(ComplexData<N_rank>& rhs) {
    Data<STD_complex,N_rank>::reference(rhs);
  }
  template<int M_rank>
  void reference(const RectView<Array<STD_complex,M_rank>&,STD_complex,M_rank>& rhs) {
    Data<STD_complex,N_rank>::detach_fmap();
    Array<STD_complex,N_rank>::reference(rhs);
  }
#else
  template<class T_expr>
  inline ComplexData<N_rank>& operator = (BZ_ETPARM(_bz_ArrayExpr<T_expr>) expr) {
    ComplexData<N_rank>::reference(ComplexData<N_rank>(expr)); // forwarding to constructor
    return *this;
  }
#endif
 
 
 
  // resolve ambiguities, appears to be obsolete with blitz-0.9 and later versions
//  Array<STD_complex,N_rank> operator + (const ComplexData<N_rank>& b) const {return Array<STD_complex,N_rank>(Array<STD_complex,N_rank>(*this)+Array<STD_complex,N_rank>(b));}
//  Array<STD_complex,N_rank> operator - (const ComplexData<N_rank>& b) const {return Array<STD_complex,N_rank>(Array<STD_complex,N_rank>(*this)-Array<STD_complex,N_rank>(b));}
//  Array<STD_complex,N_rank> operator * (const ComplexData<N_rank>& b) const {return Array<STD_complex,N_rank>(Array<STD_complex,N_rank>(*this)*Array<STD_complex,N_rank>(b));}
//  Array<STD_complex,N_rank> operator / (const ComplexData<N_rank>& b) const {return Array<STD_complex,N_rank>(Array<STD_complex,N_rank>(*this)/Array<STD_complex,N_rank>(b));}
 
 
  void fft(bool forward=true, bool cyclic_shift=true);
 
  void partial_fft(const TinyVector<bool,N_rank>& do_fft, bool forward=true, bool cyclic_shift=true);
 
 
  void modulate_offset(const TinyVector<float,N_rank>& rel_offset);
 
  void shift_subpixel(const TinyVector<float,N_rank>& shiftvec);
 
 
  Data<float,N_rank> phasemap() const;
 
 
};
 
 
 
 
template <int N_rank>
void ComplexData<N_rank>::fft(bool forward, bool cyclic_shift) {
  Log<OdinData> odinlog("ComplexData","fft");
 
  TinyVector<bool,N_rank> do_fft=true;
  ComplexData<N_rank>::partial_fft(do_fft, forward, cyclic_shift);
}
 
 
 
// call FFT of GSL
 
struct GslFftWorkSpace; // forward declaration
 
class GslFft {
 public:
  GslFft(int length);
  ~GslFft();
  void fft1d(double* complex_data, bool forward_fft);
 
 private:
  GslFftWorkSpace* ws;
};
 
 
template <int N_rank>
void ComplexData<N_rank>::partial_fft(const TinyVector<bool,N_rank>& do_fft, bool forward, bool cyclic_shift) {
  Log<OdinData> odinlog("ComplexData","partial_fft");
 
  TinyVector<int,N_rank> myshape(Array<STD_complex,N_rank>::shape());
 
  TinyVector<int,N_rank> cyclshift=(ComplexData<N_rank>::shape())/2;
 
  // Shift (back) prior to FFT to get flat phasemap
  if(cyclic_shift) {
    for(int irank=0; irank<N_rank; irank++) {
      if(do_fft(irank)) ComplexData<N_rank>::shift(irank,-cyclshift(irank));
    }
  }
 
  TinyVector<int,N_rank> indexvec;
  for(int irank=0; irank<N_rank; irank++) {
    if(do_fft(irank)) {
 
      // This holds the the shape of the subspace which is orthogonal to direction 'irank'
      TinyVector<int,N_rank> ortho_shape(myshape);
      ortho_shape(irank)=1;
 
      int oneline_size=myshape(irank);
      double *complex_data= new double[2*oneline_size];
 
      GslFft gslfft(oneline_size);
 
      // The total number of elements in the orthogonal subspace
      unsigned long n_ortho=product(ortho_shape);
      ODINLOG(odinlog,normalDebug) << "oneline_size/irank/n_ortho=" << oneline_size << "/" << irank << "/" << n_ortho << STD_endl;
 
      for(unsigned long iortho=0; iortho<n_ortho; iortho++) {
        indexvec=index2extent<N_rank>(ortho_shape,iortho);
 
        for(int j=0; j<oneline_size; j++) {
          indexvec(irank)=j;
          complex_data[2*j]=  (*this)(indexvec).real();
          complex_data[2*j+1]=(*this)(indexvec).imag();
        }
 
        gslfft.fft1d(complex_data, forward);
 
        for(int j=0; j<oneline_size; j++) {
          indexvec(irank)=j;
          float scale=1.0/sqrt(double(oneline_size));
          (*this)(indexvec)=scale*STD_complex(complex_data[2*j],complex_data[2*j+1]);
        }
 
      }
      delete[] complex_data;
    }
  }
 
  // Shift after FFT
  if(cyclic_shift) {
    for(int irank=0; irank<N_rank; irank++) {
      if(do_fft(irank)) {
        ODINLOG(odinlog,normalDebug) << "Cyclic shift in dim=" << irank << STD_endl;
        ComplexData<N_rank>::shift(irank,cyclshift(irank)); // scope reuired by GCC 4.7
      }
    }
  }
 
 
}
 
 
 
 
template <int N_rank>
void ComplexData<N_rank>::modulate_offset(const TinyVector<float,N_rank>& rel_offset) {
  Log<OdinData> odinlog("ComplexData","modulate_offset");
 
  TinyVector<int, N_rank> index;
  TinyVector<float, N_rank> findex;
  for(unsigned int i=0; i<Array<STD_complex,N_rank>::numElements(); i++) {
    index=Data<STD_complex,N_rank>::create_index(i);
    for(int irank=0; irank<N_rank; irank++) findex(irank)=index(irank);
    (*this)(index)*=expc(float2imag(-2.0*PII*sum(rel_offset*findex)));
  }
}
 
 
template <int N_rank>
void ComplexData<N_rank>::shift_subpixel(const TinyVector<float,N_rank>& shiftvec) {
  fft(true);
 
  unsigned int totalsize=Array<STD_complex,N_rank>::numElements();
  TinyVector<int, N_rank> index;
 
  for(unsigned int i=0; i<totalsize; i++) {
    index=Data<STD_complex,N_rank>::create_index(i);
    float im=0.0;
    for(int idim=0; idim<N_rank; idim++) {
      im+=-2.0*PII*(shiftvec(idim)*float(index(idim))/Array<STD_complex,N_rank>::extent(idim));
    }
    STD_complex phase(0,im);
    (*this)(index)*=exp(phase);
  }
 
 
  fft(false);
}
 
 
 
 
 
template <int N_rank>
Data<float,N_rank> ComplexData<N_rank>::phasemap() const {
 
  const ComplexData<N_rank>& indata=*this;
 
  TinyVector<int,N_rank> myshape(indata.shape());
  Data<float,N_rank> result(myshape);
 
  TinyVector<int,N_rank> ortho_shape(myshape);
  ortho_shape(N_rank-1)=1;
 
  int nlast=myshape(N_rank-1);
  Data<float,1> phasevec(nlast);
  Data<float,1> unwrapped(nlast);
 
  TinyVector<int,N_rank> indexvec;
 
  for(int iortho=0; iortho<product(ortho_shape); iortho++) {
    indexvec=index2extent<N_rank>(ortho_shape,iortho);
    int ilast;
    for(ilast=0; ilast<nlast; ilast++) {
      indexvec(N_rank-1)=ilast;
      phasevec(ilast)=phase(indata(indexvec));
    }
    unwrapped=unwrap_phase(phasevec, nlast/2);
    for(ilast=0; ilast<nlast; ilast++) {
      indexvec(N_rank-1)=ilast;
      result(indexvec)=unwrapped(ilast);
    }
  }
 
  return result;
}
 
 
 
#endif