source: Sophya/trunk/SophyaLib/Samba/localmap.cc@ 518

#include "localmap.h"
#include "nbmath.h"
#include <complex>
#include "piocmplx.h"

#include <string.h>
#include <iostream.h>

#ifdef __CXX_PRAGMA_TEMPLATES__
#pragma define_template LocalMap<double>
#pragma define_template LocalMap<float>
#pragma define_template LocalMap< complex<double> >
#pragma define_template LocalMap< complex<float> >
#pragma define_template FIO_LocalMap<double>
#pragma define_template FIO_LocalMap<float>
#pragma define_template FIO_LocalMap< complex<float> >
#pragma define_template FIO_LocalMap< complex<double> >
#endif
#if defined(ANSI_TEMPLATES) || defined(GNU_TEMPLATES)
template class LocalMap<double>;
template class LocalMap<float>;
template class LocalMap< complex<double> >;
template class LocalMap< complex<float> >;
template class FIO_LocalMap<double>;
template class FIO_LocalMap<float>;
template class FIO_LocalMap< complex<float> >;
template class FIO_LocalMap< complex<double> >;
#endif

//*****************************************************************************
//++ 
// Class        LocalMap
//
// include      localmap.h nbmath.h
//
//    A local map of a region of the sky, in cartesian coordinates.
//    It has an origin in (theta0, phi0), mapped to pixel(x0, y0)
//    (x0, y0 might be outside of this local map)
//    default value of (x0, y0) is middle of the map, center of 
//    pixel(nx/2, ny/2)
//
//    la carte est consideree comme un tableau a deux indices i et j, i etant 
//    indice de colonne et j indice de ligne. La carte est supposee resider 
//    dans un plan tangent, dont le point de tangence est repere (x0,y0) dans 
//    la carte et (theta0, phi0) sur la sphere celeste. L'extension de la 
//    carte est definie par les valeurs de deux angles couverts respectivement
//    par la totalite des pixels en x de la carte et la totalite des pixels 
//    en y. (SetSize()).
//    On considere un "plan de reference" : plan tangent a la sphere celeste 
//    aux angles theta=Pi/2 et phi=0. Dans ce plan L'origine des coordonnees 
//    est le point de tangence. L'axe Ox est la tangente parallele a 
//    l'equateur, dirige vers les phi croissants, l'axe Oy est parallele 
//    au meridien, dirige vers le pole nord.
//    De maniere interne a la classe une carte est definie dans ce plan de 
//    reference et transportee  jusqu'au point (theta0, phi0) de sorte que 
//    les axes restent paralleles aux meridiens et paralleles. L'utilisateur 
//    peut definir sa carte selon un repere en rotation par rapport au repere 
//    de reference (par l'angle entre le parallele et l'axe Ox souhaite -- 
//    methode SetOrigin(...))

//
//-- 
//++
//
// Links        Parents
//
//    PixelMap
//    
//-- 
//++
//
// Links        Descendants
//
//     
//-- 
//++
// Titre        Constructeurs
//--
//++
template<class T>
LocalMap<T>::LocalMap()
//
// Constructeur par defaut 
//--
{
  InitNul();
}

//++
template<class T>
LocalMap<T>::LocalMap(int nx, int ny) : nSzX_(nx), nSzY_(ny)
//
// Constructeur
//--
{
  InitNul();
  nPix_= nx*ny;
  pixels_.ReSize(nPix_);
  pixels_.Reset();
}

//++
template<class T>
LocalMap<T>::LocalMap(const LocalMap<T>& lm)
//
// Constructeur de recopie
//--
{
  cout<<" LocalMap:: Appel du constructeur de recopie " << endl;

  if(lm.mInfo_) mInfo_= new DVList(*lm.mInfo_);
  nSzX_= lm.nSzX_;
  nSzY_= lm.nSzY_;
  nPix_= lm.nPix_;
  originFlag_= lm.originFlag_;
  extensFlag_= lm.extensFlag_;
  x0_= lm.x0_;
  y0_= lm.y0_;
  theta0_= lm.theta0_;
  phi0_= lm.phi0_;
  angle_= lm.angle_;
  cos_angle_= lm.cos_angle_;
  sin_angle_= lm.sin_angle_;
  angleX_= lm.angleX_;
  angleY_= lm.angleY_;
  tgAngleX_= lm.tgAngleX_;
  tgAngleY_= lm.tgAngleY_;
  pixels_= lm.pixels_;
}

