BALL  1.4.2
 All Classes Namespaces Files Functions Variables Typedefs Enumerations Enumerator Properties Friends Macros Groups Pages
FFT3D.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_TFFT3D_H
6 #define BALL_MATHS_TFFT3D_H
7 
8 #ifndef BALL_COMMON_EXCEPTION_H
9 # include <BALL/COMMON/exception.h>
10 #endif
11 
12 
13 #ifndef BALL_DATATYPE_REGULARDATA3D_H
14 # include <BALL/DATATYPE/regularData3D.h>
15 #endif
16 
17 //#ifndef BALL_MATHS_VECTOR2_H
18 //# include <BALL/MATHS/vector3.h>
19 //#endif
20 
21 #include <BALL/MATHS/fftwCommon.h>
22 #include <cmath>
23 #include <complex>
24 #include <fftw3.h>
25 
26 namespace BALL
27 {
38  template <typename ComplexTraits>
39  class TFFT3D
40  : public TRegularData3D<std::complex<typename ComplexTraits::ComplexPrecision> >
41  {
42  public:
43 
44  typedef std::complex<typename ComplexTraits::ComplexPrecision> Complex;
45  typedef TRegularData3D<std::complex<typename ComplexTraits::ComplexPrecision> > ComplexVector;
46 
47  BALL_CREATE(TFFT3D)
48 
49 
52 
53 
54  TFFT3D();
55 
57  TFFT3D(const TFFT3D &data);
58 
72  // AR: ldn is not any longer the binary logarithm but the absolute number of grid points
73  TFFT3D(Size ldnX, Size ldnY, Size ldnZ, double stepPhysX=1., double stepPhysY=1., double stepPhysZ=1., Vector3 origin=Vector3(0.,0.,0), bool inFourierSpace=false);
74 
76  virtual ~TFFT3D();
77 
79 
83 
85  const TFFT3D& operator = (const TFFT3D& fft_3d);
86 
89  virtual void clear();
90 
93  virtual void destroy();
94 
96 
100 
103  bool operator == (const TFFT3D& fft3d) const;
105 
106  // @name Accessors
107 
109 
111  void doFFT();
112 
115  void doiFFT();
116 
122  bool translate(const Vector3& trans_origin);
123 
129  bool setPhysStepWidth(double new_width_x, double new_width_y, double new_width_z);
130 
133  double getPhysStepWidthX() const;
134 
137  double getPhysStepWidthY() const;
138 
141  double getPhysStepWidthZ() const;
142 
145  double getFourierStepWidthX() const;
146 
149  double getFourierStepWidthY() const;
150 
153  double getFourierStepWidthZ() const;
154 
157  double getPhysSpaceMinX() const;
158 
161  double getPhysSpaceMinY() const
162 ;
163 
166  double getPhysSpaceMinZ() const;
167 
170  double getPhysSpaceMaxX() const;
171 
174  double getPhysSpaceMaxY() const;
175 
178  double getPhysSpaceMaxZ() const;
179 
182  double getFourierSpaceMinX() const;
183 
186  double getFourierSpaceMinY() const;
187 
190  double getFourierSpaceMinZ() const;
191 
194  double getFourierSpaceMaxX() const;
195 
198  double getFourierSpaceMaxY() const;
199 
202  double getFourierSpaceMaxZ() const;
203 
209  Size getMaxXIndex() const;
210 
216  Size getMaxYIndex() const;
217 
223  Size getMaxZIndex() const;
224 
228  Size getNumberOfInverseTransforms() const;
229 
232  Vector3 getGridCoordinates(Position position) const;
233 
240  Complex getData(const Vector3& pos) const;
241 
249  Complex getInterpolatedValue(const Vector3& pos) const;
250 
258  void setData(const Vector3& pos, Complex val);
259 
265  Complex& operator[](const Vector3& pos);
266 
272  const Complex& operator[](const Vector3& pos) const;
273 
278  Complex& operator[](const Position& pos)
279  {
280  return TRegularData3D<Complex>::operator [] (pos);
281  }
282 
287  const Complex& operator[](const Position& pos) const
288  {
289  return TRegularData3D<Complex>::operator [] (pos);
290  }
291 
292  // AR:
293  void setNumberOfFFTTransforms(Size num)
294  {
295  numPhysToFourier_ = num;
296  }
297 
298  // AR:
299  void setNumberOfiFFTTransforms(Size num)
300  {
301  numFourierToPhys_ = num;
302  }
303 
304 
309  Complex phase(const Vector3& pos) const;
311 
315 
319  bool isInFourierSpace() const;
321 
322  protected:
323  Size lengthX_, lengthY_, lengthZ_;
324  bool inFourierSpace_;
325  Size numPhysToFourier_;
326  Size numFourierToPhys_;
327  Vector3 origin_;
328  double stepPhysX_, stepPhysY_, stepPhysZ_;
329  double stepFourierX_, stepFourierY_, stepFourierZ_;
330  Vector3 minPhys_, maxPhys_;
331  Vector3 minFourier_, maxFourier_;
332 
333 
334  // AR: new version for FFTW3
335  typename ComplexTraits::FftwPlan planForward_;
336  typename ComplexTraits::FftwPlan planBackward_;
337 
338  // AR: to control plan calculation with new fftw3
339  Size dataLength_;
340  Complex *dataAdress_;
341  bool planCalculated_;
342 
343  };
344 
347  typedef TFFT3D<BALL_FFTW_DEFAULT_TRAITS> FFT3D;
348 
351  template <typename ComplexTraits>
352  const TRegularData3D<typename TFFT3D<ComplexTraits>::Complex>& operator <<
353  (TRegularData3D<typename TFFT3D<ComplexTraits>::Complex>& to, const TFFT3D<ComplexTraits>& from);
354 
359  template <typename ComplexTraits>
360  const RegularData3D& operator << (RegularData3D& to, const TFFT3D<ComplexTraits>& from);
361 
362 
363  template <typename ComplexTraits>
364  TFFT3D<ComplexTraits>::TFFT3D()
365  : TRegularData3D<Complex>(),
366  dataLength_(0),
367  dataAdress_(0),
368  planCalculated_(false)
369  {
370  }
371 
372  template <typename ComplexTraits>
373  bool TFFT3D<ComplexTraits>::operator == (const TFFT3D& fft3D) const
374  {
375 
376  // AR: test whether data_.size() == fft3D.data_.size()
377  // instead of testing 3 lengths. Better for vector handling.
378 
379  if (lengthX_ == fft3D.lengthX_ &&
380  lengthY_ == fft3D.lengthY_ &&
381  lengthZ_ == fft3D.lengthZ_ &&
382  origin_ == fft3D.origin_ &&
383  stepPhysX_ == fft3D.stepPhysX_ &&
384  stepPhysY_ == fft3D.stepPhysY_ &&
385  stepPhysZ_ == fft3D.stepPhysZ_ &&
386  stepFourierX_ == fft3D.stepFourierX_ &&
387  stepFourierY_ == fft3D.stepFourierY_ &&
388  stepFourierZ_ == fft3D.stepFourierZ_ &&
389  minPhys_ == fft3D.minPhys_ &&
390  maxPhys_ == fft3D.maxPhys_ &&
391  minFourier_ == fft3D.minFourier_ &&
392  maxFourier_ == fft3D.maxFourier_ &&
393  numPhysToFourier_ == fft3D.numPhysToFourier_ &&
394  numFourierToPhys_ == fft3D.numFourierToPhys_ &&
395  planCalculated_ == fft3D.planCalculated_)
396  {
397  Vector3 min = inFourierSpace_ ? minFourier_ : minPhys_;
398  Vector3 max = inFourierSpace_ ? maxFourier_ : maxPhys_;
399  double stepX = inFourierSpace_ ? stepFourierX_ : stepPhysX_;
400  double stepY = inFourierSpace_ ? stepFourierY_ : stepPhysY_;
401  double stepZ = inFourierSpace_ ? stepFourierZ_ : stepPhysZ_;
402 
403  for (double posX=min.x; posX<=max.x; posX+=stepX)
404  {
405  for (double posY=min.y; posY<=max.y; posY+=stepY)
406  {
407  for (double posZ=min.z; posZ<=max.z; posZ+=stepZ)
408  {
409  if (getData(Vector3(posX,posY,posZ)) != fft3D.getData(Vector3(posX,posY,posZ)))
410  {
411  return false;
412  }
413  }
414  }
415  }
416 
417  return true;
418  }
419 
420  return false;
421  }
422 
423  template <typename ComplexTraits>
424  bool TFFT3D<ComplexTraits>::translate(const Vector3& trans_origin)
425  {
426  Position internalOriginX = (Position) Maths::rint(trans_origin.x*stepPhysX_);
427  Position internalOriginY = (Position) Maths::rint(trans_origin.y*stepPhysY_);
428  Position internalOriginZ = (Position) Maths::rint(trans_origin.z*stepPhysZ_);
429 
430  if ((internalOriginX <= lengthX_) && (internalOriginY <= lengthY_) && (internalOriginZ <= lengthZ_))
431  {
432  origin_.x = trans_origin.x;
433  origin_.y = trans_origin.y;
434  origin_.z = trans_origin.z;
435 
436  minPhys_ = Vector3(-origin_.x,-origin_.y,-origin_.z);
437  maxPhys_ = Vector3(((lengthX_-1)*stepPhysX_)-origin_.x,((lengthY_-1)*stepPhysY_)-origin_.y,((lengthZ_-1)*stepPhysZ_)-origin_.z);
438  minFourier_ = Vector3(-(lengthX_/2.-1)*stepFourierX_,-(lengthY_/2.-1)*stepFourierY_,-(lengthZ_/2.-1)*stepFourierZ_);
439  maxFourier_ = Vector3((lengthX_/2.)*stepFourierX_,(lengthY_/2.)*stepFourierY_,(lengthZ_/2.)*stepFourierZ_);
440 
441  return true;
442  }
443  else
444  {
445  return false;
446  }
447  }
448 
449  template <typename ComplexTraits>
450  bool TFFT3D<ComplexTraits>::setPhysStepWidth(double new_width_x, double new_width_y, double new_width_z)
451  {
452  if ((new_width_x <= 0) || (new_width_y <= 0) || (new_width_z <= 0))
453  {
454  return false;
455  }
456  else
457  {
458  stepPhysX_ = new_width_x;
459  stepPhysY_ = new_width_y;
460  stepPhysZ_ = new_width_z;
461  stepFourierX_ = 2.*M_PI/(stepPhysX_*lengthX_);
462  stepFourierY_ = 2.*M_PI/(stepPhysY_*lengthY_);
463  stepFourierZ_ = 2.*M_PI/(stepPhysZ_*lengthZ_);
464 
465  minPhys_ = Vector3(-origin_.x,-origin_.y,-origin_.z);
466  maxPhys_ = Vector3(((lengthX_-1)*stepPhysX_)-origin_.x,((lengthY_-1)*stepPhysY_)-origin_.y,((lengthZ_-1)*stepPhysZ_)-origin_.z);
467  minFourier_ = Vector3(-(lengthX_/2.-1)*stepFourierX_,-(lengthY_/2.-1)*stepFourierY_,-(lengthZ_/2.-1)*stepFourierZ_);
468  maxFourier_ = Vector3((lengthX_/2.)*stepFourierX_,(lengthY_/2.)*stepFourierY_,(lengthZ_/2.)*stepFourierZ_);
469 
470  return true;
471  }
472  }
473 
474  template <typename ComplexTraits>
475  double TFFT3D<ComplexTraits>::getPhysStepWidthX() const
476  {
477  return stepPhysX_;
478  }
479 
480  template <typename ComplexTraits>
481  double TFFT3D<ComplexTraits>::getPhysStepWidthY() const
482  {
483  return stepPhysY_;
484  }
485 
486  template <typename ComplexTraits>
487  double TFFT3D<ComplexTraits>::getPhysStepWidthZ() const
488  {
489  return stepPhysZ_;
490  }
491 
492  template <typename ComplexTraits>
493  double TFFT3D<ComplexTraits>::getFourierStepWidthX() const
494  {
495  return stepFourierX_;
496  }
497 
498  template <typename ComplexTraits>
499  double TFFT3D<ComplexTraits>::getFourierStepWidthY() const
500  {
501  return stepFourierY_;
502  }
503 
504  template <typename ComplexTraits>
505  double TFFT3D<ComplexTraits>::getFourierStepWidthZ() const
506  {
507  return stepFourierZ_;
508  }
509 
510  template <typename ComplexTraits>
511  double TFFT3D<ComplexTraits>::getPhysSpaceMinX() const
512  {
513  return minPhys_.x;
514  }
515 
516  template <typename ComplexTraits>
517  double TFFT3D<ComplexTraits>::getPhysSpaceMinY() const
518  {
519  return minPhys_.y;
520  }
521 
522  template <typename ComplexTraits>
523  double TFFT3D<ComplexTraits>::getPhysSpaceMinZ() const
524  {
525  return minPhys_.z;
526  }
527 
528  template <typename ComplexTraits>
529  double TFFT3D<ComplexTraits>::getPhysSpaceMaxX() const
530  {
531  return maxPhys_.x;
532  }
533 
534  template <typename ComplexTraits>
535  double TFFT3D<ComplexTraits>::getPhysSpaceMaxY() const
536  {
537  return maxPhys_.y;
538  }
539 
540  template <typename ComplexTraits>
541  double TFFT3D<ComplexTraits>::getPhysSpaceMaxZ() const
542  {
543  return maxPhys_.z;
544  }
545 
546  template <typename ComplexTraits>
547  double TFFT3D<ComplexTraits>::getFourierSpaceMinX() const
548 
549  {
550  return minFourier_.x;
551  }
552 
553  template <typename ComplexTraits>
554  double TFFT3D<ComplexTraits>::getFourierSpaceMinY() const
555  {
556  return minFourier_.y;
557  }
558 
559  template <typename ComplexTraits>
560  double TFFT3D<ComplexTraits>::getFourierSpaceMinZ() const
561  {
562  return minFourier_.z;
563  }
564 
565  template <typename ComplexTraits>
566  double TFFT3D<ComplexTraits>::getFourierSpaceMaxX() const
567  {
568  return maxFourier_.x;
569  }
570 
571  template <typename ComplexTraits>
572  double TFFT3D<ComplexTraits>::getFourierSpaceMaxY() const
573  {
574  return maxFourier_.y;
575  }
576 
577  template <typename ComplexTraits>
578  double TFFT3D<ComplexTraits>::getFourierSpaceMaxZ() const
579  {
580  return maxFourier_.z;
581  }
582 
583  template <typename ComplexTraits>
584  Size TFFT3D<ComplexTraits>::getMaxXIndex() const
585  {
586  return (lengthX_ - 1);
587  }
588 
589  template <typename ComplexTraits>
590  Size TFFT3D<ComplexTraits>::getMaxYIndex() const
591  {
592  return (lengthY_ - 1);
593  }
594 
595  template <typename ComplexTraits>
596  Size TFFT3D<ComplexTraits>::getMaxZIndex() const
597  {
598  return (lengthZ_ - 1);
599  }
600 
601  template <typename ComplexTraits>
602  Size TFFT3D<ComplexTraits>::getNumberOfInverseTransforms() const
603  {
604  return numFourierToPhys_;
605  }
606 
607  template <typename ComplexTraits>
608  Vector3 TFFT3D<ComplexTraits>::getGridCoordinates(Position position) const
609  {
610  if (!inFourierSpace_)
611  {
612  if (position >= ComplexVector::size())
613  {
614  throw Exception::OutOfGrid(__FILE__, __LINE__);
615  }
616 
617  Vector3 r;
618  Position x, y, z;
619 
620  z = position % lengthZ_;
621  y = (position % (lengthY_ * lengthZ_)) / lengthZ_;
622  x = position / (lengthY_ * lengthZ_);
623 
624  r.set(-origin_.x + (float)x * stepPhysX_,
625  -origin_.y + (float)y * stepPhysY_,
626  -origin_.z + (float)z * stepPhysZ_);
627 
628  return r;
629  }
630  else
631  {
632  if (position >= ComplexVector::size())
633  {
634  throw Exception::OutOfGrid(__FILE__, __LINE__);
635  }
636 
637  Vector3 r;
638  Index x, y, z;
639 
640  z = position % lengthZ_;
641  y = (position % (lengthY_ * lengthZ_)) / lengthZ_;
642  x = position / (lengthY_ * lengthZ_);
643 
644  if (x>=lengthX_/2.)
645  {
646  x-=lengthX_;
647  }
648 
649  if (y>=lengthY_/2.)
650  {
651  y-=lengthY_;
652  }
653 
654  if (z>=lengthZ_/2.)
655  {
656  z-=lengthZ_;
657  }
658 
659  r.set((float)x * stepFourierX_,
660  (float)y * stepFourierY_,
661  (float)z * stepFourierZ_);
662 
663  return r;
664  }
665  }
666 
667  template <typename ComplexTraits>
668  typename TFFT3D<ComplexTraits>::Complex TFFT3D<ComplexTraits>::getData(const Vector3& pos) const
669  {
670  Complex result;
671  double normalization=1.;
672 
673  if (!inFourierSpace_)
674  {
675  result = (*this)[pos];
676  normalization=1./((float)pow((float)(lengthX_*lengthY_*lengthZ_),(int)numFourierToPhys_));
677  }
678  else
679  {
680  // AR:
681  //old: result = (*this)[pos];
682  result = (*this)[pos]*phase(pos);
683 
684  //normalization=1./pow(sqrt(2.*M_PI),3)*(stepPhysX_*stepPhysY_*stepPhysZ_)/((float)pow((float)(lengthX_*lengthY_*lengthZ_),(int)numFourierToPhys_));
685 
686  normalization=1./pow(sqrt(2.*M_PI),3)/((float)pow((float)(lengthX_*lengthY_*lengthZ_),(int)numFourierToPhys_));
687  //normalization=1./(sqrt(2.*M_PI))*(stepPhysX_*stepPhysY_*stepPhysZ_)/((float)pow((float)(lengthX_*lengthY_*lengthZ_),(int)numFourierToPhys_));
688  }
689 
690  result *= normalization;
691 
692  return result;
693  }
694 
695  template <typename ComplexTraits>
696  typename TFFT3D<ComplexTraits>::Complex TFFT3D<ComplexTraits>::getInterpolatedValue(const Vector3& pos) const
697  {
698  Complex result;
699 
700  Vector3 min = inFourierSpace_ ? minFourier_ : minPhys_;
701  Vector3 max = inFourierSpace_ ? maxFourier_ : maxPhys_;
702  double stepX = inFourierSpace_ ? stepFourierX_ : stepPhysX_;
703  double stepY = inFourierSpace_ ? stepFourierY_ : stepPhysY_;
704  double stepZ = inFourierSpace_ ? stepFourierZ_ : stepPhysZ_;
705 
706  if ( (pos.x < min.x) || (pos.y < min.y) || (pos.z < min.z)
707  || (pos.x > max.x) || (pos.y > max.y) || (pos.z > max.z) )
708  {
709  throw Exception::OutOfGrid(__FILE__, __LINE__);
710  }
711 
712  Vector3 h(pos.x - min.x, pos.y - min.y, pos.z - min.z);
713  double modX = fmod((double)h.x,stepX);
714  double modY = fmod((double)h.y,stepY);
715  double modZ = fmod((double)h.z,stepZ);
716 
717  if (modX==0 && modY==0 && modZ==0) // we are on the grid
718  {
719  return getData(pos);
720  }
721 
722  double beforeX = floor(h.x/stepX)*stepX+min.x;
723  double beforeY = floor(h.y/stepY)*stepY+min.y;
724  double beforeZ = floor(h.z/stepZ)*stepZ+min.z;
725  double afterX = ceil(h.x/stepX)*stepX+min.x;
726  double afterY = ceil(h.y/stepY)*stepY+min.y;
727  double afterZ = ceil(h.z/stepZ)*stepZ+min.z;
728 
729  double tx = (pos.x - beforeX)/stepX;
730  double ty = (pos.y - beforeY)/stepY;
731  double tz = (pos.z - beforeZ)/stepZ;
732 
733  result = getData(Vector3(beforeX,beforeY,beforeZ))*(typename ComplexTraits::ComplexPrecision)((1.-tx)*(1.-ty)*(1.-tz));
734  result += getData(Vector3(afterX, beforeY,beforeZ))*(typename ComplexTraits::ComplexPrecision)( tx *(1.-ty)*(1.-tz));
735  result += getData(Vector3(beforeX,afterY, beforeZ))*(typename ComplexTraits::ComplexPrecision)((1.-tx)* ty *(1.-tz));
736  result += getData(Vector3(beforeX,beforeY,afterZ ))*(typename ComplexTraits::ComplexPrecision)((1.-tx)*(1.-ty)* tz );
737  result += getData(Vector3(afterX, afterY, beforeZ))*(typename ComplexTraits::ComplexPrecision)( tx * ty *(1.-tz));
738  result += getData(Vector3(afterX, beforeY,afterZ ))*(typename ComplexTraits::ComplexPrecision)( tx *(1.-ty)* tz );
739  result += getData(Vector3(beforeX,afterY, afterZ ))*(typename ComplexTraits::ComplexPrecision)((1.-tx)* ty * tz );
740  result += getData(Vector3(afterX, afterY, afterZ ))*(typename ComplexTraits::ComplexPrecision)( tx * ty * tz );
741 
742  return result;
743  }
744 
745  template <typename ComplexTraits>
746  void TFFT3D<ComplexTraits>::setData(const Vector3& pos, Complex val)
747  {
748  Complex dummy;
749 
750  if (!inFourierSpace_)
751  {
752  dummy = Complex(val.real()*((typename ComplexTraits::ComplexPrecision)pow((typename ComplexTraits::ComplexPrecision)(lengthX_*lengthY_*lengthZ_),(int)numFourierToPhys_)),
753  val.imag()*((typename ComplexTraits::ComplexPrecision)pow((typename ComplexTraits::ComplexPrecision)(lengthX_*lengthY_*lengthZ_),(int)numFourierToPhys_)));
754 
755  (*this)[pos]=dummy;
756  }
757  else
758  {
759  val*=phase(pos)*(typename ComplexTraits::ComplexPrecision)((pow(sqrt(2*M_PI),3)/(stepPhysX_*stepPhysY_*stepPhysZ_)))
760  *((typename ComplexTraits::ComplexPrecision)pow((typename ComplexTraits::ComplexPrecision)(lengthX_*lengthY_*lengthZ_),(int)numFourierToPhys_));
761  /*val*=phase(pos)*(typename ComplexTraits::ComplexPrecision)((sqrt(2*M_PI)/(stepPhysX_*stepPhysY_*stepPhysZ_)))
762  *((typename ComplexTraits::ComplexPrecision)pow((typename ComplexTraits::ComplexPrecision)(lengthX_*lengthY_*lengthZ_),(int)numFourierToPhys_));*/
763 
764  dummy = val;
765 
766  (*this)[pos]=dummy;
767  }
768  }
769 
770  template <typename ComplexTraits>
771  typename TFFT3D<ComplexTraits>::Complex& TFFT3D<ComplexTraits>::operator[](const Vector3& pos)
772  {
773  Index internalPos;
774 
775  if (!inFourierSpace_)
776  {
777  Index i, j, k;
778 
779  i = (Index) Maths::rint((pos.x+origin_.x)/stepPhysX_);
780  j = (Index) Maths::rint((pos.y+origin_.y)/stepPhysY_);
781  k = (Index) Maths::rint((pos.z+origin_.z)/stepPhysZ_);
782 
783  internalPos = (k + (j + i*lengthY_)*lengthZ_);
784 
785  /*(Index) rint( (pos.z+origin_.z)/stepPhysZ_
786  + ( (pos.y+origin_.y)/stepPhysY_
787  + (pos.x+origin_.x)/stepPhysX_*lengthY_
788  ) * lengthZ_
789  ); */
790  }
791  else
792  {
793  Index i, j, k;
794 
795  i = (Index) Maths::rint(pos.x/stepFourierX_);
796  j = (Index) Maths::rint(pos.y/stepFourierY_);
797  k = (Index) Maths::rint(pos.z/stepFourierZ_);
798 
799  if (i<0)
800  {
801  i+=lengthX_;
802  }
803 
804  if (j<0)
805  {
806  j+=lengthY_;
807  }
808 
809  if (k<0)
810  {
811  k+=lengthZ_;
812  }
813 
814  internalPos = (k + (j + i*lengthY_)*lengthZ_);
815  }
816 
817  if ((internalPos < 0) || (internalPos>=(Index) (lengthX_*lengthY_*lengthZ_)))
818  {
819  throw Exception::OutOfGrid(__FILE__, __LINE__);
820  }
821 
822  return TRegularData3D<Complex>::operator[]((Position)internalPos);
823  }
824 
825  template <typename ComplexTraits>
826  const typename TFFT3D<ComplexTraits>::Complex& TFFT3D<ComplexTraits>::operator[](const Vector3& pos) const
827  {
828  Index internalPos;
829 
830  if (!inFourierSpace_)
831  {
832  Index i, j, k;
833 
834  i = (Index) Maths::rint((pos.x+origin_.x)/stepPhysX_);
835  j = (Index) Maths::rint((pos.y+origin_.y)/stepPhysY_);
836  k = (Index) Maths::rint((pos.z+origin_.z)/stepPhysZ_);
837 
838  internalPos = (k + (j + i*lengthY_)*lengthZ_);
839 
840  /*(Index) rint( (pos.z+origin_.z)/stepPhysZ_
841  + ( (pos.y+origin_.y)/stepPhysY_
842  + (pos.x+origin_.x)/stepPhysX_*lengthY_
843  ) * lengthZ_
844  ); */
845  }
846  else
847  {
848  Index i, j, k;
849 
850  i = (Index) Maths::rint(pos.x/stepFourierX_);
851  j = (Index) Maths::rint(pos.y/stepFourierY_);
852  k = (Index) Maths::rint(pos.z/stepFourierZ_);
853 
854  if (i<0)
855  {
856  i+=lengthX_;
857  }
858 
859  if (j<0)
860  {
861  j+=lengthY_;
862  }
863 
864  if (k<0)
865  {
866  k+=lengthZ_;
867  }
868 
869  internalPos = (k + (j + i*lengthY_)*lengthZ_);
870  }
871 
872  if ((internalPos < 0) || (internalPos>=(Index) (lengthX_*lengthY_*lengthZ_)))
873  {
874  throw Exception::OutOfGrid(__FILE__, __LINE__);
875  }
876 
877  return TRegularData3D<Complex>::operator[]((Position)internalPos);
878  }
879 
880  /*Complex& TFFT3D<ComplexTraits>::operator[](const Position& pos)
881  {
882  return operator [] (pos);
883  }
884 
885  const Complex& TFFT3D<ComplexTraits>::operator[](const Position& pos) const
886  {
887  return operator [] (pos);
888  }
889 */
890  template <typename ComplexTraits>
891  typename TFFT3D<ComplexTraits>::Complex TFFT3D<ComplexTraits>::phase(const Vector3& pos) const
892  {
893 
894  // AR: old version: -2.*M_PI...
895  double phase = 2.*M_PI*( (Maths::rint(pos.x/stepFourierX_))*(Maths::rint(origin_.x/stepPhysX_))
896  /lengthX_
897  + (Maths::rint(pos.y/stepFourierY_))*(Maths::rint(origin_.y/stepPhysY_))
898  /lengthY_
899  + (Maths::rint(pos.z/stepFourierZ_))*(Maths::rint(origin_.z/stepPhysZ_))
900  /lengthZ_ );
901 
902 
903  Complex result = Complex(cos(phase), sin(phase));
904 
905  return result;
906 
907  /*double phase = -2.*M_PI*( (rint(pos.x/stepFourierX_))*(rint(origin_.x/stepPhysX_))
908  /lengthX_
909  + (Maths::rint(pos.y/stepFourierY_))*(Maths::rint(origin_.y/stepPhysY_))
910  /lengthY_
911  + (Maths::rint(pos.z/stepFourierZ_))*(Maths::rint(origin_.z/stepPhysZ_))
912  /lengthZ_ );
913 
914 
915  Complex result = Complex(cos(phase), sin(phase));
916 
917  return result;*/
918  }
919 
920  template <typename ComplexTraits>
921  bool TFFT3D<ComplexTraits>::isInFourierSpace() const
922  {
923  return inFourierSpace_;
924  }
925 
926  template <typename ComplexTraits>
927  const TRegularData3D<typename TFFT3D<ComplexTraits>::Complex>& operator<<
928  (TRegularData3D<typename TFFT3D<ComplexTraits>::Complex>& to, const TFFT3D<ComplexTraits>& from)
929  {
930  // first decide if the TFFT3D data is in Fourier space.
931  if (!from.isInFourierSpace())
932  {
933  // create a new grid
934  Size lengthX = from.getMaxXIndex()+1;
935  Size lengthY = from.getMaxYIndex()+1;
936  Size lengthZ = from.getMaxZIndex()+1;
937 
938  TRegularData3D<typename TFFT3D<ComplexTraits>::Complex> newGrid(TRegularData3D<typename TFFT3D<ComplexTraits>::Complex>::IndexType(lengthX, lengthY, lengthZ),
939  Vector3(from.getPhysSpaceMinX(), from.getPhysSpaceMinY(), from.getPhysSpaceMinZ()),
940  Vector3(from.getPhysSpaceMaxX(), from.getPhysSpaceMaxY(), from.getPhysSpaceMaxZ()));
941 
942  // and fill it
943  double normalization=1./(pow((float)(lengthX*lengthY*lengthZ),(int)from.getNumberOfInverseTransforms()));
944  typename TFFT3D<ComplexTraits>::Complex dataIn;
945  typename TFFT3D<ComplexTraits>::Complex dataOut;
946 
947  for (Position i = 0; i < from.size(); i++)
948  {
949  Position x, y, z;
950 
951  z = i % lengthZ;
952  y = (i % (lengthY * lengthZ)) / lengthZ;
953  x = i / (lengthY * lengthZ);
954 
955  dataIn = from[i];
956  dataOut = dataIn;
957 
958  newGrid[x + (y + z*lengthY)*lengthZ] = dataOut*(typename ComplexTraits::ComplexPrecision)normalization;
959  }
960 
961  to = newGrid;
962 
963  return to;
964  }
965  else
966  {
967  // we are in Fourier space
968 
969  // create a new grid
970  Size lengthX = from.getMaxXIndex()+1;
971  Size lengthY = from.getMaxYIndex()+1;
972  Size lengthZ = from.getMaxZIndex()+1;
973  //float stepPhysX = from.getPhysStepWidthX();
974  //float stepPhysY = from.getPhysStepWidthY();
975  //float stepPhysZ = from.getPhysStepWidthZ();
976  float stepFourierX = from.getFourierStepWidthX();
977  float stepFourierY = from.getFourierStepWidthY();
978  float stepFourierZ = from.getFourierStepWidthZ();
979 
980 
981 
982  TRegularData3D<typename TFFT3D<ComplexTraits>::Complex> newGrid(TRegularData3D<typename TFFT3D<ComplexTraits>::Complex>::IndexType(lengthX, lengthY, lengthZ),
983  Vector3(from.getFourierSpaceMinX(),
984  from.getFourierSpaceMinY(),
985  from.getFourierSpaceMinZ()),
986  Vector3(from.getFourierSpaceMaxX(),
987  from.getFourierSpaceMaxY(),
988  from.getFourierSpaceMaxZ()));
989 
990  // and fill it
991  // AR: old double normalization=1./(sqrt(2.*M_PI))*(stepPhysX*stepPhysY*stepPhysZ)/(pow((float)(lengthX*lengthY*lengthZ),from.getNumberOfInverseTransforms()));
992  double normalization=1./pow(sqrt(2.*M_PI),3)/(pow((float)(lengthX*lengthY*lengthZ),(int)from.getNumberOfInverseTransforms()));
993 
994 
995  Index x, y, z;
996  Vector3 r;
997  typename TFFT3D<ComplexTraits>::Complex dataIn;
998  typename TFFT3D<ComplexTraits>::Complex dataOut;
999 
1000  for (Position i = 0; i < from.size(); i++)
1001  {
1002  z = i % lengthZ;
1003  y = (i % (lengthY * lengthZ)) / lengthZ;
1004  x = i / (lengthY * lengthZ);
1005 
1006  if (x>lengthX/2.)
1007  {
1008  x-=lengthX;
1009  }
1010 
1011  if (y>lengthY/2.)
1012  {
1013  y-=lengthY;
1014  }
1015 
1016  if (z>lengthZ/2.)
1017  {
1018  z-=lengthZ;
1019  }
1020 
1021  r.set((float)x * stepFourierX,
1022  (float)y * stepFourierY,
1023  (float)z * stepFourierZ);
1024 
1025  dataIn = from[i];
1026  dataOut = dataIn;
1027 
1028  newGrid[x + (y + z*lengthY)*lengthZ] = dataOut*(typename ComplexTraits::ComplexPrecision)normalization*from.phase(r);
1029  }
1030 
1031  to = newGrid;
1032 
1033  return to;
1034  }
1035  }
1036 
1037  template <typename ComplexTraits>
1038  const RegularData3D& operator << (RegularData3D& to, const TFFT3D<ComplexTraits>& from)
1039  {
1040  // first decide if the TFFT3D data is in Fourier space.
1041  if (!from.isInFourierSpace())
1042  {
1043  // create a new grid
1044  Size lengthX = from.getMaxXIndex()+1;
1045  Size lengthY = from.getMaxYIndex()+1;
1046  Size lengthZ = from.getMaxZIndex()+1;
1047 
1048  RegularData3D newGrid(RegularData3D::IndexType(lengthX, lengthY, lengthZ), Vector3(from.getPhysSpaceMinX(), from.getPhysSpaceMinY(), from.getPhysSpaceMinZ()),
1049 Vector3(from.getPhysSpaceMaxX(), from.getPhysSpaceMaxY(), from.getPhysSpaceMaxZ()));
1050 
1051  // and fill it
1052  double normalization = 1./(pow((float)(lengthX*lengthY*lengthZ),(int)from.getNumberOfInverseTransforms()));
1053  typename TFFT3D<ComplexTraits>::Complex dataIn;
1054  typename TFFT3D<ComplexTraits>::Complex dataOut;
1055 
1056  for (Position i = 0; i < from.size(); i++)
1057  {
1058  Position x, y, z;
1059 
1060  z = i % lengthZ;
1061  y = (i % (lengthY * lengthZ)) / lengthZ;
1062  x = i / (lengthY * lengthZ);
1063 
1064  dataIn = from[i];
1065  dataOut = dataIn;
1066 
1067  newGrid[x + (y + z*lengthY)*lengthZ] = dataOut.real()*normalization;
1068  }
1069 
1070  to = newGrid;
1071 
1072  return to;
1073  }
1074  else
1075  {
1076  // we are in Fourier space
1077 
1078  // create a new grid
1079  Size lengthX = from.getMaxXIndex()+1;
1080  Size lengthY = from.getMaxYIndex()+1;
1081  Size lengthZ = from.getMaxZIndex()+1;
1082  //float stepPhysX = from.getPhysStepWidthX();
1083  //float stepPhysY = from.getPhysStepWidthY();
1084  //float stepPhysZ = from.getPhysStepWidthZ();
1085  float stepFourierX = from.getFourierStepWidthX();
1086  float stepFourierY = from.getFourierStepWidthY();
1087  float stepFourierZ = from.getFourierStepWidthZ();
1088 
1089 
1090 
1091  RegularData3D newGrid(RegularData3D::IndexType(lengthX, lengthY, lengthZ), Vector3(from.getFourierSpaceMinX(), from.getFourierSpaceMinY(), from.getFourierSpaceMinZ()), Vector3(from.getFourierSpaceMaxX(), from.getFourierSpaceMaxY(), from.getFourierSpaceMaxZ()));
1092 
1093  // and fill it
1094  // AR: old version double normalization=1./(sqrt(2.*M_PI))*(stepPhysX*stepPhysY*stepPhysZ)/(pow((float)(lengthX*lengthY*lengthZ),from.getNumberOfInverseTransforms()));
1095  double normalization=1./pow(sqrt(2.*M_PI),3)/(pow((float)(lengthX*lengthY*lengthZ),(int)from.getNumberOfInverseTransforms()));
1096 
1097  Index x, y, z;
1098  signed int xp, yp, zp;
1099  Vector3 r;
1100  typename TFFT3D<ComplexTraits>::Complex dataIn;
1101  typename TFFT3D<ComplexTraits>::Complex dataOut;
1102 
1103  for (Position i = 0; i < from.size(); i++)
1104  {
1105  z = i % lengthZ;
1106  y = (i % (lengthY * lengthZ)) / lengthZ;
1107  x = i / (lengthY * lengthZ);
1108 
1109  xp = x;
1110  yp = y;
1111  zp = z;
1112 
1113  if (xp>=lengthX/2.)
1114  {
1115  xp-=(int)lengthX;
1116  }
1117  if (yp>=lengthY/2.)
1118  {
1119  yp-=(int)lengthY;
1120  }
1121  if (zp>=lengthZ/2.)
1122  {
1123  zp-=(int)lengthZ;
1124  }
1125 
1126  if (x>=lengthX/2.)
1127  {
1128  x-=(int)(lengthX/2.);
1129  }
1130  else
1131  {
1132  x+=(int)(lengthX/2.);
1133  }
1134 
1135  if (y>=lengthY/2.)
1136  {
1137  y-=(int)(lengthY/2.);
1138  }
1139  else
1140  {
1141  y+=(int)(lengthY/2.);
1142  }
1143 
1144  if (z>=lengthZ/2.)
1145  {
1146  z-=(int)(lengthZ/2.);
1147  }
1148  else
1149  {
1150  z+=(int)(lengthZ/2.);
1151  }
1152 
1153  r.set((float)xp * stepFourierX,
1154  (float)yp * stepFourierY,
1155  (float)zp * stepFourierZ);
1156 
1157  dataIn = from[i];
1158  dataOut = dataIn;
1159 
1160  newGrid[x + (y + z*lengthY)*lengthZ] = (dataOut*(typename ComplexTraits::ComplexPrecision)normalization*from.phase(r)).real();
1161  }
1162 
1163  to = newGrid;
1164 
1165  return to;
1166  }
1167  }
1168 }
1169 
1170 #endif // BALL_MATHS_TFFT3D_H
Generated by   doxygen 1.8.3.1