libodin-dev 1.8.8-2ubuntu1 / usr / include / odindata / complexdata.h

/***************************************************************************
                          complexdata.h  -  description
                             -------------------
    begin                : Fri Apr 6 2001
    copyright            : (C) 2000-2014 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>


/**
  * @addtogroup odindata
  * @{
  */



///////////////////////////////////////////////////////


/**
  * This template class holds multidimensional complex data, most suitable for NMR signal.
  * In addition to functionality  provided by the Data Array class,
  * it offers advanced numerical routines (FFT) and block-wise
  * processing of large data files.
  */
template <int N_rank>
class ComplexData : public Data<STD_complex,N_rank> {

 public:

/**
  * Default constructor
  */
  ComplexData() {}


/**
  * Constructs a complex array with the given dimensionality
  */
  ComplexData(const TinyVector<int,N_rank>& dimvec) : Data<STD_complex,N_rank>(dimvec) {(*this)=STD_complex(0.0);}


/**
  * Constructor with a variable list of arguments to specify the extend
  * in each dimension, the number of args must match N_rank
  */
  // 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)  {}


/**
  * Copy constructor
  */
  ComplexData(const ComplexData<N_rank>& cd) : Data<STD_complex,N_rank>(cd) {}


/**
  * Copy constructor from Data
  */
  ComplexData(const Data<STD_complex,N_rank>& a) : Data<STD_complex,N_rank>(a) {}


/**
  * Constructor from expression template
  */
  template<class T_expr>
  ComplexData(BZ_ETPARM(_bz_ArrayExpr<T_expr>) expr) : Data<STD_complex,N_rank>(expr) {}


/**
  * Assignment operator
  */
  ComplexData<N_rank>& operator = (const ComplexData<N_rank>& d) {Data<STD_complex,N_rank>::operator=(d); return *this;}


/**
  * Assignment operator from Blitz::Array
  */
  ComplexData<N_rank>& operator = (const Array<STD_complex,N_rank>& a) {Array<STD_complex,N_rank>::operator=(a); return *this;}


/**
  * Fills all elements with 'val'
  */
  ComplexData<N_rank>& operator = (const STD_complex& val) {Data<STD_complex,N_rank>::operator=(val); return *this;}


/**
  * Assignment operator for expression templates
  */
  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;
  }

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


/**
  * Performs an in-place FFT with the following properties:
  * - forward:  Whether to do a forward (true) or backward (false) FFT
  * - cyclic_shift: Whether to shift data cyclically so that zero frequency
  *   is in the middle of the array after FFT
  */
  void fft(bool forward=true, bool cyclic_shift=true);

/**
  * Performs an in-place FFT over a partial number of dimensions with the following properties:
  * - do_fft: A vector where the index should be set to 'true' if the FFT should be performed for the corresponding dimension
  * - forward:  Whether to do a forward (true) or backward (false) FFT
  * - cyclic_shift: Whether to shift data cyclically so that zero frequency
  *   is in the middle of the array after FFT
  */
  void partial_fft(const TinyVector<bool,N_rank>& do_fft, bool forward=true, bool cyclic_shift=true);


/**
  * Modulate a phase gradient onto the array so that after FFT, the data will be shifted
  * in each dimension by 'rel_offset' which is given as a fraction relative to the full size,
  * i.e. a value of 0.5 shifts the array by half its size.
  *
  */
  void modulate_offset(const TinyVector<float,N_rank>& rel_offset);

/**
  * Shift the data by the number of pixels given in 'shiftvec' for each dimension using the FFT,
  * i.e. this is possible even with fractions.
  */
  void shift_subpixel(const TinyVector<float,N_rank>& shiftvec);


/**
  * Returns the phase of the complex array, whereby the phase is unwrapped in the last dimension
  */
  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;
  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 {

  Range all=Range::all();

  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