1 /************************************************************************/
2 /* */
3 /* Copyright 2009-2010 by Ullrich Koethe and Janis Fehr */
4 /* */
5 /* This file is part of the VIGRA computer vision library. */
6 /* The VIGRA Website is */
7 /* http://hci.iwr.uni-heidelberg.de/vigra/ */
8 /* Please direct questions, bug reports, and contributions to */
9 /* ullrich.koethe@iwr.uni-heidelberg.de or */
10 /* vigra@informatik.uni-hamburg.de */
11 /* */
12 /* Permission is hereby granted, free of charge, to any person */
13 /* obtaining a copy of this software and associated documentation */
14 /* files (the "Software"), to deal in the Software without */
15 /* restriction, including without limitation the rights to use, */
16 /* copy, modify, merge, publish, distribute, sublicense, and/or */
17 /* sell copies of the Software, and to permit persons to whom the */
18 /* Software is furnished to do so, subject to the following */
19 /* conditions: */
20 /* */
21 /* The above copyright notice and this permission notice shall be */
22 /* included in all copies or substantial portions of the */
23 /* Software. */
24 /* */
25 /* THE SOFTWARE IS PROVIDED "AS IS", WITHOUT WARRANTY OF ANY KIND */
26 /* EXPRESS OR IMPLIED, INCLUDING BUT NOT LIMITED TO THE WARRANTIES */
27 /* OF MERCHANTABILITY, FITNESS FOR A PARTICULAR PURPOSE AND */
28 /* NONINFRINGEMENT. IN NO EVENT SHALL THE AUTHORS OR COPYRIGHT */
29 /* HOLDERS BE LIABLE FOR ANY CLAIM, DAMAGES OR OTHER LIABILITY, */
30 /* WHETHER IN AN ACTION OF CONTRACT, TORT OR OTHERWISE, ARISING */
31 /* FROM, OUT OF OR IN CONNECTION WITH THE SOFTWARE OR THE USE OR */
32 /* OTHER DEALINGS IN THE SOFTWARE. */
33 /* */
34 /************************************************************************/
35
36 #ifndef VIGRA_WIGNER_MATRIX_HXX
37 #define VIGRA_WIGNER_MATRIX_HXX
38
39 #include <complex>
40 #include "config.hxx"
41 #include "error.hxx"
42 #include "utilities.hxx"
43 #include "mathutil.hxx"
44 #include "array_vector.hxx"
45 #include "matrix.hxx"
46 #include "tinyvector.hxx"
47 #include "quaternion.hxx"
48 #include "clebsch-gordan.hxx"
49
50 namespace vigra {
51
52 /**
53 * \class WignerMatrix
54 * \brief computation of Wigner D matrix + rotation functions
55 * in SH,VH and R�
56 *
57 * All rotations in Euler zyz' convention
58 *
59 * WARNING: not thread safe! use a new instance of WignerMatrix
60 * for each thread!!!
61 */
62 template <class Real>
63 class WignerMatrix
64 {
65 public:
66
67 // FIXME: should we rather use FFTWComplex?
68 typedef std::complex<Real> Complex;
69 typedef ArrayVector<ArrayVector<ArrayVector<Complex> > > NestedArray;
70
71 /** \brief constructor
72 *
73 * \param l_max maximum expansion band (used to pre-compute
74 * the D matrix)
75 *
76 */
77 WignerMatrix(int l_max);
78
79 /** \brief Compute D with fixed theta = pi/2, phi=0, psi=0.
80 *
81 * \param band expansion band
82 *
83 FIXME: compute_D(l, 0.0, M_PI / 2.0, 0.0) creates the transposed matrix!
84 */
85 void compute_D(int band);
86
87 /** \brief Compute D for arbitrary rotations.
88 *
89 * \param l expansion band
90 * \param phi rotation angle
91 * \param theta rotation angle
92 * \param psi rotation angle
93 *
94 */
95 void compute_D(int l, Real phi, Real theta, Real psi);
96
97 /** \brief Get the (n,m) entry of D.
98 *
99 * \param l expansion band
100 * \param n
101 * \param m
102 */
get_D(int l,int n,int m) const103 Complex get_D(int l, int n, int m) const
104 {
105 if (l>0)
106 {
107 std::string message = std::string("WignerMatrix::get_D(): index out of bounds: l=");
108 message << l << " l_max=" << D.size() << " m=" << m << " n=" << n << "\n";
109
110 vigra_precondition(l < D.size() && m+l <= 2*l+1 &&
111 n+l <= 2*l+1 && m+l >= 0 && n+l >= 0,
112 message.c_str());
113 return D[l](n+l, m+l);
114 }
115 else
116 {
117 return Complex(Real(1.0));
118 }
119 }
120
121 /** \brief Return the rotation matrix D for the lth band.
122 *
123 * \param l expansion band
124 */
get_D(int l) const125 Matrix<Complex> const & get_D(int l) const
126 {
127 std::string message = std::string("WignerMatrix::get_D(): index out of bounds: l=");
128 message << l << " l_max=" << l_max << "\n";
129
130 vigra_precondition(l > 0 && l <= l_max, message.c_str());
131 return D[l];
132 }
133
134 /** \brief Rotate in PH.
135 *
136 * \param PH input PH expansion
137 * \param phi rotation angle
138 * \param theta rotation angle
139 * \param psi rotation angle
140 *
141 * \retval PHresult PH expansion
142 */
143 void rotatePH(NestedArray const & PH, Real phi, Real theta, Real psi,
144 NestedArray & PHresult);
145
146
147 private:
148 // FIXME: is function is not called (and cannot be called from outside the class)
149 TinyVector<double,3>
rot(TinyVector<double,3> const & vec,TinyVector<double,3> const & axis,double angle)150 rot(TinyVector<double,3> const & vec, TinyVector<double,3> const & axis, double angle)
151 {
152 typedef Quaternion<double> Q;
153 Q qr = Q::createRotation(angle, axis),
154 qv(0.0, vec);
155 return (qr * qv * conj(qr)).v();
156 }
157
158 int l_max;
159 int cg1cnt;
160 ArrayVector<double> CGcoeff;
161 ArrayVector<Matrix<Complex> > D;
162 };
163
164 template <class Real>
WignerMatrix(int band)165 WignerMatrix<Real>::WignerMatrix(int band)
166 : l_max(0),
167 cg1cnt(0)
168 {
169 //precompute clebschGordan coeffs
170 for (int l = 2; l <= band+1; l++)
171 {
172 for(int m = -l; m <= l ; m++)
173 {
174 for(int n = -l; n <= l ; n++)
175 {
176 for (int m2 = -1; m2 <= 1; m2++)
177 {
178 for (int n2 = -1; n2 <= 1; n2++)
179 {
180 int m1 = m-m2;
181 int n1 = n-n2;
182 if (m1 > -l && m1 < l && n1 > -l && n1 < l)
183 {
184 CGcoeff.push_back((clebschGordan(l-1,m1,1,m2,l,m))*(clebschGordan(l-1,n1,1,n2,l,n)));
185 }
186 }
187 }
188 }
189 }
190 }
191 }
192
193 template <class Real>
194 void
compute_D(int l,Real phi,Real theta,Real psi)195 WignerMatrix<Real>::compute_D(int l, Real phi, Real theta, Real psi)
196 {
197 double s = std::sin(theta);
198 double c = std::cos(theta);
199
200 Complex i(0.0, 1.0);
201 Complex eiphi = std::exp(i*phi);
202 Complex emiphi = std::exp(-i*phi);
203 Complex eipsi = std::exp(i*psi);
204 Complex emipsi = std::exp(-i*psi);
205
206 if (D.size() < (std::size_t)(l+1))
207 D.resize(l+1);
208 D[1].reshape(MultiArrayShape<2>::type(3,3));
209
210 D[1](0,0) = emipsi * Complex(Real(0.5*(1.0+c))) * emiphi;
211 D[1](0,1) = Complex(Real(-s/M_SQRT2)) * emiphi;
212 D[1](0,2) = eipsi * Complex(Real(0.5*(1.0-c))) * emiphi;
213 D[1](1,0) = emipsi * Complex(Real(s/M_SQRT2));
214 D[1](1,1) = Complex(Real(c));
215 D[1](1,2) = eipsi * Complex(Real(-s/M_SQRT2));
216 D[1](2,0) = emipsi * Complex(Real(0.5*(1.0-c))) * eiphi;
217 D[1](2,1) = Complex(Real(s/M_SQRT2)) * eiphi;
218 D[1](2,2) = eipsi * Complex(Real(0.5*(1.0+c))) * eiphi;
219
220 l_max = 1;
221 cg1cnt = 0;
222 if(l > 1)
223 compute_D( l);
224 }
225
226
227 template <class Real>
228 void
compute_D(int l)229 WignerMatrix<Real>::compute_D(int l)
230 {
231 if (D.size() < (std::size_t)(l+1) )
232 {
233 D.resize(l+1);
234 l_max = 0;
235 }
236
237 if (l==1)
238 {
239 //precompute D0 =1 and D1 = (90 degree rot)
240 // FIXME: signs are inconsistent with above explicit formula for
241 // theta = pi/2, phi=0, psi=0 (sine terms should be negated)
242 D[1].reshape(MultiArrayShape<2>::type(3,3));
243 D[1](0,0) = Real(0.5);
244 D[1](0,1) = Real(0.5*M_SQRT2);
245 D[1](0,2) = Real(0.5);
246 D[1](1,0) = Real(-0.5*M_SQRT2);
247 D[1](1,1) = Real(0.0);
248 D[1](1,2) = Real(0.5*M_SQRT2);
249 D[1](2,0) = Real(0.5);
250 D[1](2,1) = Real(-0.5*M_SQRT2);
251 D[1](2,2) = Real(0.5);
252 l_max = 1;
253 cg1cnt = 0;
254 }
255 else
256 {
257 //compute D2-Dl_max recursive
258 if (l>l_max+1)
259 {
260 compute_D(l-1);
261 }
262
263 D[l].reshape(MultiArrayShape<2>::type(2*l+1,2*l+1));
264 D[l].init(Real(0.0));
265
266 for(int m = -l; m <= l ; m++)
267 {
268 for(int n = -l; n <= l ; n++)
269 {
270 for (int m2 = -1; m2 <= 1; m2++)
271 {
272 for (int n2 = -1; n2 <= 1; n2++)
273 {
274 int m1 = m-m2;
275 int n1 = n-n2;
276 if ((m1 > -l) && (m1 < l) && (n1 > -l) && (n1 < l))
277 {
278 D[l](m+l,n+l) += D[1](m2+1,n2+1) * D[l-1](m1+l-1,n1+l-1) * Real(CGcoeff[cg1cnt++]);
279 }
280 }
281 }
282 }
283 }
284
285 l_max = l;
286 }
287 }
288
289 template <class Real>
290 void
rotatePH(NestedArray const & PH,Real phi,Real theta,Real psi,NestedArray & PHresult)291 WignerMatrix<Real>::rotatePH(NestedArray const & PH, Real phi, Real theta, Real psi,
292 NestedArray & PHresult)
293 {
294 int band = PH[1].size()-1;
295 compute_D(band, phi, theta, psi);
296
297 PHresult.resize(PH.size());
298
299 for(int n=1; n<=band; n++)
300 {
301 PHresult[n].resize(band+1);
302 for (int l=0; l<=band; l++)
303 {
304 PHresult[n][l].resize(2*band+1);
305 for(int m=-l; m<=l; m++)
306 {
307 Complex tmp = 0;
308 for (int h=-l; h<=l; h++)
309 {
310 tmp += get_D(l,h,m) * PH[n][l][h+l];
311 }
312
313 PHresult[n][l][m+l] = tmp;
314 }
315 }
316 }
317 }
318
319 } // namespace vigra
320
321 #endif // VIGRA_WIGNER_MATRIX_HXX
322