complexdata.h
1 /***************************************************************************
2  complexdata.h - description
3  -------------------
4  begin : Fri Apr 6 2001
5  copyright : (C) 2000-2015 by Thies Jochimsen & Michael von Mengershausen
6  email : thies@jochimsen.de mengers@cns.mpg.de
7  ***************************************************************************/
8 
9 /***************************************************************************
10  * *
11  * This program is free software; you can redistribute it and/or modify *
12  * it under the terms of the GNU General Public License as published by *
13  * the Free Software Foundation; either version 2 of the License, or *
14  * (at your option) any later version. *
15  * *
16  ***************************************************************************/
17 
18 #ifndef COMPLEXDATA_H
19 #define COMPLEXDATA_H
20 
21 
22 #include <odindata/data.h>
23 #include <odindata/utils.h>
24 
25 
33 
35 
42 template <int N_rank>
43 class ComplexData : public Data<STD_complex,N_rank> {
44 
45  public:
46 
51 
52 
56  ComplexData(const TinyVector<int,N_rank>& dimvec) : Data<STD_complex,N_rank>(dimvec) {(*this)=STD_complex(0.0);}
57 
58 
63  // resolve ambiguities due to other constructors
64  ComplexData(int extent1) : Data<STD_complex,N_rank>(extent1) {}
65  ComplexData(int extent1, int extent2) : Data<STD_complex,N_rank>(extent1,extent2) {}
66  ComplexData(int extent1, int extent2,int extent3) : Data<STD_complex,N_rank>(extent1,extent2,extent3) {}
67  ComplexData(int extent1, int extent2,int extent3,int extent4) : Data<STD_complex,N_rank>(extent1,extent2,extent3,extent4) {}
68  ComplexData(int extent1, int extent2,int extent3,int extent4,int extent5) : Data<STD_complex,N_rank>(extent1,extent2,extent3,extent4,extent5) {}
69  ComplexData(int extent1, int extent2,int extent3,int extent4,int extent5,int extent6) : Data<STD_complex,N_rank>(extent1,extent2,extent3,extent4,extent5,extent6) {}
70  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) {}
71  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) {}
72  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) {}
73 
74 
78  ComplexData(const ComplexData<N_rank>& cd) : Data<STD_complex,N_rank>(cd) {}
79 
80 
84  ComplexData(const Data<STD_complex,N_rank>& a) : Data<STD_complex,N_rank>(a) {}
85 
86 
90  template<class T_expr>
91  ComplexData(BZ_ETPARM(_bz_ArrayExpr<T_expr>) expr) : Data<STD_complex,N_rank>(expr) {}
92 
93 
98 
99 
103  ComplexData<N_rank>& operator = (const Array<STD_complex,N_rank>& a) {Array<STD_complex,N_rank>::operator=(a); return *this;}
104 
105 
109  ComplexData<N_rank>& operator = (const STD_complex& val) {Data<STD_complex,N_rank>::operator=(val); return *this;}
110 
111 
115  template<class T_expr>
116  inline ComplexData<N_rank>& operator = (BZ_ETPARM(_bz_ArrayExpr<T_expr>) expr) {
117  ComplexData<N_rank>::reference(ComplexData<N_rank>(expr)); // forwarding to constructor
118  return *this;
119  }
120 
121  // resolve ambiguities, appears to be obsolete with blitz-0.9 and later versions
122 // 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));}
123 // 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));}
124 // 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));}
125 // 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));}
126 
127 
134  void fft(bool forward=true, bool cyclic_shift=true);
135 
143  void partial_fft(const TinyVector<bool,N_rank>& do_fft, bool forward=true, bool cyclic_shift=true);
144 
145 
152  void modulate_offset(const TinyVector<float,N_rank>& rel_offset);
153 
158  void shift_subpixel(const TinyVector<float,N_rank>& shiftvec);
159 
160 
165 
166 
167 };
168 
169 
174 
176 
177 template <int N_rank>
178 void ComplexData<N_rank>::fft(bool forward, bool cyclic_shift) {
179  Log<OdinData> odinlog("ComplexData","fft");
180 
181  TinyVector<bool,N_rank> do_fft=true;
182  ComplexData<N_rank>::partial_fft(do_fft, forward, cyclic_shift);
183 }
184 
185 
187 
188 // call FFT of GSL
189 
190 struct GslFftWorkSpace; // forward declaration
191 
192 class GslFft {
193  public:
194  GslFft(int length);
195  ~GslFft();
196  void fft1d(double* complex_data, bool forward_fft);
197 
198  private:
199  GslFftWorkSpace* ws;
200 };
201 
203 
204 template <int N_rank>
205 void ComplexData<N_rank>::partial_fft(const TinyVector<bool,N_rank>& do_fft, bool forward, bool cyclic_shift) {
206  Log<OdinData> odinlog("ComplexData","partial_fft");
207 
208  TinyVector<int,N_rank> myshape(Array<STD_complex,N_rank>::shape());
209 
210  TinyVector<int,N_rank> cyclshift=(ComplexData<N_rank>::shape())/2;
211 
212  // Shift (back) prior to FFT to get flat phasemap
213  if(cyclic_shift) {
214  for(int irank=0; irank<N_rank; irank++) {
215  if(do_fft(irank)) ComplexData<N_rank>::shift(irank,-cyclshift(irank));
216  }
217  }
218 
219  TinyVector<int,N_rank> indexvec;
220  for(int irank=0; irank<N_rank; irank++) {
221  if(do_fft(irank)) {
222 
223  // This holds the the shape of the subspace which is orthogonal to direction 'irank'
224  TinyVector<int,N_rank> ortho_shape(myshape);
225  ortho_shape(irank)=1;
226 
227  int oneline_size=myshape(irank);
228  double *complex_data= new double[2*oneline_size];
229 
230  GslFft gslfft(oneline_size);
231 
232  // The total number of elements in the orthogonal subspace
233  unsigned long n_ortho=product(ortho_shape);
234  ODINLOG(odinlog,normalDebug) << "oneline_size/irank/n_ortho=" << oneline_size << "/" << irank << "/" << n_ortho << STD_endl;
235 
236  for(unsigned long iortho=0; iortho<n_ortho; iortho++) {
237  indexvec=index2extent<N_rank>(ortho_shape,iortho);
238 
239  for(int j=0; j<oneline_size; j++) {
240  indexvec(irank)=j;
241  complex_data[2*j]= (*this)(indexvec).real();
242  complex_data[2*j+1]=(*this)(indexvec).imag();
243  }
244 
245  gslfft.fft1d(complex_data, forward);
246 
247  for(int j=0; j<oneline_size; j++) {
248  indexvec(irank)=j;
249  float scale=1.0/sqrt(double(oneline_size));
250  (*this)(indexvec)=scale*STD_complex(complex_data[2*j],complex_data[2*j+1]);
251  }
252 
253  }
254  delete[] complex_data;
255  }
256  }
257 
258  // Shift after FFT
259  if(cyclic_shift) {
260  for(int irank=0; irank<N_rank; irank++) {
261  if(do_fft(irank)) {
262  ODINLOG(odinlog,normalDebug) << "Cyclic shift in dim=" << irank << STD_endl;
263  ComplexData<N_rank>::shift(irank,cyclshift(irank)); // scope reuired by GCC 4.7
264  }
265  }
266  }
267 
268 
269 }
270 
271 
273 
274 
275 template <int N_rank>
276 void ComplexData<N_rank>::modulate_offset(const TinyVector<float,N_rank>& rel_offset) {
277  Log<OdinData> odinlog("ComplexData","modulate_offset");
278 
279  TinyVector<int, N_rank> index;
280  for(unsigned int i=0; i<Array<STD_complex,N_rank>::numElements(); i++) {
282  (*this)(index)*=exp(float2imag(-2.0*PII*sum(rel_offset*index)));
283  }
284 }
285 
287 
288 template <int N_rank>
289 void ComplexData<N_rank>::shift_subpixel(const TinyVector<float,N_rank>& shiftvec) {
290  fft(true);
291 
292  unsigned int totalsize=Array<STD_complex,N_rank>::numElements();
293  TinyVector<int, N_rank> index;
294 
295  for(unsigned int i=0; i<totalsize; i++) {
297  float im=0.0;
298  for(int idim=0; idim<N_rank; idim++) {
299  im+=-2.0*PII*(shiftvec(idim)*float(index(idim))/Array<STD_complex,N_rank>::extent(idim));
300  }
301  STD_complex phase(0,im);
302  (*this)(index)*=exp(phase);
303  }
304 
305 
306  fft(false);
307 }
308 
309 
311 
312 
313 
314 template <int N_rank>
316 
317  const ComplexData<N_rank>& indata=*this;
318 
319  TinyVector<int,N_rank> myshape(indata.shape());
320  Data<float,N_rank> result(myshape);
321 
322  TinyVector<int,N_rank> ortho_shape(myshape);
323  ortho_shape(N_rank-1)=1;
324 
325  int nlast=myshape(N_rank-1);
326  Data<float,1> phasevec(nlast);
327  Data<float,1> unwrapped(nlast);
328 
329  TinyVector<int,N_rank> indexvec;
330 
331  for(int iortho=0; iortho<product(ortho_shape); iortho++) {
332  indexvec=index2extent<N_rank>(ortho_shape,iortho);
333  int ilast;
334  for(ilast=0; ilast<nlast; ilast++) {
335  indexvec(N_rank-1)=ilast;
336  phasevec(ilast)=phase(indata(indexvec));
337  }
338  unwrapped=unwrap_phase(phasevec, nlast/2);
339  for(ilast=0; ilast<nlast; ilast++) {
340  indexvec(N_rank-1)=ilast;
341  result(indexvec)=unwrapped(ilast);
342  }
343  }
344 
345  return result;
346 }
347 
348 
349 
350 #endif
TinyVector< int, N_rank > create_index(unsigned long index) const
void shift_subpixel(const TinyVector< float, N_rank > &shiftvec)
Definition: complexdata.h:288
void shift(unsigned int shift_dim, int shift)
fvector imag(const cvector &cv)
fvector real(const cvector &cv)
ComplexData(const Data< STD_complex, N_rank > &a)
Definition: complexdata.h:84
void modulate_offset(const TinyVector< float, N_rank > &rel_offset)
Definition: complexdata.h:275
Data< float, N_rank > phasemap() const
Definition: complexdata.h:314
Definition: tjlog.h:218
ComplexData(BZ_ETPARM(_bz_ArrayExpr< T_expr >) expr)
Definition: complexdata.h:91
ComplexData(const TinyVector< int, N_rank > &dimvec)
Definition: complexdata.h:56
ComplexData(int extent1)
Definition: complexdata.h:64
ComplexData< N_rank > & operator=(const ComplexData< N_rank > &d)
Definition: complexdata.h:97
fvector phase(const cvector &cv)
Data< float, 1 > unwrap_phase(const Data< float, 1 > &phase, int startindex=0)
Data< T, N_rank > & operator=(const Data< T, N_rank > &d)
void partial_fft(const TinyVector< bool, N_rank > &do_fft, bool forward=true, bool cyclic_shift=true)
Definition: complexdata.h:204
ComplexData(const ComplexData< N_rank > &cd)
Definition: complexdata.h:78
void fft(bool forward=true, bool cyclic_shift=true)
Definition: complexdata.h:177