eic-opticks/ana/tangential.cc

0001 /**
0002 tangential.cc
0003 ===============
0004
0005 1. Creates tangential frame transforms for a point on a sphere
0006    configured via envvar RTP  radius-theta-phi where theta and
0007    phi are expressed in units of pi.
0008
0009 2. Performs test conversions between conventional and tangential frames
0010    which are dumped to stdout
0011
0012 **/
0013
0014 #include <cstdlib>
0015 #include <vector>
0016 #include <array>
0017 #include <iostream>
0018 #include <iomanip>
0019 #include <string>
0020 #include <sstream>
0021
0022 #include <glm/glm.hpp>
0023 #include <glm/gtx/string_cast.hpp>
0024 #include <glm/gtc/type_ptr.hpp>
0025
0026
0027 template<typename T>
0028 struct Tangential
0029 {
0030     static const T pi ;
0031     static void GetEVector( std::vector<T>& vec, const char* ekey, const char* fallback );
0032     static glm::tmat4x4<T> Rotation( T theta, T phi ) ;
0033     static glm::tmat4x4<T> IRotation( T theta, T phi ) ;
0034     static glm::tmat4x4<T> Translation( T radius,  T theta, T phi ) ;
0035     static glm::tmat4x4<T> ITranslation( T radius, T theta, T phi ) ;
0036
0037     static void CartesianToSpherical( glm::tvec3<T>& radius_theta_phi, const glm::tvec4<T>& xyzw );
0038
0039     Tangential( T radius, T theta, T phi );
0040
0041     void conventional_to_tangential( glm::tvec4<T>& rtpw , const glm::tvec4<T>& xyzw );
0042     void tangential_to_conventional( glm::tvec4<T>& xyzw , const glm::tvec4<T>& rtpw );
0043     void dump() const ;
0044
0045     glm::tmat4x4<T> rot ;
0046     glm::tmat4x4<T> iro ;
0047     glm::tmat4x4<T> tra ;
0048     glm::tmat4x4<T> itr ;
0049
0050     glm::tmat4x4<T> rot_itr ;
0051     glm::tmat4x4<T> tra_iro ;
0052 };
0053
0054 template<typename T>
0055 const T Tangential<T>::pi = glm::pi<T>() ;
0056
0057
0058 template<typename T>
0059 void Tangential<T>::CartesianToSpherical( glm::tvec3<T>& radius_theta_phi, const glm::tvec4<T>& xyzw )
0060 {
0061     const T x = xyzw.x ;
0062     const T y = xyzw.y ;
0063     const T z = xyzw.z ;
0064     const T radius =  sqrt( x*x + y*y + z*z ) ;
0065     const T zero(0.) ;
0066
0067     const T theta = radius == zero ? zero : acos( z/radius );
0068     const T phi = atan2(y, x) ;
0069
0070     radius_theta_phi.x = radius ;
0071     radius_theta_phi.y = theta ;
0072     radius_theta_phi.z = phi ;
0073 }
0074
0075
0076
0077 /**
0078 Tangential::Rotation
0079 -----------------------
0080
0081 Spherical unit vectors (ru tu pu) related to cartesian unit vectors (xu yu zu)
0082 via orthogonal rotation matrix::
0083
0084       | ru |     |   sin(theta)cos(phi)    sin(theta)sin(phi)      cos(theta)   |  | xu |
0085       |    |     |                                                              |  |    |
0086       | tu | =   |  cos(theta)cos(phi)    cos(theta)sin(phi)     -sin(theta)    |  | yu |
0087       |    |     |                                                              |  |    |
0088       | pu |     |  -sin(phi)                 cos(phi)              0           |  | zu |
0089
0090 **/
0091
0092 template<typename T>
0093 glm::tmat4x4<T> Tangential<T>::Rotation( T theta, T phi )
0094 {
0095     std::array<T, 16> _rot = {{
0096         sin(theta)*cos(phi),   sin(theta)*sin(phi),   cos(theta),  0. ,
0097         cos(theta)*cos(phi),   cos(theta)*sin(phi),  -sin(theta),  0. ,
0098         -sin(phi),             cos(phi),              0.,          0. ,
0099          0.,                   0.,                    0.,          1.
0100       }} ;
0101     return glm::make_mat4x4<T>(_rot.data()) ;
0102 }
0103
0104 /**
0105 Tangential::IRotation
0106 -----------------------
0107
0108 Cartesian unit vectors (xu yu zu) in terms of spherical unit vectors related by
0109 the inverse of the above transform which is its transpose::
0110
0111
0112       | xu |     |   sin(theta)cos(phi)    cos(theta)cos(phi)      -sin(phi)    |  | ru |
0113       |    |     |                                                              |  |    |
0114       | yu | =   |  sin(theta)sin(phi)    cos(theta)sin(phi)        cos(phi)    |  | tu |
0115       |    |     |                                                              |  |    |
0116       | zu |     |   cos(theta)               -sin(theta)              0        |  | pu |
0117
0118 **/
0119
0120 template<typename T>
0121 glm::tmat4x4<T> Tangential<T>::IRotation(  T theta, T phi )
0122 {
0123     std::array<T, 16> _iro = {{
0124         sin(theta)*cos(phi),   cos(theta)*cos(phi),  -sin(phi),     0. ,
0125         sin(theta)*sin(phi),   cos(theta)*sin(phi),  cos(phi),      0. ,
0126         cos(theta),            -sin(theta),          0.,            0. ,
0127          0.,                   0.,                   0.,            1.
0128       }} ;
0129     return glm::make_mat4x4<T>(_iro.data()) ;
0130 }
0131
0132 template<typename T>
0133 glm::tmat4x4<T> Tangential<T>::Translation( T radius, T theta, T phi )
0134 {
0135     std::array<T, 16> _tra = {{
0136          1.,                           0.,                          0.,                  0. ,
0137          0.,                           1.,                          0.,                  0. ,
0138          0.,                           0.,                          1.,                  0. ,
0139          radius*sin(theta)*cos(phi),   radius*sin(theta)*sin(phi),  radius*cos(theta),   1.
0140       }} ;
0141     return glm::make_mat4x4<T>(_tra.data()) ;
0142 }
0143
0144 template<typename T>
0145 glm::tmat4x4<T> Tangential<T>::ITranslation( T radius, T theta, T phi )
0146 {
0147     std::array<T, 16> _itr = {{
0148          1.,                           0.,                          0.,                  0. ,
0149          0.,                           1.,                          0.,                  0. ,
0150          0.,                           0.,                          1.,                  0. ,
0151          -radius*sin(theta)*cos(phi),  -radius*sin(theta)*sin(phi), -radius*cos(theta),  1.
0152       }} ;
0153     return glm::make_mat4x4<T>(_itr.data()) ;
0154 }
0155
0156
0157 template<typename T>
0158 Tangential<T>::Tangential(T radius, T theta, T phi )
0159     :
0160     rot(Rotation(pi*theta,pi*phi)),
0161     iro(IRotation(pi*theta,pi*phi)),
0162     tra(Translation(radius,pi*theta,pi*phi)),
0163     itr(ITranslation(radius,pi*theta,pi*phi)),
0164     rot_itr(rot * itr),
0165     tra_iro(tra * iro)
0166 {
0167 }
0168
0169
0170 template<typename T>
0171 void Tangential<T>::tangential_to_conventional( glm::tvec4<T>& xyzw , const glm::tvec4<T>& rtpw )
0172 {
0173     //xyzw = rot * itr * rtpw ;
0174     xyzw = rot_itr * rtpw ;
0175 }
0176
0177 template<typename T>
0178 void Tangential<T>::conventional_to_tangential( glm::tvec4<T>& rtpw , const glm::tvec4<T>& xyzw )
0179 {
0180     //rtpw = tra * iro * xyzw ;
0181     rtpw = tra_iro * xyzw ;
0182 }
0183
0184
0185 template<typename T>
0186 void Tangential<T>::dump() const
0187 {
0188     std::cout  << " rot " << glm::to_string(rot) << std::endl ;
0189     std::cout  << " iro " << glm::to_string(iro) << std::endl ;
0190     std::cout  << " tra " << glm::to_string(tra) << std::endl ;
0191     std::cout  << " itr " << glm::to_string(itr) << std::endl ;
0192 }
0193
0194 template<typename T>
0195 void Tangential<T>::GetEVector( std::vector<T>& vec, const char* ekey, const char* fallback )  // static
0196 {
0197     const char* sval = getenv(ekey);
0198     std::stringstream ss(sval ? sval : fallback);
0199     std::string s ;
0200     while(getline(ss, s, ','))
0201     {
0202         T val(0) ;
0203         std::stringstream vv(s.c_str());
0204         vv >> val ;
0205         vec.push_back(val);
0206     }
0207 }
0208
0209 int main(int argc, char** argv)
0210 {
0211     std::vector<double> rtp ;
0212     Tangential<double>::GetEVector(rtp, "RTP", "20,0.25,0.0") ;
0213     for(unsigned i=0 ; i < rtp.size() ; i++) std::cout << std::setw(4) << i << " : " << rtp[i] << std::endl ;
0214     assert( rtp.size() == 3 );
0215     double radius = rtp[0] ;
0216     double theta  = rtp[1] ;
0217     double phi    = rtp[2] ;
0218
0219     Tangential<double> ta(radius, theta, phi );
0220     ta.dump();
0221
0222     glm::tvec4<double> rtpw(0., 0., 0., 1.) ;
0223     glm::tvec4<double> xyzw(0., 0., 0., 1.) ;
0224     glm::tvec4<double> rtpw2(0., 0., 0., 1.) ;
0225
0226     for(int t=-3 ; t <= 3 ; t++)
0227     for(int p=-3 ; p <= 3 ; p++)
0228     {
0229         rtpw.y = double(t) ;
0230         rtpw.z = double(p) ;
0231
0232         ta.tangential_to_conventional(xyzw, rtpw );
0233         ta.conventional_to_tangential(rtpw2, xyzw );
0234
0235         std::cout
0236             << "rtpw"
0237             << "(" << std::setw(10) << std::fixed << std::setprecision(4) << rtpw.x
0238             << " " << std::setw(10) << std::fixed << std::setprecision(4) << rtpw.y
0239             << " " << std::setw(10) << std::fixed << std::setprecision(4) << rtpw.z
0240             << " " << std::setw(10) << std::fixed << std::setprecision(4) << rtpw.w
0241             << ")"
0242             << " xyzw "
0243             << "(" << std::setw(10) << std::fixed << std::setprecision(4) << xyzw.x
0244             << " " << std::setw(10) << std::fixed << std::setprecision(4) << xyzw.y
0245             << " " << std::setw(10) << std::fixed << std::setprecision(4) << xyzw.z
0246             << " " << std::setw(10) << std::fixed << std::setprecision(4) << xyzw.w
0247             << ")"
0248             << "rtpw2"
0249             << "(" << std::setw(10) << std::fixed << std::setprecision(4) << rtpw2.x
0250             << " " << std::setw(10) << std::fixed << std::setprecision(4) << rtpw2.y
0251             << " " << std::setw(10) << std::fixed << std::setprecision(4) << rtpw2.z
0252             << " " << std::setw(10) << std::fixed << std::setprecision(4) << rtpw2.w
0253             << ")"
0254             << std::endl
0255             ;
0256     }
0257
0258     return 0 ;
0259 }