eic-opticks/sysrap/stran.h

0001 #pragma once
0002 /**
0003 stran.h : Tran
0004 ================
0005
0006 * FOR CLARITY AM MOVING SINGLE TRANFORM FUNCTIONALITY TO stra.h
0007
0008 Transform handler that creates inverse transforms at every stage
0009 loosely based on NPY nmat4triple_
0010
0011 Aim to avoid having to take the inverse eg with glm::inverse
0012 which inevitably risks numerical issues.
0013
0014 **/
0015
0016 #include <iostream>
0017 #include <iomanip>
0018 #include <string>
0019 #include <sstream>
0020 #include <array>
0021 #include <vector>
0022 #include <csignal>
0023 #include <cstdlib>
0024
0025 #include "scuda.h"
0026 #include "sqat4.h"
0027 #include "stra.h"
0028 #include "NP.hh"
0029
0030 #include "glm/glm.hpp"
0031 #include "glm/gtx/string_cast.hpp"
0032 #include <glm/gtx/transform.hpp>
0033 #include <glm/gtc/type_ptr.hpp>
0034
0035
0036
0037 template<typename D, typename S>
0038 void TranConvert(glm::tmat4x4<D>& dst, const glm::tmat4x4<S>& src )  // static
0039 {
0040     const S* src_ptr = glm::value_ptr(src);
0041     D* dst_ptr = glm::value_ptr(dst);
0042     for(int i=0 ; i < 16 ; i++) dst_ptr[i] = D(src_ptr[i]) ;
0043 }
0044
0045
0046
0047 template<typename T>
0048 struct Tran
0049 {
0050     glm::tmat4x4<T> t ;  // transform
0051     glm::tmat4x4<T> v ;  // inverse
0052     glm::tmat4x4<T> i ;  // identity
0053
0054
0055     static constexpr const T EPSILON = 1e-6 ;
0056     static constexpr const bool VERBOSE = false ;
0057
0058     static const Tran<T>* make_translate( const T tx, const T ty, const T tz, const T sc);
0059     static const Tran<T>* make_translate( const T tx, const T ty, const T tz);
0060     static const Tran<T>* make_identity();
0061     static       Tran<T>* make_identity_();
0062     static const Tran<T>* make_scale(     const T sx, const T sy, const T sz);
0063     static const Tran<T>* make_rotate(    const T ax, const T ay, const T az, const T angle_deg);
0064
0065     static const Tran<T>* make_rotate_a2b( const glm::tvec3<T>& a, const glm::tvec3<T>& b, bool flip );
0066     static const Tran<T>* make_rotate_a2b( const T ax, const T ay, const T az, const T bx, const T by, const T bz, bool flip );
0067
0068     static const Tran<T>* product(const Tran<T>* a, const Tran<T>* b, bool reverse);
0069     static const Tran<T>* product(const Tran<T>* a, const Tran<T>* b, const Tran<T>* c, bool reverse);
0070     static const Tran<T>* product(const std::vector<const Tran<T>*>& tt, bool reverse );
0071
0072     static Tran<T>* ConvertToTran( const qat4* q, T epsilon=EPSILON );
0073     static Tran<T>* ConvertFromQat(const qat4* q_, T epsilon=EPSILON );
0074
0075     static Tran<T>* ConvertFromData(const T* data);
0076
0077     static const qat4* Invert( const qat4* q, T epsilon=EPSILON );
0078     static Tran<T>* FromPair( const qat4* t, const qat4* v, T epsilon=EPSILON ); // WIDENS from float
0079     static glm::tmat4x4<T> MatFromQat( const qat4* q );
0080     static glm::tmat4x4<T> MatFromData(const T* data );
0081
0082     static qat4*    ConvertFrom(const glm::tmat4x4<T>& tr );
0083
0084     // ctor
0085     Tran( const T* transform, const T* inverse ) ;
0086     Tran( const glm::tmat4x4<T>& transform, const glm::tmat4x4<T>& inverse ) ;
0087
0088     T    maxdiff_from_identity(char mat='t') const ;
0089     bool is_identity(char mat='t', T epsilon=1e-6) const ;
0090     std::string brief(bool only_tlate=false, char mat='t', unsigned wid=6, unsigned prec=1) const ;
0091     std::string desc() const ;
0092     bool checkIsIdentity(char mat='i', const char* caller="caller", T epsilon=EPSILON);
0093
0094     void write(T* dst, unsigned num_values=3*4*4) const ;
0095
0096
0097     NP* as_array()const ;
0098     void save_(const char* path) const ;
0099     void save(const char* dir, const char* name="stran.npy") const ;
0100
0101     static void Apply(        const glm::tmat4x4<T>& tr, T* p0    , T w    , unsigned count, unsigned stride, unsigned offset, bool normalize );
0102     static void ApplyToFloat( const glm::tmat4x4<T>& tr, float* p0, float w, unsigned count, unsigned stride, unsigned offset, bool normalize );
0103
0104     void photon_transform( NP* ph, bool normalize, bool inverse ) const ;
0105
0106     static NP* PhotonTransform( const NP* ph, bool normalize, const Tran<T>* tr, bool inverse );
0107     static void AddTransform( T* ttk, const char* opt, const glm::tvec3<T>& a, const glm::tvec3<T>& b, const glm::tvec3<T>& c );
0108
0109     const T* tdata() const ;
0110     const T* vdata() const ;
0111
0112 };
0113
0114
0115 /*
0116 template<typename T>
0117 inline std::ostream& operator<< (std::ostream& out, const glm::tmat4x4<T>& m  )
0118 {
0119     int prec = 4 ;
0120     int wid = 10 ;
0121     bool flip = false ;
0122     for(int i=0 ; i < 4 ; i++)
0123     {
0124         for(int j=0 ; j < 4 ; j++) out << std::setprecision(prec) << std::fixed << std::setw(wid) << ( flip ? m[j][i] : m[i][j] ) << " " ;
0125         out << std::endl ;
0126     }
0127     return out ;
0128 }
0129 */
0130
0131 template<typename T>
0132 inline std::ostream& operator<< (std::ostream& out, const Tran<T>& tr)
0133 {
0134     out
0135        << std::endl
0136        << "tr.t"
0137        << std::endl
0138        <<  tr.t
0139        << std::endl
0140        << "tr.v"
0141        << std::endl
0142        <<  tr.v
0143        << "tr.i"
0144        << std::endl
0145        <<  tr.i
0146        << std::endl
0147        ;
0148     return out;
0149 }
0150
0151
0152 template<typename T>
0153 inline const Tran<T>* Tran<T>::make_translate( const T tx, const T ty, const T tz, const T sc)
0154 {
0155     glm::tvec3<T> tlate(tx*sc,ty*sc,tz*sc);
0156     glm::tmat4x4<T> t = glm::translate(glm::tmat4x4<T>(1.),   tlate ) ;
0157     glm::tmat4x4<T> v = glm::translate(glm::tmat4x4<T>(1.),  -tlate ) ;
0158     return new Tran<T>(t, v);
0159 }
0160 template<typename T>
0161 inline const Tran<T>* Tran<T>::make_translate( const T tx, const T ty, const T tz)
0162 {
0163     glm::tvec3<T> tlate(tx,ty,tz);
0164     glm::tmat4x4<T> t = glm::translate(glm::tmat4x4<T>(1.),   tlate ) ;
0165     glm::tmat4x4<T> v = glm::translate(glm::tmat4x4<T>(1.),  -tlate ) ;
0166     return new Tran<T>(t, v);
0167 }
0168
0169 template<typename T>
0170 inline const Tran<T>* Tran<T>::make_identity()
0171 {
0172     glm::tmat4x4<T> t(1.) ;
0173     glm::tmat4x4<T> v(1.) ;
0174     return new Tran<T>(t, v);
0175 }
0176 template<typename T>
0177 inline Tran<T>* Tran<T>::make_identity_()
0178 {
0179     glm::tmat4x4<T> t(1.) ;
0180     glm::tmat4x4<T> v(1.) ;
0181     return new Tran<T>(t, v);
0182 }
0183
0184
0185 template<typename T>
0186 inline const Tran<T>* Tran<T>::make_scale( const T sx, const T sy, const T sz)
0187 {
0188     glm::tvec3<T> scal(sx,sy,sz);
0189     glm::tvec3<T> isca(1./sx,1./sy,1./sz);
0190     glm::tmat4x4<T> t = glm::scale(glm::tmat4x4<T>(1.),   scal ) ;
0191     glm::tmat4x4<T> v = glm::scale(glm::tmat4x4<T>(1.),   isca ) ;
0192     return new Tran<T>(t, v);
0193 }
0194
0195 template<typename T>
0196 inline const Tran<T>* Tran<T>::make_rotate( const T ax, const T ay, const T az, const T angle_deg)
0197 {
0198     T angle_rad = glm::pi<T>()*angle_deg/T(180.) ;
0199     glm::tvec3<T> axis(ax,ay,az);
0200     glm::tmat4x4<T> t = glm::rotate(glm::tmat4x4<T>(1.),  angle_rad, axis ) ;
0201     glm::tmat4x4<T> v = glm::rotate(glm::tmat4x4<T>(1.), -angle_rad, axis ) ;
0202     return new Tran<T>(t, v);
0203 }
0204
0205
0206 template<typename T>
0207 inline const Tran<T>* Tran<T>::make_rotate_a2b(const T ax, const T ay, const T az, const T bx, const T by, const T bz, bool flip )
0208 {
0209     glm::tvec3<T> a(ax,ay,az);
0210     glm::tvec3<T> b(bx,by,bz);
0211     return make_rotate_a2b(a,b, flip);
0212 }
0213
0214 template<typename T>
0215 inline const Tran<T>* Tran<T>::make_rotate_a2b(const glm::tvec3<T>& a, const glm::tvec3<T>& b, bool flip )
0216 {
0217     glm::tmat4x4<T> t = stra<T>::RotateA2B(a,b,flip) ;
0218     glm::tmat4x4<T> v = stra<T>::RotateA2B(b,a,flip) ;
0219     return new Tran<T>(t, v);
0220 }
0221
0222
0223
0224
0225
0226
0227 template<typename T>
0228 inline const Tran<T>* Tran<T>::product(const Tran<T>* a, const Tran<T>* b, bool reverse)
0229 {
0230     std::vector<const Tran<T>*> tt ;
0231     tt.push_back(a);
0232     tt.push_back(b);
0233     return Tran<T>::product( tt, reverse );
0234 }
0235
0236 template<typename T>
0237 inline const Tran<T>* Tran<T>::product(const Tran<T>* a, const Tran<T>* b, const Tran<T>* c, bool reverse)
0238 {
0239     std::vector<const Tran<T>*> tt ;
0240     tt.push_back(a);
0241     tt.push_back(b);
0242     tt.push_back(c);
0243     return Tran<T>::product( tt, reverse );
0244 }
0245
0246
0247
0248
0249
0250 /**
0251 Tran::product
0252 --------------
0253
0254 Tran houses transforms paired with their inverse transforms, the product
0255 of the transforms and opposite order product of the inverse transforms
0256 is done using inclusive indices to access Tran from the vector::
0257
0258    i: 0 -> ntt - 1      ascending
0259    j: ntt - 1 -> 0      descending (from last transform down to first)
0260
0261 Use *reverse=true* when the transforms are in reverse heirarchical order, ie when
0262 they have been collected by starting from the leaf node and then following parent
0263 links back up to the root node.
0264
0265 When combining s:scale r:rotation and t:translate the typical ordering is s-r-t
0266 because wish to scale and orient about a nearby (local frame) origin before
0267 translating into position.
0268
0269
0270 **/
0271 template<typename T>
0272 inline const Tran<T>* Tran<T>::product(const std::vector<const Tran<T>*>& tt, bool reverse )
0273 {
0274     unsigned ntt = tt.size();
0275     if(ntt==0) return NULL ;
0276     if(ntt==1) return tt[0] ;
0277
0278     glm::tmat4x4<T> t(T(1)) ;
0279     glm::tmat4x4<T> v(T(1)) ;
0280
0281     for(unsigned i=0,j=ntt-1 ; i < ntt ; i++,j-- )
0282     {
0283         const Tran<T>* ii = tt[reverse ? j : i] ;  // reverse:true starts from the last (ie root node)
0284         const Tran<T>* jj = tt[reverse ? i : j] ;  // reverse:true starts from the first (ie leaf node)
0285
0286         if( ii != nullptr )  t *= ii->t ;
0287         if( jj != nullptr )  v *= jj->v ;  // inverse-transform product in opposite order
0288     }
0289     return new Tran<T>(t, v) ;
0290 }
0291
0292
0293
0294
0295
0296 template<typename T>
0297 inline Tran<T>::Tran( const T* transform, const T* inverse )
0298     :
0299     t(glm::make_mat4x4<T>(transform)),
0300     v(glm::make_mat4x4<T>(inverse)),
0301     i(t*v)
0302 {
0303 }
0304
0305 template<typename T>
0306 inline Tran<T>::Tran( const glm::tmat4x4<T>& transform, const glm::tmat4x4<T>& inverse )
0307     :
0308     t(transform),
0309     v(inverse),
0310     i(transform*inverse)
0311 {
0312 }
0313
0314
0315 template<typename T>
0316 inline const T*  Tran<T>::tdata() const
0317 {
0318     return glm::value_ptr(t) ;
0319 }
0320
0321 template<typename T>
0322 inline const T*  Tran<T>::vdata() const
0323 {
0324     return glm::value_ptr(v) ;
0325 }
0326
0327
0328 template<typename T>
0329 inline T Tran<T>::maxdiff_from_identity(char mat) const
0330 {
0331     const glm::tmat4x4<T>& m = mat == 't' ? t : ( mat == 'v' ? v : i ) ;
0332     return stra<double>::Maxdiff_from_Identity(m);
0333 }
0334
0335 template<typename T>
0336 inline bool Tran<T>::is_identity(char mat, T epsilon) const
0337 {
0338     T mxdif = maxdiff_from_identity(mat) ;
0339     return mxdif < epsilon ;
0340 }
0341
0342
0343 template<typename T>
0344 inline std::string Tran<T>::brief(bool only_tlate, char mat, unsigned wid, unsigned prec) const
0345 {
0346     std::stringstream ss ;
0347     ss << mat << ":" ;
0348     if(is_identity(mat))
0349     {
0350        ss << "identity" ;
0351     }
0352     else
0353     {
0354         const glm::tmat4x4<T>& m = mat == 't' ? t : ( mat == 'v' ? v : i ) ;
0355         int j0 = only_tlate ? 3 : 0 ;
0356         for(int j=j0 ; j < 4 ; j++ )
0357         {
0358             ss << "[" ;
0359             for(int k=0 ; k < 4 ; k++ ) ss << std::setw(wid) << std::fixed << std::setprecision(prec) << m[j][k] << " " ;
0360             ss << "]" ;
0361         }
0362     }
0363     std::string s = ss.str() ;
0364     return s ;
0365 }
0366
0367 template<typename T>
0368 inline std::string Tran<T>::desc() const
0369 {
0370     bool only_tlate = false ;
0371     std::stringstream ss ;
0372     ss << brief(only_tlate, 't' ) << std::endl ;
0373     ss << brief(only_tlate, 'v' ) << std::endl ;
0374     ss << brief(only_tlate, 'i' ) << std::endl ;
0375     std::string s = ss.str() ;
0376     return s ;
0377 }
0378
0379 template<typename T>
0380 bool Tran<T>::checkIsIdentity(char mat, const char* caller, T epsilon)
0381 {
0382     bool ok = is_identity('i', epsilon);
0383     if(!ok)
0384     {
0385         T mxdif = maxdiff_from_identity('i');
0386         std::cerr
0387             << "stran.h : Tran::checkIsIdentity FAIL : "
0388             << " caller " << caller
0389             << " epsilon " << epsilon
0390             << " mxdif_from_identity " << mxdif
0391             << std::endl
0392             ;
0393
0394         if(getenv("stran_checkIsIdentity_SIGINT")) std::raise(SIGINT);
0395     }
0396     return ok ;
0397 }
0398
0399 template<typename T>
0400 Tran<T>* Tran<T>::ConvertToTran(const qat4* q_, T epsilon )
0401 {
0402     qat4 q(q_->cdata());
0403     q.clearIdentity();
0404
0405     glm::tmat4x4<T> tran = MatFromQat(&q) ;
0406     glm::tmat4x4<T> itra = glm::inverse(tran) ;
0407     Tran<T>* tr = new Tran<T>(tran, itra) ;
0408     bool ok = tr->checkIsIdentity('i', "ConvertToTran");
0409     if(!ok) std::cerr << "stran.h Tran::ConvertToTran checkIsIdentity FAIL " << std::endl ;
0410
0411     return tr ;
0412 }
0413
0414 /**
0415 Tran<T>::ConvertFromQat
0416 -------------------------
0417
0418 Caller is owner of the Tran
0419
0420 **/
0421
0422 template<typename T>
0423 Tran<T>* Tran<T>::ConvertFromQat(const qat4* q_, T epsilon )
0424 {
0425     qat4 q(q_->cdata());
0426     q.clearIdentity();
0427
0428     glm::tmat4x4<T> tran = MatFromQat(&q) ;
0429     glm::tmat4x4<T> itra = glm::inverse(tran) ;
0430     Tran<T>* tr = new Tran<T>(tran, itra) ;
0431     bool ok = tr->checkIsIdentity('i', "ConvertFromQat");
0432     if(!ok) std::cerr << "stran.h Tran::ConvertFromQat checkIsIdentity FAIL " << std::endl ;
0433     return tr ;
0434 }
0435
0436
0437
0438
0439 template<typename T>
0440 Tran<T>* Tran<T>::ConvertFromData(const T* data )
0441 {
0442     glm::tmat4x4<T> tran = MatFromData(data) ;
0443     glm::tmat4x4<T> itra = glm::inverse(tran) ;
0444     Tran<T>* tr = new Tran<T>(tran, itra) ;
0445     bool ok = tr->checkIsIdentity('i', "ConvertFromData");
0446     if(!ok) std::cerr << "stran.h Tran::ConvertFromData checkIsIdentity FAIL " << std::endl ;
0447     return tr ;
0448 }
0449
0450
0451
0452
0453
0454 template<typename T>
0455 const qat4* Tran<T>::Invert( const qat4* q, T epsilon )
0456 {
0457
0458     int ins_idx,  gas_idx, sensor_identifier, sensor_index ;
0459     q->getIdentity(ins_idx,  gas_idx, sensor_identifier, sensor_index );
0460
0461     Tran<T>* tr = ConvertToTran(q) ;
0462
0463     qat4* v = ConvertFrom(tr->v);
0464     v->setIdentity(ins_idx, gas_idx, sensor_identifier, sensor_index ) ;
0465
0466     return v ;
0467 }
0468
0469
0470 template<typename T>
0471 Tran<T>* Tran<T>::FromPair(const qat4* t, const qat4* v, T epsilon ) // static
0472 {
0473     glm::tmat4x4<T> tran = MatFromQat(t) ;
0474     glm::tmat4x4<T> itra = MatFromQat(v) ;
0475     Tran<T>* tr = new Tran<T>(tran, itra) ;
0476     bool ok = tr->checkIsIdentity('i', "FromPair", epsilon );
0477     if(!ok)
0478     {
0479          std::cerr << "stran.h Tran::FromPair checkIsIdentity FAIL " << std::endl ;
0480          const char* path = "/tmp/stran_FromPair_checkIsIdentity_FAIL.npy" ;
0481          std::cerr << "stran.h save to path " << path << std::endl ;
0482          tr->save_(path);
0483     }
0484     return tr ;
0485 }
0486
0487 template<typename T>
0488 glm::tmat4x4<T> Tran<T>::MatFromQat(const qat4* q )  // static
0489 {
0490     const float* q_data = q->cdata();
0491     glm::tmat4x4<T> tran(1.);
0492     T* tran_ptr = glm::value_ptr(tran) ;
0493     for(int i=0 ; i < 16 ; i++) tran_ptr[i] = T(q_data[i]) ;
0494     return tran ;
0495 }
0496
0497 template<typename T>
0498 glm::tmat4x4<T> Tran<T>::MatFromData(const T* data)  // static
0499 {
0500     glm::tmat4x4<T> tran(1.);
0501     T* tran_ptr = glm::value_ptr(tran) ;
0502     for(int i=0 ; i < 16 ; i++) tran_ptr[i] = data[i] ;
0503     return tran ;
0504 }
0505
0506
0507
0508
0509
0510
0511
0512
0513
0514
0515
0516
0517
0518
0519
0520
0521
0522
0523 /**
0524 Tran::ConvertFrom
0525 -------------------
0526
0527 With T=double will narrow to floats within qat4
0528
0529 **/
0530
0531 template<typename T>
0532 qat4* Tran<T>::ConvertFrom(const glm::tmat4x4<T>& transform )
0533 {
0534     const T* ptr = glm::value_ptr(transform) ;
0535
0536     float ff[16] ;
0537
0538     for(int i=0 ; i < 16 ; i++ ) ff[i] = float(ptr[i]) ;
0539
0540     return new qat4(ff) ;
0541 }
0542
0543 template<typename T>
0544 void Tran<T>::write(T* dst, unsigned num_values) const
0545 {
0546     unsigned matrix_values = 4*4 ;
0547     assert( num_values == 3*matrix_values );
0548
0549     unsigned matrix_bytes = matrix_values*sizeof(T) ;
0550     char* dst_bytes = (char*)dst ;
0551
0552     memcpy( dst_bytes + 0*matrix_bytes , (char*)glm::value_ptr(t), matrix_bytes );
0553     memcpy( dst_bytes + 1*matrix_bytes , (char*)glm::value_ptr(v), matrix_bytes );
0554     memcpy( dst_bytes + 2*matrix_bytes , (char*)glm::value_ptr(i), matrix_bytes );
0555 }
0556
0557 template<typename T>
0558 NP* Tran<T>::as_array()const
0559 {
0560     NP* a = NP::Make<T>(3, 4, 4 );
0561     unsigned num_values = 3*4*4 ;
0562     write( a->values<T>(), num_values ) ;
0563     return a ;
0564 }
0565
0566 template<typename T>
0567 void Tran<T>::save(const char* dir, const char* name) const
0568 {
0569     NP* a = as_array();
0570     a->save(dir, name);
0571 }
0572 template<typename T>
0573 void Tran<T>::save_(const char* path) const
0574 {
0575     NP* a = as_array();
0576     a->save(path);
0577 }
0578
0579
0580
0581
0582
0583
0584 /**
0585 Tran::Apply
0586 -------------
0587
0588 Applies the transform *count* times using 4th component w
0589 (w is usually 1. for transforming as a position or 0. for transforming as a direction).
0590
0591 **/
0592
0593 template<typename T>
0594 void Tran<T>::Apply( const glm::tmat4x4<T>& tr,  T* p0, T w, unsigned count, unsigned stride, unsigned offset, bool normalize )
0595 {
0596     for(unsigned i=0 ; i < count ; i++)
0597     {
0598         T* a_ = p0 + i*stride + offset ;
0599
0600         glm::tvec4<T> a(a_[0],a_[1],a_[2],w);
0601         glm::tvec4<T> ta = tr * a ;
0602         glm::tvec3<T> ta3(ta);
0603         glm::tvec3<T> nta3 = normalize ? glm::normalize(ta3) : ta3  ;
0604         T* ta_ = glm::value_ptr( nta3 ) ;
0605
0606         for(unsigned j=0 ; j < 3 ; j++) a_[j] = ta_[j] ;
0607
0608         //std::cout << " apply: a " << glm::to_string( a ) << std::endl ;
0609         //std::cout << " apply: ta= " << glm::to_string( ta ) << std::endl ;
0610     }
0611 }
0612
0613 /**
0614 Tran::ApplyToFloat : inplace application of transform
0615 ----------------------------------------------------
0616
0617 Unlike the above Tran::apply this performs the multiplication
0618 in the T precision even though the input is in float precision.
0619 The is pragmatic approach, as it will often be preferable to use a
0620 double precision transform and single precision array.
0621
0622 1. widening copy pos/mom/pol into tvec4 using w=1./0./0. for pos/mom/pol
0623 2. left-multiply the "t" transform with the pos/mom/pol, tvec3 normalize when chosen
0624 3. copy back into src array
0625
0626 **/
0627
0628
0629 template<typename T>
0630 void Tran<T>::ApplyToFloat( const glm::tmat4x4<T>& tr, float* p0, float w, unsigned count, unsigned stride, unsigned offset, bool normalize )
0631 {
0632     for(unsigned i=0 ; i < count ; i++)
0633     {
0634         float* a_ = p0 + i*stride + offset ;
0635
0636         // 1. widening copy pos/mom/pol into tvec4 using w=1./0./0. for pos/mom/pol
0637         glm::tvec4<T> a(0.,0.,0.,0.)  ;
0638         T* aa = glm::value_ptr(a) ;
0639         for(unsigned j=0 ; j < 3 ; j++) aa[j] = T(a_[j]) ;  // potentially widen
0640         aa[3] = T(w) ;
0641
0642         // 2. left-multiply the "t" transform with the pos/mom/pol, tvec3 normalize when chosen
0643         glm::tvec4<T> ta = tr * a ;
0644         glm::tvec3<T> ta3(ta);
0645         glm::tvec3<T> nta3 = normalize ? glm::normalize(ta3) : ta3  ;
0646         T* ta_ = glm::value_ptr( nta3 ) ;
0647
0648         // 3. copy back into src array
0649         for(unsigned j=0 ; j < 3 ; j++) a_[j] = float(ta_[j]) ;   // potentially narrow
0650     }
0651 }
0652
0653
0654 template<typename T>
0655 void Tran<T>::photon_transform( NP* ph, bool normalize, bool inverse ) const
0656 {
0657     T one(1.);
0658     T zero(0.);
0659
0660     const glm::tmat4x4<T>& tv = inverse ? v : t ;  // transform
0661
0662     assert( ph->has_shape(-1,4,4) );
0663     size_t count  = ph->shape[0] ;
0664     size_t stride = 4*4 ;
0665
0666     if( ph->ebyte == sizeof(T) )
0667     {
0668         T* p0 = ph->values<T>();
0669         Apply( tv, p0, one,  count, stride, 0, false );  // transform pos as position
0670         Apply( tv, p0, zero, count, stride, 4, normalize );  // transform mom as direction
0671         Apply( tv, p0, zero, count, stride, 8, normalize );  // transform pol as direction
0672     }
0673     else if( ph->ebyte == 4 && sizeof(T) == 8 )
0674     {
0675         float* p0 = ph->values<float>();
0676         ApplyToFloat( tv, p0, one,  count, stride, 0, false );  // transform pos as position
0677         ApplyToFloat( tv, p0, zero, count, stride, 4, normalize );  // transform mom as direction
0678         ApplyToFloat( tv, p0, zero, count, stride, 8, normalize );  // transform pol as direction
0679     }
0680 }
0681
0682 /**
0683 Tran::PhotonTransform
0684 -----------------------
0685
0686 Test::
0687
0688    cd ~/opticks/CSG/tests
0689    ./CSGFoundry_getFrame_Test.sh
0690
0691 Note that the returned transformed photon array is always in double precision.
0692
0693 **/
0694
0695 template<typename T>
0696 NP* Tran<T>::PhotonTransform( const NP* ph, bool normalize, const Tran<T>* tr, bool inverse ) // static
0697 {
0698     NP* b = NP::MakeWideIfNarrow(ph) ;
0699     tr->photon_transform(b, normalize, inverse );
0700     assert( b->ebyte == 8 );
0701     return b ;
0702 }
0703
0704
0705 /**
0706 Tran::AddTransform
0707 ---------------------
0708
0709 * The option string determines the meaning of the a,b,c vectors
0710 * lowercase/UPPERCASE controls flip:true/FALSE
0711
0712 +--------+-------------+------------------------------+-------------------+
0713 | opt    |  flip       |  impl                        |  note             |
0714 +========+=============+==============================+===================+
0715 | TR/tr  | false/true  |  stra::Place(a,b,c,flip)     |                   |
0716 +--------+-------------+------------------------------+-------------------+
0717 | R/r    | false/true  |  stra::RotateA2B(a,b, flip)  | c is ignored      |
0718 +--------+-------------+------------------------------+-------------------+
0719 | T/t    | false/true  |  stra::Translate(c, flip)    | a,b ignored       |
0720 +--------+-------------+------------------------------+-------------------+
0721 | D/d    | false/true  |  stra::Dupe(a,b,c,flip)      |                   |
0722 +--------+-------------+------------------------------+-------------------+
0723
0724 TODO: move this into stra.h
0725
0726 **/
0727
0728 template<typename T>
0729 inline void Tran<T>::AddTransform( T* ttk, const char* opt, const glm::tvec3<T>& a, const glm::tvec3<T>& b, const glm::tvec3<T>& c )
0730 {
0731     for(unsigned i=0 ; i < 16 ; i++) ttk[i] = T(0.) ; // start zeroed
0732
0733     if(strcmp(opt,"TR") == 0 || strcmp(opt,"tr") == 0 )
0734     {
0735         bool flip = strcmp(opt,"tr") == 0 ;
0736         glm::tmat4x4<T> tr = stra<T>::Place( a, b, c, flip );
0737         T* src = glm::value_ptr(tr) ;
0738         //std::cout << Desc(src, 16) << std::endl ;
0739         memcpy( ttk, src, 16*sizeof(T) );
0740     }
0741     else if(strcmp(opt,"R") == 0 || strcmp(opt,"r") == 0)
0742     {
0743         bool flip = strcmp(opt,"r") == 0 ;
0744         glm::tmat4x4<T> tr = stra<T>::RotateA2B( a, b, flip );
0745         T* src = glm::value_ptr(tr) ;
0746         memcpy( ttk, src, 16*sizeof(T) );
0747     }
0748     else if(strcmp(opt,"T") == 0 || strcmp(opt,"t") == 0)
0749     {
0750         bool flip = strcmp(opt,"t") == 0 ;
0751         glm::tmat4x4<T> tr = stra<T>::Translate(c, flip );
0752         T* src = glm::value_ptr(tr) ;
0753         memcpy( ttk, src, 16*sizeof(T) );
0754     }
0755     else if(strcmp(opt,"D")== 0 || strcmp(opt, "d") == 0)
0756     {
0757         bool flip = strcmp(opt,"d") == 0 ;
0758         glm::tmat4x4<T> tr = stra<T>::Dupe(a, b, c, flip );
0759         T* src = glm::value_ptr(tr) ;
0760         memcpy( ttk, src, 16*sizeof(T) );
0761     }
0762     else
0763     {
0764         std::cout << "Tran::AddTransform : ERROR opt is not handled [" << opt << "]" << std::endl ;
0765     }
0766 }
0767
0768
0769
0770
0771 template struct Tran<float> ;
0772 template struct Tran<double> ;
0773
0774
0775