//++
template<class T>
LocalMap<T>::~LocalMap()
//
// Destructeur
//--
{
  InitNul();
}

//++
template<class T>
void LocalMap<T>::InitNul()
//
// Initialise à zéro certaines variables internes
//--
{ 
  originFlag_= false;
  extensFlag_= false;
  cos_angle_= 1.0;
  sin_angle_= 0.0;
  pixels_.Reset();
}

//++
template<class T>
void LocalMap<T>::ReSize(int nx, int ny) 
//
// Redimensionne l'espace de stokage des pixels
//--
{
  InitNul();
  nSzX_ = nx;
  nSzY_ = ny;
  nPix_= nx*ny;
  pixels_.ReSize(nPix_);
  pixels_.Reset();
}

//++
// Titre        Methodes
//--


//++
template<class T>
int LocalMap<T>::NbPixels() const
//
//    Retourne le nombre de pixels du découpage
//--
{
  return(nPix_);
}

//++
template<class T>
T& LocalMap<T>::PixVal(int k)
//
//    Retourne la valeur du contenu du pixel d'indice k
//--
{
  if((k < 0) || (k >= nPix_)) 
    {
      cout << " LocalMap::PIxVal : exceptions  a mettre en place" <<endl;
      //  THROW(out_of_range("LocalMap::PIxVal Pixel index out of range"));
      THROW(rangeCheckErr);
      //throw "LocalMap::PIxVal Pixel index out of range";
    }
  return(pixels_(k));
}

//++

template<class T>
T const& LocalMap<T>::PixVal(int k) const
//
//    Retourne la valeur du contenu du pixel d'indice k
//--
{
  if((k < 0) || (k >= nPix_)) 
    {
      cout << " LocalMap::PIxVal : exceptions  a mettre en place" <<endl;
      //    THROW(out_of_range("LocalMap::PIxVal Pixel index out of range"));
      
      throw "LocalMap::PIxVal Pixel index out of range"; 
    } 
  return *(pixels_.Data()+k);
}

//++
template<class T>
bool LocalMap<T>::ContainsSph(double theta, double phi) const
//--
{
  return(true);  // $CHECK$  A MODIFIER - Reza 26/10/99
}

//++
template<class T>
int LocalMap<T>::PixIndexSph(double theta,double phi) const
//
//    Retourne l'indice du pixel vers lequel pointe une direction définie par 
//    ses coordonnées sphériques
//--
{
  int i,j;
  if(!(originFlag_) || !(extensFlag_)) 
    {
      cout << " LocalMap: correspondance carte-sphere non etablie" << endl;
      exit(0);
    }

  // theta et phi en coordonnees relatives (on se ramene a une situation par rapport au plan de reference)
  double theta_aux= theta;
  double phi_aux  = phi;
  UserToReference(theta_aux, phi_aux);

  // coordonnees dans le plan local en unites de pixels
  double x,y;
  AngleProjToPix(theta_aux,phi_aux, x, y);

  double xmin= -x0_-0.5;
  double xmax= xmin+nSzX_;
  if((x > xmax) || (x < xmin)) return(-1);
  double xcurrent= xmin;
  for(i = 0; i < nSzX_; i++ ) 
    {
      xcurrent += 1.;
      if( x < xcurrent ) break;
    }
  double ymin= -y0_-0.5;
  double ymax= ymin+nSzY_;
  if((y >  ymax) || (y < ymin)) return(-1);
  double ycurrent= ymin;
  for(j = 0; j < nSzY_; j++ ) 
    {
      ycurrent += 1.;
      if( y < ycurrent ) break;
    }
  return (j*nSzX_+i);
}

//++
template<class T>
void LocalMap<T>::PixThetaPhi(int k,double& theta,double& phi) const
//
//    Retourne les coordonnées (theta,phi) du milieu du pixel d'indice k
//--
{
  if(!(originFlag_) || !(extensFlag_)) 
    {
      cout << " LocalMap: correspondance carte-sphere non etablie" << endl;
      exit(0);
    }

  int i,j;
  Getij(k,i,j);

  double X= double(i-x0_);
  double Y= double(j-y0_);
  // situation de ce pixel dans le plan de reference
  double x= X*cos_angle_-Y*sin_angle_;
  double y= X*sin_angle_+Y* cos_angle_;
  // projection sur la sphere
  PixProjToAngle(x, y, theta, phi);
  // retour au plan utilisateur
  ReferenceToUser(theta, phi);
}
 
