BALL
1.4.1
|
00001 // -*- Mode: C++; tab-width: 2; -*- 00002 // vi: set ts=2: 00003 // 00004 00005 #ifndef BALL_MATHS_TFFT3D_H 00006 #define BALL_MATHS_TFFT3D_H 00007 00008 #ifndef BALL_COMMON_EXCEPTION_H 00009 # include <BALL/COMMON/exception.h> 00010 #endif 00011 00012 00013 #ifndef BALL_DATATYPE_REGULARDATA3D_H 00014 # include <BALL/DATATYPE/regularData3D.h> 00015 #endif 00016 00017 //#ifndef BALL_MATHS_VECTOR2_H 00018 //# include <BALL/MATHS/vector3.h> 00019 //#endif 00020 00021 #include <BALL/MATHS/fftwCommon.h> 00022 #include <math.h> 00023 #include <complex> 00024 #include <fftw3.h> 00025 00026 namespace BALL 00027 { 00038 template <typename ComplexTraits> 00039 class TFFT3D 00040 : public TRegularData3D<std::complex<typename ComplexTraits::ComplexPrecision> > 00041 { 00042 public: 00043 00044 typedef std::complex<typename ComplexTraits::ComplexPrecision> Complex; 00045 typedef TRegularData3D<std::complex<typename ComplexTraits::ComplexPrecision> > ComplexVector; 00046 00047 BALL_CREATE(TFFT3D) 00048 00049 00052 00053 00054 TFFT3D(); 00055 00057 TFFT3D(const TFFT3D &data); 00058 00072 // AR: ldn is not any longer the binary logarithm but the absolute number of grid points 00073 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); 00074 00076 virtual ~TFFT3D(); 00077 00079 00083 00085 const TFFT3D& operator = (const TFFT3D& fft_3d); 00086 00089 virtual void clear(); 00090 00093 virtual void destroy(); 00094 00096 00100 00103 bool operator == (const TFFT3D& fft3d) const; 00105 00106 // @name Accessors 00107 00109 00111 void doFFT(); 00112 00115 void doiFFT(); 00116 00122 bool translate(const Vector3& trans_origin); 00123 00129 bool setPhysStepWidth(double new_width_x, double new_width_y, double new_width_z); 00130 00133 double getPhysStepWidthX() const; 00134 00137 double getPhysStepWidthY() const; 00138 00141 double getPhysStepWidthZ() const; 00142 00145 double getFourierStepWidthX() const; 00146 00149 double getFourierStepWidthY() const; 00150 00153 double getFourierStepWidthZ() const; 00154 00157 double getPhysSpaceMinX() const; 00158 00161 double getPhysSpaceMinY() const 00162 ; 00163 00166 double getPhysSpaceMinZ() const; 00167 00170 double getPhysSpaceMaxX() const; 00171 00174 double getPhysSpaceMaxY() const; 00175 00178 double getPhysSpaceMaxZ() const; 00179 00182 double getFourierSpaceMinX() const; 00183 00186 double getFourierSpaceMinY() const; 00187 00190 double getFourierSpaceMinZ() const; 00191 00194 double getFourierSpaceMaxX() const; 00195 00198 double getFourierSpaceMaxY() const; 00199 00202 double getFourierSpaceMaxZ() const; 00203 00209 Size getMaxXIndex() const; 00210 00216 Size getMaxYIndex() const; 00217 00223 Size getMaxZIndex() const; 00224 00228 Size getNumberOfInverseTransforms() const; 00229 00232 Vector3 getGridCoordinates(Position position) const; 00233 00240 Complex getData(const Vector3& pos) const; 00241 00249 Complex getInterpolatedValue(const Vector3& pos) const; 00250 00258 void setData(const Vector3& pos, Complex val); 00259 00265 Complex& operator[](const Vector3& pos); 00266 00272 const Complex& operator[](const Vector3& pos) const; 00273 00278 Complex& operator[](const Position& pos) 00279 { 00280 return TRegularData3D<Complex>::operator [] (pos); 00281 } 00282 00287 const Complex& operator[](const Position& pos) const 00288 { 00289 return TRegularData3D<Complex>::operator [] (pos); 00290 } 00291 00292 // AR: 00293 void setNumberOfFFTTransforms(Size num) 00294 { 00295 numPhysToFourier_ = num; 00296 } 00297 00298 // AR: 00299 void setNumberOfiFFTTransforms(Size num) 00300 { 00301 numFourierToPhys_ = num; 00302 } 00303 00304 00309 Complex phase(const Vector3& pos) const; 00311 00315 00319 bool isInFourierSpace() const; 00321 00322 protected: 00323 Size lengthX_, lengthY_, lengthZ_; 00324 bool inFourierSpace_; 00325 Size numPhysToFourier_; 00326 Size numFourierToPhys_; 00327 Vector3 origin_; 00328 double stepPhysX_, stepPhysY_, stepPhysZ_; 00329 double stepFourierX_, stepFourierY_, stepFourierZ_; 00330 Vector3 minPhys_, maxPhys_; 00331 Vector3 minFourier_, maxFourier_; 00332 00333 00334 // AR: new version for FFTW3 00335 typename ComplexTraits::FftwPlan planForward_; 00336 typename ComplexTraits::FftwPlan planBackward_; 00337 00338 // AR: to control plan calculation with new fftw3 00339 Size dataLength_; 00340 Complex *dataAdress_; 00341 bool planCalculated_; 00342 00343 }; 00344 00347 typedef TFFT3D<BALL_FFTW_DEFAULT_TRAITS> FFT3D; 00348 00351 template <typename ComplexTraits> 00352 const TRegularData3D<typename TFFT3D<ComplexTraits>::Complex>& operator << 00353 (TRegularData3D<typename TFFT3D<ComplexTraits>::Complex>& to, const TFFT3D<ComplexTraits>& from); 00354 00359 template <typename ComplexTraits> 00360 const RegularData3D& operator << (RegularData3D& to, const TFFT3D<ComplexTraits>& from); 00361 00362 00363 template <typename ComplexTraits> 00364 TFFT3D<ComplexTraits>::TFFT3D() 00365 : TRegularData3D<Complex>(), 00366 dataLength_(0), 00367 dataAdress_(0), 00368 planCalculated_(false) 00369 { 00370 } 00371 00372 template <typename ComplexTraits> 00373 bool TFFT3D<ComplexTraits>::operator == (const TFFT3D& fft3D) const 00374 { 00375 00376 // AR: test whether data_.size() == fft3D.data_.size() 00377 // instead of testing 3 lengths. Better for vector handling. 00378 00379 if (lengthX_ == fft3D.lengthX_ && 00380 lengthY_ == fft3D.lengthY_ && 00381 lengthZ_ == fft3D.lengthZ_ && 00382 origin_ == fft3D.origin_ && 00383 stepPhysX_ == fft3D.stepPhysX_ && 00384 stepPhysY_ == fft3D.stepPhysY_ && 00385 stepPhysZ_ == fft3D.stepPhysZ_ && 00386 stepFourierX_ == fft3D.stepFourierX_ && 00387 stepFourierY_ == fft3D.stepFourierY_ && 00388 stepFourierZ_ == fft3D.stepFourierZ_ && 00389 minPhys_ == fft3D.minPhys_ && 00390 maxPhys_ == fft3D.maxPhys_ && 00391 minFourier_ == fft3D.minFourier_ && 00392 maxFourier_ == fft3D.maxFourier_ && 00393 numPhysToFourier_ == fft3D.numPhysToFourier_ && 00394 numFourierToPhys_ == fft3D.numFourierToPhys_ && 00395 planCalculated_ == fft3D.planCalculated_) 00396 { 00397 Vector3 min = inFourierSpace_ ? minFourier_ : minPhys_; 00398 Vector3 max = inFourierSpace_ ? maxFourier_ : maxPhys_; 00399 double stepX = inFourierSpace_ ? stepFourierX_ : stepPhysX_; 00400 double stepY = inFourierSpace_ ? stepFourierY_ : stepPhysY_; 00401 double stepZ = inFourierSpace_ ? stepFourierZ_ : stepPhysZ_; 00402 00403 for (double posX=min.x; posX<=max.x; posX+=stepX) 00404 { 00405 for (double posY=min.y; posY<=max.y; posY+=stepY) 00406 { 00407 for (double posZ=min.z; posZ<=max.z; posZ+=stepZ) 00408 { 00409 if (getData(Vector3(posX,posY,posZ)) != fft3D.getData(Vector3(posX,posY,posZ))) 00410 { 00411 return false; 00412 } 00413 } 00414 } 00415 } 00416 00417 return true; 00418 } 00419 00420 return false; 00421 } 00422 00423 template <typename ComplexTraits> 00424 bool TFFT3D<ComplexTraits>::translate(const Vector3& trans_origin) 00425 { 00426 Position internalOriginX = (Position) Maths::rint(trans_origin.x*stepPhysX_); 00427 Position internalOriginY = (Position) Maths::rint(trans_origin.y*stepPhysY_); 00428 Position internalOriginZ = (Position) Maths::rint(trans_origin.z*stepPhysZ_); 00429 00430 if ((internalOriginX <= lengthX_) && (internalOriginY <= lengthY_) && (internalOriginZ <= lengthZ_)) 00431 { 00432 origin_.x = trans_origin.x; 00433 origin_.y = trans_origin.y; 00434 origin_.z = trans_origin.z; 00435 00436 minPhys_ = Vector3(-origin_.x,-origin_.y,-origin_.z); 00437 maxPhys_ = Vector3(((lengthX_-1)*stepPhysX_)-origin_.x,((lengthY_-1)*stepPhysY_)-origin_.y,((lengthZ_-1)*stepPhysZ_)-origin_.z); 00438 minFourier_ = Vector3(-(lengthX_/2.-1)*stepFourierX_,-(lengthY_/2.-1)*stepFourierY_,-(lengthZ_/2.-1)*stepFourierZ_); 00439 maxFourier_ = Vector3((lengthX_/2.)*stepFourierX_,(lengthY_/2.)*stepFourierY_,(lengthZ_/2.)*stepFourierZ_); 00440 00441 return true; 00442 } 00443 else 00444 { 00445 return false; 00446 } 00447 } 00448 00449 template <typename ComplexTraits> 00450 bool TFFT3D<ComplexTraits>::setPhysStepWidth(double new_width_x, double new_width_y, double new_width_z) 00451 { 00452 if ((new_width_x <= 0) || (new_width_y <= 0) || (new_width_z <= 0)) 00453 { 00454 return false; 00455 } 00456 else 00457 { 00458 stepPhysX_ = new_width_x; 00459 stepPhysY_ = new_width_y; 00460 stepPhysZ_ = new_width_z; 00461 stepFourierX_ = 2.*M_PI/(stepPhysX_*lengthX_); 00462 stepFourierY_ = 2.*M_PI/(stepPhysY_*lengthY_); 00463 stepFourierZ_ = 2.*M_PI/(stepPhysZ_*lengthZ_); 00464 00465 minPhys_ = Vector3(-origin_.x,-origin_.y,-origin_.z); 00466 maxPhys_ = Vector3(((lengthX_-1)*stepPhysX_)-origin_.x,((lengthY_-1)*stepPhysY_)-origin_.y,((lengthZ_-1)*stepPhysZ_)-origin_.z); 00467 minFourier_ = Vector3(-(lengthX_/2.-1)*stepFourierX_,-(lengthY_/2.-1)*stepFourierY_,-(lengthZ_/2.-1)*stepFourierZ_); 00468 maxFourier_ = Vector3((lengthX_/2.)*stepFourierX_,(lengthY_/2.)*stepFourierY_,(lengthZ_/2.)*stepFourierZ_); 00469 00470 return true; 00471 } 00472 } 00473 00474 template <typename ComplexTraits> 00475 double TFFT3D<ComplexTraits>::getPhysStepWidthX() const 00476 { 00477 return stepPhysX_; 00478 } 00479 00480 template <typename ComplexTraits> 00481 double TFFT3D<ComplexTraits>::getPhysStepWidthY() const 00482 { 00483 return stepPhysY_; 00484 } 00485 00486 template <typename ComplexTraits> 00487 double TFFT3D<ComplexTraits>::getPhysStepWidthZ() const 00488 { 00489 return stepPhysZ_; 00490 } 00491 00492 template <typename ComplexTraits> 00493 double TFFT3D<ComplexTraits>::getFourierStepWidthX() const 00494 { 00495 return stepFourierX_; 00496 } 00497 00498 template <typename ComplexTraits> 00499 double TFFT3D<ComplexTraits>::getFourierStepWidthY() const 00500 { 00501 return stepFourierY_; 00502 } 00503 00504 template <typename ComplexTraits> 00505 double TFFT3D<ComplexTraits>::getFourierStepWidthZ() const 00506 { 00507 return stepFourierZ_; 00508 } 00509 00510 template <typename ComplexTraits> 00511 double TFFT3D<ComplexTraits>::getPhysSpaceMinX() const 00512 { 00513 return minPhys_.x; 00514 } 00515 00516 template <typename ComplexTraits> 00517 double TFFT3D<ComplexTraits>::getPhysSpaceMinY() const 00518 { 00519 return minPhys_.y; 00520 } 00521 00522 template <typename ComplexTraits> 00523 double TFFT3D<ComplexTraits>::getPhysSpaceMinZ() const 00524 { 00525 return minPhys_.z; 00526 } 00527 00528 template <typename ComplexTraits> 00529 double TFFT3D<ComplexTraits>::getPhysSpaceMaxX() const 00530 { 00531 return maxPhys_.x; 00532 } 00533 00534 template <typename ComplexTraits> 00535 double TFFT3D<ComplexTraits>::getPhysSpaceMaxY() const 00536 { 00537 return maxPhys_.y; 00538 } 00539 00540 template <typename ComplexTraits> 00541 double TFFT3D<ComplexTraits>::getPhysSpaceMaxZ() const 00542 { 00543 return maxPhys_.z; 00544 } 00545 00546 template <typename ComplexTraits> 00547 double TFFT3D<ComplexTraits>::getFourierSpaceMinX() const 00548 00549 { 00550 return minFourier_.x; 00551 } 00552 00553 template <typename ComplexTraits> 00554 double TFFT3D<ComplexTraits>::getFourierSpaceMinY() const 00555 { 00556 return minFourier_.y; 00557 } 00558 00559 template <typename ComplexTraits> 00560 double TFFT3D<ComplexTraits>::getFourierSpaceMinZ() const 00561 { 00562 return minFourier_.z; 00563 } 00564 00565 template <typename ComplexTraits> 00566 double TFFT3D<ComplexTraits>::getFourierSpaceMaxX() const 00567 { 00568 return maxFourier_.x; 00569 } 00570 00571 template <typename ComplexTraits> 00572 double TFFT3D<ComplexTraits>::getFourierSpaceMaxY() const 00573 { 00574 return maxFourier_.y; 00575 } 00576 00577 template <typename ComplexTraits> 00578 double TFFT3D<ComplexTraits>::getFourierSpaceMaxZ() const 00579 { 00580 return maxFourier_.z; 00581 } 00582 00583 template <typename ComplexTraits> 00584 Size TFFT3D<ComplexTraits>::getMaxXIndex() const 00585 { 00586 return (lengthX_ - 1); 00587 } 00588 00589 template <typename ComplexTraits> 00590 Size TFFT3D<ComplexTraits>::getMaxYIndex() const 00591 { 00592 return (lengthY_ - 1); 00593 } 00594 00595 template <typename ComplexTraits> 00596 Size TFFT3D<ComplexTraits>::getMaxZIndex() const 00597 { 00598 return (lengthZ_ - 1); 00599 } 00600 00601 template <typename ComplexTraits> 00602 Size TFFT3D<ComplexTraits>::getNumberOfInverseTransforms() const 00603 { 00604 return numFourierToPhys_; 00605 } 00606 00607 template <typename ComplexTraits> 00608 Vector3 TFFT3D<ComplexTraits>::getGridCoordinates(Position position) const 00609 { 00610 if (!inFourierSpace_) 00611 { 00612 if (position >= ComplexVector::size()) 00613 { 00614 throw Exception::OutOfGrid(__FILE__, __LINE__); 00615 } 00616 00617 Vector3 r; 00618 Position x, y, z; 00619 00620 z = position % lengthZ_; 00621 y = (position % (lengthY_ * lengthZ_)) / lengthZ_; 00622 x = position / (lengthY_ * lengthZ_); 00623 00624 r.set(-origin_.x + (float)x * stepPhysX_, 00625 -origin_.y + (float)y * stepPhysY_, 00626 -origin_.z + (float)z * stepPhysZ_); 00627 00628 return r; 00629 } 00630 else 00631 { 00632 if (position >= ComplexVector::size()) 00633 { 00634 throw Exception::OutOfGrid(__FILE__, __LINE__); 00635 } 00636 00637 Vector3 r; 00638 Index x, y, z; 00639 00640 z = position % lengthZ_; 00641 y = (position % (lengthY_ * lengthZ_)) / lengthZ_; 00642 x = position / (lengthY_ * lengthZ_); 00643 00644 if (x>=lengthX_/2.) 00645 { 00646 x-=lengthX_; 00647 } 00648 00649 if (y>=lengthY_/2.) 00650 { 00651 y-=lengthY_; 00652 } 00653 00654 if (z>=lengthZ_/2.) 00655 { 00656 z-=lengthZ_; 00657 } 00658 00659 r.set((float)x * stepFourierX_, 00660 (float)y * stepFourierY_, 00661 (float)z * stepFourierZ_); 00662 00663 return r; 00664 } 00665 } 00666 00667 template <typename ComplexTraits> 00668 typename TFFT3D<ComplexTraits>::Complex TFFT3D<ComplexTraits>::getData(const Vector3& pos) const 00669 { 00670 Complex result; 00671 double normalization=1.; 00672 00673 if (!inFourierSpace_) 00674 { 00675 result = (*this)[pos]; 00676 normalization=1./((float)pow((float)(lengthX_*lengthY_*lengthZ_),(int)numFourierToPhys_)); 00677 } 00678 else 00679 { 00680 // AR: 00681 //old: result = (*this)[pos]; 00682 result = (*this)[pos]*phase(pos); 00683 00684 //normalization=1./pow(sqrt(2.*M_PI),3)*(stepPhysX_*stepPhysY_*stepPhysZ_)/((float)pow((float)(lengthX_*lengthY_*lengthZ_),(int)numFourierToPhys_)); 00685 00686 normalization=1./pow(sqrt(2.*M_PI),3)/((float)pow((float)(lengthX_*lengthY_*lengthZ_),(int)numFourierToPhys_)); 00687 //normalization=1./(sqrt(2.*M_PI))*(stepPhysX_*stepPhysY_*stepPhysZ_)/((float)pow((float)(lengthX_*lengthY_*lengthZ_),(int)numFourierToPhys_)); 00688 } 00689 00690 result *= normalization; 00691 00692 return result; 00693 } 00694 00695 template <typename ComplexTraits> 00696 typename TFFT3D<ComplexTraits>::Complex TFFT3D<ComplexTraits>::getInterpolatedValue(const Vector3& pos) const 00697 { 00698 Complex result; 00699 00700 Vector3 min = inFourierSpace_ ? minFourier_ : minPhys_; 00701 Vector3 max = inFourierSpace_ ? maxFourier_ : maxPhys_; 00702 double stepX = inFourierSpace_ ? stepFourierX_ : stepPhysX_; 00703 double stepY = inFourierSpace_ ? stepFourierY_ : stepPhysY_; 00704 double stepZ = inFourierSpace_ ? stepFourierZ_ : stepPhysZ_; 00705 00706 if ( (pos.x < min.x) || (pos.y < min.y) || (pos.z < min.z) 00707 || (pos.x > max.x) || (pos.y > max.y) || (pos.z > max.z) ) 00708 { 00709 throw Exception::OutOfGrid(__FILE__, __LINE__); 00710 } 00711 00712 Vector3 h(pos.x - min.x, pos.y - min.y, pos.z - min.z); 00713 double modX = fmod((double)h.x,stepX); 00714 double modY = fmod((double)h.y,stepY); 00715 double modZ = fmod((double)h.z,stepZ); 00716 00717 if (modX==0 && modY==0 && modZ==0) // we are on the grid 00718 { 00719 return getData(pos); 00720 } 00721 00722 double beforeX = floor(h.x/stepX)*stepX+min.x; 00723 double beforeY = floor(h.y/stepY)*stepY+min.y; 00724 double beforeZ = floor(h.z/stepZ)*stepZ+min.z; 00725 double afterX = ceil(h.x/stepX)*stepX+min.x; 00726 double afterY = ceil(h.y/stepY)*stepY+min.y; 00727 double afterZ = ceil(h.z/stepZ)*stepZ+min.z; 00728 00729 double tx = (pos.x - beforeX)/stepX; 00730 double ty = (pos.y - beforeY)/stepY; 00731 double tz = (pos.z - beforeZ)/stepZ; 00732 00733 result = getData(Vector3(beforeX,beforeY,beforeZ))*(typename ComplexTraits::ComplexPrecision)((1.-tx)*(1.-ty)*(1.-tz)); 00734 result += getData(Vector3(afterX, beforeY,beforeZ))*(typename ComplexTraits::ComplexPrecision)( tx *(1.-ty)*(1.-tz)); 00735 result += getData(Vector3(beforeX,afterY, beforeZ))*(typename ComplexTraits::ComplexPrecision)((1.-tx)* ty *(1.-tz)); 00736 result += getData(Vector3(beforeX,beforeY,afterZ ))*(typename ComplexTraits::ComplexPrecision)((1.-tx)*(1.-ty)* tz ); 00737 result += getData(Vector3(afterX, afterY, beforeZ))*(typename ComplexTraits::ComplexPrecision)( tx * ty *(1.-tz)); 00738 result += getData(Vector3(afterX, beforeY,afterZ ))*(typename ComplexTraits::ComplexPrecision)( tx *(1.-ty)* tz ); 00739 result += getData(Vector3(beforeX,afterY, afterZ ))*(typename ComplexTraits::ComplexPrecision)((1.-tx)* ty * tz ); 00740 result += getData(Vector3(afterX, afterY, afterZ ))*(typename ComplexTraits::ComplexPrecision)( tx * ty * tz ); 00741 00742 return result; 00743 } 00744 00745 template <typename ComplexTraits> 00746 void TFFT3D<ComplexTraits>::setData(const Vector3& pos, Complex val) 00747 { 00748 Complex dummy; 00749 00750 if (!inFourierSpace_) 00751 { 00752 dummy = Complex(val.real()*((typename ComplexTraits::ComplexPrecision)pow((typename ComplexTraits::ComplexPrecision)(lengthX_*lengthY_*lengthZ_),(int)numFourierToPhys_)), 00753 val.imag()*((typename ComplexTraits::ComplexPrecision)pow((typename ComplexTraits::ComplexPrecision)(lengthX_*lengthY_*lengthZ_),(int)numFourierToPhys_))); 00754 00755 (*this)[pos]=dummy; 00756 } 00757 else 00758 { 00759 val*=phase(pos)*(typename ComplexTraits::ComplexPrecision)((pow(sqrt(2*M_PI),3)/(stepPhysX_*stepPhysY_*stepPhysZ_))) 00760 *((typename ComplexTraits::ComplexPrecision)pow((typename ComplexTraits::ComplexPrecision)(lengthX_*lengthY_*lengthZ_),(int)numFourierToPhys_)); 00761 /*val*=phase(pos)*(typename ComplexTraits::ComplexPrecision)((sqrt(2*M_PI)/(stepPhysX_*stepPhysY_*stepPhysZ_))) 00762 *((typename ComplexTraits::ComplexPrecision)pow((typename ComplexTraits::ComplexPrecision)(lengthX_*lengthY_*lengthZ_),(int)numFourierToPhys_));*/ 00763 00764 dummy = val; 00765 00766 (*this)[pos]=dummy; 00767 } 00768 } 00769 00770 template <typename ComplexTraits> 00771 typename TFFT3D<ComplexTraits>::Complex& TFFT3D<ComplexTraits>::operator[](const Vector3& pos) 00772 { 00773 Index internalPos; 00774 00775 if (!inFourierSpace_) 00776 { 00777 Index i, j, k; 00778 00779 i = (Index) Maths::rint((pos.x+origin_.x)/stepPhysX_); 00780 j = (Index) Maths::rint((pos.y+origin_.y)/stepPhysY_); 00781 k = (Index) Maths::rint((pos.z+origin_.z)/stepPhysZ_); 00782 00783 internalPos = (k + (j + i*lengthY_)*lengthZ_); 00784 00785 /*(Index) rint( (pos.z+origin_.z)/stepPhysZ_ 00786 + ( (pos.y+origin_.y)/stepPhysY_ 00787 + (pos.x+origin_.x)/stepPhysX_*lengthY_ 00788 ) * lengthZ_ 00789 ); */ 00790 } 00791 else 00792 { 00793 Index i, j, k; 00794 00795 i = (Index) Maths::rint(pos.x/stepFourierX_); 00796 j = (Index) Maths::rint(pos.y/stepFourierY_); 00797 k = (Index) Maths::rint(pos.z/stepFourierZ_); 00798 00799 if (i<0) 00800 { 00801 i+=lengthX_; 00802 } 00803 00804 if (j<0) 00805 { 00806 j+=lengthY_; 00807 } 00808 00809 if (k<0) 00810 { 00811 k+=lengthZ_; 00812 } 00813 00814 internalPos = (k + (j + i*lengthY_)*lengthZ_); 00815 } 00816 00817 if ((internalPos < 0) || (internalPos>=(Index) (lengthX_*lengthY_*lengthZ_))) 00818 { 00819 throw Exception::OutOfGrid(__FILE__, __LINE__); 00820 } 00821 00822 return TRegularData3D<Complex>::operator[]((Position)internalPos); 00823 } 00824 00825 template <typename ComplexTraits> 00826 const typename TFFT3D<ComplexTraits>::Complex& TFFT3D<ComplexTraits>::operator[](const Vector3& pos) const 00827 { 00828 Index internalPos; 00829 00830 if (!inFourierSpace_) 00831 { 00832 Index i, j, k; 00833 00834 i = (Index) Maths::rint((pos.x+origin_.x)/stepPhysX_); 00835 j = (Index) Maths::rint((pos.y+origin_.y)/stepPhysY_); 00836 k = (Index) Maths::rint((pos.z+origin_.z)/stepPhysZ_); 00837 00838 internalPos = (k + (j + i*lengthY_)*lengthZ_); 00839 00840 /*(Index) rint( (pos.z+origin_.z)/stepPhysZ_ 00841 + ( (pos.y+origin_.y)/stepPhysY_ 00842 + (pos.x+origin_.x)/stepPhysX_*lengthY_ 00843 ) * lengthZ_ 00844 ); */ 00845 } 00846 else 00847 { 00848 Index i, j, k; 00849 00850 i = (Index) Maths::rint(pos.x/stepFourierX_); 00851 j = (Index) Maths::rint(pos.y/stepFourierY_); 00852 k = (Index) Maths::rint(pos.z/stepFourierZ_); 00853 00854 if (i<0) 00855 { 00856 i+=lengthX_; 00857 } 00858 00859 if (j<0) 00860 { 00861 j+=lengthY_; 00862 } 00863 00864 if (k<0) 00865 { 00866 k+=lengthZ_; 00867 } 00868 00869 internalPos = (k + (j + i*lengthY_)*lengthZ_); 00870 } 00871 00872 if ((internalPos < 0) || (internalPos>=(Index) (lengthX_*lengthY_*lengthZ_))) 00873 { 00874 throw Exception::OutOfGrid(__FILE__, __LINE__); 00875 } 00876 00877 return TRegularData3D<Complex>::operator[]((Position)internalPos); 00878 } 00879 00880 /*Complex& TFFT3D<ComplexTraits>::operator[](const Position& pos) 00881 { 00882 return operator [] (pos); 00883 } 00884 00885 const Complex& TFFT3D<ComplexTraits>::operator[](const Position& pos) const 00886 { 00887 return operator [] (pos); 00888 } 00889 */ 00890 template <typename ComplexTraits> 00891 typename TFFT3D<ComplexTraits>::Complex TFFT3D<ComplexTraits>::phase(const Vector3& pos) const 00892 { 00893 00894 // AR: old version: -2.*M_PI... 00895 double phase = 2.*M_PI*( (Maths::rint(pos.x/stepFourierX_))*(Maths::rint(origin_.x/stepPhysX_)) 00896 /lengthX_ 00897 + (Maths::rint(pos.y/stepFourierY_))*(Maths::rint(origin_.y/stepPhysY_)) 00898 /lengthY_ 00899 + (Maths::rint(pos.z/stepFourierZ_))*(Maths::rint(origin_.z/stepPhysZ_)) 00900 /lengthZ_ ); 00901 00902 00903 Complex result = Complex(cos(phase), sin(phase)); 00904 00905 return result; 00906 00907 /*double phase = -2.*M_PI*( (rint(pos.x/stepFourierX_))*(rint(origin_.x/stepPhysX_)) 00908 /lengthX_ 00909 + (Maths::rint(pos.y/stepFourierY_))*(Maths::rint(origin_.y/stepPhysY_)) 00910 /lengthY_ 00911 + (Maths::rint(pos.z/stepFourierZ_))*(Maths::rint(origin_.z/stepPhysZ_)) 00912 /lengthZ_ ); 00913 00914 00915 Complex result = Complex(cos(phase), sin(phase)); 00916 00917 return result;*/ 00918 } 00919 00920 template <typename ComplexTraits> 00921 bool TFFT3D<ComplexTraits>::isInFourierSpace() const 00922 { 00923 return inFourierSpace_; 00924 } 00925 00926 template <typename ComplexTraits> 00927 const TRegularData3D<typename TFFT3D<ComplexTraits>::Complex>& operator<< 00928 (TRegularData3D<typename TFFT3D<ComplexTraits>::Complex>& to, const TFFT3D<ComplexTraits>& from) 00929 { 00930 // first decide if the TFFT3D data is in Fourier space. 00931 if (!from.isInFourierSpace()) 00932 { 00933 // create a new grid 00934 Size lengthX = from.getMaxXIndex()+1; 00935 Size lengthY = from.getMaxYIndex()+1; 00936 Size lengthZ = from.getMaxZIndex()+1; 00937 00938 TRegularData3D<typename TFFT3D<ComplexTraits>::Complex> newGrid(TRegularData3D<typename TFFT3D<ComplexTraits>::Complex>::IndexType(lengthX, lengthY, lengthZ), 00939 Vector3(from.getPhysSpaceMinX(), from.getPhysSpaceMinY(), from.getPhysSpaceMinZ()), 00940 Vector3(from.getPhysSpaceMaxX(), from.getPhysSpaceMaxY(), from.getPhysSpaceMaxZ())); 00941 00942 // and fill it 00943 double normalization=1./(pow((float)(lengthX*lengthY*lengthZ),(int)from.getNumberOfInverseTransforms())); 00944 typename TFFT3D<ComplexTraits>::Complex dataIn; 00945 typename TFFT3D<ComplexTraits>::Complex dataOut; 00946 00947 for (Position i = 0; i < from.size(); i++) 00948 { 00949 Position x, y, z; 00950 00951 z = i % lengthZ; 00952 y = (i % (lengthY * lengthZ)) / lengthZ; 00953 x = i / (lengthY * lengthZ); 00954 00955 dataIn = from[i]; 00956 dataOut = dataIn; 00957 00958 newGrid[x + (y + z*lengthY)*lengthZ] = dataOut*(typename ComplexTraits::ComplexPrecision)normalization; 00959 } 00960 00961 to = newGrid; 00962 00963 return to; 00964 } 00965 else 00966 { 00967 // we are in Fourier space 00968 00969 // create a new grid 00970 Size lengthX = from.getMaxXIndex()+1; 00971 Size lengthY = from.getMaxYIndex()+1; 00972 Size lengthZ = from.getMaxZIndex()+1; 00973 //float stepPhysX = from.getPhysStepWidthX(); 00974 //float stepPhysY = from.getPhysStepWidthY(); 00975 //float stepPhysZ = from.getPhysStepWidthZ(); 00976 float stepFourierX = from.getFourierStepWidthX(); 00977 float stepFourierY = from.getFourierStepWidthY(); 00978 float stepFourierZ = from.getFourierStepWidthZ(); 00979 00980 00981 00982 TRegularData3D<typename TFFT3D<ComplexTraits>::Complex> newGrid(TRegularData3D<typename TFFT3D<ComplexTraits>::Complex>::IndexType(lengthX, lengthY, lengthZ), 00983 Vector3(from.getFourierSpaceMinX(), 00984 from.getFourierSpaceMinY(), 00985 from.getFourierSpaceMinZ()), 00986 Vector3(from.getFourierSpaceMaxX(), 00987 from.getFourierSpaceMaxY(), 00988 from.getFourierSpaceMaxZ())); 00989 00990 // and fill it 00991 // AR: old double normalization=1./(sqrt(2.*M_PI))*(stepPhysX*stepPhysY*stepPhysZ)/(pow((float)(lengthX*lengthY*lengthZ),from.getNumberOfInverseTransforms())); 00992 double normalization=1./pow(sqrt(2.*M_PI),3)/(pow((float)(lengthX*lengthY*lengthZ),(int)from.getNumberOfInverseTransforms())); 00993 00994 00995 Index x, y, z; 00996 Vector3 r; 00997 typename TFFT3D<ComplexTraits>::Complex dataIn; 00998 typename TFFT3D<ComplexTraits>::Complex dataOut; 00999 01000 for (Position i = 0; i < from.size(); i++) 01001 { 01002 z = i % lengthZ; 01003 y = (i % (lengthY * lengthZ)) / lengthZ; 01004 x = i / (lengthY * lengthZ); 01005 01006 if (x>lengthX/2.) 01007 { 01008 x-=lengthX; 01009 } 01010 01011 if (y>lengthY/2.) 01012 { 01013 y-=lengthY; 01014 } 01015 01016 if (z>lengthZ/2.) 01017 { 01018 z-=lengthZ; 01019 } 01020 01021 r.set((float)x * stepFourierX, 01022 (float)y * stepFourierY, 01023 (float)z * stepFourierZ); 01024 01025 dataIn = from[i]; 01026 dataOut = dataIn; 01027 01028 newGrid[x + (y + z*lengthY)*lengthZ] = dataOut*(typename ComplexTraits::ComplexPrecision)normalization*from.phase(r); 01029 } 01030 01031 to = newGrid; 01032 01033 return to; 01034 } 01035 } 01036 01037 template <typename ComplexTraits> 01038 const RegularData3D& operator << (RegularData3D& to, const TFFT3D<ComplexTraits>& from) 01039 { 01040 // first decide if the TFFT3D data is in Fourier space. 01041 if (!from.isInFourierSpace()) 01042 { 01043 // create a new grid 01044 Size lengthX = from.getMaxXIndex()+1; 01045 Size lengthY = from.getMaxYIndex()+1; 01046 Size lengthZ = from.getMaxZIndex()+1; 01047 01048 RegularData3D newGrid(RegularData3D::IndexType(lengthX, lengthY, lengthZ), Vector3(from.getPhysSpaceMinX(), from.getPhysSpaceMinY(), from.getPhysSpaceMinZ()), 01049 Vector3(from.getPhysSpaceMaxX(), from.getPhysSpaceMaxY(), from.getPhysSpaceMaxZ())); 01050 01051 // and fill it 01052 double normalization = 1./(pow((float)(lengthX*lengthY*lengthZ),(int)from.getNumberOfInverseTransforms())); 01053 typename TFFT3D<ComplexTraits>::Complex dataIn; 01054 typename TFFT3D<ComplexTraits>::Complex dataOut; 01055 01056 for (Position i = 0; i < from.size(); i++) 01057 { 01058 Position x, y, z; 01059 01060 z = i % lengthZ; 01061 y = (i % (lengthY * lengthZ)) / lengthZ; 01062 x = i / (lengthY * lengthZ); 01063 01064 dataIn = from[i]; 01065 dataOut = dataIn; 01066 01067 newGrid[x + (y + z*lengthY)*lengthZ] = dataOut.real()*normalization; 01068 } 01069 01070 to = newGrid; 01071 01072 return to; 01073 } 01074 else 01075 { 01076 // we are in Fourier space 01077 01078 // create a new grid 01079 Size lengthX = from.getMaxXIndex()+1; 01080 Size lengthY = from.getMaxYIndex()+1; 01081 Size lengthZ = from.getMaxZIndex()+1; 01082 //float stepPhysX = from.getPhysStepWidthX(); 01083 //float stepPhysY = from.getPhysStepWidthY(); 01084 //float stepPhysZ = from.getPhysStepWidthZ(); 01085 float stepFourierX = from.getFourierStepWidthX(); 01086 float stepFourierY = from.getFourierStepWidthY(); 01087 float stepFourierZ = from.getFourierStepWidthZ(); 01088 01089 01090 01091 RegularData3D newGrid(RegularData3D::IndexType(lengthX, lengthY, lengthZ), Vector3(from.getFourierSpaceMinX(), from.getFourierSpaceMinY(), from.getFourierSpaceMinZ()), Vector3(from.getFourierSpaceMaxX(), from.getFourierSpaceMaxY(), from.getFourierSpaceMaxZ())); 01092 01093 // and fill it 01094 // AR: old version double normalization=1./(sqrt(2.*M_PI))*(stepPhysX*stepPhysY*stepPhysZ)/(pow((float)(lengthX*lengthY*lengthZ),from.getNumberOfInverseTransforms())); 01095 double normalization=1./pow(sqrt(2.*M_PI),3)/(pow((float)(lengthX*lengthY*lengthZ),(int)from.getNumberOfInverseTransforms())); 01096 01097 Index x, y, z; 01098 signed int xp, yp, zp; 01099 Vector3 r; 01100 typename TFFT3D<ComplexTraits>::Complex dataIn; 01101 typename TFFT3D<ComplexTraits>::Complex dataOut; 01102 01103 for (Position i = 0; i < from.size(); i++) 01104 { 01105 z = i % lengthZ; 01106 y = (i % (lengthY * lengthZ)) / lengthZ; 01107 x = i / (lengthY * lengthZ); 01108 01109 xp = x; 01110 yp = y; 01111 zp = z; 01112 01113 if (xp>=lengthX/2.) 01114 { 01115 xp-=(int)lengthX; 01116 } 01117 if (yp>=lengthY/2.) 01118 { 01119 yp-=(int)lengthY; 01120 } 01121 if (zp>=lengthZ/2.) 01122 { 01123 zp-=(int)lengthZ; 01124 } 01125 01126 if (x>=lengthX/2.) 01127 { 01128 x-=(int)(lengthX/2.); 01129 } 01130 else 01131 { 01132 x+=(int)(lengthX/2.); 01133 } 01134 01135 if (y>=lengthY/2.) 01136 { 01137 y-=(int)(lengthY/2.); 01138 } 01139 else 01140 { 01141 y+=(int)(lengthY/2.); 01142 } 01143 01144 if (z>=lengthZ/2.) 01145 { 01146 z-=(int)(lengthZ/2.); 01147 } 01148 else 01149 { 01150 z+=(int)(lengthZ/2.); 01151 } 01152 01153 r.set((float)xp * stepFourierX, 01154 (float)yp * stepFourierY, 01155 (float)zp * stepFourierZ); 01156 01157 dataIn = from[i]; 01158 dataOut = dataIn; 01159 01160 newGrid[x + (y + z*lengthY)*lengthZ] = (dataOut*(typename ComplexTraits::ComplexPrecision)normalization*from.phase(r)).real(); 01161 } 01162 01163 to = newGrid; 01164 01165 return to; 01166 } 01167 } 01168 } 01169 01170 #endif // BALL_MATHS_TFFT3D_H