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

wigner-matrix.hxx VIGRA

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³ * * All rotations in Euler zyz' convention * * WARNING: not thread safe! use a new instance of WignerMatrix * for each thread!!! */ template <class Real> class WignerMatrix { public: // FIXME: should we rather use FFTWComplex? typedef std::complex<Real> Complex; typedef ArrayVector<ArrayVector<ArrayVector<Complex> > > NestedArray; /** \brief constructor * * \param l_max maximum expansion band (used to pre-compute * the D matrix) * */ WignerMatrix(int l_max); /** \brief Compute D with fixed theta = pi/2, phi=0, psi=0. * * \param band expansion band * FIXME: compute_D(l, 0.0, M_PI / 2.0, 0.0) creates the transposed matrix! */ void compute_D(int band); /** \brief Compute D for arbitrary rotations. * * \param l expansion band * \param phi rotation angle * \param theta rotation angle * \param psi rotation angle * */ void compute_D(int l, Real phi, Real theta, Real psi); /** \brief Get the (n,m) entry of D. * * \param l expansion band * \param n * \param m */ Complex get_D(int l, int n, int m) const { if (l>0) { std::string message = std::string("WignerMatrix::get_D(): index out of bounds: l="); message << l << " l_max=" << D.size() << " m=" << m << " n=" << n << "\n"; vigra_precondition(l < D.size() && m+l <= 2*l+1 && n+l <= 2*l+1 && m+l >= 0 && n+l >= 0, message.c_str()); return D[l](n+l, m+l); } else { return Complex(Real(1.0)); } } /** \brief Return the rotation matrix D for the lth band. * * \param l expansion band */ Matrix<Complex> const & get_D(int l) const { std::string message = std::string("WignerMatrix::get_D(): index out of bounds: l="); message << l << " l_max=" << l_max << "\n"; vigra_precondition(l > 0 && l <= l_max, message.c_str()); return D[l]; } /** \brief Rotate in PH. * * \param PH input PH expansion * \param phi rotation angle * \param theta rotation angle * \param psi rotation angle * * \retval PHresult PH expansion */ void rotatePH(NestedArray const & PH, Real phi, Real theta, Real psi, NestedArray & PHresult); private: // FIXME: is function is not called (and cannot be called from outside the class) TinyVector<double,3> rot(TinyVector<double,3> const & vec, TinyVector<double,3> const & axis, double angle) { typedef Quaternion<double> Q; Q qr = Q::createRotation(angle, axis), qv(0.0, vec); return (qr * qv * conj(qr)).v(); } int l_max; int cg1cnt; ArrayVector<double> CGcoeff; ArrayVector<Matrix<Complex> > D; }; template <class Real> WignerMatrix<Real>::WignerMatrix(int band) : l_max(0), cg1cnt(0) { //precompute clebschGordan coeffs for (int l = 2; l <= band+1; l++) { for(int m = -l; m <= l ; m++) { for(int n = -l; n <= l ; n++) { for (int m2 = -1; m2 <= 1; m2++) { for (int n2 = -1; n2 <= 1; n2++) { int m1 = m-m2; int n1 = n-n2; if (m1 > -l && m1 < l && n1 > -l && n1 < l) { CGcoeff.push_back((clebschGordan(l-1,m1,1,m2,l,m))*(clebschGordan(l-1,n1,1,n2,l,n))); } } } } } } } template <class Real> void WignerMatrix<Real>::compute_D(int l, Real phi, Real theta, Real psi) { double s = std::sin(theta); double c = std::cos(theta); Complex i(0.0, 1.0); Complex eiphi = std::exp(i*phi); Complex emiphi = std::exp(-i*phi); Complex eipsi = std::exp(i*psi); Complex emipsi = std::exp(-i*psi); if (D.size() < (std::size_t)(l+1)) D.resize(l+1); D[1].reshape(MultiArrayShape<2>::type(3,3)); D[1](0,0) = emipsi * Complex(Real(0.5*(1.0+c))) * emiphi; D[1](0,1) = Complex(Real(-s/M_SQRT2)) * emiphi; D[1](0,2) = eipsi * Complex(Real(0.5*(1.0-c))) * emiphi; D[1](1,0) = emipsi * Complex(Real(s/M_SQRT2)); D[1](1,1) = Complex(Real(c)); D[1](1,2) = eipsi * Complex(Real(-s/M_SQRT2)); D[1](2,0) = emipsi * Complex(Real(0.5*(1.0-c))) * eiphi; D[1](2,1) = Complex(Real(s/M_SQRT2)) * eiphi; D[1](2,2) = eipsi * Complex(Real(0.5*(1.0+c))) * eiphi; l_max = 1; cg1cnt = 0; if(l > 1) compute_D( l); } template <class Real> void WignerMatrix<Real>::compute_D(int l) { if (D.size() < (std::size_t)(l+1) ) { D.resize(l+1); l_max = 0; } if (l==1) { //precompute D0 =1 and D1 = (90 degree rot) // FIXME: signs are inconsistent with above explicit formula for // theta = pi/2, phi=0, psi=0 (sine terms should be negated) D[1].reshape(MultiArrayShape<2>::type(3,3)); D[1](0,0) = Real(0.5); D[1](0,1) = Real(0.5*M_SQRT2); D[1](0,2) = Real(0.5); D[1](1,0) = Real(-0.5*M_SQRT2); D[1](1,1) = Real(0.0); D[1](1,2) = Real(0.5*M_SQRT2); D[1](2,0) = Real(0.5); D[1](2,1) = Real(-0.5*M_SQRT2); D[1](2,2) = Real(0.5); l_max = 1; cg1cnt = 0; } else { //compute D2-Dl_max recursive if (l>l_max+1) { compute_D(l-1); } D[l].reshape(MultiArrayShape<2>::type(2*l+1,2*l+1)); D[l].init(Real(0.0)); for(int m = -l; m <= l ; m++) { for(int n = -l; n <= l ; n++) { for (int m2 = -1; m2 <= 1; m2++) { for (int n2 = -1; n2 <= 1; n2++) { int m1 = m-m2; int n1 = n-n2; if ((m1 > -l) && (m1 < l) && (n1 > -l) && (n1 < l)) { D[l](m+l,n+l) += D[1](m2+1,n2+1) * D[l-1](m1+l-1,n1+l-1) * Real(CGcoeff[cg1cnt++]); } } } } } l_max = l; } } template <class Real> void WignerMatrix<Real>::rotatePH(NestedArray const & PH, Real phi, Real theta, Real psi, NestedArray & PHresult) { int band = PH[1].size()-1; compute_D(band, phi, theta, psi); PHresult.resize(PH.size()); for(int n=1; n<=band; n++) { PHresult[n].resize(band+1); for (int l=0; l<=band; l++) { PHresult[n][l].resize(2*band+1); for(int m=-l; m<=l; m++) { Complex tmp = 0; for (int h=-l; h<=l; h++) { tmp += get_D(l,h,m) * PH[n][l][h+l]; } PHresult[n][l][m+l] = tmp; } } } } } // namespace vigra #endif // VIGRA_WIGNER_MATRIX_HXX
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  */
103  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  */
125  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>
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>
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
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
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
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
vigra::Quaternion
Definition: quaternion.hxx:61
vigra::conj
FFTWComplex< R > conj(const FFTWComplex< R > &a)
complex conjugate
Definition: fftw3.hxx:1030
vigra::WignerMatrix
computation of Wigner D matrix + rotation functions in SH,VH and R³
Definition: wigner-matrix.hxx:63
vigra::linalg::sin
linalg::TemporaryMatrix< T > sin(MultiArrayView< 2, T, C > const &v)
vigra::linalg::exp
linalg::TemporaryMatrix< T > exp(MultiArrayView< 2, T, C > const &v)
vigra::WignerMatrix::rotatePH
void rotatePH(NestedArray const &PH, Real phi, Real theta, Real psi, NestedArray &PHresult)
Rotate in PH.
Definition: wigner-matrix.hxx:291
vigra::ArrayVector
Definition: array_vector.hxx:58
vigra::MultiArrayShape
Definition: multi_fwd.hxx:63
vigra::WignerMatrix::compute_D
void compute_D(int band)
Compute D with fixed theta = pi/2, phi=0, psi=0.
Definition: wigner-matrix.hxx:229
vigra::TinyVector
Class for fixed size vectors.This class contains an array of size SIZE of the specified VALUETYPE...
Definition: accessor.hxx:940
vigra::WignerMatrix::get_D
Complex get_D(int l, int n, int m) const
Get the (n,m) entry of D.
Definition: wigner-matrix.hxx:103
vigra::ArrayVectorView::size
size_type size() const
Definition: array_vector.hxx:358
vigra::WignerMatrix::WignerMatrix
WignerMatrix(int l_max)
constructor
Definition: wigner-matrix.hxx:165
vigra::linalg::cos
linalg::TemporaryMatrix< T > cos(MultiArrayView< 2, T, C > const &v)
vigra::WignerMatrix::get_D
Matrix< Complex > const & get_D(int l) const
Return the rotation matrix D for the lth band.
Definition: wigner-matrix.hxx:125

© 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.11.1 (Fri May 19 2017)