//++
template<class T>
double LocalMap<T>::PixSolAngle(int k) const
//
// Pixel Solid angle  (steradians)
// All the pixels have not necessarly the same size in (theta, phi)
// because of the projection scheme which is not yet fixed.   
//--
{
  int i,j;
  Getij(k,i,j);
  double X= double(i-x0_);
  double Y= double(j-y0_);
  double XR= X+double(i)*0.5;
  double XL= X-double(i)*0.5;
  double YU= Y+double(j)*0.5;
  double YL= Y-double(j)*0.5;

  // situation  dans le plan de reference
  double x0= XL*cos_angle_-YL*sin_angle_;
  double y0= XL*sin_angle_+YL*cos_angle_;
  double xa= XR*cos_angle_-YL*sin_angle_;
  double ya= XR*sin_angle_+YL*cos_angle_;
  double xb= XL*cos_angle_-YU*sin_angle_;
  double yb= XL*sin_angle_+YU*cos_angle_;

  // projection sur la sphere
  double thet0,phi0,theta,phia,thetb,phib;
  PixProjToAngle(x0, y0, thet0, phi0);
  PixProjToAngle(xa, ya, theta, phia);
  PixProjToAngle(xb, yb, thetb, phib);

  // angle solide
  double sol= fabs((xa-x0)*(yb-y0)-(xb-x0)*(ya-y0));
  return sol;
}

//++
template<class T>
void LocalMap<T>::SetOrigin(double theta0,double phi0,double angle)
//
// Fixe la repere de reference ( angles en degres) 
//--
{
  theta0_= theta0;
  phi0_  = phi0;
  angle_ = angle;
  x0_= nSzX_/2;
  y0_= nSzY_/2;
  cos_angle_= cos(angle*Pi/180.);
  sin_angle_= sin(angle*Pi/180.);
  originFlag_= true; 
  cout << " LocalMap:: set origin 1 done" << endl;
}
 
//++
template<class T>
void LocalMap<T>::SetOrigin(double theta0,double phi0,int x0,int y0,double angle)
//
// Fixe le repere de reference (angles en degres)  
//--
{
  theta0_= theta0;
  phi0_  = phi0;
  angle_ = angle;
  x0_= x0;
  y0_= y0;
  cos_angle_= cos(angle*Pi/180.);
  sin_angle_= sin(angle*Pi/180.);
  originFlag_= true;
  cout << " LocalMap:: set origin 2 done" << endl;
} 

//++
template<class T>
void LocalMap<T>::SetSize(double angleX,double angleY)
//
// Fixe l'extension de la carte (angles en degres)
//--
{
  angleX_= angleX;
  angleY_= angleY;

  // tangente de la moitie de l'ouverture angulaire totale
  tgAngleX_= tan(0.5*angleX_*Pi/180.);
  tgAngleY_= tan(0.5*angleY_*Pi/180.);

  extensFlag_= true;
  cout << " LocalMap:: set extension done" << endl;
}

//++
template<class T>
void LocalMap<T>::Project(SphericalMap<T>& sphere) const
//
// Projection to spherical map   
//--
{
  for(int m = 0; m < nPix_; m++) 
    {
      double theta,phi;
      PixThetaPhi(m,theta,phi);
      sphere(theta,phi)= pixels_(m);
      // cout << "theta " << theta << " phi " << phi << " valeur " << sphere(theta,phi)<< endl;
    }
}

//++
template<class T>
void LocalMap<T>::Getij(int k,int& i,int& j) const
//
//--
{
  i= (k+1)%nSzX_-1;
  if(i == -1) i= nSzX_-1;
  j= (k-i+2)/nSzX_;
}

//++
template<class T>
void  LocalMap<T>::ReferenceToUser(double& theta,double& phi) const
//
// --     
{
  if(theta > Pi || theta < 0. || phi < 0. || phi >= 2*Pi) 
    {
      throw "LocalMap::ReferenceToUser (theta,phi) out of range"; 
    } 

  theta= theta0_*Pi/180.+theta-Pi*0.5;
  if(theta < 0.) 
    {
      theta= -theta;
      phi += Pi;
    } 
  else 
    {
      if(theta > Pi) 
        {
          theta= 2.*Pi-theta;
          phi += Pi;
        } 
    }

  phi= phi0_*Pi/180.+phi;
  while(phi >= 2.*Pi) phi-= 2.*Pi;

  if(theta > Pi || theta < 0. || phi < 0. || phi >= 2*Pi) 
    {
      cout <<  " LocalMap::ReferenceToUser : erreur bizarre dans le transfert a la carte utilisateur " << endl;
      cout << " theta= " << theta << " phi= " << phi << endl;
      exit(0);
    }
}

