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

vigra/wigner-matrix.hxx |

    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.9.0 (Tue Nov 6 2012)