[ VIGRA Homepage | Function Index | Class Index | Namespaces | File List | Main Page ]

vigra/wigner-matrix.hxx VIGRA

00001 /************************************************************************/
00002 /*                                                                      */
00003 /*        Copyright 2009-2010 by Ullrich Koethe and Janis Fehr          */
00004 /*                                                                      */
00005 /*    This file is part of the VIGRA computer vision library.           */
00006 /*    The VIGRA Website is                                              */
00007 /*        http://hci.iwr.uni-heidelberg.de/vigra/                       */
00008 /*    Please direct questions, bug reports, and contributions to        */
00009 /*        ullrich.koethe@iwr.uni-heidelberg.de    or                    */
00010 /*        vigra@informatik.uni-hamburg.de                               */
00011 /*                                                                      */
00012 /*    Permission is hereby granted, free of charge, to any person       */
00013 /*    obtaining a copy of this software and associated documentation    */
00014 /*    files (the "Software"), to deal in the Software without           */
00015 /*    restriction, including without limitation the rights to use,      */
00016 /*    copy, modify, merge, publish, distribute, sublicense, and/or      */
00017 /*    sell copies of the Software, and to permit persons to whom the    */
00018 /*    Software is furnished to do so, subject to the following          */
00019 /*    conditions:                                                       */
00020 /*                                                                      */
00021 /*    The above copyright notice and this permission notice shall be    */
00022 /*    included in all copies or substantial portions of the             */
00023 /*    Software.                                                         */
00024 /*                                                                      */
00025 /*    THE SOFTWARE IS PROVIDED "AS IS", WITHOUT WARRANTY OF ANY KIND    */
00026 /*    EXPRESS OR IMPLIED, INCLUDING BUT NOT LIMITED TO THE WARRANTIES   */
00027 /*    OF MERCHANTABILITY, FITNESS FOR A PARTICULAR PURPOSE AND          */
00028 /*    NONINFRINGEMENT. IN NO EVENT SHALL THE AUTHORS OR COPYRIGHT       */
00029 /*    HOLDERS BE LIABLE FOR ANY CLAIM, DAMAGES OR OTHER LIABILITY,      */
00030 /*    WHETHER IN AN ACTION OF CONTRACT, TORT OR OTHERWISE, ARISING      */
00031 /*    FROM, OUT OF OR IN CONNECTION WITH THE SOFTWARE OR THE USE OR     */
00032 /*    OTHER DEALINGS IN THE SOFTWARE.                                   */
00033 /*                                                                      */
00034 /************************************************************************/
00035 
00036 #ifndef VIGRA_WIGNER_MATRIX_HXX
00037 #define VIGRA_WIGNER_MATRIX_HXX
00038 
00039 #include <complex>
00040 #include "config.hxx"
00041 #include "error.hxx"
00042 #include "utilities.hxx"
00043 #include "mathutil.hxx"
00044 #include "array_vector.hxx"
00045 #include "matrix.hxx"
00046 #include "tinyvector.hxx"
00047 #include "quaternion.hxx"
00048 #include "clebsch-gordan.hxx"
00049 
00050 namespace vigra {
00051 
00052     /*!
00053      *  \class WignerMatrix 
00054      *  \brief computation of Wigner D matrix + rotation functions 
00055      *         in SH,VH and R³
00056      *
00057      * All rotations in Euler zyz' convention   
00058      *
00059      * WARNING: not thread safe! use a new instance of WignerMatrix
00060      * for each thread!!!
00061      */
00062 template <class Real>
00063 class WignerMatrix
00064 {   
00065   public:
00066   
00067     // FIXME: should we rather use FFTWComplex?
00068     typedef std::complex<Real> Complex;
00069     typedef ArrayVector<ArrayVector<ArrayVector<Complex> > > NestedArray;
00070     
00071          /*! \brief constructor
00072           * 
00073           * \param l_max    maximum expansion band (used to pre-compute 
00074           *         the D matrix)
00075           *
00076           */
00077     WignerMatrix(int l_max);
00078     
00079          /*! \brief Compute D with fixed theta = pi/2, phi=0, psi=0.
00080           *
00081           * \param band expansion band
00082           *
00083              FIXME: compute_D(l, 0.0, M_PI / 2.0, 0.0) creates the transposed matrix!
00084           */
00085     void compute_D(int band);
00086 
00087          /*! \brief Compute D for arbitrary rotations.
00088           *
00089           * \param l        expansion band
00090           * \param phi  rotation angle  
00091           * \param theta    rotation angle
00092           * \param psi  rotation angle
00093           *
00094           */
00095     void compute_D(int l, Real phi, Real theta, Real psi);
00096 
00097          /*! \brief  Get the (n,m) entry of D.
00098           *
00099           * \param l        expansion band
00100           * \param n        
00101           * \param m    
00102           */
00103     Complex get_D(int l, int n, int m) const
00104     {
00105         if (l>0)
00106         {
00107             std::string message = std::string("WignerMatrix::get_D(): index out of bounds: l=");
00108             message << l << " l_max=" << D.size() << " m=" << m << " n=" << n << "\n";
00109                                   
00110             vigra_precondition(l < D.size() && m+l <= 2*l+1 &&
00111                                n+l <= 2*l+1 && m+l >= 0 && n+l >= 0,
00112                                message.c_str());
00113             return D[l](n+l, m+l);
00114         }
00115         else 
00116         {
00117             return Complex(Real(1.0));
00118         }
00119     }
00120 
00121          /*! \brief Return the rotation matrix D for the lth band.
00122           *
00123           * \param l        expansion band
00124           */
00125     Matrix<Complex> const & get_D(int l) const
00126     {
00127         std::string message = std::string("WignerMatrix::get_D(): index out of bounds: l=");
00128         message << l << " l_max=" << l_max << "\n";
00129                               
00130         vigra_precondition(l > 0 && l <= l_max, message.c_str());
00131         return D[l];
00132     }
00133 
00134          /*! \brief Rotate in PH.
00135           *
00136           * \param PH       input PH expansion 
00137           * \param phi      rotation angle
00138           * \param theta    rotation angle
00139           * \param psi      rotation angle
00140           *
00141           * \retval PHresult PH expansion   
00142           */
00143     void rotatePH(NestedArray const & PH, Real phi, Real theta, Real psi,
00144                   NestedArray & PHresult);
00145 
00146 
00147   private:
00148     // FIXME: is function is not called (and cannot be called from outside the class)
00149     TinyVector<double,3> 
00150     rot(TinyVector<double,3> const & vec, TinyVector<double,3> const & axis, double angle)
00151     {
00152         typedef Quaternion<double> Q;
00153         Q qr = Q::createRotation(angle, axis),
00154           qv(0.0, vec);
00155         return (qr * qv * conj(qr)).v();
00156     }
00157 
00158     int l_max;
00159     int cg1cnt;
00160     ArrayVector<double> CGcoeff;
00161     ArrayVector<Matrix<Complex> > D;
00162 };
00163 
00164 template <class Real>
00165 WignerMatrix<Real>::WignerMatrix(int band)
00166 : l_max(0),
00167   cg1cnt(0)
00168 {
00169     //precompute clebschGordan coeffs
00170     for (int l = 2; l <= band+1; l++)
00171     {
00172         for(int m = -l; m <= l ; m++)
00173         {
00174             for(int n = -l; n <= l ; n++)
00175             {
00176                 for (int m2 = -1; m2 <= 1; m2++)
00177                 {
00178                     for (int n2 = -1; n2 <= 1; n2++)
00179                     {
00180                         int m1 = m-m2;
00181                         int n1 = n-n2;
00182                         if (m1 > -l && m1 < l && n1 > -l && n1 < l)
00183                         {
00184                             CGcoeff.push_back((clebschGordan(l-1,m1,1,m2,l,m))*(clebschGordan(l-1,n1,1,n2,l,n)));
00185                         }
00186                     }
00187                 }
00188             }
00189         }
00190     }
00191 }
00192 
00193 template <class Real>
00194 void
00195 WignerMatrix<Real>::compute_D(int l, Real phi, Real theta, Real psi)
00196 {
00197     double s = std::sin(theta);
00198     double c = std::cos(theta);
00199     
00200     Complex i(0.0, 1.0);
00201     Complex eiphi = std::exp(i*phi);
00202     Complex emiphi = std::exp(-i*phi);
00203     Complex eipsi = std::exp(i*psi);
00204     Complex emipsi = std::exp(-i*psi);
00205     
00206     if (D.size() < (std::size_t)(l+1)) 
00207         D.resize(l+1);
00208     D[1].reshape(MultiArrayShape<2>::type(3,3));
00209     
00210     D[1](0,0) = emipsi * Complex(Real(0.5*(1.0+c))) * emiphi;
00211     D[1](0,1) = Complex(Real(-s/M_SQRT2)) * emiphi;
00212     D[1](0,2) = eipsi * Complex(Real(0.5*(1.0-c))) * emiphi;
00213     D[1](1,0) = emipsi * Complex(Real(s/M_SQRT2));
00214     D[1](1,1) = Complex(Real(c));
00215     D[1](1,2) = eipsi * Complex(Real(-s/M_SQRT2));
00216     D[1](2,0) = emipsi * Complex(Real(0.5*(1.0-c))) * eiphi; 
00217     D[1](2,1) = Complex(Real(s/M_SQRT2)) * eiphi;
00218     D[1](2,2) = eipsi * Complex(Real(0.5*(1.0+c))) * eiphi;
00219 
00220     l_max = 1;
00221     cg1cnt = 0;
00222     if(l > 1)
00223         compute_D( l);
00224 }
00225 
00226 
00227 template <class Real>
00228 void
00229 WignerMatrix<Real>::compute_D(int l)
00230 {
00231     if (D.size() < (std::size_t)(l+1) ) 
00232     {
00233         D.resize(l+1);
00234         l_max = 0;
00235     }
00236 
00237     if (l==1)
00238     {
00239         //precompute D0 =1 and D1 = (90 degree rot)
00240         // FIXME: signs are inconsistent with above explicit formula for 
00241         //        theta = pi/2, phi=0, psi=0 (sine terms should be negated)
00242         D[1].reshape(MultiArrayShape<2>::type(3,3));
00243         D[1](0,0) = Real(0.5);
00244         D[1](0,1) = Real(0.5*M_SQRT2);
00245         D[1](0,2) = Real(0.5);
00246         D[1](1,0) = Real(-0.5*M_SQRT2);
00247         D[1](1,1) = Real(0.0);
00248         D[1](1,2) = Real(0.5*M_SQRT2);
00249         D[1](2,0) = Real(0.5);
00250         D[1](2,1) = Real(-0.5*M_SQRT2);
00251         D[1](2,2) = Real(0.5);
00252         l_max = 1;
00253         cg1cnt = 0;
00254     }
00255     else
00256     {
00257         //compute D2-Dl_max recursive
00258         if (l>l_max+1)
00259         {
00260             compute_D(l-1);
00261         }
00262 
00263         D[l].reshape(MultiArrayShape<2>::type(2*l+1,2*l+1));
00264         D[l].init(Real(0.0));      
00265         
00266         for(int m = -l; m <= l ; m++)
00267         {
00268             for(int n = -l; n <= l ; n++)
00269             {
00270                 for (int m2 = -1; m2 <= 1; m2++)
00271                 {
00272                     for (int n2 = -1; n2 <= 1; n2++)
00273                     {
00274                         int m1 = m-m2;
00275                         int n1 = n-n2;
00276                         if ((m1 > -l) && (m1 < l) && (n1 > -l) && (n1 < l))
00277                         {
00278                             D[l](m+l,n+l) += D[1](m2+1,n2+1) * D[l-1](m1+l-1,n1+l-1) * Real(CGcoeff[cg1cnt++]);
00279                         }
00280                     }
00281                 }
00282             }
00283         }
00284 
00285         l_max = l;
00286     }
00287 }
00288 
00289 template <class Real>
00290 void 
00291 WignerMatrix<Real>::rotatePH(NestedArray const & PH, Real phi, Real theta, Real psi,
00292                              NestedArray & PHresult)
00293 {
00294     int band = PH[1].size()-1;
00295     compute_D(band, phi, theta, psi);
00296 
00297     PHresult.resize(PH.size());
00298 
00299     for(int n=1; n<=band; n++)
00300     {
00301         PHresult[n].resize(band+1);
00302         for (int l=0; l<=band; l++)
00303         {
00304             PHresult[n][l].resize(2*band+1);
00305             for(int m=-l; m<=l; m++)
00306             {
00307                 Complex tmp = 0;
00308                 for (int h=-l; h<=l; h++)
00309                 {
00310                     tmp += get_D(l,h,m) * PH[n][l][h+l];
00311                 }
00312 
00313                 PHresult[n][l][m+l] = tmp;
00314             }
00315         }
00316     }
00317 }
00318 
00319 } // namespace vigra 
00320 
00321 #endif // VIGRA_WIGNER_MATRIX_HXX

© Ullrich Köthe (ullrich.koethe@iwr.uni-heidelberg.de)
Heidelberg Collaboratory for Image Processing, University of Heidelberg, Germany

html generated using doxygen and Python
vigra 1.8.0 (20 Sep 2011)