BALL 1.5.0
Loading...
Searching...
No Matches
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
26namespace 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>
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.*Constants::PI/(stepPhysX_*lengthX_);
462 stepFourierY_ = 2.*Constants::PI/(stepPhysY_*lengthY_);
463 stepFourierZ_ = 2.*Constants::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>
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.*Constants::PI),3)*(stepPhysX_*stepPhysY_*stepPhysZ_)/((float)pow((float)(lengthX_*lengthY_*lengthZ_),(int)numFourierToPhys_));
685
686 normalization=1./pow(sqrt(2.*Constants::PI),3)/((float)pow((float)(lengthX_*lengthY_*lengthZ_),(int)numFourierToPhys_));
687 //normalization=1./(sqrt(2.*Constants::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*Constants::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*Constants::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.*Constants::PI...
895 double phase = 2.*Constants::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.*Constants::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.*Constants::PI))*(stepPhysX*stepPhysY*stepPhysZ)/(pow((float)(lengthX*lengthY*lengthZ),from.getNumberOfInverseTransforms()));
992 double normalization=1./pow(sqrt(2.*Constants::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()),
1049Vector3(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.*Constants::PI))*(stepPhysX*stepPhysY*stepPhysZ)/(pow((float)(lengthX*lengthY*lengthZ),from.getNumberOfInverseTransforms()));
1095 double normalization=1./pow(sqrt(2.*Constants::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
BALL::TFFT3D::Complex
std::complex< typename ComplexTraits::ComplexPrecision > Complex
Definition FFT3D.h:44
BALL::TFFT3D::setNumberOfiFFTTransforms
void setNumberOfiFFTTransforms(Size num)
Definition FFT3D.h:299
BALL::TFFT3D< BALL_FFTW_DEFAULT_TRAITS >::getGridCoordinates
Vector3 getGridCoordinates(Position position) const
BALL::TFFT3D< BALL_FFTW_DEFAULT_TRAITS >::planForward_
BALL_FFTW_DEFAULT_TRAITS::FftwPlan planForward_
Definition FFT3D.h:335
BALL::TFFT3D::setNumberOfFFTTransforms
void setNumberOfFFTTransforms(Size num)
Definition FFT3D.h:293
BALL::TFFT3D::operator[]
Complex & operator[](const Vector3 &pos)
Definition FFT3D.h:771
BALL::TFFT3D< BALL_FFTW_DEFAULT_TRAITS >::setData
void setData(const Vector3 &pos, Complex val)
BALL::TFFT3D::operator==
bool operator==(const TFFT3D &fft3d) const
Definition FFT3D.h:373
BALL::TFFT3D< BALL_FFTW_DEFAULT_TRAITS >::planBackward_
BALL_FFTW_DEFAULT_TRAITS::FftwPlan planBackward_
Definition FFT3D.h:336
BALL::TFFT3D::isInFourierSpace
bool isInFourierSpace() const
Definition FFT3D.h:921
BALL::TFFT3D< BALL_FFTW_DEFAULT_TRAITS >::setPhysStepWidth
bool setPhysStepWidth(double new_width_x, double new_width_y, double new_width_z)
BALL::TFFT3D::phase
Complex phase(const Vector3 &pos) const
Definition FFT3D.h:891
BALL::TFFT3D< BALL_FFTW_DEFAULT_TRAITS >::translate
bool translate(const Vector3 &trans_origin)
BALL::TFFT3D::operator[]
const Complex & operator[](const Position &pos) const
Definition FFT3D.h:287
BALL::TFFT3D::ComplexVector
TRegularData3D< std::complex< typename ComplexTraits::ComplexPrecision > > ComplexVector
Definition FFT3D.h:45
BALL::TFFT3D< BALL_FFTW_DEFAULT_TRAITS >::getData
Complex getData(const Vector3 &pos) const
BALL::TFFT3D< BALL_FFTW_DEFAULT_TRAITS >::getInterpolatedValue
Complex getInterpolatedValue(const Vector3 &pos) const
BALL::TFFT3D::TFFT3D
TFFT3D()
Default constructor.
Definition FFT3D.h:364
BALL::TRegularData3D< std::complex< ComplexTraits::ComplexPrecision > >::getData
const vector< std::complex< ComplexTraits::ComplexPrecision > > & getData() const
Definition regularData3D.h:877
BALL::TRegularData3D::operator[]
const ValueType & operator[](const IndexType &index) const
Definition regularData3D.h:226
BALL::TVector3::set
void set(const T *ptr)
Definition vector3.h:615
BALL_CREATE
#define BALL_CREATE(name)
Definition create.h:62
BALL::operator<<
BALL_EXPORT std::ostream & operator<<(std::ostream &os, const Exception::GeneralException &e)
BALL::Constants::PI
BALL_EXTERN_VARIABLE const double PI
PI.
Definition constants.h:35
BALL::Maths::rint
double rint(double x)
round to integral value in floating-point format
Definition MATHS/common.h:327
BALL::Size
BALL_SIZE_TYPE Size
Definition COMMON/global.h:114
BALL::FFT3D
TFFT3D< BALL_FFTW_DEFAULT_TRAITS > FFT3D
Definition FFT3D.h:347
BALL::RegularData3D
TRegularData3D< float > RegularData3D
Definition regularData3D.h:515
BALL::Position
BALL_SIZE_TYPE Position
Definition COMMON/global.h:140
BALL::Vector3
TVector3< float > Vector3
Definition vector3.h:1084
BALL::Complex
std::complex< BALL_COMPLEX_PRECISION > Complex
Definition complex.h:21
BALL::Index
BALL_INDEX_TYPE Index
Definition COMMON/global.h:105
Generated by doxygen 1.15.0