eic-opticks/u4/U4Transform.h

0001 #pragma once
0002
0003 #include <array>
0004 #include <csignal>
0005
0006 #include "G4VPhysicalVolume.hh"
0007 #include "G4Transform3D.hh"
0008 #include "G4RotationMatrix.hh"
0009 #include "G4BooleanSolid.hh"
0010 #include "G4MultiUnion.hh"
0011
0012 #include <glm/glm.hpp>
0013 #include <glm/gtc/type_ptr.hpp>
0014 #include "glm/gtx/string_cast.hpp"
0015
0016 struct U4Transform
0017 {
0018     static void Read(     const glm::tmat4x4<double>& dst, const double* src );
0019     static void ReadWide( const glm::tmat4x4<double>& dst, const float*  src );
0020
0021     static void WriteTransform(      double* dst, const glm::tmat4x4<double>& src );
0022     static void WriteObjectTransform(double* dst, const G4VPhysicalVolume* const pv) ;  // MOST USED BY OPTICKS
0023     static void WriteFrameTransform( double* dst, const G4VPhysicalVolume* const pv) ;
0024
0025     static void GetObjectTransform(glm::tmat4x4<double>& tr, const G4VPhysicalVolume* const pv) ; // MOST USED BY OPTICKS
0026     static void GetFrameTransform( glm::tmat4x4<double>& tr, const G4VPhysicalVolume* const pv) ;
0027     static void GetDispTransform(  glm::tmat4x4<double>& tr, const G4DisplacedSolid* disp );
0028     static void GetMultiUnionItemTransform(  glm::tmat4x4<double>& tr, const G4MultiUnion* const muni, unsigned item );
0029
0030     static void GetScaleTransform( glm::tmat4x4<double>& tr, double sx, double sy, double sz );
0031
0032     template<typename T>
0033     static void Convert(glm::tmat4x4<T>& dst,  const std::array<T,16>& src );
0034
0035     template<typename T>
0036     static void Convert_RotateThenTranslate(glm::tmat4x4<T>& dst,  const G4Transform3D& src );
0037
0038     template<typename T>
0039     static void Convert_RotateThenTranslate(glm::tmat4x4<T>& d, const G4RotationMatrix& rot, const G4ThreeVector& tla, bool f );
0040
0041     template<typename T>
0042     static void Convert_TranslateThenRotate(glm::tmat4x4<T>& dst,  const G4Transform3D& src );
0043
0044     template<typename T>
0045     static unsigned Check(const std::array<T,16>& a);
0046 };
0047
0048
0049 // HMM: no G4 types here, this should be elsewhere
0050 inline void U4Transform::WriteTransform( double* dst, const glm::tmat4x4<double>& src )
0051 {
0052     memcpy( dst, glm::value_ptr(src), sizeof(double)*16 );
0053 }
0054 inline void U4Transform::WriteObjectTransform(double* dst, const G4VPhysicalVolume* const pv)
0055 {
0056     glm::tmat4x4<double> tr(1.);
0057     GetObjectTransform(tr, pv );
0058     WriteTransform(dst, tr );
0059 }
0060 inline void U4Transform::WriteFrameTransform(double* dst, const G4VPhysicalVolume* const pv)
0061 {
0062     glm::tmat4x4<double> tr(1.);
0063     GetFrameTransform(tr, pv );
0064     WriteTransform(dst, tr );
0065 }
0066 inline void U4Transform::GetObjectTransform(glm::tmat4x4<double>& tr, const G4VPhysicalVolume* const pv)
0067 {
0068    // preferred for interop with glm/Opticks : obj relative to mother
0069     G4RotationMatrix rot = pv->GetObjectRotationValue() ;
0070     G4ThreeVector    tla = pv->GetObjectTranslation() ;
0071     G4Transform3D    tra(rot,tla);
0072     Convert_RotateThenTranslate(tr, tra);
0073 }
0074 inline void U4Transform::GetFrameTransform(glm::tmat4x4<double>& tr, const G4VPhysicalVolume* const pv)
0075 {
0076     const G4RotationMatrix* rotp = pv->GetFrameRotation() ;
0077     G4ThreeVector    tla = pv->GetFrameTranslation() ;
0078     G4Transform3D    tra(rotp ? *rotp : G4RotationMatrix(),tla);
0079     Convert_TranslateThenRotate(tr, tra );
0080 }
0081
0082 /**
0083 U4Transform::GetDispTransform
0084 ------------------------------
0085
0086 It looks a bit fishy using GetFrameRotation and GetObjectTranslation,
0087 but looking at the impl those are both from fDirectTransform.
0088
0089 g4-cls G4DisplacedSolid::
0090
0091     240 G4RotationMatrix G4DisplacedSolid::GetFrameRotation() const
0092     241 {
0093     242   G4RotationMatrix InvRotation= fDirectTransform->NetRotation();
0094     243   return InvRotation;
0095     244 }
0096
0097     281 G4ThreeVector  G4DisplacedSolid::GetObjectTranslation() const
0098     282 {
0099     283   return fDirectTransform->NetTranslation();
0100     284 }
0101
0102
0103
0104 **/
0105
0106 inline void U4Transform::GetDispTransform(glm::tmat4x4<double>& tr, const G4DisplacedSolid* disp )
0107 {
0108     assert(disp) ;
0109     G4RotationMatrix rot = disp->GetFrameRotation();
0110     G4ThreeVector    tla = disp->GetObjectTranslation();
0111
0112     //G4RotationMatrix rot = disp->GetObjectRotation();
0113     //G4ThreeVector    tla = disp->GetFrameTranslation();
0114
0115     Convert_RotateThenTranslate(tr, rot, tla, true );
0116 }
0117
0118
0119 /**
0120 U4Transform::GetMultiUnionItemTransform
0121 -----------------------------------------
0122
0123 Gave very incorrect bbox for WaterDistributorPartIIIUnion (a multiunion with
0124 rotations and translations on the items) whilst using Convert_TranslateThenRotate.
0125 After switched to Convert_RotateThenTranslate the bbox looks consistent with
0126 expectations from SMesh.
0127
0128 **/
0129
0130
0131 inline void U4Transform::GetMultiUnionItemTransform(  glm::tmat4x4<double>& xf, const G4MultiUnion* const muni, unsigned item )
0132 {
0133     [[maybe_unused]] unsigned num_item = muni->GetNumberOfSolids() ;
0134     assert( item < num_item );
0135     const G4Transform3D& tr = muni->GetTransformation(item) ;
0136
0137     Convert_RotateThenTranslate(xf,  tr );
0138 }
0139
0140
0141
0142 inline void U4Transform::GetScaleTransform( glm::tmat4x4<double>& tr, double sx, double sy, double sz )
0143 {
0144     glm::tvec3<double> sc(sx, sy, sz);
0145     tr = glm::scale(glm::tmat4x4<double>(1.), sc) ;
0146 }
0147
0148
0149
0150
0151
0152 template<typename T>
0153 inline void U4Transform::Convert(glm::tmat4x4<T>& d,  const std::array<T,16>& s ) // static
0154 {
0155     unsigned n = Check(s);
0156     bool n_expect = n == 0 ;
0157     assert( n_expect );
0158     if(!n_expect) std::raise(SIGINT);
0159
0160     memcpy( glm::value_ptr(d), s.data(), sizeof(T)*16 );
0161 }
0162
0163 /**
0164
0165 U4Transform::Convert_RotateThenTranslate
0166 -------------------------------------------
0167
0168 The canonical form of 4x4 transform with the translation visible
0169 in the last row assumes that the rotation is done first followed by
0170 the translation. This ordering is appropriate for model2world "m2w" transforms
0171 from "GetObjectTransform".
0172
0173 BUT that order is not appropriate for world2model : in that case need the converse.
0174 Translation first followed by rotation.
0175
0176 **/
0177
0178 template<typename T>
0179 inline void U4Transform::Convert_RotateThenTranslate(glm::tmat4x4<T>& d,  const G4Transform3D& s ) // static
0180 {
0181     T zero(0.);
0182     T one(1.);
0183
0184     std::array<T, 16> a = {{
0185              s.xx(), s.yx(), s.zx(), zero ,
0186              s.xy(), s.yy(), s.zy(), zero ,
0187              s.xz(), s.yz(), s.zz(), zero ,
0188              s.dx(), s.dy(), s.dz(), one   }} ;
0189     Convert(d, a);
0190 }
0191
0192
0193 /**
0194 U4Transform::Convert_RotateThenTranslate
0195 ------------------------------------------
0196
0197 Caution getting the correct transpose is always problematic...
0198
0199 Cannot rely on G4 streamers presenting in a way that matches
0200 Opticks/glm presentation of matrices.
0201
0202 G4AffineTransform::NetRotation
0203
0204 **/
0205
0206 template<typename T>
0207 inline void U4Transform::Convert_RotateThenTranslate(glm::tmat4x4<T>& d, const G4RotationMatrix& r, const G4ThreeVector& t, bool f ) // static
0208 {
0209     T zero(0.);
0210     T one(1.);
0211
0212     if( f == false )
0213     {
0214         std::array<T, 16> a = {{
0215                  r.xx(), r.yx(), r.zx(), zero ,
0216                  r.xy(), r.yy(), r.zy(), zero ,
0217                  r.xz(), r.yz(), r.zz(), zero ,
0218                  t.x(),  t.y(),  t.z(),  one   }} ;
0219
0220         Convert(d, a);
0221     }
0222     else
0223     {
0224         std::array<T, 16> a = {{
0225                  r.xx(), r.xy(), r.xz(), zero ,
0226                  r.yx(), r.yy(), r.yz(), zero ,
0227                  r.zx(), r.zy(), r.zz(), zero ,
0228                  t.x(),  t.y(),  t.z(),  one   }} ;
0229
0230         Convert(d, a);
0231     }
0232
0233 }
0234
0235
0236
0237
0238
0239
0240 /**
0241 U4Transform::Convert_TranslateThenRotate
0242 -------------------------------------------
0243
0244 See ana/translate_rotate.py for sympy demo::
0245
0246     rxx⋅tx + rxy⋅ty + rxz⋅tz
0247
0248     ryx⋅tx + ryy⋅ty + ryz⋅tz
0249
0250     rzx⋅tx + rzy⋅ty + rzz⋅tz
0251
0252 **/
0253
0254 template<typename T>
0255 inline void U4Transform::Convert_TranslateThenRotate(glm::tmat4x4<T>& d,  const G4Transform3D& s ) // static
0256 {
0257     T rxx = s.xx() ;
0258     T rxy = s.xy() ;
0259     T rxz = s.xz() ;
0260
0261     T ryx = s.yx() ;
0262     T ryy = s.yy() ;
0263     T ryz = s.yz() ;
0264
0265     T rzx = s.zx() ;
0266     T rzy = s.zy() ;
0267     T rzz = s.zz() ;
0268
0269     T tx = s.dx() ;
0270     T ty = s.dy() ;
0271     T tz = s.dz() ;
0272
0273     T RTx =  rxx*tx + rxy*ty + rxz*tz  ;
0274     T RTy =  ryx*tx + ryy*ty + ryz*tz  ;
0275     T RTz =  rzx*tx + rzy*ty + rzz*tz  ;
0276
0277     T zero(0.);
0278     T one(1.);
0279     std::array<T, 16> a = {{
0280              rxx   , ryx  , rzx  , zero ,
0281              rxy   , ryy  , rzy  , zero ,
0282              rxz   , ryz  , rzz  , zero ,
0283              RTx   , RTy  , RTz  , one   }} ;
0284
0285     unsigned n = Check(a);
0286     bool n_expect = n == 0 ;
0287     assert( n_expect );
0288     if(!n_expect) std::raise(SIGINT);
0289     memcpy( glm::value_ptr(d), a.data(), sizeof(T)*16 );
0290 }
0291
0292
0293
0294
0295
0296
0297 template<typename T>
0298 inline unsigned U4Transform::Check(const std::array<T,16>& a) // static
0299 {
0300     unsigned num_nan(0);
0301     unsigned num_inf(0);
0302     for(unsigned i=0 ; i < 16 ; i++) if(std::isnan(a[i])) num_nan++ ;
0303     for(unsigned i=0 ; i < 16 ; i++) if(std::isinf(a[i])) num_inf++ ;
0304     return num_nan + num_inf ;
0305 }
0306