BALL  1.4.2
 All Classes Namespaces Files Functions Variables Typedefs Enumerations Enumerator Properties Friends Macros Groups Pages
FFT2D.h
Go to the documentation of this file.
1 // -*- Mode: C++; tab-width: 2; -*-
2 // vi: set ts=2:
3 //
4 
5 #ifndef BALL_MATHS_TFFT2D_H
6 #define BALL_MATHS_TFFT2D_H
7 
8 #ifndef BALL_COMMON_EXCEPTION_H
9 # include <BALL/COMMON/exception.h>
10 #endif
11 
12 #ifndef BALL_DATATYPE_REGULARDATA2D_H
13 # include <BALL/DATATYPE/regularData2D.h>
14 #endif
15 
16 #ifndef BALL_MATHS_VECTOR2_H
17 # include <BALL/MATHS/vector2.h>
18 #endif
19 
20 #include <cmath>
21 #include <complex>
22 #include <fftw3.h>
23 
24 #include <BALL/MATHS/fftwCommon.h>
25 
26 
27 namespace BALL
28 {
40  template <typename ComplexTraits>
41  class TFFT2D
42  : public TRegularData2D<std::complex<typename ComplexTraits::ComplexPrecision> >
43  {
44  public:
45 
46  typedef std::complex<typename ComplexTraits::ComplexPrecision> Complex;
47  typedef TRegularData2D<std::complex<typename ComplexTraits::ComplexPrecision> > ComplexVector;
48  typedef typename TRegularData2D<std::complex<typename ComplexTraits::ComplexPrecision> >::IndexType IndexType;
49 
50  BALL_CREATE(TFFT2D)
51 
52 
55 
56 
57  TFFT2D();
58 
60  TFFT2D(const TFFT2D &data);
61 
71  // ldn is not any longer the binary logarithm but the absolute number of grid points
72  TFFT2D(Size nX, Size nY, double stepPhysX=1., double stepPhysY=1., Vector2 origin=Vector2(0.,0.), bool inFourierSpace=false);
73 
75  virtual ~TFFT2D();
76 
78 
82 
84  const TFFT2D& operator = (const TFFT2D& fft_2d);
85 
88  virtual void clear();
89 
92  virtual void destroy();
93 
95 
99 
102  bool operator == (const TFFT2D& fft_2d) const;
104 
105  // @name Accessors
106 
109  void doFFT();
110 
113  void doiFFT();
114 
120  bool translate(const Vector2& trans_origin);
121 
127  bool setPhysStepWidth(double new_width_x, double new_width_y);
128 
131  double getPhysStepWidthX() const;
132 
135  double getPhysStepWidthY() const;
136 
139  double getFourierStepWidthX() const;
140 
143  double getFourierStepWidthY() const;
144 
147  double getPhysSpaceMinX() const;
148 
151  double getPhysSpaceMinY() const;
152 
155  double getPhysSpaceMaxX() const;
156 
159  double getPhysSpaceMaxY() const;
160 
163  double getFourierSpaceMinX() const;
164 
167  double getFourierSpaceMinY() const;
168 
171  double getFourierSpaceMaxX() const;
172 
175  double getFourierSpaceMaxY() const;
176 
182  Size getMaxXIndex() const;
183 
189  Size getMaxYIndex() const;
190 
194  Size getNumberOfInverseTransforms() const;
195 
198  Vector2 getGridCoordinates(Position position) const;
199 
206  Complex getData(const Vector2& pos) const;
207 
215  Complex getInterpolatedValue(const Vector2& pos) const;
216 
224  void setData(const Vector2& pos, Complex val);
225 
231  Complex& operator[](const Vector2& pos);
232 
238  const Complex& operator[](const Vector2& pos) const;
239 
244  const Complex& operator [] (const IndexType& index) const { return TRegularData2D<Complex>::operator [](index); }
245 
250  Complex& operator [] (const IndexType& index) { return TRegularData2D<Complex>::operator [](index); }
251 
256  Complex& operator[](const Position& pos)
257  {
258  return TRegularData2D<Complex>::operator [] (pos);
259  }
260 
265  const Complex& operator[](const Position& pos) const
266  {
267  return TRegularData2D<Complex>::operator [] (pos);
268  }
269 
270  //
271  void setNumberOfFFTTransforms(Size num)
272  {
273  numPhysToFourier_ = num;
274  }
275 
276  //
277  void setNumberOfiFFTTransforms(Size num)
278  {
279  numFourierToPhys_ = num;
280  }
281 
286  Complex phase(const Vector2& pos) const;
287 
291  bool isInFourierSpace() const;
292 
293  protected:
294  Size lengthX_, lengthY_;
295  bool inFourierSpace_;
296  Size numPhysToFourier_;
297  Size numFourierToPhys_;
298  Vector2 origin_;
299  double stepPhysX_, stepPhysY_;
300  double stepFourierX_, stepFourierY_;
301  Vector2 minPhys_, maxPhys_;
302  Vector2 minFourier_, maxFourier_;
303 
304  // new version for FFTW3
305  typename ComplexTraits::FftwPlan planForward_;
306  typename ComplexTraits::FftwPlan planBackward_;
307 
308 
309  // to control plan calculation with new fftw3
310  Size dataLength_;
311  Complex *dataAdress_;
312  bool planCalculated_;
313 
314  };
315 
318  typedef TFFT2D<BALL_FFTW_DEFAULT_TRAITS> FFT2D;
319 
322  template <typename ComplexTraits>
323  const TRegularData2D<typename TFFT2D<ComplexTraits>::Complex>& operator<<
324  (TRegularData2D<typename TFFT2D<ComplexTraits>::Complex>& to, const TFFT2D<ComplexTraits>& from);
325 
330  template <typename ComplexTraits>
331  const RegularData2D& operator << (RegularData2D& to, const TFFT2D<ComplexTraits>& from);
332 
333  template <typename ComplexTraits>
334  TFFT2D<ComplexTraits>::TFFT2D()
335  : TRegularData2D<Complex>(),
336  dataLength_(0),
337  dataAdress_(0),
338  planCalculated_(false)
339  {
340  }
341 
342  template <typename ComplexTraits>
343  bool TFFT2D<ComplexTraits>::operator == (const TFFT2D<ComplexTraits>& fft2D) const
344  {
345  // test whether data_.size() == fft2D.data_.size()
346  // instead of testing 2 lengths. Better for vector handling.
347 
348  if (lengthX_ == fft2D.lengthX_ &&
349  lengthY_ == fft2D.lengthY_ &&
350  origin_ == fft2D.origin_ &&
351  stepPhysX_ == fft2D.stepPhysX_ &&
352  stepPhysY_ == fft2D.stepPhysY_ &&
353  stepFourierX_ == fft2D.stepFourierX_ &&
354  stepFourierY_ == fft2D.stepFourierY_ &&
355  minPhys_ == fft2D.minPhys_ &&
356  maxPhys_ == fft2D.maxPhys_ &&
357  minFourier_ == fft2D.minFourier_ &&
358  maxFourier_ == fft2D.maxFourier_ &&
359  numPhysToFourier_ == fft2D.numPhysToFourier_ &&
360  numFourierToPhys_ == fft2D.numFourierToPhys_)
361  {
362  Vector2 min = inFourierSpace_ ? minFourier_ : minPhys_;
363  Vector2 max = inFourierSpace_ ? maxFourier_ : maxPhys_;
364  double stepX = inFourierSpace_ ? stepFourierX_ : stepPhysX_;
365  double stepY = inFourierSpace_ ? stepFourierY_ : stepPhysY_;
366 
367  for (double posX=min.x; posX<=max.x; posX+=stepX)
368  {
369  for (double posY=min.y; posY<=max.y; posY+=stepY)
370  {
371  if (getData(Vector2(posX,posY)) != fft2D.getData(Vector2(posX,posY)))
372  {
373  return false;
374  }
375  }
376  }
377 
378  return true;
379  }
380 
381  return false;
382  }
383 
384  template <typename ComplexTraits>
385  bool TFFT2D<ComplexTraits>::translate(const Vector2& trans_origin)
386  {
387  Position internalOriginX = (Position) Maths::rint(trans_origin.x*stepPhysX_);
388  Position internalOriginY = (Position) Maths::rint(trans_origin.y*stepPhysY_);
389 
390  if ((internalOriginX <= lengthX_) && (internalOriginY <= lengthY_))
391  {
392  origin_.x = trans_origin.x;
393  origin_.y = trans_origin.y;
394 
395  minPhys_ = Vector2(-origin_.x,-origin_.y);
396  maxPhys_ = Vector2(((lengthX_-1)*stepPhysX_)-origin_.x,((lengthY_-1)*stepPhysY_)-origin_.y);
397  minFourier_ = Vector2(-(lengthX_/2.-1)*stepFourierX_,-(lengthY_/2.-1)*stepFourierY_);
398  maxFourier_ = Vector2((lengthX_/2.)*stepFourierX_,(lengthY_/2.)*stepFourierY_);
399 
400  return true;
401  }
402  else
403  {
404  return false;
405  }
406  }
407 
408  template <typename ComplexTraits>
409  bool TFFT2D<ComplexTraits>::setPhysStepWidth(double new_width_x, double new_width_y)
410  {
411  if ((new_width_x <= 0) || (new_width_y <= 0))
412  {
413  return false;
414  }
415  else
416  {
417  stepPhysX_ = new_width_x;
418  stepPhysY_ = new_width_y;
419  stepFourierX_ = 2.*M_PI/(stepPhysX_*lengthX_);
420  stepFourierY_ = 2.*M_PI/(stepPhysY_*lengthY_);
421 
422  minPhys_ = Vector2(-origin_.x,-origin_.y);
423  maxPhys_ = Vector2(((lengthX_-1)*stepPhysX_)-origin_.x,((lengthY_-1)*stepPhysY_)-origin_.y);
424  minFourier_ = Vector2(-(lengthX_/2.-1)*stepFourierX_,-(lengthY_/2.-1)*stepFourierY_);
425  maxFourier_ = Vector2((lengthX_/2.)*stepFourierX_,(lengthY_/2.)*stepFourierY_);
426 
427  return true;
428  }
429  }
430 
431  template <typename ComplexTraits>
432  double TFFT2D<ComplexTraits>::getPhysStepWidthX() const
433  {
434  return stepPhysX_;
435  }
436 
437  template <typename ComplexTraits>
438  double TFFT2D<ComplexTraits>::getPhysStepWidthY() const
439  {
440  return stepPhysY_;
441  }
442 
443  template <typename ComplexTraits>
444  double TFFT2D<ComplexTraits>::getFourierStepWidthX() const
445  {
446  return stepFourierX_;
447  }
448 
449  template <typename ComplexTraits>
450  double TFFT2D<ComplexTraits>::getFourierStepWidthY() const
451  {
452  return stepFourierY_;
453  }
454 
455  template <typename ComplexTraits>
456  double TFFT2D<ComplexTraits>::getPhysSpaceMinX() const
457  {
458  return minPhys_.x;
459  }
460 
461  template <typename ComplexTraits>
462  double TFFT2D<ComplexTraits>::getPhysSpaceMinY() const
463  {
464  return minPhys_.y;
465  }
466 
467  template <typename ComplexTraits>
468  double TFFT2D<ComplexTraits>::getPhysSpaceMaxX() const
469  {
470  return maxPhys_.x;
471  }
472 
473  template <typename ComplexTraits>
474  double TFFT2D<ComplexTraits>::getPhysSpaceMaxY() const
475  {
476  return maxPhys_.y;
477  }
478 
479  template <typename ComplexTraits>
480  double TFFT2D<ComplexTraits>::getFourierSpaceMinX() const
481  {
482  return minFourier_.x;
483  }
484 
485  template <typename ComplexTraits>
486  double TFFT2D<ComplexTraits>::getFourierSpaceMinY() const
487  {
488  return minFourier_.y;
489  }
490 
491  template <typename ComplexTraits>
492  double TFFT2D<ComplexTraits>::getFourierSpaceMaxX() const
493  {
494  return maxFourier_.x;
495  }
496 
497  template <typename ComplexTraits>
498  double TFFT2D<ComplexTraits>::getFourierSpaceMaxY() const
499  {
500  return maxFourier_.y;
501  }
502 
503  template <typename ComplexTraits>
504  Size TFFT2D<ComplexTraits>::getMaxXIndex() const
505  {
506  return (lengthX_ - 1);
507  }
508 
509  template <typename ComplexTraits>
510  Size TFFT2D<ComplexTraits>::getMaxYIndex() const
511  {
512  return (lengthY_ - 1);
513  }
514 
515  template <typename ComplexTraits>
516  Size TFFT2D<ComplexTraits>::getNumberOfInverseTransforms() const
517  {
518  return numFourierToPhys_;
519  }
520 
521  // new for compatibility with FFT3D
522  template <typename ComplexTraits>
523  Vector2 TFFT2D<ComplexTraits>::getGridCoordinates(Position position) const
524  {
525  if (!inFourierSpace_)
526  {
527  if (position >= ComplexVector::size())
528  {
529  throw Exception::OutOfGrid(__FILE__, __LINE__);
530  }
531 
532  Vector2 r;
533  Position x, y;
534 
535 
536  // AR: ??????
537  y = position % lengthY_;
538  x = position / lengthY_;
539 
540  r.set(-origin_.x + (float)x * stepPhysX_,
541  -origin_.y + (float)y * stepPhysY_);
542 
543  return r;
544  }
545  else
546  {
547  if (position >= ComplexVector::size())
548  {
549  throw Exception::OutOfGrid(__FILE__, __LINE__);
550  }
551 
552  Vector2 r;
553  Index x, y;
554 
555  // AR: ??????
556  y = position % lengthY_;
557  x = position / lengthY_;
558 
559  if (x>=lengthX_/2.)
560  {
561  x-=lengthX_;
562  }
563 
564  if (y>=lengthY_/2.)
565  {
566  y-=lengthY_;
567  }
568 
569  r.set((float)x * stepFourierX_,
570  (float)y * stepFourierY_);
571 
572  return r;
573  }
574  }
575 
576 
577 
578  template <typename ComplexTraits>
579  typename TFFT2D<ComplexTraits>::Complex TFFT2D<ComplexTraits>::getData(const Vector2& pos) const
580  {
581  Complex result;
582  double normalization=1.;
583 
584  if (!inFourierSpace_)
585  {
586  result = (*this)[pos];
587  normalization=1./((float)pow((float)(lengthX_*lengthY_),(int)numFourierToPhys_));
588  }
589  else
590  {
591  result = (*this)[pos] * phase(pos);
592  normalization=1./(2.*M_PI)*(stepPhysX_*stepPhysY_)/((float)pow((float)(lengthX_*lengthY_),(int)numFourierToPhys_));
593  }
594 
595  result *= normalization;
596 
597  return result;
598  }
599 
600  template <typename ComplexTraits>
601  typename TFFT2D<ComplexTraits>::Complex TFFT2D<ComplexTraits>::getInterpolatedValue(const Vector2& pos) const
602  {
603  Complex result;
604 
605  Vector2 min = inFourierSpace_ ? minFourier_ : minPhys_;
606  Vector2 max = inFourierSpace_ ? maxFourier_ : maxPhys_;
607  double stepX = inFourierSpace_ ? stepFourierX_ : stepPhysX_;
608  double stepY = inFourierSpace_ ? stepFourierY_ : stepPhysY_;
609 
610  if ( (pos.x < min.x) || (pos.y < min.y)
611  || (pos.x > max.x) || (pos.y > max.y) )
612  {
613  throw Exception::OutOfGrid(__FILE__, __LINE__);
614  }
615 
616  Vector2 h(pos.x - min.x, pos.y - min.y);
617  double modX = fmod((double)h.x,stepX);
618  double modY = fmod((double)h.y,stepY);
619 
620  if (modX==0 && modY ==0) // we are on the grid
621  {
622  return getData(pos);
623  }
624 
625  double beforeX = floor(h.x/stepX)*stepX+ min.x;
626  double beforeY = floor(h.y/stepY)*stepY+ min.y;
627  double afterX = ceil(h.x/stepX)*stepX+ min.x;
628  double afterY = ceil(h.y/stepY)*stepY+ min.y;
629 
630  double tx = (pos.x - beforeX)/stepX;
631  double ty = (pos.y - beforeY)/stepY;
632 
633  result = getData(Vector2(beforeX,beforeY))*(typename ComplexTraits::ComplexPrecision)((1.-tx)*(1.-ty));
634  result += getData(Vector2(afterX, beforeY))*(typename ComplexTraits::ComplexPrecision)( tx *(1.-ty));
635  result += getData(Vector2(beforeX,afterY ))*(typename ComplexTraits::ComplexPrecision)((1.-tx)* ty );
636  result += getData(Vector2(afterX, afterY ))*(typename ComplexTraits::ComplexPrecision)( tx * ty );
637 
638  return result;
639  }
640 
641  template <typename ComplexTraits>
642  void TFFT2D<ComplexTraits>::setData(const Vector2& pos, Complex val)
643  {
644  Complex dummy;
645 
646  if (!inFourierSpace_)
647  {
648  dummy = Complex(val.real()*((float)pow((float)(lengthX_*lengthY_),(int)numFourierToPhys_)),
649  val.imag()*((float)pow((float)(lengthX_*lengthY_),(int)numFourierToPhys_)));
650 
651  (*this)[pos]=dummy;
652  }
653  else
654  {
655  val*=phase(pos)*(typename ComplexTraits::ComplexPrecision)((2*M_PI/(stepPhysX_*stepPhysY_)))
656  *(typename ComplexTraits::ComplexPrecision)pow((typename ComplexTraits::ComplexPrecision)(lengthX_*lengthY_),(int)numFourierToPhys_);
657 
658  dummy = val;
659 
660  (*this)[pos]=dummy;
661  }
662  }
663 
664  template <typename ComplexTraits>
665  typename TFFT2D<ComplexTraits>::Complex& TFFT2D<ComplexTraits>::operator[](const Vector2& pos)
666  {
667  Index internalPos;
668 
669  if (!inFourierSpace_)
670  {
671  Index i, j;
672 
673  i = (Index) Maths::rint((pos.x+origin_.x)/stepPhysX_);
674  j = (Index) Maths::rint((pos.y+origin_.y)/stepPhysY_);
675 
676  internalPos = j + i*lengthY_;
677  }
678  else
679  {
680  Index i, j;
681 
682  i = (Index) Maths::rint(pos.x/stepFourierX_);
683  j = (Index) Maths::rint(pos.y/stepFourierY_);
684 
685  if (i<0)
686  {
687  i+=lengthX_;
688  }
689 
690  if (j<0)
691  {
692  j+=lengthY_;
693  }
694 
695  internalPos = (j + i*lengthY_);
696  }
697 
698  if ((internalPos < 0) || (internalPos>=(Index) (lengthX_*lengthY_)))
699  {
700  throw Exception::OutOfGrid(__FILE__, __LINE__);
701  }
702 
703  return operator [] (internalPos);
704  }
705 
706  template <typename ComplexTraits>
707  const typename TFFT2D<ComplexTraits>::Complex& TFFT2D<ComplexTraits>::operator[](const Vector2& pos) const
708  {
709  Index internalPos;
710 
711  if (!inFourierSpace_)
712  {
713  Index i, j;
714 
715  i = (Index) Maths::rint((pos.x+origin_.x)/stepPhysX_);
716  j = (Index) Maths::rint((pos.y+origin_.y)/stepPhysY_);
717 
718  internalPos = j + i*lengthY_;
719  }
720  else
721  {
722  Index i, j;
723 
724  i = (Index) Maths::rint(pos.x/stepFourierX_);
725  j = (Index) Maths::rint(pos.y/stepFourierY_);
726 
727  if (i<0)
728  {
729  i+=lengthX_;
730  }
731 
732  if (j<0)
733  {
734  j+=lengthY_;
735  }
736 
737  internalPos = (j + i*lengthY_);
738  }
739 
740  if ((internalPos < 0) || (internalPos>=(Index) (lengthX_*lengthY_)))
741  {
742  throw Exception::OutOfGrid(__FILE__, __LINE__);
743  }
744 
745  return operator [] (internalPos);
746  }
747 
748  template <typename ComplexTraits>
749  typename TFFT2D<ComplexTraits>::Complex TFFT2D<ComplexTraits>::phase(const Vector2& pos) const
750  {
751  double phase = 2.*M_PI*( Maths::rint(pos.x/stepFourierX_)*Maths::rint(origin_.x/stepPhysX_)
752  /lengthX_
753  + Maths::rint(pos.y/stepFourierY_)*Maths::rint(origin_.y/stepPhysY_)
754  /lengthY_ );
755 
756  Complex result = Complex(cos(phase), sin(phase));
757 
758  return result;
759  }
760 
761  template <typename ComplexTraits>
762  bool TFFT2D<ComplexTraits>::isInFourierSpace() const
763  {
764  return inFourierSpace_;
765  }
766 
767  template <typename ComplexTraits>
768  const TRegularData2D<typename TFFT2D<ComplexTraits>::Complex>& operator <<
769  (TRegularData2D<typename TFFT2D<ComplexTraits>::Complex>& to, const TFFT2D<ComplexTraits>& from)
770  {
771  // first decide if the FFT3D data is in Fourier space.
772  if (!from.isInFourierSpace())
773  {
774  // create a new grid
775  Size lengthX = from.getMaxXIndex()+1;
776  Size lengthY = from.getMaxYIndex()+1;
777 
778  TRegularData2D<typename TFFT2D<ComplexTraits>::Complex> newGrid(TRegularData2D<typename TFFT2D<ComplexTraits>::Complex>::IndexType(lengthX, lengthY),
779  Vector2(from.getPhysSpaceMinX(), from.getPhysSpaceMinY()),
780  Vector2(from.getPhysSpaceMaxX(), from.getPhysSpaceMaxY()));
781 
782  // and fill it
783  double normalization=1./(pow((float)(lengthX*lengthY),from.getNumberOfInverseTransforms()));
784  typename TFFT2D<ComplexTraits>::Complex dataIn;
785  typename TFFT2D<ComplexTraits>::Complex dataOut;
786 
787  for (Position i = 0; i < from.size(); i++)
788  {
789  Position x, y;
790 
791  y = i % lengthY;
792  x = i / lengthY;
793 
794  dataIn = from[i];
795  dataOut = dataIn;
796 
797  newGrid[x + y*lengthY] = dataOut*(typename ComplexTraits::ComplexPrecision)normalization;
798  }
799 
800  to = newGrid;
801 
802  return to;
803  }
804  else
805  {
806  // we are in Fourier space
807 
808  // create a new grid
809  Size lengthX = from.getMaxXIndex()+1;
810  Size lengthY = from.getMaxYIndex()+1;
811  //float stepPhysX = from.getPhysStepWidthX();
812  //float stepPhysY = from.getPhysStepWidthY();
813  float stepFourierX = from.getFourierStepWidthX();
814  float stepFourierY = from.getFourierStepWidthY();
815 
816 
817 
818  TRegularData2D<typename TFFT2D<ComplexTraits>::Complex> newGrid(TRegularData2D<typename TFFT2D<ComplexTraits>::Complex>::IndexType(lengthX, lengthY),
819  Vector2(from.getFourierSpaceMinX(),
820  from.getFourierSpaceMinY()),
821  Vector2(from.getFourierSpaceMaxX(),
822  from.getFourierSpaceMaxY()));
823 
824  // and fill it
825  // AR: old double normalization=1./(sqrt(2.*M_PI))*(stepPhysX*stepPhysY*stepPhysZ)/(pow((float)(lengthX*lengthY*lengthZ),from.getNumberOfInverseTransforms()));
826  double normalization=1./(2.*M_PI)/(pow((float)(lengthX*lengthY),from.getNumberOfInverseTransforms()));
827 
828 
829  Index x, y;
830  Vector2 r;
831  typename TFFT2D<ComplexTraits>::Complex dataIn;
832  typename TFFT2D<ComplexTraits>::Complex dataOut;
833 
834  for (Position i = 0; i < from.size(); i++)
835  {
836  y = i % lengthY;
837  x = i / lengthY;
838 
839  if (x>lengthX/2.)
840  {
841  x-=lengthX;
842  }
843 
844  if (y>lengthY/2.)
845  {
846  y-=lengthY;
847  }
848 
849  r.set((float)x * stepFourierX,
850  (float)y * stepFourierY);
851 
852  dataIn = from[i];
853  dataOut = dataIn;
854 
855  newGrid[x + y*lengthY] = dataOut*(typename ComplexTraits::ComplexPrecision)normalization*from.phase(r);
856  }
857 
858  to = newGrid;
859 
860  return to;
861  }
862  }
863 
864  template <typename ComplexTraits>
865  const RegularData2D& operator << (RegularData2D& to, const TFFT2D<ComplexTraits>& from)
866  {
867  // first decide if the FFT2D data is in Fourier space.
868  if (!from.isInFourierSpace())
869  {
870  // create a new grid
871  Size lengthX = from.getMaxXIndex()+1;
872  Size lengthY = from.getMaxYIndex()+1;
873 
874  RegularData2D newGrid(RegularData2D::IndexType(lengthX, lengthY),
875  Vector2(from.getPhysSpaceMinX(),
876  from.getPhysSpaceMinY()),
877  Vector2(from.getPhysSpaceMaxX(),
878  from.getPhysSpaceMaxY()));
879 
880  // and fill it
881  double normalization = 1./(pow((float)(lengthX*lengthY),from.getNumberOfInverseTransforms()));
882  typename TFFT2D<ComplexTraits>::Complex dataIn;
883  typename TFFT2D<ComplexTraits>::Complex dataOut;
884 
885  typename TFFT2D<ComplexTraits>::IndexType current_index;
886  typename RegularData2D::IndexType regdat_index;
887  for (current_index.x = 0; current_index.x < lengthX; current_index.x++)
888  {
889  for (current_index.y = 0; current_index.y < lengthY; current_index.y++)
890  {
891  regdat_index.x = current_index.x;
892  regdat_index.y = current_index.y;
893 
894  dataIn = from[current_index];
895  dataOut = dataIn;
896 
897  newGrid[regdat_index] = dataOut.real()*normalization;
898  }
899  }
900 
901  to = newGrid;
902 
903  return to;
904  }
905  else
906  {
907  // we are in Fourier space
908 
909  // create a new grid
910  Size lengthX = from.getMaxXIndex()+1;
911  Size lengthY = from.getMaxYIndex()+1;
912  //float stepPhysX = from.getPhysStepWidthX();
913  //float stepPhysY = from.getPhysStepWidthY();
914  float stepFourierX = from.getFourierStepWidthX();
915  float stepFourierY = from.getFourierStepWidthY();
916 
917  RegularData2D newGrid(RegularData2D::IndexType(lengthX, lengthY), Vector2(from.getFourierSpaceMinX(), from.getFourierSpaceMinY()), Vector2(from.getFourierSpaceMaxX(), from.getFourierSpaceMaxY()));
918 
919  // and fill it
920  // AR: old version double normalization=1./(sqrt(2.*M_PI))*(stepPhysX*stepPhysY*stepPhysZ)/(pow((float)(lengthX*lengthY*lengthZ),from.getNumberOfInverseTransforms()));
921  double normalization=1./(2.*M_PI)/(pow((float)(lengthX*lengthY),from.getNumberOfInverseTransforms()));
922 
923  Index x, y;
924  signed int xp, yp;
925  Vector2 r;
926  typename TFFT2D<ComplexTraits>::Complex dataIn;
927  typename TFFT2D<ComplexTraits>::Complex dataOut;
928 
929  RegularData2D::IndexType current_index;
930  for (Position i = 0; i < from.size(); i++)
931  {
932  y = i % lengthY;
933  x = i / lengthY;
934 
935  xp = x;
936  yp = y;
937 
938  if (xp>=lengthX/2.)
939  {
940  xp-=(int)lengthX;
941  }
942  if (yp>=lengthY/2.)
943  {
944  yp-=(int)lengthY;
945  }
946 
947  if (x>=lengthX/2.)
948  {
949  x-=(int)(lengthX/2.);
950  }
951  else
952  {
953  x+=(int)(lengthX/2.);
954  }
955 
956  if (y>=lengthY/2.)
957  {
958  y-=(int)(lengthY/2.);
959  }
960  else
961  {
962  y+=(int)(lengthY/2.);
963  }
964 
965 
966  r.set((float)xp * stepFourierX,
967  (float)yp * stepFourierY);
968 
969  dataIn = from[i];
970  dataOut = dataIn;
971 
972  current_index.x = x;
973  current_index.y = y;
974 
975  newGrid[current_index] = (dataOut*(typename ComplexTraits::ComplexPrecision)normalization*from.phase(r)).real();
976  }
977 
978  to = newGrid;
979 
980  return to;
981  }
982  }
983 
984 
985 
986 }
987 
988 #endif // BALL_MATHS_TFFT2D_H
Generated by   doxygen 1.8.3.1