//++
template<class T>
void  LocalMap<T>::UserToReference(double& theta,double& phi) const
//
// --     
{
  if(theta > Pi || theta < 0. || phi < 0. || phi >= 2*Pi) 
    {
      cout<<" LocalMap::UserToReference: exceptions a mettre en place" <<endl;
      // THROW(out_of_range("LocalMap::PIxVal Pixel index out of range"));
      throw "LocalMap::UserToReference (theta,phi) out of range"; 
    }

  double phi1= phi-phi0_*Pi/180.;
  if(phi1 < 0.) phi1+= 2.*Pi;

  double theta1= theta-theta0_*Pi/180.+Pi*0.5;
  if(theta1 < 0.) 
    {
      theta= -theta1;
      phi1+= Pi;
    } 
  else
    { 
      if(theta1 > Pi) 
        {
          theta= 2.*Pi-theta1;
          phi1+= Pi;
        }
    }

  while(phi1 >= 2.*Pi) phi1-= 2.*Pi;
  phi= phi1;
  if(theta > Pi || theta < 0. || phi < 0. || phi >= 2*Pi) 
    {
      cout <<  " LocalMap::UserToReference : erreur bizarre dans le transfert a la carte de reference " << endl;
      cout << " theta= " << theta << " phi= " << phi << endl;
      exit(0);
    }
}

//++
template<class T>
void LocalMap<T>::PixProjToAngle(double x,double y,double& theta,double& phi) const
//
// (x,y) representent les coordonnees en unites de pixels d'un point DANS LE PLAN DE REFERENCE. 
// On recupere (theta,phi) par rapport au repere "absolu" theta=pi/2 et phi=0.
//--
{
  theta= Pi*0.5-atan(2.*y*tgAngleY_/(double)nSzY_);
  phi= atan2(2.*x*tgAngleX_,(double)nSzX_);
  if(phi < 0.) phi += DeuxPi;
}

//++
template<class T>
void LocalMap<T>::AngleProjToPix(double theta,double phi,double& x,double& y) const
//
// (theta,phi) par rapport au repere "absolu" theta=pi/2,phi=0. On recupere 
// (i,j)  DANS LE PLAN DE REFERENCE. 
//--
{
  if(phi >= Pi) phi-= DeuxPi;
  //  y=0.5*mSzY_*cot(theta)/tgAngleY_;  $CHECK-REZA-04/99$
  y= 0.5*nSzY_/tan(theta)/tgAngleY_;  // ? cot = 1/tan ?  
  x= 0.5*nSzX_*tan(phi)/tgAngleX_;
}

template<class T>
void LocalMap<T>::print(ostream& os) const
{
  os<<" SzX= "<<nSzX_<<", SzY= "<<nSzY_<<", NPix= "<<nPix_<<endl;
  if(LocalMap_isDone())
    {
      os<<" theta0= "<<theta0_<<", phi0= "<<phi0_<<", angle= "<<angle_<<endl;
      os<<" x0= "<<x0_<<", y0= "<<y0_<<endl;
      os<<" cos= "<<cos_angle_<<", & sin= "<<sin_angle_<<endl;
      os<<" angleX= "<<angleX_<<", angleY= "<<angleY_<<endl;
      os<<" tg(angleX)= "<<tgAngleX_<<", tg(angleY)= "<<tgAngleY_<<endl;
    }

  os << " contenu de pixels : ";
  for(int i=0; i < nPix_; i++)
    {
      if(i%5 == 0) os << endl;
      os <<  pixels_(i) <<", ";
    }
  os << endl;
}

//*******************************************************************
// class FIO_LocalMap<T>
//  Les objets delegues pour la gestion de persistance 
//*******************************************************************

template <class T>
FIO_LocalMap<T>::FIO_LocalMap()
{
  dobj= new LocalMap<T>;
  ownobj= true;
}

