include/Geant4/G4AffineTransform.icc

0001 //
0002 // ********************************************************************
0003 // * License and Disclaimer                                           *
0004 // *                                                                  *
0005 // * The  Geant4 software  is  copyright of the Copyright Holders  of *
0006 // * the Geant4 Collaboration.  It is provided  under  the terms  and *
0007 // * conditions of the Geant4 Software License,  included in the file *
0008 // * LICENSE and available at  http://cern.ch/geant4/license .  These *
0009 // * include a list of copyright holders.                             *
0010 // *                                                                  *
0011 // * Neither the authors of this software system, nor their employing *
0012 // * institutes,nor the agencies providing financial support for this *
0013 // * work  make  any representation or  warranty, express or implied, *
0014 // * regarding  this  software system or assume any liability for its *
0015 // * use.  Please see the license in the file  LICENSE  and URL above *
0016 // * for the full disclaimer and the limitation of liability.         *
0017 // *                                                                  *
0018 // * This  code  implementation is the result of  the  scientific and *
0019 // * technical work of the GEANT4 collaboration.                      *
0020 // * By using,  copying,  modifying or  distributing the software (or *
0021 // * any work based  on the software)  you  agree  to acknowledge its *
0022 // * use  in  resulting  scientific  publications,  and indicate your *
0023 // * acceptance of all terms of the Geant4 Software license.          *
0024 // ********************************************************************
0025 //
0026 // G4AffineTransformation Inline implementation
0027 //
0028 // Author: Paul R C Kent (CERN), 06.08.1996 - Initial version
0029 //         E.Tcherniaev (CERN), 19.09.1996, 06.05.2018 - Revised
0030 // --------------------------------------------------------------------
0031
0032 inline G4AffineTransform::G4AffineTransform()
0033  : rxx(1),rxy(0),rxz(0),
0034    ryx(0),ryy(1),ryz(0),
0035    rzx(0),rzy(0),rzz(1),
0036    tx(0),ty(0),tz(0)
0037 {}
0038
0039 inline G4AffineTransform::G4AffineTransform(const G4ThreeVector& tlate)
0040  : rxx(1),rxy(0),rxz(0),
0041    ryx(0),ryy(1),ryz(0),
0042    rzx(0),rzy(0),rzz(1),
0043    tx(tlate.x()),ty(tlate.y()),tz(tlate.z())
0044 {}
0045
0046 inline G4AffineTransform::G4AffineTransform(const G4RotationMatrix& rot)
0047  : rxx(rot.xx()),rxy(rot.xy()),rxz(rot.xz()),
0048    ryx(rot.yx()),ryy(rot.yy()),ryz(rot.yz()),
0049    rzx(rot.zx()),rzy(rot.zy()),rzz(rot.zz()),
0050    tx(0),ty(0),tz(0)
0051 {}
0052
0053 inline G4AffineTransform::G4AffineTransform( const G4RotationMatrix& rot,
0054                                              const G4ThreeVector& tlate )
0055  : rxx(rot.xx()),rxy(rot.xy()),rxz(rot.xz()),
0056    ryx(rot.yx()),ryy(rot.yy()),ryz(rot.yz()),
0057    rzx(rot.zx()),rzy(rot.zy()),rzz(rot.zz()),
0058    tx(tlate.x()),ty(tlate.y()),tz(tlate.z())
0059 {}
0060
0061 inline G4AffineTransform::G4AffineTransform( const G4RotationMatrix* rot,
0062                                              const G4ThreeVector& tlate )
0063  : tx(tlate.x()),ty(tlate.y()),tz(tlate.z())
0064 {
0065   if (rot != nullptr)
0066   {
0067     rxx=rot->xx();rxy=rot->xy();rxz=rot->xz();
0068     ryx=rot->yx();ryy=rot->yy();ryz=rot->yz();
0069     rzx=rot->zx();rzy=rot->zy();rzz=rot->zz();
0070   }
0071   else
0072   {
0073     rxx=1; rxy=0; rxz=0;
0074     ryx=0; ryy=1; ryz=0;
0075     rzx=0; rzy=0; rzz=1;
0076   }
0077 }
0078
0079 inline
0080 G4AffineTransform::
0081 G4AffineTransform( const G4double prxx,const G4double prxy,const G4double prxz,
0082                    const G4double pryx,const G4double pryy,const G4double pryz,
0083                    const G4double przx,const G4double przy,const G4double przz,
0084                    const G4double ptx,const G4double pty,const G4double ptz )
0085  : rxx(prxx),rxy(prxy),rxz(prxz),
0086    ryx(pryx),ryy(pryy),ryz(pryz),
0087    rzx(przx),rzy(przy),rzz(przz),
0088    tx(ptx),ty(pty),tz(ptz)
0089 {}
0090
0091 inline G4AffineTransform&
0092 G4AffineTransform::operator = (const G4AffineTransform& rhs)
0093 {
0094   if (this == &rhs)  { return *this; }
0095   rxx = rhs.rxx; rxy = rhs.rxy; rxz = rhs.rxz;
0096   ryx = rhs.ryx; ryy = rhs.ryy; ryz = rhs.ryz;
0097   rzx = rhs.rzx; rzy = rhs.rzy; rzz = rhs.rzz;
0098   tx = rhs.tx; ty = rhs.ty; tz = rhs.tz;
0099   return *this;
0100 }
0101
0102 inline G4AffineTransform
0103 G4AffineTransform::operator * (const G4AffineTransform& tf) const
0104 {
0105   return G4AffineTransform(rxx*tf.rxx+rxy*tf.ryx+rxz*tf.rzx,
0106                            rxx*tf.rxy+rxy*tf.ryy+rxz*tf.rzy,
0107                            rxx*tf.rxz+rxy*tf.ryz+rxz*tf.rzz,
0108
0109                            ryx*tf.rxx+ryy*tf.ryx+ryz*tf.rzx,
0110                            ryx*tf.rxy+ryy*tf.ryy+ryz*tf.rzy,
0111                            ryx*tf.rxz+ryy*tf.ryz+ryz*tf.rzz,
0112
0113                            rzx*tf.rxx+rzy*tf.ryx+rzz*tf.rzx,
0114                            rzx*tf.rxy+rzy*tf.ryy+rzz*tf.rzy,
0115                            rzx*tf.rxz+rzy*tf.ryz+rzz*tf.rzz,
0116
0117                            tx*tf.rxx+ty*tf.ryx+tz*tf.rzx+tf.tx,
0118                            tx*tf.rxy+ty*tf.ryy+tz*tf.rzy+tf.ty,
0119                            tx*tf.rxz+ty*tf.ryz+tz*tf.rzz+tf.tz);
0120 }
0121
0122 inline G4AffineTransform&
0123 G4AffineTransform::operator *= (const G4AffineTransform& tf)
0124 {
0125   // Use temporaries for `in place' compound transform computation
0126   //
0127   G4double nrxx=rxx*tf.rxx+rxy*tf.ryx+rxz*tf.rzx;
0128   G4double nrxy=rxx*tf.rxy+rxy*tf.ryy+rxz*tf.rzy;
0129   G4double nrxz=rxx*tf.rxz+rxy*tf.ryz+rxz*tf.rzz;
0130
0131   G4double nryx=ryx*tf.rxx+ryy*tf.ryx+ryz*tf.rzx;
0132   G4double nryy=ryx*tf.rxy+ryy*tf.ryy+ryz*tf.rzy;
0133   G4double nryz=ryx*tf.rxz+ryy*tf.ryz+ryz*tf.rzz;
0134
0135   G4double nrzx=rzx*tf.rxx+rzy*tf.ryx+rzz*tf.rzx;
0136   G4double nrzy=rzx*tf.rxy+rzy*tf.ryy+rzz*tf.rzy;
0137   G4double nrzz=rzx*tf.rxz+rzy*tf.ryz+rzz*tf.rzz;
0138
0139   G4double ntx=tx*tf.rxx+ty*tf.ryx+tz*tf.rzx+tf.tx;
0140   G4double nty=tx*tf.rxy+ty*tf.ryy+tz*tf.rzy+tf.ty;
0141   G4double ntz=tx*tf.rxz+ty*tf.ryz+tz*tf.rzz+tf.tz;
0142
0143   tx=ntx; ty=nty; tz=ntz;
0144   rxx=nrxx; rxy=nrxy; rxz=nrxz;
0145   ryx=nryx; ryy=nryy; ryz=nryz;
0146   rzx=nrzx; rzy=nrzy; rzz=nrzz;
0147
0148   return *this;
0149 }
0150
0151 inline G4AffineTransform&
0152 G4AffineTransform::Product( const G4AffineTransform& tf1,
0153                             const G4AffineTransform& tf2 )
0154 {
0155   rxx = tf1.rxx*tf2.rxx + tf1.rxy*tf2.ryx + tf1.rxz*tf2.rzx;
0156   rxy = tf1.rxx*tf2.rxy + tf1.rxy*tf2.ryy + tf1.rxz*tf2.rzy;
0157   rxz = tf1.rxx*tf2.rxz + tf1.rxy*tf2.ryz + tf1.rxz*tf2.rzz;
0158
0159   ryx = tf1.ryx*tf2.rxx + tf1.ryy*tf2.ryx + tf1.ryz*tf2.rzx;
0160   ryy = tf1.ryx*tf2.rxy + tf1.ryy*tf2.ryy + tf1.ryz*tf2.rzy;
0161   ryz = tf1.ryx*tf2.rxz + tf1.ryy*tf2.ryz + tf1.ryz*tf2.rzz;
0162
0163   rzx = tf1.rzx*tf2.rxx + tf1.rzy*tf2.ryx + tf1.rzz*tf2.rzx;
0164   rzy = tf1.rzx*tf2.rxy + tf1.rzy*tf2.ryy + tf1.rzz*tf2.rzy;
0165   rzz = tf1.rzx*tf2.rxz + tf1.rzy*tf2.ryz + tf1.rzz*tf2.rzz;
0166
0167   tx = tf1.tx*tf2.rxx + tf1.ty*tf2.ryx + tf1.tz*tf2.rzx + tf2.tx;
0168   ty = tf1.tx*tf2.rxy + tf1.ty*tf2.ryy + tf1.tz*tf2.rzy + tf2.ty;
0169   tz = tf1.tx*tf2.rxz + tf1.ty*tf2.ryz + tf1.tz*tf2.rzz + tf2.tz;
0170
0171   return *this;
0172 }
0173
0174 inline G4AffineTransform&
0175 G4AffineTransform::InverseProduct( const G4AffineTransform& tf1,
0176                                    const G4AffineTransform& tf2 )
0177 {
0178   if ((tf2.rxx + tf2.ryy + tf2.rzz) == 3.)
0179   {
0180     rxx = tf1.rxx;
0181     rxy = tf1.rxy;
0182     rxz = tf1.rxz;
0183
0184     ryx = tf1.ryx;
0185     ryy = tf1.ryy;
0186     ryz = tf1.ryz;
0187
0188     rzx = tf1.rzx;
0189     rzy = tf1.rzy;
0190     rzz = tf1.rzz;
0191
0192     tx = tf1.tx - tf2.tx;
0193     ty = tf1.ty - tf2.ty;
0194     tz = tf1.tz - tf2.tz;
0195   }
0196   else
0197   {
0198     G4double tf1rxx = tf1.rxx, tf1rxy = tf1.rxy, tf1rxz = tf1.rxz;
0199     rxx = tf1rxx*tf2.rxx + tf1rxy*tf2.rxy + tf1rxz*tf2.rxz;
0200     rxy = tf1rxx*tf2.ryx + tf1rxy*tf2.ryy + tf1rxz*tf2.ryz;
0201     rxz = tf1rxx*tf2.rzx + tf1rxy*tf2.rzy + tf1rxz*tf2.rzz;
0202
0203     G4double tf1ryx = tf1.ryx, tf1ryy = tf1.ryy, tf1ryz = tf1.ryz;
0204     ryx = tf1ryx*tf2.rxx + tf1ryy*tf2.rxy + tf1ryz*tf2.rxz;
0205     ryy = tf1ryx*tf2.ryx + tf1ryy*tf2.ryy + tf1ryz*tf2.ryz;
0206     ryz = tf1ryx*tf2.rzx + tf1ryy*tf2.rzy + tf1ryz*tf2.rzz;
0207
0208     G4double tf1rzx = tf1.rzx, tf1rzy = tf1.rzy, tf1rzz = tf1.rzz;
0209     rzx = tf1rzx*tf2.rxx + tf1rzy*tf2.rxy + tf1rzz*tf2.rxz;
0210     rzy = tf1rzx*tf2.ryx + tf1rzy*tf2.ryy + tf1rzz*tf2.ryz;
0211     rzz = tf1rzx*tf2.rzx + tf1rzy*tf2.rzy + tf1rzz*tf2.rzz;
0212
0213     G4double tf1_2tx = tf1.tx - tf2.tx;
0214     G4double tf1_2ty = tf1.ty - tf2.ty;
0215     G4double tf1_2tz = tf1.tz - tf2.tz;
0216     tx = tf1_2tx*tf2.rxx + tf1_2ty*tf2.rxy + tf1_2tz*tf2.rxz;
0217     ty = tf1_2tx*tf2.ryx + tf1_2ty*tf2.ryy + tf1_2tz*tf2.ryz;
0218     tz = tf1_2tx*tf2.rzx + tf1_2ty*tf2.rzy + tf1_2tz*tf2.rzz;
0219   }
0220   return *this;
0221 }
0222
0223 inline G4ThreeVector
0224 G4AffineTransform::TransformPoint(const G4ThreeVector& vec) const
0225 {
0226   G4double vecx = vec.x(), vecy = vec.y(), vecz = vec.z();
0227   return { vecx*rxx + vecy*ryx + vecz*rzx + tx,
0228            vecx*rxy + vecy*ryy + vecz*rzy + ty,
0229            vecx*rxz + vecy*ryz + vecz*rzz + tz };
0230 }
0231
0232 inline G4ThreeVector
0233 G4AffineTransform::InverseTransformPoint(const G4ThreeVector& vec) const
0234 {
0235   G4double vecx = vec.x()-tx, vecy = vec.y()-ty, vecz = vec.z()-tz;
0236   return { vecx*rxx + vecy*rxy + vecz*rxz,
0237            vecx*ryx + vecy*ryy + vecz*ryz,
0238            vecx*rzx + vecy*rzy + vecz*rzz };
0239 }
0240
0241 inline G4ThreeVector
0242 G4AffineTransform::TransformAxis(const G4ThreeVector& axis) const
0243 {
0244   G4double axisx = axis.x(), axisy = axis.y(), axisz = axis.z();
0245   return { axisx*rxx + axisy*ryx + axisz*rzx,
0246            axisx*rxy + axisy*ryy + axisz*rzy,
0247            axisx*rxz + axisy*ryz + axisz*rzz };
0248 }
0249
0250 inline G4ThreeVector
0251 G4AffineTransform::InverseTransformAxis(const G4ThreeVector& axis) const
0252 {
0253   G4double axisx = axis.x(), axisy = axis.y(), axisz = axis.z();
0254   return { axisx*rxx + axisy*rxy + axisz*rxz,
0255            axisx*ryx + axisy*ryy + axisz*ryz,
0256            axisx*rzx + axisy*rzy + axisz*rzz };
0257 }
0258
0259 inline
0260 void G4AffineTransform::ApplyPointTransform(G4ThreeVector& vec) const
0261 {
0262   G4double vecx = vec.x(), vecy = vec.y(), vecz = vec.z();
0263   vec.setX( vecx*rxx + vecy*ryx + vecz*rzx + tx );
0264   vec.setY( vecx*rxy + vecy*ryy + vecz*rzy + ty );
0265   vec.setZ( vecx*rxz + vecy*ryz + vecz*rzz + tz );
0266 }
0267
0268 inline
0269 void G4AffineTransform::ApplyAxisTransform(G4ThreeVector& axis) const
0270 {
0271   G4double axisx = axis.x(), axisy = axis.y(), axisz = axis.z();
0272   axis.setX( axisx*rxx + axisy*ryx + axisz*rzx );
0273   axis.setY( axisx*rxy + axisy*ryy + axisz*rzy );
0274   axis.setZ( axisx*rxz + axisy*ryz + axisz*rzz );
0275 }
0276
0277 inline
0278 G4AffineTransform G4AffineTransform::Inverse() const
0279 {
0280   G4double ttx = -tx, tty = -ty, ttz = -tz;
0281   return G4AffineTransform( rxx, ryx, rzx,
0282                             rxy, ryy, rzy,
0283                             rxz, ryz, rzz,
0284                             ttx*rxx + tty*rxy + ttz*rxz,
0285                             ttx*ryx + tty*ryy + ttz*ryz,
0286                             ttx*rzx + tty*rzy + ttz*rzz );
0287 }
0288
0289 inline
0290 G4AffineTransform& G4AffineTransform::Invert()
0291 {
0292   G4double ttx = -tx, tty = -ty, ttz = -tz;
0293            tx = ttx*rxx + tty*rxy + ttz*rxz;
0294            ty = ttx*ryx + tty*ryy + ttz*ryz;
0295            tz = ttx*rzx + tty*rzy + ttz*rzz;
0296
0297   G4double tmp1=ryx; ryx=rxy; rxy=tmp1;
0298   G4double tmp2=rzx; rzx=rxz; rxz=tmp2;
0299   G4double tmp3=rzy; rzy=ryz; ryz=tmp3;
0300
0301   return *this;
0302 }
0303
0304 inline
0305 G4AffineTransform& G4AffineTransform::operator +=(const G4ThreeVector& tlate)
0306 {
0307   tx += tlate.x();
0308   ty += tlate.y();
0309   tz += tlate.z();
0310
0311   return *this;
0312 }
0313
0314 inline
0315 G4AffineTransform& G4AffineTransform::operator -=(const G4ThreeVector& tlate)
0316 {
0317   tx -= tlate.x();
0318   ty -= tlate.y();
0319   tz -= tlate.z();
0320
0321   return *this;
0322 }
0323
0324 inline
0325 G4bool G4AffineTransform::operator == (const G4AffineTransform& tf) const
0326 {
0327   return tx==tf.tx&&ty==tf.ty&&tz==tf.tz&&
0328          rxx==tf.rxx&&rxy==tf.rxy&&rxz==tf.rxz&&
0329          ryx==tf.ryx&&ryy==tf.ryy&&ryz==tf.ryz&&
0330          rzx==tf.rzx&&rzy==tf.rzy&&rzz==tf.rzz;
0331 }
0332
0333 inline
0334 G4bool G4AffineTransform::operator != (const G4AffineTransform& tf) const
0335 {
0336   return tx!=tf.tx||ty!=tf.ty||tz!=tf.tz||
0337          rxx!=tf.rxx||rxy!=tf.rxy||rxz!=tf.rxz||
0338          ryx!=tf.ryx||ryy!=tf.ryy||ryz!=tf.ryz||
0339          rzx!=tf.rzx||rzy!=tf.rzy||rzz!=tf.rzz;
0340 }
0341
0342 inline
0343 G4double G4AffineTransform::operator [] (const G4int n) const
0344 {
0345   G4double v = 0.0;
0346   switch(n)
0347   {
0348     case 0:
0349             v=rxx;
0350             break;
0351     case 1:
0352             v=rxy;
0353             break;
0354     case 2:
0355             v=rxz;
0356             break;
0357     case 4:
0358             v=ryx;
0359             break;
0360     case 5:
0361             v=ryy;
0362             break;
0363     case 6:
0364             v=ryz;
0365             break;
0366     case 8:
0367             v=rzx;
0368             break;
0369     case 9:
0370             v=rzy;
0371             break;
0372     case 10:
0373             v=rzz;
0374             break;
0375     case 12:
0376             v=tx;
0377             break;
0378     case 13:
0379             v=ty;
0380             break;
0381     case 14:
0382             v=tz;
0383             break;
0384     case 3:
0385     case 7:
0386     case 11:
0387             break;
0388     case 15:
0389             v=1.0;
0390             break;
0391   }
0392   return v;
0393 }
0394
0395 inline
0396 G4bool G4AffineTransform::IsRotated() const
0397 {
0398   return rxx!=1.0 || ryy!=1.0 || rzz!=1.0;
0399 }
0400
0401 inline
0402 G4bool G4AffineTransform::IsTranslated() const
0403 {
0404   return (tx != 0.0) || (ty != 0.0) || (tz != 0.0);
0405 }
0406
0407 inline G4RotationMatrix G4AffineTransform::NetRotation() const
0408 {
0409   return G4Rep3x3(rxx,rxy,rxz,
0410                   ryx,ryy,ryz,
0411                   rzx,rzy,rzz);
0412 }
0413
0414 inline G4RotationMatrix G4AffineTransform::InverseNetRotation() const
0415 {
0416   return G4Rep3x3(rxx,ryx,rzx,
0417                   rxy,ryy,rzy,
0418                   rxz,ryz,rzz);
0419 }
0420
0421 inline
0422 G4ThreeVector G4AffineTransform::NetTranslation() const
0423 {
0424   return {tx,ty,tz};
0425 }
0426
0427 inline
0428 G4ThreeVector G4AffineTransform::InverseNetTranslation() const
0429 {
0430   G4double ttx = -tx, tty = -ty, ttz = -tz;
0431   G4double invtx = ttx*rxx + tty*rxy + ttz*rxz;
0432   G4double invty = ttx*ryx + tty*ryy + ttz*ryz;
0433   G4double invtz = ttx*rzx + tty*rzy + ttz*rzz;
0434   return {invtx,invty,invtz};
0435 }
0436
0437 inline
0438 void G4AffineTransform::SetNetRotation(const G4RotationMatrix& rot)
0439 {
0440   rxx=rot.xx();
0441   rxy=rot.xy();
0442   rxz=rot.xz();
0443   ryx=rot.yx();
0444   ryy=rot.yy();
0445   ryz=rot.yz();
0446   rzx=rot.zx();
0447   rzy=rot.zy();
0448   rzz=rot.zz();
0449 }
0450
0451 inline
0452 void G4AffineTransform::SetNetTranslation(const G4ThreeVector& tlate)
0453 {
0454   tx=tlate.x();
0455   ty=tlate.y();
0456   tz=tlate.z();
0457 }
0458
0459 inline
0460 G4AffineTransform::operator G4Transform3D () const
0461 {
0462   return G4Transform3D(NetRotation().inverse(),NetTranslation());
0463 }
0464
0465 inline
0466 std::ostream& operator << (std::ostream& os, const G4AffineTransform& transf)
0467 {
0468   std::streamsize oldPrec = os.precision(6);
0469   G4double DeviationTolerance = 1.0e-05;
0470
0471   G4double diagDeviation = 0.0;
0472   G4double offdDeviationUL = 0.0;
0473   G4double offdDeviationDR = 0.0;
0474   G4double offdDeviation = 0.0;
0475
0476   os << "  Transformation: " << G4endl;
0477
0478   // double  a = std::max ( 1, 2, 3 ) ;
0479
0480   G4bool UnitTr = ! transf.IsRotated();
0481   diagDeviation = std::max( std::fabs( transf[0] - 1.0 ) ,   //  |rxx - 1|
0482                             std::fabs( transf[5] - 1.0 ) );  //  |ryy - 1|
0483   diagDeviation = std::max( diagDeviation,
0484                             std::fabs( transf[10] - 1.0 ) ); //  |rzz - 1|
0485
0486   offdDeviationUL = std::max( std::fabs( transf[1] ) ,       //  |rxy|
0487                               std::fabs( transf[2] ) );      //  |rxz|
0488   offdDeviationUL = std::max( offdDeviationUL,
0489                               std::fabs( transf[4] ) );      //  |ryx|
0490
0491   offdDeviationDR = std::max( std::fabs( transf[6] ) ,       //  |ryz|
0492                               std::fabs( transf[8] ) );      //  |rzx|
0493   offdDeviationDR = std::max( offdDeviationDR,
0494                               std::fabs( transf[9] ) );      //  |rzy|
0495   offdDeviation = std::max( offdDeviationUL, offdDeviationDR );
0496
0497   if( UnitTr || std::max(diagDeviation, offdDeviation) < DeviationTolerance )
0498   {
0499      os << "    UNIT  Rotation " << G4endl;
0500   }
0501   else
0502   {
0503      os << "rx/x,y,z: "
0504         << transf[0]  << " " << transf[1] << " " << transf[2] << G4endl
0505         << "ry/x,y,z: "
0506         << transf[4]  << " " << transf[5] << " " << transf[6] << G4endl
0507         << "rz/x,y,z: "
0508         << transf[8]  << " " << transf[9] << " " << transf[10] << G4endl;
0509   }
0510
0511   os << "tr/x,y,z: " << transf[12]  << " " << transf[13] << " " << transf[14]
0512      << G4endl;
0513
0514   os.precision(oldPrec);
0515
0516   return os;
0517 }