[ 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) |
html generated using doxygen and Python
|