template <class T>
FIO_LocalMap<T>::FIO_LocalMap(string const& filename)
{
  dobj= new LocalMap<T>;
  ownobj= true;
  Read(filename);
}

template <class T>
FIO_LocalMap<T>::FIO_LocalMap(const LocalMap<T>& obj)
{
  dobj= new LocalMap<T>(obj);
  ownobj= true;
}

template <class T>
FIO_LocalMap<T>::FIO_LocalMap(LocalMap<T>* obj)
{
  dobj= obj;
  ownobj= false;
}

template <class T>
FIO_LocalMap<T>::~FIO_LocalMap()
{
  if (ownobj && dobj) delete dobj;
}

template <class T>
AnyDataObj* FIO_LocalMap<T>::DataObj()
{
  return(dobj);
}

template <class T>
void FIO_LocalMap<T>::ReadSelf(PInPersist& is)
{
  cout << " FIO_LocalMap:: ReadSelf " << endl;

  if(dobj == NULL) 
    {
      dobj= new LocalMap<T>;
    }

// Let's Read the SphereCoordSys object  -- ATTENTIOn - $CHECK$
  SphereCoordSys* cs = dynamic_cast<SphereCoordSys*>(is.ReadObject());
  dobj->SetCoordSys(cs);

  // Pour savoir s'il y avait un DVList Info associe
  char strg[256];
  is.GetLine(strg, 255);
  bool hadinfo= false;
  if(strncmp(strg+strlen(strg)-7, "HasInfo", 7) == 0)  hadinfo= true;
  if(hadinfo) 
    {    // Lecture eventuelle du DVList Info
      is >> dobj->Info(); 
    }

  int nSzX;
  is.GetI4(nSzX);
  dobj->setSize_x(nSzX);

  int nSzY;
  is.GetI4(nSzY);
  dobj->setSize_y(nSzY);

  int nPix;
  is.GetI4(nPix);
  dobj->setNbPixels(nPix);

  string ss("local mapping is done");
  string sso;
  is.GetStr(sso);
  if(sso == ss)
    {
      cout<<" ReadSelf:: local mapping"<<endl;
      int x0, y0;
      double theta, phi, angle;
      is.GetI4(x0);
      is.GetI4(y0);
      is.GetR8(theta);
      is.GetR8(phi);
      is.GetR8(angle);
      dobj->SetOrigin(theta, phi, x0, y0, angle);

      double angleX, angleY;
      is.GetR8(angleX);
      is.GetR8(angleY);
      dobj->SetSize(angleX, angleY);
    }

  T* pixels= new T[nPix];
  PIOSReadArray(is, pixels, nPix);
  dobj->setDataBlock(pixels, nPix);
  delete [] pixels;
}

template <class T>
void FIO_LocalMap<T>::WriteSelf(POutPersist& os) const
{
  cout << " FIO_LocalMap:: WriteSelf " << endl;

  if(dobj == NULL) 
    {
      cout << " WriteSelf:: dobj= null " << endl;
      return;
    }

// Let's write the SphereCoordSys object
  dobj->GetCoordSys()->Write(os);
  
  char strg[256];
  int nSzX= dobj->Size_x();
  int nSzY= dobj->Size_y();
  int nPix= dobj->NbPixels();
 
  if(dobj->ptrInfo()) 
    {
      sprintf(strg,"LocalMap: NPixX=%6d  NPixY=%9d HasInfo",nSzX,nSzY);
      os.PutLine(strg);
      os << dobj->Info();
    }
  else 
    { 
      sprintf(strg,"LocalMap: NPixX=%6d  NPixY=%9d ",nSzX,nSzY);
      os.PutLine(strg);  
    }

  os.PutI4(nSzX);
  os.PutI4(nSzY);
  os.PutI4(nPix);

  if(dobj->LocalMap_isDone())
    {
      string ss("local mapping is done");
      os.PutStr(ss);
      int x0, y0;
      double theta, phi, angle;
      dobj->Origin(theta, phi, x0, y0, angle);
      os.PutI4(x0);
      os.PutI4(y0);
      os.PutR8(theta);
      os.PutR8(phi);
      os.PutR8(angle);

      double angleX, angleY;
      dobj->Aperture(angleX, angleY);
      os.PutR8(angleX);
      os.PutR8(angleY);
    }
  else
    {
      string ss("no local mapping");
      os.PutStr(ss);
    }

  PIOSWriteArray(os,(dobj->getDataBlock())->Data(), nPix);
}