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