complexdata.h

/***************************************************************************
                          complexdata.h  -  description
                             -------------------
    begin                : Fri Apr 6 2001
    copyright            : (C) 2001 by Thies Jochimsen & Michael von Mengershausen
    email                : jochimse@cns.mpg.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, const STD_complex& val=0) : Data<STD_complex,N_rank>(dimvec) {(*this)=val;}


  // 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) {}


  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;}


  template<class T_expr>
  inline ComplexData<N_rank>& operator = (BZ_ETPARM(_bz_ArrayExpr<T_expr>) expr) {
    typedef _bz_typename T_expr::T_numtype T_numtype;
    evaluate(expr, _bz_update<T_numtype, _bz_typename T_expr::T_numtype>());
    return *this;
  }


  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:
  int length_cache;
  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;
        shift(irank,cyclshift(irank));
      }
    }
  }


}




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;
  for(int i=0; i<Array<STD_complex,N_rank>::numElements(); i++) {
    index=Data<STD_complex,N_rank>::create_index(i);
    (*this)(index)*=exp(float2imag(-2.0*PII*sum(rel_offset*index)));
  }
}


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 {
  Data<float,N_rank> phasemap_result(phase(*this));
  unwrapPhase1d(phasemap_result);
  return phasemap_result;
}



#endif

---