Changeset 470 in Sophya for trunk/SophyaLib/Samba/localmap.cc


Ignore:
Timestamp:
Oct 15, 1999, 5:43:30 PM (26 years ago)
Author:
ansari
Message:

versions templatees, NdataBlocks etc. 15-OCT-99-GLM

File:
1 edited

Legend:

Unmodified
Added
Removed

  • trunk/SophyaLib/Samba/localmap.cc

    r276 r470  
    11#include "localmap.h"
     2#include "nbmath.h"
     3#include <complex>
     4#include "piocmplx.h"
     5
     6#include <string.h>
    27#include <iostream.h>
    3 #include "nbmath.h"
     8
     9#ifdef __CXX_PRAGMA_TEMPLATES__
     10#pragma define_template LocalMap<double>
     11#pragma define_template LocalMap<float>
     12#pragma define_template LocalMap< complex<double> >
     13#pragma define_template LocalMap< complex<float> >
     14#pragma define_template FIO_LocalMap<double>
     15#pragma define_template FIO_LocalMap<float>
     16#pragma define_template FIO_LocalMap< complex<float> >
     17#pragma define_template FIO_LocalMap< complex<double> >
     18#endif
     19#if defined(ANSI_TEMPLATES) || defined(GNU_TEMPLATES)
     20template class LocalMap<double>;
     21template class LocalMap<float>;
     22template class LocalMap< complex<double> >;
     23template class LocalMap< complex<float> >;
     24template class FIO_LocalMap<double>;
     25template class FIO_LocalMap<float>;
     26template class FIO_LocalMap< complex<float> >;
     27template class FIO_LocalMap< complex<double> >;
     28#endif
     29
    430//*****************************************************************************
    531//++
    … …  
    5278//--
    5379//++
    54 
    55 LocalMap::LocalMap()
    56 
    57 //--
    58 
    59 {InitNul();}
    60 
    61 //++
    62 
    63 LocalMap::LocalMap(int_4 nx, int_4 ny) : mSzX_(nx), mSzY_(ny)
    64 
     80template<class T>
     81LocalMap<T>::LocalMap()
     82//
     83// Constructeur par defaut
    6584//--
    6685{
    6786  InitNul();
    68   mNPix_=nx*ny;
    69   mPix_ = new r_8[mNPix_];
    70 }
    71 
    72 
    73 //++
    74 // Titre        Destructeur
    75 //--
    76 //++
    77 LocalMap::~LocalMap()
    78 
    79 //--
    80 {
    81   Clear();
    82 }
    83 
    84 
    85 
    86 void        LocalMap::WriteSelf(POutPersist& s) const
    87 {
    88   char strg[256];
    89   if (mInfo_) {sprintf(strg, "LocalMap: NPixX=%6d  NPixY=%9d HasInfo",  mSzX_, mSzY_);} else  {
    90     sprintf(strg, "LocalMap: NPixX=%6d  NPixY=%9d ",  mSzX_, mSzY_);
    91   }         
    92   s.PutLine(strg);
    93   if (mInfo_)  s << (*mInfo_);
    94   s.PutI4(mSzX_);
    95   s.PutI4(mSzY_);
    96   s.PutI4(mNPix_);
    97   s.PutI4(x0_);
    98   s.PutI4(y0_);
    99   s.PutI4(originFlag_);
    100   s.PutI4(SzFlag_);
    101   s.PutR4(cos_angle_);
    102   s.PutR4(sin_angle_);
    103   s.PutR8(theta0_);
    104   s.PutR8( phi0_);
    105   s.PutR8( tgAngleX_);
    106   s.PutR8(tgAngleY_);
    107   s.PutR8s(mPix_, mNPix_);
    108 }
    109 void        LocalMap::ReadSelf(PInPersist& s)
    110 {
    111   Clear();
    112   char strg[256];
    113   s.GetLine(strg, 255);
    114 // Pour savoir s'il y avait un DVList Info associe
    115   bool hadinfo = false;
    116   if (strncmp(strg+strlen(strg)-7, "HasInfo", 7) == 0)  hadinfo = true;
    117   if (hadinfo) {    // Lecture eventuelle du DVList Info
    118     if (mInfo_ == NULL)  mInfo_ = new DVList;
    119     s >> (*mInfo_);
    120   }
    121   s.GetI4(mSzX_);
    122   s.GetI4(mSzY_);
    123   s.GetI4(mNPix_);
    124   s.GetI4(x0_);
    125   s.GetI4(y0_);
    126   s.GetI4(originFlag_);
    127   s.GetI4(SzFlag_);
    128   s.GetR4(cos_angle_);
    129   s.GetR4(sin_angle_);
    130   s.GetR8(theta0_);
    131   s.GetR8( phi0_);
    132   s.GetR8( tgAngleX_);
    133   s.GetR8(tgAngleY_);
    134   mPix_ = new r_8[mNPix_+1];
    135   s.GetR8s(mPix_, mNPix_);
     87}
     88
     89//++
     90template<class T>
     91LocalMap<T>::LocalMap(int_4 nx, int_4 ny) : nSzX_(nx), nSzY_(ny)
     92//
     93// Constructeur
     94//--
     95{
     96  InitNul();
     97  nPix_= nx*ny;
     98  pixels_.ReSize(nPix_);
     99  pixels_.Reset();
     100}
     101
     102//++
     103template<class T>
     104LocalMap<T>::LocalMap(const LocalMap<T>& lm)
     105//
     106// Constructeur de recopie
     107//--
     108{
     109  cout<<" LocalMap:: Appel du constructeur de recopie " << endl;
     110
     111  if(lm.mInfo_) mInfo_= new DVList(*lm.mInfo_);
     112  nSzX_= lm.nSzX_;
     113  nSzY_= lm.nSzY_;
     114  nPix_= lm.nPix_;
     115  originFlag_= lm.originFlag_;
     116  extensFlag_= lm.extensFlag_;
     117  x0_= lm.x0_;
     118  y0_= lm.y0_;
     119  theta0_= lm.theta0_;
     120  phi0_= lm.phi0_;
     121  angle_= lm.angle_;
     122  cos_angle_= lm.cos_angle_;
     123  sin_angle_= lm.sin_angle_;
     124  angleX_= lm.angleX_;
     125  angleY_= lm.angleY_;
     126  tgAngleX_= lm.tgAngleX_;
     127  tgAngleY_= lm.tgAngleY_;
     128  pixels_= lm.pixels_;
     129}
     130
     131//++
     132template<class T>
     133LocalMap<T>::~LocalMap()
     134//
     135// Destructeur
     136//--
     137{
     138  InitNul();
     139}
     140
     141//++
     142template<class T>
     143void LocalMap<T>::InitNul()
     144//
     145// Initialise à zéro certaines variables internes
     146//--
     147{
     148  originFlag_= false;
     149  extensFlag_= false;
     150  cos_angle_= 1.0;
     151  sin_angle_= 0.0;
     152  pixels_.Reset();
     153}
     154
     155//++
     156template<class T>
     157void LocalMap<T>::ReSize(int_4 nx, int_4 ny)
     158//
     159// Redimensionne l'espace de stokage des pixels
     160//--
     161{
     162  InitNul();
     163  nSzX_ = nx;
     164  nSzY_ = ny;
     165  nPix_= nx*ny;
     166  pixels_.ReSize(nPix_);
     167  pixels_.Reset();
    136168}
    137169
    … …  
    142174
    143175//++
    144 int_4 LocalMap::NbPixels() const
    145 
     176template<class T>
     177int_4 LocalMap<T>::NbPixels() const
     178//
    146179//    Retourne le nombre de pixels du découpage
    147180//--
    148181{
    149   return(mNPix_);
    150 }
    151 
    152 //++
    153 r_8& LocalMap::PixVal(int_4 k)
    154 
     182  return(nPix_);
     183}
     184
     185//++
     186template<class T>
     187T& LocalMap<T>::PixVal(int_4 k)
     188//
    155189//    Retourne la valeur du contenu du pixel d'indice k
    156190//--
    157191{
    158   if ( (k<0) || (k >= mNPix_) ) {
    159     cout << " LocalMap::PIxVal : exceptions  a mettre en place" <<endl;
    160     //  THROW(out_of_range("LocalMap::PIxVal Pixel index out of range"));
    161   THROW(rangeCheckErr);
    162 
    163     //throw "LocalMap::PIxVal Pixel index out of range";
    164   }
    165   //return (*(pixel_->Data()+k));
    166   return(mPix_[k]);
    167 }
    168 
    169 //++
    170 
    171 r_8 const& LocalMap::PixVal(int_4 k) const
    172 
     192  if((k < 0) || (k >= nPix_))
     193    {
     194      cout << " LocalMap::PIxVal : exceptions  a mettre en place" <<endl;
     195      //  THROW(out_of_range("LocalMap::PIxVal Pixel index out of range"));
     196      THROW(rangeCheckErr);
     197      //throw "LocalMap::PIxVal Pixel index out of range";
     198    }
     199  return(pixels_(k));
     200}
     201
     202//++
     203
     204template<class T>
     205T const& LocalMap<T>::PixVal(int_4 k) const
     206//
    173207//    Retourne la valeur du contenu du pixel d'indice k
    174208//--
    175209{
    176   if ( (k<0) || (k >= mNPix_) ) {
    177     cout << " LocalMap::PIxVal : exceptions  a mettre en place" <<endl;
    178   //    THROW(out_of_range("LocalMap::PIxVal Pixel index out of range"));
    179 
    180     throw "LocalMap::PIxVal Pixel index out of range";
    181   }
    182   //return *(pixel_->Data()+k);
    183   return(mPix_[k]);
    184 }
    185 
    186 //++
    187 int_4 LocalMap::PixIndexSph(float theta, float phi) const
    188 
     210  if((k < 0) || (k >= nPix_))
     211    {
     212      cout << " LocalMap::PIxVal : exceptions  a mettre en place" <<endl;
     213      //    THROW(out_of_range("LocalMap::PIxVal Pixel index out of range"));
     214     
     215      throw "LocalMap::PIxVal Pixel index out of range";
     216    }
     217  return *(pixels_.Data()+k);
     218}
     219
     220//++
     221template<class T>
     222int_4 LocalMap<T>::PixIndexSph(float theta, float phi) const
     223//
    189224//    Retourne l'indice du pixel vers lequel pointe une direction définie par
    190225//    ses coordonnées sphériques
    … …  
    192227{
    193228  int i,j;
    194   if (!(originFlag_==1) || !(SzFlag_==1) ) {
    195     cout << " LocalMap: correspondance carte-sphere non etablie" << endl;
    196     exit(0);
    197   }
    198   // theta et phi en coordonnees relatives (on se ramene a une situation
    199   // par rapport au plan de reference)
     229  if(!(originFlag_) || !(extensFlag_))
     230    {
     231      cout << " LocalMap: correspondance carte-sphere non etablie" << endl;
     232      exit(0);
     233    }
     234
     235  // theta et phi en coordonnees relatives (on se ramene a une situation par rapport au plan de reference)
    200236  float theta_aux=theta;
    201237  float phi_aux=phi;
    202238  UserToReference(theta_aux, phi_aux);
     239
    203240  // coordonnees dans le plan local en unites de pixels
    204241  float x,y;
    205   AngleProjToPix(theta_aux, phi_aux, x,y);
    206 
    207   float xmin=-x0_-0.5;
    208   float xmax=xmin+mSzX_;
    209   if( (x >  xmax) || (x < xmin ) ) return(-1);
    210   float xcurrent=xmin;
    211   for( i=0; i < mSzX_; i++ ) {
    212     xcurrent += 1.;
    213     if( x < xcurrent ) break;
    214   }
    215   float ymin=-y0_-0.5;
    216   float ymax=ymin+mSzY_;
    217   if( (y >  ymax) || (y < ymin ) ) return(-1);
    218   float ycurrent=ymin;
    219   for( j=0; j < mSzY_; j++ ) {
    220     ycurrent += 1.;
    221     if( y < ycurrent ) break;
    222   }
    223 
    224 
    225   return (j*mSzX_+i);
    226 }
    227 
    228 //++
    229 void  LocalMap::PixThetaPhi(int_4 k, float& theta, float& phi) const
    230 
     242  AngleProjToPix(theta_aux, phi_aux, x, y);
     243
     244  float xmin= -x0_-0.5;
     245  float xmax= xmin+nSzX_;
     246  if((x >  xmax) || (x < xmin)) return(-1);
     247  float xcurrent= xmin;
     248  for(i = 0; i < nSzX_; i++ )
     249    {
     250      xcurrent += 1.;
     251      if( x < xcurrent ) break;
     252    }
     253  float ymin= -y0_-0.5;
     254  float ymax= ymin+nSzY_;
     255  if((y >  ymax) || (y < ymin)) return(-1);
     256  float ycurrent= ymin;
     257  for(j = 0; j < nSzY_; j++ )
     258    {
     259      ycurrent += 1.;
     260      if( y < ycurrent ) break;
     261    }
     262  return (j*nSzX_+i);
     263}
     264
     265//++
     266template<class T>
     267void LocalMap<T>::PixThetaPhi(int_4 k, float& theta, float& phi) const
     268//
    231269//    Retourne les coordonnées (theta,phi) du milieu du pixel d'indice k
    232270//--
    233271{
    234   if (!(originFlag_==1) || !(SzFlag_==1) ) {
    235     cout << " LocalMap: correspondance carte-sphere non etablie" << endl;
    236     exit(0);
    237   }
     272  if(!(originFlag_) || !(extensFlag_))
     273    {
     274      cout << " LocalMap: correspondance carte-sphere non etablie" << endl;
     275      exit(0);
     276    }
     277
    238278  int i,j;
    239279  Getij(k,i,j);
    240   float X=float(i-x0_);
    241   float Y=float(j-y0_);
     280
     281  float X= float(i-x0_);
     282  float Y= float(j-y0_);
    242283  // situation de ce pixel dans le plan de reference
    243   float x=X*cos_angle_-Y*sin_angle_;
    244   float y=X*sin_angle_+Y* cos_angle_;
     284  float x= X*cos_angle_-Y*sin_angle_;
     285  float y= X*sin_angle_+Y* cos_angle_;
    245286  // projection sur la sphere
    246   PixProjToAngle(x, y,theta, phi);
     287  PixProjToAngle(x, y, theta, phi);
    247288  // retour au plan utilisateur
    248289  ReferenceToUser(theta, phi);
    … …  
    250291 
    251292//++
    252 r_8    LocalMap::PixSolAngle(int_4 k) const
    253 //    Pixel Solid angle  (steradians)
    254 //    All the pixels have not necessarly the same size in (theta, phi)
    255 //    because of the projection scheme which is not yet fixed.
    256 //   
     293template<class T>
     294r_8 LocalMap<T>::PixSolAngle(int_4 k) const
     295//
     296// Pixel Solid angle  (steradians)
     297// All the pixels have not necessarly the same size in (theta, phi)
     298// because of the projection scheme which is not yet fixed.   
    257299//--
    258300{
    259301  int i,j;
    260302  Getij(k,i,j);
    261   float X=float(i-x0_);
    262   float Y=float(j-y0_);
    263   float XR=X+float(i)*0.5;
    264   float XL=X-float(i)*0.5;
    265   float YU=Y+float(j)*0.5;
    266   float YL=Y-float(j)*0.5;
     303  float X= float(i-x0_);
     304  float Y= float(j-y0_);
     305  float XR= X+float(i)*0.5;
     306  float XL= X-float(i)*0.5;
     307  float YU= Y+float(j)*0.5;
     308  float YL= Y-float(j)*0.5;
    267309  // situation  dans le plan de reference
    268   float x0=XL*cos_angle_-YL*sin_angle_;
    269   float y0=XL*sin_angle_+YL* cos_angle_;
    270   float xa=XR*cos_angle_-YL*sin_angle_;
    271   float ya=XR*sin_angle_+YL* cos_angle_;
    272   float xb=XL*cos_angle_-YU*sin_angle_;
    273   float yb=XL*sin_angle_+YU* cos_angle_;
     310  float x0= XL*cos_angle_-YL*sin_angle_;
     311  float y0= XL*sin_angle_+YL*cos_angle_;
     312  float xa= XR*cos_angle_-YL*sin_angle_;
     313  float ya= XR*sin_angle_+YL*cos_angle_;
     314  float xb= XL*cos_angle_-YU*sin_angle_;
     315  float yb= XL*sin_angle_+YU*cos_angle_;
    274316  // projection sur la sphere
    275317  float tet0,phi0,teta,phia,tetb,phib;
    276   PixProjToAngle(x0, y0,tet0, phi0);
    277   PixProjToAngle(xa, ya,teta, phia);
    278   PixProjToAngle(xb, yb,tetb, phib);
     318  PixProjToAngle(x0, y0, tet0, phi0);
     319  PixProjToAngle(xa, ya, teta, phia);
     320  PixProjToAngle(xb, yb, tetb, phib);
    279321  // angle solide
    280   float sol=fabs((xa-x0)*(yb-y0)-(xb-x0)*(ya-y0));
     322  float sol= fabs((xa-x0)*(yb-y0)-(xb-x0)*(ya-y0));
    281323  return r_8(sol);
    282324}
    283325
    284 
    285 //++
    286 void  LocalMap::SetOrigin(float theta0, float phi0, float angle)
    287 
    288 //  
    289 //--
    290 {
    291   theta0_=theta0;
    292   phi0_=phi0;
    293   x0_=mSzX_/2;
    294   y0_=mSzY_/2;
    295   angle=angle*Pi/180.;
    296   cos_angle_=cos(angle);
    297   sin_angle_=sin(angle);
    298   originFlag_ =1;
    299 }
    300 
     326//++
     327template<class T>
     328void  LocalMap<T>::SetOrigin(float theta0, float phi0, float angle)
     329//
     330// Fixe la repere de reference ( angles en degres)
     331//--
     332{
     333  theta0_= (double)theta0;
     334  phi0_  = (double)phi0;
     335  angle_ = (double)angle;
     336  x0_= nSzX_/2;
     337  y0_= nSzY_/2;
     338  cos_angle_= cosdf(angle);
     339  sin_angle_= sindf(angle);
     340  originFlag_= true;
     341  cout << " LocalMap:: set origin 1 done" << endl;
     342}
    301343 
    302344//++
    303 void  LocalMap::SetOrigin(float theta0, float phi0, int_4 x0, int_4 y0, float angle)
    304 
    305 //   
    306 //--
    307 {
    308   theta0_=theta0;
    309   phi0_=phi0;
    310   x0_=x0;
    311   y0_=y0;
    312   angle=angle*Pi/180.;
    313   cos_angle_=cos(angle);
    314   sin_angle_=sin(angle);
    315   originFlag_ =1;
     345template<class T>
     346void  LocalMap<T>::SetOrigin(float theta0, float phi0, int_4 x0, int_4 y0, float angle)
     347//
     348// Fixe le repere de reference (angles en degres) 
     349//--
     350{
     351  theta0_= (double)theta0;
     352  phi0_  = (double)phi0;
     353  angle_ = (double)angle;
     354  x0_= x0;
     355  y0_= y0;
     356  cos_angle_= cosdf(angle);
     357  sin_angle_= sindf(angle);
     358  originFlag_= true;
     359  cout << " LocalMap:: set origin 2 done" << endl;
    316360}
    317361
    318362//++
    319 void  LocalMap::SetSize(float angleX, float angleY)
    320 
    321 //   
    322 //--
    323 {
    324  tgAngleX_=r_8(tan(angleX*Pi/360.));
    325  tgAngleY_=r_8(tan(angleY*Pi/360.));
    326  SzFlag_=1;
    327 }
    328 //++
    329 void  LocalMap::Project(SphericalMap& sphere) const
    330 
    331 //    Projection to spherical map   
    332 //--
    333 {
    334   for (int m=0; m<mNPix_; m++) {
    335     float theta,phi;
    336     PixThetaPhi(m,theta,phi);
    337     sphere(theta,phi)=mPix_[m];
    338     //    cout << "theta " << theta << " phi " << phi << " valeur " << sphere(theta,phi)<< endl;
    339   }
    340 }
    341 
    342 void LocalMap::InitNul()
    343 //
    344 //    initialise à zéro certaines variables internes
    345 {
    346 originFlag_=0;
    347 SzFlag_=0;
    348 cos_angle_=1.;
    349 sin_angle_=0.;
    350 mPix_=NULL;
    351 }
    352 
    353 void LocalMap::Clear()
    354 //
    355 //   
    356 {
    357   if (mPix_)   delete[] mPix_;
    358 }
    359 
    360 void    LocalMap::Getij(int k, int& i, int& j) const
    361 {
    362   i=(k+1)%mSzX_-1;
    363   if (i==-1) i=mSzX_-1;
    364   j=(k-i+2)/mSzX_;
    365 }
    366 
    367 void  LocalMap::ReferenceToUser(float &theta, float &phi) const
    368 
    369 //       
    370 {
    371   if (theta > (r_4)Pi || theta < 0.  || phi<0. || phi >=2* (r_4)Pi ) {
    372     //cout << " LocalMap::ReferenceToUser : exceptions  a mettre en place" <<endl;
    373   //    THROW(out_of_range("LocalMap::PIxVal Pixel index out of range"));
    374 
    375     throw "LocalMap::ReferenceToUser (theta,phi) out of range";
    376   }
     363template<class T>
     364void LocalMap<T>::SetSize(float angleX, float angleY)
     365//
     366// Fixe l'extension de la carte (angles en degres)
     367//--
     368{
     369  angleX_= (double)angleX;
     370  angleY_= (double)angleY;
     371
     372  //tgAngleX_= T(tan(angleX*Pi/360.));
     373  //tgAngleY_= T(tan(angleY*Pi/360.));
     374
     375
     376  // tangente de la moitie de l'ouverture angulaire totale
     377  tgAngleX_= tand(0.5*angleX_);
     378  tgAngleY_= tand(0.5*angleY_);
     379
     380  extensFlag_= true;
     381  cout << " LocalMap:: set extension done" << endl;
     382}
     383
     384//++
     385template<class T>
     386void LocalMap<T>::Project(SphericalMap<T>& sphere) const
     387//
     388// Projection to spherical map   
     389//--
     390{
     391  for(int m = 0; m < nPix_; m++)
     392    {
     393      float theta,phi;
     394      PixThetaPhi(m,theta,phi);
     395      sphere(theta,phi)= pixels_(m);
     396      // cout << "theta " << theta << " phi " << phi << " valeur " << sphere(theta,phi)<< endl;
     397    }
     398}
     399
     400//++
     401template<class T>
     402void LocalMap<T>::Getij(int k, int& i, int& j) const
     403//
     404//--
     405{
     406  i= (k+1)%nSzX_-1;
     407  if(i == -1) i= nSzX_-1;
     408  j= (k-i+2)/nSzX_;
     409}
     410
     411//++
     412template<class T>
     413void  LocalMap<T>::ReferenceToUser(float &theta, float &phi) const
     414//
     415// --     
     416{
     417  if(theta > (r_4)Pi || theta < 0.  || phi < 0. || phi >= 2*(r_4)Pi)
     418    {
     419      //cout << " LocalMap::ReferenceToUser : exceptions  a mettre en place" <<endl;
     420      //    THROW(out_of_range("LocalMap::PIxVal Pixel index out of range"));
     421      throw "LocalMap::ReferenceToUser (theta,phi) out of range";
     422    }
     423
    377424  //cout << " ReferenceToUser entree, t= " << theta << " phi= " << phi << endl;
    378   theta=(r_4)theta0_+theta-(r_4)Pi*0.5;
    379   if (theta<0.) {
    380     theta=-theta;
    381     phi+=(r_4)Pi;
    382   } else
    383     if (theta>(r_4)Pi) {
    384       theta=2.*(r_4)Pi-theta;
    385       phi+=(r_4)Pi;
     425  theta= (r_4)theta0_*Pi/180.+theta-(r_4)Pi*0.5;
     426  if(theta < 0.)
     427    {
     428      theta= -theta;
     429      phi += (r_4)Pi;
    386430    }
    387  
    388   phi=(r_4)phi0_+phi;
    389   while (phi>=2.*(r_4)Pi) phi-=2.*(r_4)Pi;
    390   if (theta > (r_4)Pi || theta < 0.  || phi<0. || phi >=2* (r_4)Pi ) {
    391     cout <<  " LocalMap::ReferenceToUser : erreur bizarre dans le transfert a la carte utilisateur " << endl;
    392     cout << " theta= " << theta << " phi= " << phi << endl;
    393     exit(0);
    394   }
     431  else
     432    {
     433      if(theta > (r_4)Pi)
     434        {
     435          theta= 2.*(r_4)Pi-theta;
     436          phi += (r_4)Pi;
     437        }
     438    }
     439
     440  phi= (r_4)phi0_*Pi/180.+phi;
     441  while(phi >= 2.*(r_4)Pi) phi-=2.*(r_4)Pi;
     442
     443  if(theta > (r_4)Pi || theta < 0.  || phi < 0. || phi >= 2*(r_4)Pi)
     444    {
     445      cout <<  " LocalMap::ReferenceToUser : erreur bizarre dans le transfert a la carte utilisateur " << endl;
     446      cout << " theta= " << theta << " phi= " << phi << endl;
     447      exit(0);
     448    }
    395449  //cout << " ReferenceToUser sortie, t= " << theta << " phi= " << phi << endl;
    396450}
    397451
    398 void  LocalMap::UserToReference(float &theta, float &phi) const
    399 
    400 //       
    401 {
    402   if (theta > (r_4)Pi || theta < 0.  || phi<0. || phi >=2* (r_4)Pi ) {
    403     cout << " LocalMap::UserToReference : exceptions  a mettre en place" <<endl;
    404   //    THROW(out_of_range("LocalMap::PIxVal Pixel index out of range"));
    405 
    406     throw "LocalMap::UserToReference (theta,phi) out of range";
    407   }
    408   float phi1=phi-(r_4)phi0_;
    409   if (phi1<0.) phi1+=2.*(r_4)Pi;
    410 
    411   float theta1=theta-(r_4)theta0_+(r_4)Pi*0.5;
    412   if (theta1<0.) {
    413     theta=-theta1;
    414     phi1+=(r_4)Pi;
    415   } else
    416     if (theta1>(r_4)Pi) {
    417       theta=2.*(r_4)Pi-theta1;
    418       phi1+=(r_4)Pi;
    419     }
    420   while (phi1>=2.*(r_4)Pi) phi1-=2.*(r_4)Pi;
    421   phi=phi1;
    422   if (theta > (r_4)Pi || theta < 0.  || phi<0. || phi >=2* (r_4)Pi ) {
    423     cout <<  " LocalMap::UserToReference : erreur bizarre dans le transfert a la carte de reference " << endl;
    424     cout << " theta= " << theta << " phi= " << phi << endl;
    425     exit(0);
    426   }
    427 }
    428 
    429 void LocalMap::PixProjToAngle(float x, float y,float &theta, float &phi) const {
    430 
    431 // (x,y) representent les coordonnees en unites de pixels d'un point DANS
    432 // LE PLAN DE REFERENCE. On recupere (theta,phi) par rapport au repere "absolu"
    433 // theta=pi/2,phi=0.
    434 
    435 theta=(r_4)Pi*0.5-atan(2*y* tgAngleY_/(float) mSzY_);
    436 phi=atan2(2*x* tgAngleX_,(float) mSzX_);
    437 if( phi< 0. ) phi+= DeuxPi;
    438 }
    439 
    440 void  LocalMap::AngleProjToPix(float theta, float phi, float& x, float& y) const {
    441 
     452//++
     453template<class T>
     454void  LocalMap<T>::UserToReference(float &theta, float &phi) const
     455//
     456// --     
     457{
     458  if(theta > (r_4)Pi || theta < 0.  || phi<0. || phi >= 2*(r_4)Pi )
     459    {
     460      cout << " LocalMap::UserToReference : exceptions  a mettre en place" <<endl;
     461      //    THROW(out_of_range("LocalMap::PIxVal Pixel index out of range"));
     462      throw "LocalMap::UserToReference (theta,phi) out of range";
     463    }
     464
     465  float phi1=phi-(r_4)phi0_*Pi/180.;
     466  if(phi1 < 0.) phi1+=2.*(r_4)Pi;
     467
     468  float theta1= theta-(r_4)theta0_*Pi/180.+(r_4)Pi*0.5;
     469  if(theta1 < 0.)
     470    {
     471      theta= -theta1;
     472      phi1+= (r_4)Pi;
     473    }
     474  else
     475    {
     476      if(theta1 > (r_4)Pi)
     477        {
     478          theta= 2.*(r_4)Pi-theta1;
     479          phi1+= (r_4)Pi;
     480        }
     481    }
     482
     483  while(phi1 >= 2.*(r_4)Pi) phi1-=2.*(r_4)Pi;
     484  phi= phi1;
     485  if(theta > (r_4)Pi || theta < 0.  || phi < 0. || phi >= 2*(r_4)Pi )
     486    {
     487      cout <<  " LocalMap::UserToReference : erreur bizarre dans le transfert a la carte de reference " << endl;
     488      cout << " theta= " << theta << " phi= " << phi << endl;
     489      exit(0);
     490    }
     491}
     492
     493//++
     494template<class T>
     495void LocalMap<T>::PixProjToAngle(float x, float y, float& theta, float& phi) const
     496//
     497// (x,y) representent les coordonnees en unites de pixels d'un point DANS  LE PLAN DE REFERENCE.
     498// On recupere (theta,phi) par rapport au repere "absolu" theta=pi/2 et phi=0.
     499//--
     500{
     501  theta= (r_4)Pi*0.5-atan(2*y*tgAngleY_/(float)nSzY_);
     502  phi= atan2(2*x*tgAngleX_,(float)nSzX_);
     503  if(phi < 0.) phi += DeuxPi;
     504}
     505
     506//++
     507template<class T>
     508void LocalMap<T>::AngleProjToPix(float theta, float phi, float& x, float& y) const
     509//
    442510// (theta,phi) par rapport au repere "absolu" theta=pi/2,phi=0. On recupere
    443511// (i,j)  DANS LE PLAN DE REFERENCE.
    444   if (phi>=(r_4)Pi) phi-= DeuxPi;
     512//--
     513{
     514  if(phi >= (r_4)Pi) phi-= DeuxPi;
    445515  //  y=0.5*mSzY_*cot(theta)/tgAngleY_;  $CHECK-REZA-04/99$
    446   y=0.5*mSzY_/tan(theta)/tgAngleY_;  // ? cot = 1/tan ? 
    447   x=0.5*mSzX_*tan(phi)/tgAngleX_;
    448 }
    449 
    450 
    451 
    452 
    453 
     516  y= 0.5*nSzY_/tan(theta)/tgAngleY_;  // ? cot = 1/tan ? 
     517  x= 0.5*nSzX_*tan(phi)/tgAngleX_;
     518}
     519
     520template<class T>
     521void LocalMap<T>::print(ostream& os) const
     522{
     523  os<<" SzX= "<<nSzX_<<", SzY= "<<nSzY_<<", NPix= "<<nPix_<<endl;
     524  if(LocalMap_isDone())
     525    {
     526      os<<" theta0= "<<theta0_<<", phi0= "<<phi0_<<", angle= "<<angle_<<endl;
     527      os<<" x0= "<<x0_<<", y0= "<<y0_<<endl;
     528      os<<" cos= "<<cos_angle_<<", & sin= "<<sin_angle_<<endl;
     529      os<<" angleX= "<<angleX_<<", angleY= "<<angleY_<<endl;
     530      os<<" tg(angleX)= "<<tgAngleX_<<", tg(angleY)= "<<tgAngleY_<<endl;
     531    }
     532
     533  os << " contenu de pixels : ";
     534  for(int i=0; i < nPix_; i++)
     535    {
     536      if(i%5 == 0) os << endl;
     537      os <<  pixels_(i) <<", ";
     538    }
     539  os << endl;
     540}
     541
     542//*******************************************************************
     543// class FIO_LocalMap<T>
     544//  Les objets delegues pour la gestion de persistance
     545//*******************************************************************
     546
     547template <class T>
     548FIO_LocalMap<T>::FIO_LocalMap()
     549{
     550  dobj= new LocalMap<T>;
     551  ownobj= true;
     552}
     553
     554template <class T>
     555FIO_LocalMap<T>::FIO_LocalMap(string const& filename)
     556{
     557  dobj= new LocalMap<T>;
     558  ownobj= true;
     559  Read(filename);
     560}
     561
     562template <class T>
     563FIO_LocalMap<T>::FIO_LocalMap(const LocalMap<T>& obj)
     564{
     565  dobj= new LocalMap<T>(obj);
     566  ownobj= true;
     567}
     568
     569template <class T>
     570FIO_LocalMap<T>::FIO_LocalMap(LocalMap<T>* obj)
     571{
     572  dobj= obj;
     573  ownobj= false;
     574}
     575
     576template <class T>
     577FIO_LocalMap<T>::~FIO_LocalMap()
     578{
     579  if (ownobj && dobj) delete dobj;
     580}
     581
     582template <class T>
     583AnyDataObj* FIO_LocalMap<T>::DataObj()
     584{
     585  return(dobj);
     586}
     587
     588template <class T>
     589void FIO_LocalMap<T>::ReadSelf(PInPersist& is)
     590{
     591  cout << " FIO_LocalMap:: ReadSelf " << endl;
     592
     593  if(dobj == NULL)
     594    {
     595      dobj= new LocalMap<T>;
     596    }
     597
     598  // Pour savoir s'il y avait un DVList Info associe
     599  char strg[256];
     600  is.GetLine(strg, 255);
     601  bool hadinfo= false;
     602  if(strncmp(strg+strlen(strg)-7, "HasInfo", 7) == 0)  hadinfo= true;
     603  if(hadinfo)
     604    {    // Lecture eventuelle du DVList Info
     605      is >> dobj->Info();
     606    }
     607
     608  int_4 nSzX;
     609  is.GetI4(nSzX);
     610  dobj->setSize_x(nSzX);
     611
     612  int_4 nSzY;
     613  is.GetI4(nSzY);
     614  dobj->setSize_y(nSzY);
     615
     616  int_4 nPix;
     617  is.GetI4(nPix);
     618  dobj->setNbPixels(nPix);
     619
     620  string ss("local mapping is done");
     621  string sso;
     622  is.GetStr(sso);
     623  if(sso == ss)
     624    {
     625      cout<<" ReadSelf:: local mapping"<<endl;
     626      int_4 x0, y0;
     627      float theta, phi, angle;
     628      is.GetI4(x0);
     629      is.GetI4(y0);
     630      is.GetR4(theta);
     631      is.GetR4(phi);
     632      is.GetR4(angle);
     633      dobj->SetOrigin(theta, phi, x0, y0, angle);
     634
     635      float angleX, angleY;
     636      is.GetR4(angleX);
     637      is.GetR4(angleY);
     638      dobj->SetSize(angleX, angleY);
     639    }
     640
     641  T* pixels= new T[nPix];
     642  PIOSReadArray(is, pixels, nPix);
     643  dobj->setDataBlock(pixels, nPix);
     644  delete [] pixels;
     645}
     646
     647template <class T>
     648void FIO_LocalMap<T>::WriteSelf(POutPersist& os) const
     649{
     650  cout << " FIO_LocalMap:: WriteSelf " << endl;
     651
     652  if(dobj == NULL)
     653    {
     654      cout << " WriteSelf:: dobj= null " << endl;
     655      return;
     656    }
     657
     658  char strg[256];
     659  int_4 nSzX= dobj->Size_x();
     660  int_4 nSzY= dobj->Size_y();
     661  int_4 nPix= dobj->NbPixels();
     662 
     663  if(dobj->ptrInfo())
     664    {
     665      sprintf(strg,"LocalMap: NPixX=%6d  NPixY=%9d HasInfo",nSzX,nSzY);
     666      os.PutLine(strg);
     667      os << dobj->Info();
     668    }
     669  else
     670    {
     671      sprintf(strg,"LocalMap: NPixX=%6d  NPixY=%9d ",nSzX,nSzY);
     672      os.PutLine(strg); 
     673    }
     674
     675  os.PutI4(nSzX);
     676  os.PutI4(nSzY);
     677  os.PutI4(nPix);
     678
     679  if(dobj->LocalMap_isDone())
     680    {
     681      string ss("local mapping is done");
     682      os.PutStr(ss);
     683      int_4 x0, y0;
     684      float theta, phi, angle;
     685      dobj->Origin(theta, phi, x0, y0, angle);
     686      os.PutI4(x0);
     687      os.PutI4(y0);
     688      os.PutR4(theta);
     689      os.PutR4(phi);
     690      os.PutR4(angle);
     691
     692      float angleX, angleY;
     693      dobj->Aperture(angleX, angleY);
     694      os.PutR4(angleX);
     695      os.PutR4(angleY);
     696    }
     697  else
     698    {
     699      string ss("no local mapping");
     700      os.PutStr(ss);
     701    }
     702
     703  PIOSWriteArray(os,(dobj->getDataBlock())->Data(), nPix);
     704}
     705
Note: See TracChangeset for help on using the changeset viewer.