BALL
1.4.1
|
00001 // -*- Mode: C++; tab-width: 2; -*- 00002 // vi: set ts=2: 00003 // 00004 00005 #ifndef BALL_MATHS_TFFT2D_H 00006 #define BALL_MATHS_TFFT2D_H 00007 00008 #ifndef BALL_COMMON_EXCEPTION_H 00009 # include <BALL/COMMON/exception.h> 00010 #endif 00011 00012 #ifndef BALL_DATATYPE_REGULARDATA2D_H 00013 # include <BALL/DATATYPE/regularData2D.h> 00014 #endif 00015 00016 #ifndef BALL_MATHS_VECTOR2_H 00017 # include <BALL/MATHS/vector2.h> 00018 #endif 00019 00020 #include <math.h> 00021 #include <complex> 00022 #include <fftw3.h> 00023 00024 #include <BALL/MATHS/fftwCommon.h> 00025 00026 00027 namespace BALL 00028 { 00040 template <typename ComplexTraits> 00041 class TFFT2D 00042 : public TRegularData2D<std::complex<typename ComplexTraits::ComplexPrecision> > 00043 { 00044 public: 00045 00046 typedef std::complex<typename ComplexTraits::ComplexPrecision> Complex; 00047 typedef TRegularData2D<std::complex<typename ComplexTraits::ComplexPrecision> > ComplexVector; 00048 typedef typename TRegularData2D<std::complex<typename ComplexTraits::ComplexPrecision> >::IndexType IndexType; 00049 00050 BALL_CREATE(TFFT2D) 00051 00052 00055 00056 00057 TFFT2D(); 00058 00060 TFFT2D(const TFFT2D &data); 00061 00071 // ldn is not any longer the binary logarithm but the absolute number of grid points 00072 TFFT2D(Size nX, Size nY, double stepPhysX=1., double stepPhysY=1., Vector2 origin=Vector2(0.,0.), bool inFourierSpace=false); 00073 00075 virtual ~TFFT2D(); 00076 00078 00082 00084 const TFFT2D& operator = (const TFFT2D& fft_2d); 00085 00088 virtual void clear(); 00089 00092 virtual void destroy(); 00093 00095 00099 00102 bool operator == (const TFFT2D& fft_2d) const; 00104 00105 // @name Accessors 00106 00109 void doFFT(); 00110 00113 void doiFFT(); 00114 00120 bool translate(const Vector2& trans_origin); 00121 00127 bool setPhysStepWidth(double new_width_x, double new_width_y); 00128 00131 double getPhysStepWidthX() const; 00132 00135 double getPhysStepWidthY() const; 00136 00139 double getFourierStepWidthX() const; 00140 00143 double getFourierStepWidthY() const; 00144 00147 double getPhysSpaceMinX() const; 00148 00151 double getPhysSpaceMinY() const; 00152 00155 double getPhysSpaceMaxX() const; 00156 00159 double getPhysSpaceMaxY() const; 00160 00163 double getFourierSpaceMinX() const; 00164 00167 double getFourierSpaceMinY() const; 00168 00171 double getFourierSpaceMaxX() const; 00172 00175 double getFourierSpaceMaxY() const; 00176 00182 Size getMaxXIndex() const; 00183 00189 Size getMaxYIndex() const; 00190 00194 Size getNumberOfInverseTransforms() const; 00195 00198 Vector2 getGridCoordinates(Position position) const; 00199 00206 Complex getData(const Vector2& pos) const; 00207 00215 Complex getInterpolatedValue(const Vector2& pos) const; 00216 00224 void setData(const Vector2& pos, Complex val); 00225 00231 Complex& operator[](const Vector2& pos); 00232 00238 const Complex& operator[](const Vector2& pos) const; 00239 00244 const Complex& operator [] (const IndexType& index) const { return TRegularData2D<Complex>::operator [](index); } 00245 00250 Complex& operator [] (const IndexType& index) { return TRegularData2D<Complex>::operator [](index); } 00251 00256 Complex& operator[](const Position& pos) 00257 { 00258 return TRegularData2D<Complex>::operator [] (pos); 00259 } 00260 00265 const Complex& operator[](const Position& pos) const 00266 { 00267 return TRegularData2D<Complex>::operator [] (pos); 00268 } 00269 00270 // 00271 void setNumberOfFFTTransforms(Size num) 00272 { 00273 numPhysToFourier_ = num; 00274 } 00275 00276 // 00277 void setNumberOfiFFTTransforms(Size num) 00278 { 00279 numFourierToPhys_ = num; 00280 } 00281 00286 Complex phase(const Vector2& pos) const; 00287 00291 bool isInFourierSpace() const; 00292 00293 protected: 00294 Size lengthX_, lengthY_; 00295 bool inFourierSpace_; 00296 Size numPhysToFourier_; 00297 Size numFourierToPhys_; 00298 Vector2 origin_; 00299 double stepPhysX_, stepPhysY_; 00300 double stepFourierX_, stepFourierY_; 00301 Vector2 minPhys_, maxPhys_; 00302 Vector2 minFourier_, maxFourier_; 00303 00304 // new version for FFTW3 00305 typename ComplexTraits::FftwPlan planForward_; 00306 typename ComplexTraits::FftwPlan planBackward_; 00307 00308 00309 // to control plan calculation with new fftw3 00310 Size dataLength_; 00311 Complex *dataAdress_; 00312 bool planCalculated_; 00313 00314 }; 00315 00318 typedef TFFT2D<BALL_FFTW_DEFAULT_TRAITS> FFT2D; 00319 00322 template <typename ComplexTraits> 00323 const TRegularData2D<typename TFFT2D<ComplexTraits>::Complex>& operator<< 00324 (TRegularData2D<typename TFFT2D<ComplexTraits>::Complex>& to, const TFFT2D<ComplexTraits>& from); 00325 00330 template <typename ComplexTraits> 00331 const RegularData2D& operator << (RegularData2D& to, const TFFT2D<ComplexTraits>& from); 00332 00333 template <typename ComplexTraits> 00334 TFFT2D<ComplexTraits>::TFFT2D() 00335 : TRegularData2D<Complex>(), 00336 dataLength_(0), 00337 dataAdress_(0), 00338 planCalculated_(false) 00339 { 00340 } 00341 00342 template <typename ComplexTraits> 00343 bool TFFT2D<ComplexTraits>::operator == (const TFFT2D<ComplexTraits>& fft2D) const 00344 { 00345 // test whether data_.size() == fft2D.data_.size() 00346 // instead of testing 2 lengths. Better for vector handling. 00347 00348 if (lengthX_ == fft2D.lengthX_ && 00349 lengthY_ == fft2D.lengthY_ && 00350 origin_ == fft2D.origin_ && 00351 stepPhysX_ == fft2D.stepPhysX_ && 00352 stepPhysY_ == fft2D.stepPhysY_ && 00353 stepFourierX_ == fft2D.stepFourierX_ && 00354 stepFourierY_ == fft2D.stepFourierY_ && 00355 minPhys_ == fft2D.minPhys_ && 00356 maxPhys_ == fft2D.maxPhys_ && 00357 minFourier_ == fft2D.minFourier_ && 00358 maxFourier_ == fft2D.maxFourier_ && 00359 numPhysToFourier_ == fft2D.numPhysToFourier_ && 00360 numFourierToPhys_ == fft2D.numFourierToPhys_) 00361 { 00362 Vector2 min = inFourierSpace_ ? minFourier_ : minPhys_; 00363 Vector2 max = inFourierSpace_ ? maxFourier_ : maxPhys_; 00364 double stepX = inFourierSpace_ ? stepFourierX_ : stepPhysX_; 00365 double stepY = inFourierSpace_ ? stepFourierY_ : stepPhysY_; 00366 00367 for (double posX=min.x; posX<=max.x; posX+=stepX) 00368 { 00369 for (double posY=min.y; posY<=max.y; posY+=stepY) 00370 { 00371 if (getData(Vector2(posX,posY)) != fft2D.getData(Vector2(posX,posY))) 00372 { 00373 return false; 00374 } 00375 } 00376 } 00377 00378 return true; 00379 } 00380 00381 return false; 00382 } 00383 00384 template <typename ComplexTraits> 00385 bool TFFT2D<ComplexTraits>::translate(const Vector2& trans_origin) 00386 { 00387 Position internalOriginX = (Position) Maths::rint(trans_origin.x*stepPhysX_); 00388 Position internalOriginY = (Position) Maths::rint(trans_origin.y*stepPhysY_); 00389 00390 if ((internalOriginX <= lengthX_) && (internalOriginY <= lengthY_)) 00391 { 00392 origin_.x = trans_origin.x; 00393 origin_.y = trans_origin.y; 00394 00395 minPhys_ = Vector2(-origin_.x,-origin_.y); 00396 maxPhys_ = Vector2(((lengthX_-1)*stepPhysX_)-origin_.x,((lengthY_-1)*stepPhysY_)-origin_.y); 00397 minFourier_ = Vector2(-(lengthX_/2.-1)*stepFourierX_,-(lengthY_/2.-1)*stepFourierY_); 00398 maxFourier_ = Vector2((lengthX_/2.)*stepFourierX_,(lengthY_/2.)*stepFourierY_); 00399 00400 return true; 00401 } 00402 else 00403 { 00404 return false; 00405 } 00406 } 00407 00408 template <typename ComplexTraits> 00409 bool TFFT2D<ComplexTraits>::setPhysStepWidth(double new_width_x, double new_width_y) 00410 { 00411 if ((new_width_x <= 0) || (new_width_y <= 0)) 00412 { 00413 return false; 00414 } 00415 else 00416 { 00417 stepPhysX_ = new_width_x; 00418 stepPhysY_ = new_width_y; 00419 stepFourierX_ = 2.*M_PI/(stepPhysX_*lengthX_); 00420 stepFourierY_ = 2.*M_PI/(stepPhysY_*lengthY_); 00421 00422 minPhys_ = Vector2(-origin_.x,-origin_.y); 00423 maxPhys_ = Vector2(((lengthX_-1)*stepPhysX_)-origin_.x,((lengthY_-1)*stepPhysY_)-origin_.y); 00424 minFourier_ = Vector2(-(lengthX_/2.-1)*stepFourierX_,-(lengthY_/2.-1)*stepFourierY_); 00425 maxFourier_ = Vector2((lengthX_/2.)*stepFourierX_,(lengthY_/2.)*stepFourierY_); 00426 00427 return true; 00428 } 00429 } 00430 00431 template <typename ComplexTraits> 00432 double TFFT2D<ComplexTraits>::getPhysStepWidthX() const 00433 { 00434 return stepPhysX_; 00435 } 00436 00437 template <typename ComplexTraits> 00438 double TFFT2D<ComplexTraits>::getPhysStepWidthY() const 00439 { 00440 return stepPhysY_; 00441 } 00442 00443 template <typename ComplexTraits> 00444 double TFFT2D<ComplexTraits>::getFourierStepWidthX() const 00445 { 00446 return stepFourierX_; 00447 } 00448 00449 template <typename ComplexTraits> 00450 double TFFT2D<ComplexTraits>::getFourierStepWidthY() const 00451 { 00452 return stepFourierY_; 00453 } 00454 00455 template <typename ComplexTraits> 00456 double TFFT2D<ComplexTraits>::getPhysSpaceMinX() const 00457 { 00458 return minPhys_.x; 00459 } 00460 00461 template <typename ComplexTraits> 00462 double TFFT2D<ComplexTraits>::getPhysSpaceMinY() const 00463 { 00464 return minPhys_.y; 00465 } 00466 00467 template <typename ComplexTraits> 00468 double TFFT2D<ComplexTraits>::getPhysSpaceMaxX() const 00469 { 00470 return maxPhys_.x; 00471 } 00472 00473 template <typename ComplexTraits> 00474 double TFFT2D<ComplexTraits>::getPhysSpaceMaxY() const 00475 { 00476 return maxPhys_.y; 00477 } 00478 00479 template <typename ComplexTraits> 00480 double TFFT2D<ComplexTraits>::getFourierSpaceMinX() const 00481 { 00482 return minFourier_.x; 00483 } 00484 00485 template <typename ComplexTraits> 00486 double TFFT2D<ComplexTraits>::getFourierSpaceMinY() const 00487 { 00488 return minFourier_.y; 00489 } 00490 00491 template <typename ComplexTraits> 00492 double TFFT2D<ComplexTraits>::getFourierSpaceMaxX() const 00493 { 00494 return maxFourier_.x; 00495 } 00496 00497 template <typename ComplexTraits> 00498 double TFFT2D<ComplexTraits>::getFourierSpaceMaxY() const 00499 { 00500 return maxFourier_.y; 00501 } 00502 00503 template <typename ComplexTraits> 00504 Size TFFT2D<ComplexTraits>::getMaxXIndex() const 00505 { 00506 return (lengthX_ - 1); 00507 } 00508 00509 template <typename ComplexTraits> 00510 Size TFFT2D<ComplexTraits>::getMaxYIndex() const 00511 { 00512 return (lengthY_ - 1); 00513 } 00514 00515 template <typename ComplexTraits> 00516 Size TFFT2D<ComplexTraits>::getNumberOfInverseTransforms() const 00517 { 00518 return numFourierToPhys_; 00519 } 00520 00521 // new for compatibility with FFT3D 00522 template <typename ComplexTraits> 00523 Vector2 TFFT2D<ComplexTraits>::getGridCoordinates(Position position) const 00524 { 00525 if (!inFourierSpace_) 00526 { 00527 if (position >= ComplexVector::size()) 00528 { 00529 throw Exception::OutOfGrid(__FILE__, __LINE__); 00530 } 00531 00532 Vector2 r; 00533 Position x, y; 00534 00535 00536 // AR: ?????? 00537 y = position % lengthY_; 00538 x = position / lengthY_; 00539 00540 r.set(-origin_.x + (float)x * stepPhysX_, 00541 -origin_.y + (float)y * stepPhysY_); 00542 00543 return r; 00544 } 00545 else 00546 { 00547 if (position >= ComplexVector::size()) 00548 { 00549 throw Exception::OutOfGrid(__FILE__, __LINE__); 00550 } 00551 00552 Vector2 r; 00553 Index x, y; 00554 00555 // AR: ?????? 00556 y = position % lengthY_; 00557 x = position / lengthY_; 00558 00559 if (x>=lengthX_/2.) 00560 { 00561 x-=lengthX_; 00562 } 00563 00564 if (y>=lengthY_/2.) 00565 { 00566 y-=lengthY_; 00567 } 00568 00569 r.set((float)x * stepFourierX_, 00570 (float)y * stepFourierY_); 00571 00572 return r; 00573 } 00574 } 00575 00576 00577 00578 template <typename ComplexTraits> 00579 typename TFFT2D<ComplexTraits>::Complex TFFT2D<ComplexTraits>::getData(const Vector2& pos) const 00580 { 00581 Complex result; 00582 double normalization=1.; 00583 00584 if (!inFourierSpace_) 00585 { 00586 result = (*this)[pos]; 00587 normalization=1./((float)pow((float)(lengthX_*lengthY_),(int)numFourierToPhys_)); 00588 } 00589 else 00590 { 00591 result = (*this)[pos] * phase(pos); 00592 normalization=1./(2.*M_PI)*(stepPhysX_*stepPhysY_)/((float)pow((float)(lengthX_*lengthY_),(int)numFourierToPhys_)); 00593 } 00594 00595 result *= normalization; 00596 00597 return result; 00598 } 00599 00600 template <typename ComplexTraits> 00601 typename TFFT2D<ComplexTraits>::Complex TFFT2D<ComplexTraits>::getInterpolatedValue(const Vector2& pos) const 00602 { 00603 Complex result; 00604 00605 Vector2 min = inFourierSpace_ ? minFourier_ : minPhys_; 00606 Vector2 max = inFourierSpace_ ? maxFourier_ : maxPhys_; 00607 double stepX = inFourierSpace_ ? stepFourierX_ : stepPhysX_; 00608 double stepY = inFourierSpace_ ? stepFourierY_ : stepPhysY_; 00609 00610 if ( (pos.x < min.x) || (pos.y < min.y) 00611 || (pos.x > max.x) || (pos.y > max.y) ) 00612 { 00613 throw Exception::OutOfGrid(__FILE__, __LINE__); 00614 } 00615 00616 Vector2 h(pos.x - min.x, pos.y - min.y); 00617 double modX = fmod((double)h.x,stepX); 00618 double modY = fmod((double)h.y,stepY); 00619 00620 if (modX==0 && modY ==0) // we are on the grid 00621 { 00622 return getData(pos); 00623 } 00624 00625 double beforeX = floor(h.x/stepX)*stepX+ min.x; 00626 double beforeY = floor(h.y/stepY)*stepY+ min.y; 00627 double afterX = ceil(h.x/stepX)*stepX+ min.x; 00628 double afterY = ceil(h.y/stepY)*stepY+ min.y; 00629 00630 double tx = (pos.x - beforeX)/stepX; 00631 double ty = (pos.y - beforeY)/stepY; 00632 00633 result = getData(Vector2(beforeX,beforeY))*(typename ComplexTraits::ComplexPrecision)((1.-tx)*(1.-ty)); 00634 result += getData(Vector2(afterX, beforeY))*(typename ComplexTraits::ComplexPrecision)( tx *(1.-ty)); 00635 result += getData(Vector2(beforeX,afterY ))*(typename ComplexTraits::ComplexPrecision)((1.-tx)* ty ); 00636 result += getData(Vector2(afterX, afterY ))*(typename ComplexTraits::ComplexPrecision)( tx * ty ); 00637 00638 return result; 00639 } 00640 00641 template <typename ComplexTraits> 00642 void TFFT2D<ComplexTraits>::setData(const Vector2& pos, Complex val) 00643 { 00644 Complex dummy; 00645 00646 if (!inFourierSpace_) 00647 { 00648 dummy = Complex(val.real()*((float)pow((float)(lengthX_*lengthY_),(int)numFourierToPhys_)), 00649 val.imag()*((float)pow((float)(lengthX_*lengthY_),(int)numFourierToPhys_))); 00650 00651 (*this)[pos]=dummy; 00652 } 00653 else 00654 { 00655 val*=phase(pos)*(typename ComplexTraits::ComplexPrecision)((2*M_PI/(stepPhysX_*stepPhysY_))) 00656 *(typename ComplexTraits::ComplexPrecision)pow((typename ComplexTraits::ComplexPrecision)(lengthX_*lengthY_),(int)numFourierToPhys_); 00657 00658 dummy = val; 00659 00660 (*this)[pos]=dummy; 00661 } 00662 } 00663 00664 template <typename ComplexTraits> 00665 typename TFFT2D<ComplexTraits>::Complex& TFFT2D<ComplexTraits>::operator[](const Vector2& pos) 00666 { 00667 Index internalPos; 00668 00669 if (!inFourierSpace_) 00670 { 00671 Index i, j; 00672 00673 i = (Index) Maths::rint((pos.x+origin_.x)/stepPhysX_); 00674 j = (Index) Maths::rint((pos.y+origin_.y)/stepPhysY_); 00675 00676 internalPos = j + i*lengthY_; 00677 } 00678 else 00679 { 00680 Index i, j; 00681 00682 i = (Index) Maths::rint(pos.x/stepFourierX_); 00683 j = (Index) Maths::rint(pos.y/stepFourierY_); 00684 00685 if (i<0) 00686 { 00687 i+=lengthX_; 00688 } 00689 00690 if (j<0) 00691 { 00692 j+=lengthY_; 00693 } 00694 00695 internalPos = (j + i*lengthY_); 00696 } 00697 00698 if ((internalPos < 0) || (internalPos>=(Index) (lengthX_*lengthY_))) 00699 { 00700 throw Exception::OutOfGrid(__FILE__, __LINE__); 00701 } 00702 00703 return operator [] (internalPos); 00704 } 00705 00706 template <typename ComplexTraits> 00707 const typename TFFT2D<ComplexTraits>::Complex& TFFT2D<ComplexTraits>::operator[](const Vector2& pos) const 00708 { 00709 Index internalPos; 00710 00711 if (!inFourierSpace_) 00712 { 00713 Index i, j; 00714 00715 i = (Index) Maths::rint((pos.x+origin_.x)/stepPhysX_); 00716 j = (Index) Maths::rint((pos.y+origin_.y)/stepPhysY_); 00717 00718 internalPos = j + i*lengthY_; 00719 } 00720 else 00721 { 00722 Index i, j; 00723 00724 i = (Index) Maths::rint(pos.x/stepFourierX_); 00725 j = (Index) Maths::rint(pos.y/stepFourierY_); 00726 00727 if (i<0) 00728 { 00729 i+=lengthX_; 00730 } 00731 00732 if (j<0) 00733 { 00734 j+=lengthY_; 00735 } 00736 00737 internalPos = (j + i*lengthY_); 00738 } 00739 00740 if ((internalPos < 0) || (internalPos>=(Index) (lengthX_*lengthY_))) 00741 { 00742 throw Exception::OutOfGrid(__FILE__, __LINE__); 00743 } 00744 00745 return operator [] (internalPos); 00746 } 00747 00748 template <typename ComplexTraits> 00749 typename TFFT2D<ComplexTraits>::Complex TFFT2D<ComplexTraits>::phase(const Vector2& pos) const 00750 { 00751 double phase = 2.*M_PI*( Maths::rint(pos.x/stepFourierX_)*Maths::rint(origin_.x/stepPhysX_) 00752 /lengthX_ 00753 + Maths::rint(pos.y/stepFourierY_)*Maths::rint(origin_.y/stepPhysY_) 00754 /lengthY_ ); 00755 00756 Complex result = Complex(cos(phase), sin(phase)); 00757 00758 return result; 00759 } 00760 00761 template <typename ComplexTraits> 00762 bool TFFT2D<ComplexTraits>::isInFourierSpace() const 00763 { 00764 return inFourierSpace_; 00765 } 00766 00767 template <typename ComplexTraits> 00768 const TRegularData2D<typename TFFT2D<ComplexTraits>::Complex>& operator << 00769 (TRegularData2D<typename TFFT2D<ComplexTraits>::Complex>& to, const TFFT2D<ComplexTraits>& from) 00770 { 00771 // first decide if the FFT3D data is in Fourier space. 00772 if (!from.isInFourierSpace()) 00773 { 00774 // create a new grid 00775 Size lengthX = from.getMaxXIndex()+1; 00776 Size lengthY = from.getMaxYIndex()+1; 00777 00778 TRegularData2D<typename TFFT2D<ComplexTraits>::Complex> newGrid(TRegularData2D<typename TFFT2D<ComplexTraits>::Complex>::IndexType(lengthX, lengthY), 00779 Vector2(from.getPhysSpaceMinX(), from.getPhysSpaceMinY()), 00780 Vector2(from.getPhysSpaceMaxX(), from.getPhysSpaceMaxY())); 00781 00782 // and fill it 00783 double normalization=1./(pow((float)(lengthX*lengthY),from.getNumberOfInverseTransforms())); 00784 typename TFFT2D<ComplexTraits>::Complex dataIn; 00785 typename TFFT2D<ComplexTraits>::Complex dataOut; 00786 00787 for (Position i = 0; i < from.size(); i++) 00788 { 00789 Position x, y; 00790 00791 y = i % lengthY; 00792 x = i / lengthY; 00793 00794 dataIn = from[i]; 00795 dataOut = dataIn; 00796 00797 newGrid[x + y*lengthY] = dataOut*(typename ComplexTraits::ComplexPrecision)normalization; 00798 } 00799 00800 to = newGrid; 00801 00802 return to; 00803 } 00804 else 00805 { 00806 // we are in Fourier space 00807 00808 // create a new grid 00809 Size lengthX = from.getMaxXIndex()+1; 00810 Size lengthY = from.getMaxYIndex()+1; 00811 //float stepPhysX = from.getPhysStepWidthX(); 00812 //float stepPhysY = from.getPhysStepWidthY(); 00813 float stepFourierX = from.getFourierStepWidthX(); 00814 float stepFourierY = from.getFourierStepWidthY(); 00815 00816 00817 00818 TRegularData2D<typename TFFT2D<ComplexTraits>::Complex> newGrid(TRegularData2D<typename TFFT2D<ComplexTraits>::Complex>::IndexType(lengthX, lengthY), 00819 Vector2(from.getFourierSpaceMinX(), 00820 from.getFourierSpaceMinY()), 00821 Vector2(from.getFourierSpaceMaxX(), 00822 from.getFourierSpaceMaxY())); 00823 00824 // and fill it 00825 // AR: old double normalization=1./(sqrt(2.*M_PI))*(stepPhysX*stepPhysY*stepPhysZ)/(pow((float)(lengthX*lengthY*lengthZ),from.getNumberOfInverseTransforms())); 00826 double normalization=1./(2.*M_PI)/(pow((float)(lengthX*lengthY),from.getNumberOfInverseTransforms())); 00827 00828 00829 Index x, y; 00830 Vector2 r; 00831 typename TFFT2D<ComplexTraits>::Complex dataIn; 00832 typename TFFT2D<ComplexTraits>::Complex dataOut; 00833 00834 for (Position i = 0; i < from.size(); i++) 00835 { 00836 y = i % lengthY; 00837 x = i / lengthY; 00838 00839 if (x>lengthX/2.) 00840 { 00841 x-=lengthX; 00842 } 00843 00844 if (y>lengthY/2.) 00845 { 00846 y-=lengthY; 00847 } 00848 00849 r.set((float)x * stepFourierX, 00850 (float)y * stepFourierY); 00851 00852 dataIn = from[i]; 00853 dataOut = dataIn; 00854 00855 newGrid[x + y*lengthY] = dataOut*(typename ComplexTraits::ComplexPrecision)normalization*from.phase(r); 00856 } 00857 00858 to = newGrid; 00859 00860 return to; 00861 } 00862 } 00863 00864 template <typename ComplexTraits> 00865 const RegularData2D& operator << (RegularData2D& to, const TFFT2D<ComplexTraits>& from) 00866 { 00867 // first decide if the FFT2D data is in Fourier space. 00868 if (!from.isInFourierSpace()) 00869 { 00870 // create a new grid 00871 Size lengthX = from.getMaxXIndex()+1; 00872 Size lengthY = from.getMaxYIndex()+1; 00873 00874 RegularData2D newGrid(RegularData2D::IndexType(lengthX, lengthY), 00875 Vector2(from.getPhysSpaceMinX(), 00876 from.getPhysSpaceMinY()), 00877 Vector2(from.getPhysSpaceMaxX(), 00878 from.getPhysSpaceMaxY())); 00879 00880 // and fill it 00881 double normalization = 1./(pow((float)(lengthX*lengthY),from.getNumberOfInverseTransforms())); 00882 typename TFFT2D<ComplexTraits>::Complex dataIn; 00883 typename TFFT2D<ComplexTraits>::Complex dataOut; 00884 00885 typename TFFT2D<ComplexTraits>::IndexType current_index; 00886 typename RegularData2D::IndexType regdat_index; 00887 for (current_index.x = 0; current_index.x < lengthX; current_index.x++) 00888 { 00889 for (current_index.y = 0; current_index.y < lengthY; current_index.y++) 00890 { 00891 regdat_index.x = current_index.x; 00892 regdat_index.y = current_index.y; 00893 00894 dataIn = from[current_index]; 00895 dataOut = dataIn; 00896 00897 newGrid[regdat_index] = dataOut.real()*normalization; 00898 } 00899 } 00900 00901 to = newGrid; 00902 00903 return to; 00904 } 00905 else 00906 { 00907 // we are in Fourier space 00908 00909 // create a new grid 00910 Size lengthX = from.getMaxXIndex()+1; 00911 Size lengthY = from.getMaxYIndex()+1; 00912 //float stepPhysX = from.getPhysStepWidthX(); 00913 //float stepPhysY = from.getPhysStepWidthY(); 00914 float stepFourierX = from.getFourierStepWidthX(); 00915 float stepFourierY = from.getFourierStepWidthY(); 00916 00917 RegularData2D newGrid(RegularData2D::IndexType(lengthX, lengthY), Vector2(from.getFourierSpaceMinX(), from.getFourierSpaceMinY()), Vector2(from.getFourierSpaceMaxX(), from.getFourierSpaceMaxY())); 00918 00919 // and fill it 00920 // AR: old version double normalization=1./(sqrt(2.*M_PI))*(stepPhysX*stepPhysY*stepPhysZ)/(pow((float)(lengthX*lengthY*lengthZ),from.getNumberOfInverseTransforms())); 00921 double normalization=1./(2.*M_PI)/(pow((float)(lengthX*lengthY),from.getNumberOfInverseTransforms())); 00922 00923 Index x, y; 00924 signed int xp, yp; 00925 Vector2 r; 00926 typename TFFT2D<ComplexTraits>::Complex dataIn; 00927 typename TFFT2D<ComplexTraits>::Complex dataOut; 00928 00929 RegularData2D::IndexType current_index; 00930 for (Position i = 0; i < from.size(); i++) 00931 { 00932 y = i % lengthY; 00933 x = i / lengthY; 00934 00935 xp = x; 00936 yp = y; 00937 00938 if (xp>=lengthX/2.) 00939 { 00940 xp-=(int)lengthX; 00941 } 00942 if (yp>=lengthY/2.) 00943 { 00944 yp-=(int)lengthY; 00945 } 00946 00947 if (x>=lengthX/2.) 00948 { 00949 x-=(int)(lengthX/2.); 00950 } 00951 else 00952 { 00953 x+=(int)(lengthX/2.); 00954 } 00955 00956 if (y>=lengthY/2.) 00957 { 00958 y-=(int)(lengthY/2.); 00959 } 00960 else 00961 { 00962 y+=(int)(lengthY/2.); 00963 } 00964 00965 00966 r.set((float)xp * stepFourierX, 00967 (float)yp * stepFourierY); 00968 00969 dataIn = from[i]; 00970 dataOut = dataIn; 00971 00972 current_index.x = x; 00973 current_index.y = y; 00974 00975 newGrid[current_index] = (dataOut*(typename ComplexTraits::ComplexPrecision)normalization*from.phase(r)).real(); 00976 } 00977 00978 to = newGrid; 00979 00980 return to; 00981 } 00982 } 00983 00984 00985 00986 } 00987 00988 #endif // BALL_MATHS_TFFT2D_H