/vol/vipdata/irtk/packages/transformation/include/irtkQuaternionTransformation.h

00001 /*=========================================================================
00002 
00003   Library   : Image Registration Toolkit (IRTK)
00004   Module    : $Id: irtkQuaternionTransformation.h 2 2008-12-23 12:40:14Z dr $
00005   Copyright : Imperial College, Department of Computing
00006               Visual Information Processing (VIP), 2008 onwards
00007   Date      : $Date: 2008-12-23 12:40:14 +0000 (Tue, 23 Dec 2008) $
00008   Version   : $Revision: 2 $
00009   Changes   : $Author: dr $
00010 
00011 =========================================================================*/
00012 
00013 #ifndef _IRTKQUATERNIONTRANSFORMATION_H
00014 
00015 #define _IRTKQUATERNIONTRANSFORMATION_H
00016 
00017 #include <irtkTransformation.h>
00018 
00019 class irtkQuaternionTransformation : public irtkTransformation
00020 {
00021 private:
00023   double _origin[3];
00024 
00026   double _translation[3];
00027 
00029   double _scale[3];
00030 
00032   double _alphaAxis[3];
00033 
00035   double _betaAxis[3];
00036 
00038   double _gammaAxis[3];
00039 
00041   double _finalAlphaAxis[3];
00042 
00044   double _finalBetaAxis[3];
00045 
00047   double _finalGammaAxis[3];
00048 
00050   double _alpha;
00051 
00053   double _beta;
00054 
00056   double _gamma;
00057 
00059   irtkMatrix _matrix;
00060 
00062   static const int NUMBER_OF_DOFS = 9;
00063 
00064 protected:
00066   void UpdateGammaAxis();
00067 
00072   void UpdateAxes();
00073 
00075   void Normalize(double* vec);
00076 
00077 public:
00079   irtkQuaternionTransformation();
00080 
00083   irtkQuaternionTransformation(double* origin, double* alphaAxis, double* betaAxis);
00084 
00086   irtkQuaternionTransformation(const irtkQuaternionTransformation& trans);
00087 
00089   irtkQuaternionTransformation& operator=(const irtkQuaternionTransformation& trans);
00090 
00092   void SetAlpha(double angle);
00093 
00095   void SetBeta(double angle);
00096 
00098   void SetGamma(double angle);
00099 
00101   double GetAlpha() const;
00102 
00104   double GetBeta() const;
00105 
00107   double GetGamma() const;
00108 
00110   void SetAxes(double* alphaAxis, double* betaAxis);
00111 
00113   void GetAxes(double* alphaAxis, double* betaAxis) const;
00114 
00116   void GetAxes(double* alphaAxis, double* betaAxis, double* gammaAxis) const;
00117 
00120   void GetRotatedAxes(double* alphaAxis, double* betaAxis) const;
00121 
00124   void GetRotatedAxes(double* alphaAxis, double* betaAxis, double* gammaAxis) const;
00125 
00127   void SetOrigin(double* origin);
00128 
00130   void GetOrigin(double* origin) const;
00131 
00133   void SetTranslation(double* translation);
00134 
00136   void GetTranslation(double* translation) const;
00137 
00139   void SetScale(double* scale);
00140 
00142   void GetScale(double* scale) const;
00143 
00147   void ResetOrigin();
00148 
00150   const irtkMatrix& GetMatrix() const;
00151 
00153   void UpdateMatrix();
00154 
00156   virtual int NumberOfDOFs() const;
00157 
00159   virtual double Get(int param) const;
00160 
00162   virtual void Put(int param, double val);
00163 
00165   virtual void Transform(double& x, double& y, double& z);
00166 
00168   virtual void GlobalTransform(double &, double &, double &);
00169 
00171   virtual void LocalTransform (double &, double &, double &);
00172 
00174   virtual void GlobalDisplacement(double &, double &, double &);
00175 
00177   virtual void LocalDisplacement(double &, double &, double &);
00178 
00180   virtual void Invert();
00181 
00183   virtual void Jacobian(double x, double y, double z, irtkMatrix& jac);
00184 
00186   virtual double Jacobian(double x, double y, double z);
00187 
00194   virtual void LocalJacobian(double x, double y, double z, irtkMatrix& jac);
00195 
00201   virtual void GlobalJacobian(double x, double y, double z, irtkMatrix& jac);
00202 
00205   static int CheckHeader(char* pFileName);
00206 
00208   virtual Bool IsIdentity();
00209 
00211   virtual irtkCifstream& Read(irtkCifstream& is);
00212 
00214   virtual istream& Import(istream& is);
00215 
00217   virtual irtkCofstream& Write(irtkCofstream& os);
00218 
00220   virtual ostream& Export(ostream& os);
00221 
00223   virtual void Print();
00224 
00226   virtual const char* NameOfClass();
00227 };
00228 
00229 inline void irtkQuaternionTransformation::UpdateGammaAxis()
00230 {
00231   _gammaAxis[0] = _alphaAxis[1]*_betaAxis[2] - _alphaAxis[2]*_betaAxis[1];
00232   _gammaAxis[1] = _alphaAxis[2]*_betaAxis[0] - _alphaAxis[0]*_betaAxis[2];
00233   _gammaAxis[2] = _alphaAxis[0]*_betaAxis[1] - _alphaAxis[1]*_betaAxis[0];
00234 }
00235 
00236 inline void irtkQuaternionTransformation::UpdateAxes()
00237 {
00238   Normalize(_alphaAxis);
00239 
00240   UpdateGammaAxis();
00241 
00242   if (_gammaAxis[0] == 0 && _gammaAxis[1] == 0 && _gammaAxis[2] == 0) {
00243     std::stringstream mesg;
00244     mesg << "Alpha axis and beta axis are collinear:" << std::endl
00245     << "_alphaAxis = (" << _alphaAxis[0] << ", " << _alphaAxis[1] << ", " << _alphaAxis[2] << ")" << std::endl
00246     << "_betaAxis = (" << _betaAxis[0] << ", " << _betaAxis[1] << ", " << _betaAxis[2] << ")";
00247   }
00248 
00249   double dp = _alphaAxis[0]*_betaAxis[0] + _alphaAxis[1]*_betaAxis[1] +
00250               _alphaAxis[2]*_betaAxis[2];
00251 
00252   _betaAxis[0] -= dp*_alphaAxis[0];
00253   _betaAxis[1] -= dp*_alphaAxis[1];
00254   _betaAxis[2] -= dp*_alphaAxis[2];
00255 
00256   Normalize(_betaAxis);
00257   Normalize(_gammaAxis);
00258 }
00259 
00260 inline void irtkQuaternionTransformation::Normalize(double* vec)
00261 {
00262   double len = sqrt(vec[0]*vec[0] + vec[1]*vec[1] + vec[2]*vec[2]);
00263 
00264   if (len == 0)
00265     throw irtkException("Tried to normalize vector of length 0.", __FILE__,
00266                         __LINE__);
00267 
00268   vec[0] /= len;
00269   vec[1] /= len;
00270   vec[2] /= len;
00271 }
00272 
00273 inline void irtkQuaternionTransformation::SetAlpha(double angle)
00274 {
00275   _alpha = angle*M_PI/180;
00276   UpdateMatrix();
00277 }
00278 
00279 inline void irtkQuaternionTransformation::SetBeta(double angle)
00280 {
00281   _beta = angle*M_PI/180;
00282   UpdateMatrix();
00283 }
00284 
00285 inline void irtkQuaternionTransformation::SetGamma(double angle)
00286 {
00287   _gamma = angle*M_PI/180;
00288   UpdateMatrix();
00289 }
00290 
00291 inline double irtkQuaternionTransformation::GetAlpha() const
00292 {
00293   return _alpha*180/M_PI;
00294 }
00295 
00296 inline double irtkQuaternionTransformation::GetBeta() const
00297 {
00298   return _beta*180/M_PI;
00299 }
00300 
00301 inline double irtkQuaternionTransformation::GetGamma() const
00302 {
00303   return _gamma*180/M_PI;
00304 }
00305 
00306 inline void irtkQuaternionTransformation::SetAxes(double* alphaAxis, double* betaAxis)
00307 {
00308   for (int i = 0; i < 3; i++) {
00309     _alphaAxis[i] = alphaAxis[i];
00310     _betaAxis[i] = betaAxis[i];
00311   }
00312 
00313   UpdateAxes();
00314   UpdateMatrix();
00315 }
00316 
00317 inline void irtkQuaternionTransformation::GetAxes(double* alphaAxis, double* betaAxis) const
00318 {
00319   for (int i = 0; i < 3; i++) {
00320     alphaAxis[i] = _alphaAxis[i];
00321     betaAxis[i] = _betaAxis[i];
00322   }
00323 }
00324 
00325 inline void irtkQuaternionTransformation::GetAxes(double* alphaAxis, double* betaAxis, double* gammaAxis) const
00326 {
00327   for (int i = 0; i < 3; i++) {
00328     alphaAxis[i] = _alphaAxis[i];
00329     betaAxis[i] = _betaAxis[i];
00330     gammaAxis[i] = _gammaAxis[i];
00331   }
00332 }
00333 
00334 inline void irtkQuaternionTransformation::GetRotatedAxes(double* alphaAxis, double* betaAxis) const
00335 {
00336   for (int i = 0; i < 3; i++) {
00337     alphaAxis[i] = _finalAlphaAxis[i];
00338     betaAxis[i] = _finalBetaAxis[i];
00339   }
00340 }
00341 
00342 inline void irtkQuaternionTransformation::GetRotatedAxes(double* alphaAxis, double* betaAxis, double* gammaAxis) const
00343 {
00344   for (int i = 0; i < 3; i++) {
00345     alphaAxis[i] = _finalAlphaAxis[i];
00346     betaAxis[i] = _finalBetaAxis[i];
00347     gammaAxis[i] = _finalGammaAxis[i];
00348   }
00349 }
00350 
00351 inline void irtkQuaternionTransformation::SetOrigin(double* origin)
00352 {
00353   _origin[0] = origin[0];
00354   _origin[1] = origin[1];
00355   _origin[2] = origin[2];
00356 
00357   UpdateMatrix();
00358 }
00359 
00360 inline void irtkQuaternionTransformation::GetOrigin(double* origin) const
00361 {
00362   for (int i = 0; i < 3; i++)
00363     origin[i] = _origin[i];
00364 }
00365 
00366 inline void irtkQuaternionTransformation::SetTranslation(double* translation)
00367 {
00368   _translation[0] = translation[0];
00369   _translation[1] = translation[1];
00370   _translation[2] = translation[2];
00371 
00372   UpdateMatrix();
00373 }
00374 
00375 inline void irtkQuaternionTransformation::GetTranslation(double* translation) const
00376 {
00377   translation[0] = _translation[0];
00378   translation[1] = _translation[1];
00379   translation[2] = _translation[2];
00380 }
00381 
00382 inline void irtkQuaternionTransformation::SetScale(double* scale)
00383 {
00384   _scale[0] = scale[0];
00385   _scale[1] = scale[1];
00386   _scale[2] = scale[2];
00387 }
00388 
00389 inline void irtkQuaternionTransformation::GetScale(double* scale) const
00390 {
00391   scale[0] = _scale[0];
00392   scale[1] = _scale[1];
00393   scale[2] = _scale[2];
00394 }
00395 
00396 inline const irtkMatrix& irtkQuaternionTransformation::GetMatrix() const
00397 {
00398   return _matrix;
00399 }
00400 
00401 inline int irtkQuaternionTransformation::NumberOfDOFs() const
00402 {
00403   return NUMBER_OF_DOFS;
00404 }
00405 
00406 inline void irtkQuaternionTransformation::GlobalDisplacement(double &x, double &y, double &z)
00407 {
00408   cerr << "irtkQuaternionTransformation::GlobalDisplacement: No implemented yet" << endl;
00409 }
00410 
00411 inline void irtkQuaternionTransformation::LocalDisplacement(double &x, double &y, double &z)
00412 {
00413   cerr << "irtkQuaternionTransformation::LocalDisplacement: No implemented yet" << endl;
00414 }
00415 
00416 inline void irtkQuaternionTransformation::GlobalTransform(double &x, double &y, double &z)
00417 {
00418   cerr << "irtkQuaternionTransformation::GlobalTransform: No implemented yet" << endl;
00419 }
00420 
00421 inline void irtkQuaternionTransformation::LocalTransform(double &x, double &y, double &z)
00422 {
00423   cerr << "irtkQuaternionTransformation::LocalTransform: No implemented yet" << endl;
00424 }
00425 
00426 inline double irtkQuaternionTransformation::Get(int param) const
00427 {
00428   std::stringstream mesg;
00429 
00430   switch (param) {
00431   case 0:
00432     return _alpha*180/M_PI;
00433     break;
00434 
00435   case 1:
00436     return _beta*180/M_PI;
00437     break;
00438 
00439   case 2:
00440     return _gamma*180/M_PI;
00441     break;
00442 
00443   case 3:
00444     return _translation[0];
00445     break;
00446 
00447   case 4:
00448     return _translation[1];
00449     break;
00450 
00451   case 5:
00452     return _translation[2];
00453     break;
00454 
00455   case 6:
00456     return _scale[0];
00457     break;
00458 
00459   case 7:
00460     return _scale[1];
00461     break;
00462 
00463   case 8:
00464     return _scale[2];
00465     break;
00466 
00467   default:
00468     mesg << "Invalid parameter: " << param;
00469     throw irtkException(mesg.str(), __FILE__, __LINE__);
00470     break;
00471   }
00472 }
00473 
00474 inline void irtkQuaternionTransformation::Put(int param, double val)
00475 {
00476   std::stringstream mesg;
00477 
00478   switch (param) {
00479   case 0:
00480     _alpha = val*M_PI/180;
00481     break;
00482 
00483   case 1:
00484     _beta = val*M_PI/180;
00485     break;
00486 
00487   case 2:
00488     _gamma = val*M_PI/180;
00489     break;
00490 
00491   case 3:
00492     _translation[0] = val;
00493     break;
00494 
00495   case 4:
00496     _translation[1] = val;
00497     break;
00498 
00499   case 5:
00500     _translation[2] = val;
00501     break;
00502 
00503   case 6:
00504     _scale[0] = val;
00505     break;
00506 
00507   case 7:
00508     _scale[1] = val;
00509     break;
00510 
00511   case 8:
00512     _scale[2] = val;
00513     break;
00514 
00515   default:
00516     mesg << "Invalid parameter: " << param;
00517     throw irtkException(mesg.str(), __FILE__, __LINE__);
00518     break;
00519   }
00520 
00521   UpdateMatrix();
00522 }
00523 
00524 inline const char* irtkQuaternionTransformation::NameOfClass()
00525 {
00526   return "irtkQuaternionTransformation";
00527 }
00528 
00529 #endif // __IRTKQUATERNIONTRANSFORMATION_H