source: Sophya/trunk/Poubelle/DPC:FitsIOServer/Samba/localmap.cc@ 1203

#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)
//
//    A local map is a 2 dimensional array, with i as column index and j
//    as row index. The map is supposed to lie on a plan tangent to the 
//    celestial sphere in a point whose coordinates are (x0,y0) on the local
//    map and (theta0, phi0) on the sphere. The range of the map is defined
//    by two values of angles covered respectively by all the pixels in 
//    x direction and all the pixels in y direction (SetSize()).
//
//    A "reference plane" is considered : this plane is tangent to the 
//    celestial sphere in a point with angles theta=Pi/2 and phi=0. This 
//    point is the origine of coordinates is of the reference plane. The 
//    x-axis is the tangent parallel to the equatorial line and oriented 
//    toward the increasing phi's ; the y-axis is parallel to the meridian 
//    line and oriented toward the north pole.
//    
//    Internally, a map is first defined within this reference plane and
//    tranported until the point (theta0, phi0) in such a way that both 
//    axes are kept parallel to meridian and parallel lines of the sphere.
//    The user can define its own map with axes rotated with respect to
//    reference axes (this rotation is characterized by angle between 
//    the local parallel line and the wanted x-axis-- see method 
//    SetOrigin(...))
//    
//

//
//-- 
//++
//
// Links        Parents
//
//    PixelMap
//    
//-- 
//++
//
// Links        Childs
//
//     
//-- 
//++
// Titre        Constructors
//--
//++
template<class T>
LocalMap<T>::LocalMap()
//
//     
//--
{
  InitNul();
  pixels_.Reset();
}

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

//++
template<class T>
LocalMap<T>::LocalMap(const LocalMap<T>& lm, bool share)
  : pixels_(lm.pixels_, share)

//
//    copy constructor
//--
{
  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_;
}

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



//++
// Titre        Public Methods
//--

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

//++
template<class T>
int LocalMap<T>::NbPixels() const
//
//    Return number of pixels
//--
{
  return(nPix_);
}

//++
template<class T>
T& LocalMap<T>::PixVal(int k)
//
//    Return value of pixel with index 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
//
//   const version of previous method
//--
{
  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
{
// $CHECK$  Reza 11/11/99 - A modifier
return(true);
}

//++
template<class T>
int LocalMap<T>::PixIndexSph(double theta,double phi) const
//
//    Return index of the pixel with spherical coordinates (theta,phi)
//--
{
  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
//
//    Return (theta, phi) coordinates of pixel with index 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>
T LocalMap<T>::SetPixels(T v)
{
pixels_.Reset(v);
return(v);
}
 
//++
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)
//
//    set the referential of the map (angles in degrees)
//    (default x0=siz_x/2,  y0=siz_y/2) 
//--
{
  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)
//
//    set the referential of the map (angles in degrees)
//--
{
  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)
//
//    angle range of tthe map (angles in degrees)
//--
{
  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 a 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;
    }
}
// Titre        Private Methods
//++
template<class T>
void LocalMap<T>::InitNul()
//
//    set some attributes to zero
//--
{ 
  originFlag_= false;
  extensFlag_= false;
  cos_angle_= 1.0;
  sin_angle_= 0.0;
//  pixels_.Reset(); Pas de reset par InitNul (en cas de share) - Reza 20/11/99 $CHECK$
}

//++
template<class T>
void LocalMap<T>::Getij(int k,int& i,int& j) const
//
//    Return 2 indices corresponding to the pixel number k
//--
{
  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
//
//    Transform a pair of coordinates (theta, phi) given in 
//    reference coordinates into map coordinates
//--     
{
  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
//
//    Transform a pair of coordinates (theta, phi) given in 
//    map coordinates into reference coordinates
//--     
{
  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
//
//  
//    Given coordinates in pixel units in the REFERENCE PLANE, return
//    (theta, phi) in "absolute" referential theta=pi/2 ,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
//
//    Given coordinates  (theta, phi) in "absolute" referential 
//    theta=pi/2 ,phi=0  return pixel indices  (i,j) in the REFERENCE PLANE.
//--
{
  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;
}
//++
// Titre        class FIO_LocalMap
//    Delegated objects for persitance management
//--

//*******************************************************************
// 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>;
  dobj->DataBlock().SetTemp(true);
  ownobj= true;
  Read(filename);
}

//++
template <class T>
FIO_LocalMap<T>::FIO_LocalMap(const LocalMap<T>& obj)
//
//--
{
  dobj= new LocalMap<T>(obj, true);
  dobj->DataBlock().SetTemp(true);
  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)
//
//--
{

  if(dobj == NULL) 
    {
      dobj= new LocalMap<T>;
      dobj->DataBlock().SetTemp(true);
      ownobj= true;
    }

  // 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);
    }

// On lit le DataBlock;
  is >> dobj->DataBlock();
}

//++
template <class T>
void FIO_LocalMap<T>::WriteSelf(POutPersist& os) const
//
//--
{
  if(dobj == NULL) 
    {
      cout << " FIO_LocalMap::WriteSelf:: dobj= null " << endl;
      return;
    }

  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);
    }

// On ecrit le dataBlock 
  os << dobj->DataBlock();
}