Changeset 564 in Sophya for trunk/SophyaExt

Timestamp:

Nov 10, 1999, 12:14:20 PM (26 years ago)

Author:

ansari

Message:

ajout methode load pour spheregorski

Location:

trunk/SophyaExt/FitsIOServer

Files:

• fitsioserver.cc (modified) (3 diffs)
• fitsioserver.h (modified) (4 diffs)

trunk/SophyaExt/FitsIOServer/fitsioserver.cc

@@ -291 +291 @@
   for (int j = 0; j < sph.NbPixels(); j++) sph(j)= (double)r_8tab_[j];
 }
+
+
+void FitsIoServer::load(SphereGorski<float>& sph, char flnm[], int hdunum) 
+{
+  int npixels=0;
+  int nside=0;
+
+  FITS_tab_typ_ = TFLOAT;
+
+  DVList dvl;
+  planck_read_bntbl(flnm, hdunum, npixels, dvl);
+  //dvl.Print();
+  nside= dvl.GetI("NSIDE");
+  const char* ordering= dvl.GetS("ORDERING").c_str();
+  char* ring = "'RING";
+  if (  strncmp(ordering, ring,3) != 0) 
+
+    //  if (ordering!="RING    ")
+    {
+      cout << " numerotation non RING" << endl;
+    }
+  // number of pixels in the sphere
+  if (sph.NbPixels() != npixels)
+    {
+      cout << " found " << npixels << " pixels" << endl;
+      cout << " expected " << sph.NbPixels() << endl;
+      if (nside <= 0 ) 
+        {
+          cout<<" FITSIOSERVER: no resolution parameter on fits file "<<endl;
+          exit(0);
+        }
+      if (nside != sph.SizeIndex())
+        { 
+          sph.Resize(nside);
+          cout << " resizing the sphere to nside=  " << nside  << endl;
+        }
+      else
+        {
+          cout << " FITSIOSERVER : same resolution, surprising ! " << endl;
+          exit(0);
+        }
+    }
+  for (int j = 0; j < sph.NbPixels(); j++) sph(j)= (float)r_4tab_[j];
+}
+
 void FitsIoServer::load(SphericalMap<float>& sph, char flnm[]) 
 {
@@ -1106 +1151 @@
 
 
+void FitsIoServer::planck_read_bntbl(char flnm[], int hdunum, int& npixels, DVList& dvl)
+{
+  int status=0;
+  int nkeys,keypos;
+  int hdutype;
+  int tfields;
+  int datype;
+  long lastpix;
+  long repeat, width;
+  long nrows;
+  long extend;
+  char* comment=NULL;
+
+  // pointer to the FITS file, defined in fitsio.h
+  fitsfile *fptr;
+  // initialize status before calling fitsio routines
+  fits_open_file(&fptr, flnm, READONLY, &status);
+  if( status )  printerror( status );
+  fits_read_key_lng(fptr, "EXTEND", &extend, comment, &status);
+  if( status )  printerror( status );
+  if (extend!=1) 
+    {
+      cout << "FitsIoServer:: le fichier fits ne contient pas d'extension binary table" << endl;
+      return;
+    }
+  fits_movabs_hdu(fptr, hdunum,&hdutype,&status);
+  if( status )  printerror( status );
+  if (hdutype!=BINARY_TBL) 
+    {
+      cout << "FitsIoServer:: this HDU is not a binary table " << endl;
+      exit(status);
+    }
+   char xtension[FLEN_VALUE];
+   fits_read_key_str(fptr,"XTENSION",xtension,NULL,&status);
+   if( status )  printerror( status );
+   cout << "xtension= " << xtension << endl;
+
+  char binta[] = "BINTABLE";
+  if (  strncmp(xtension, binta,8) != 0) 
+    //   if (xtension !="BINTABLE")
+     {
+      cout << "FitsIoServer::  not a binary table " << endl;
+      exit(status);
+     }
+   fits_get_hdrpos(fptr,&nkeys,&keypos,&status);
+   if( status )  printerror( status );
+   //cout << " nombre de mots-cles : " << nkeys << endl;
+  fits_get_num_cols(fptr, &tfields, &status);
+  cout << "tfields= " << tfields << endl;
+  if (tfields != 1)
+    {
+      cout << "FitsIoServer:: il y a plus d'une colonne" << endl;
+      return;
+    }
+  fits_get_num_rows(fptr, &nrows, &status);
+  //cout << "nblignes= " << nrows << endl;
+  fits_get_coltype(fptr, 1, &datype, &repeat, &width, &status);
+  if( status )  printerror( status );
+  //cout << " type de donnees : " << datype << endl;
+  //cout << " repeat : " << repeat << endl;
+  //cout << " width : " << width << endl;
+  fits_read_key_lng(fptr, "LASTPIX", &lastpix, comment, &status);
+  if( status )  printerror( status );
+
+  long nelements= nrows*repeat;
+  if (nelements!=lastpix+1) 
+    {
+      cout << " erreur sur longueur du vecteur " << endl;
+      cout << " nelements= " << nelements << " lastpix+1=" << lastpix+1 << endl;
+    }
+  npixels=nelements;
+  int anynull;
+  r_8 dnullval=0.;
+  r_4 fnullval=0.;
+  int_4 inullval=0;
+  // on laisse a fits le soin de convertir le type du tableau lu vers
+  // le type suppose par l'utilisateur de fitsioserver 
+  // 
+  switch ( FITS_tab_typ_) 
+    {
+    case TDOUBLE :
+      if (datype !=  TDOUBLE) 
+        {
+          cout << " FitsIOServer : the data type on fits file " << flnm << " is not double, "
+               << " conversion to double will be achieved by cfitsio lib " << endl; 
+        }
+      if (r_8tab_ != NULL) delete [] r_8tab_;
+      r_8tab_=new r_8[ npixels];
+      fits_read_col(fptr, TDOUBLE, 1, 1, 1, nelements,  &dnullval, 
+                    r_8tab_, 
+                    &anynull, &status);
+      if( status )  printerror( status );
+      break;
+    case TFLOAT :
+      if (datype !=  TFLOAT) 
+        {
+
+          cout << " FitsIOServer : the data type on fits file " << flnm << " is not float, "
+               << " conversion to float will be achieved by cfitsio lib " << endl; 
+        }
+      if (r_4tab_ != NULL) delete [] r_4tab_;
+      r_4tab_=new r_4[nelements];
+      fits_read_col(fptr, TFLOAT, 1, 1, 1, nelements,  &fnullval, 
+                    r_4tab_, &anynull, &status);
+      if( status )  printerror( status );
+      break;
+
+
+    case TINT :
+      if (datype !=  TLONG) 
+        {
+          cout << " FitsIOServer : the data type on fits file " << flnm << " is not long, "
+               << " conversion to long will be achieved by cfitsio lib " << endl; 
+        }
+      if (i_4tab_ != NULL) delete [] i_4tab_;
+      i_4tab_=new int_4[nelements];
+      fits_read_col(fptr, TLONG, 1, 1, 1, nelements,  &inullval, 
+                    i_4tab_, &anynull, &status);
+      //fits_read_img(fptr, TINT, 1, nelements,  &inullval, i_4tab_, 
+      //                    &anynull, &status);
+  //if( status )  printerror( status );
+      break;
+
+
+    default :
+      cout << " FitsIOServer::read_img : type non traite: " << FITS_tab_typ_  << endl;
+      break;
+    }
+  char card[FLEN_CARD];  
+  char keyname[LEN_KEYWORD]= "";
+  char strval[FLEN_VALUE];
+  char comment1[FLEN_COMMENT];
+  char dtype;
+  //char bidon[LEN_KEYWORD];
+  char comkey[] = "COMMENT";
+
+   for(int j = 1; j <= nkeys; j++)
+     {
+       //       fits_read_record(fptr, j, card, &status);
+       // strncpy(keyname,card,LEN_KEYWORD-1);
+       // cout << " bidon= " << keyname << endl;
+       // if ( strncmp(keyname,comkey ,LEN_KEYWORD-1) != 0)
+       fits_read_keyn(fptr,j,card,strval,comment1,&status);
+       strncpy(keyname,card,LEN_KEYWORD-1);
+
+        if ( strncmp(keyname,comkey ,LEN_KEYWORD-1) != 0)
+          {
+            fits_get_keytype(strval,&dtype,&status);
+            strip (keyname, 'B',' '); 
+            strip(strval, 'B',' ');
+            switch( dtype ) 
+              {
+              case 'C':
+                dvl[keyname]= strval;     
+                break;
+              case 'I':
+                int ival;
+                ctoi(strval,&ival);
+                dvl[keyname]= (int_4)ival;
+                break;
+              case 'F':
+                double dval;
+                ctof(strval,&dval);
+                dvl[keyname]= dval;
+                break;
+              default :
+                cout << " FitsIOServer : mot-cle bizerre " << endl;
+                break;
+              }
+          }
+     }
+
+   // close the file
+   status=0;
+   fits_close_file(fptr, &status);
+   if( status )  printerror( status );
+}
+
 
 // projects a SphericalMap<double>, according sinus-method, and saves onto 
@@ -1180 +1403 @@
 }
 
+// projects a SphericalMap<double>, according Mollweide-method, and saves onto 
+// a FITS-file
+void FitsIoServer::Mollweide_picture_projection(SphericalMap<float>& sph, char filename[])
+{
+  // le code de cete methode duplique celui de la precedente, la seule 
+  //difference etant le type de sphere en entree. Ces methodes de projection
+  // sont provisoires, et ne servent que pour les tests. C est pourquoi je 
+  // ne me suis pas casse la tete, pour l instant
+    
+  long naxes[2]={600, 300};
+  float* map = new float[ 600*300 ];
+  int npixels=  naxes[0]*naxes[1];
+
+   cout << " image FITS en projection MOLLWEIDE" << endl;
+  // table will have npixels rows
+   for(int j=0; j < npixels; j++) map[j]=0.;
+   for(int j=0; j<naxes[1]; j++) 
+     {
+       double yd = (j+0.5)/naxes[1]-0.5;
+       double facteur=2.*Pi/sin(acos(yd*2));
+       double theta = (0.5-yd)*Pi;
+       for(int i=0; i<naxes[0]; i++)  
+         {
+           double xa = (i+0.5)/naxes[0]-0.5;
+           double phi =  xa*facteur+Pi;
+           float th=float(theta);
+           float ff=float(phi);
+           if (phi<2*Pi && phi>= 0) 
+             {
+               map[j*naxes[0]+i] = sph.PixValSph(th, ff);
+             }
+         }
+     }
+   
+   write_picture(naxes, map,  filename);
+   delete [] map;
+   
+}
+
+
+
 // saves a (LocalMap<double> on a FITS-file in order to be visualized
 // (for tests)

trunk/SophyaExt/FitsIOServer/fitsioserver.h

@@ -5 +5 @@
 #include "tvector.h"
 #include "sphericalmap.h"
+#include "spheregorski.h"
 #include "localmap.h"
 #include "ntuple.h"
@@ -33 +34 @@
  void load(ImageR4& DpcImg,char flnm[]);
  void load(ImageI4& DpcImg,char flnm[]);
+ void load(SphereGorski<float>& sph, char flnm[], int hdunum); 
  void save(TMatrix<double>& mat, char flnm[]);
  void save(NTuple& ntpl,char flnm[]);
@@ -42 +44 @@
  void sinus_picture_projection(SphericalMap<double>& sph, char flnm[]);
  void sinus_picture_projection(SphericalMap<float>& sph, char flnm[]);
+ void Mollweide_picture_projection(SphericalMap<float>& sph, char flnm[]);
  void picture(LocalMap<double>& lcm, char flnm[]);
 
@@ -56 +59 @@
 
  void planck_read_img(char flnm[],long& naxis,int& n1,int& n2,int& n3,DVList& dvl);
+ void planck_read_bntbl(char flnm[], int hdunum, int& npixels, DVList& dvl);
  void write_picture(long* naxes,float* map,char* filename) const;