Changeset 689 in Sophya for trunk

Timestamp:

Dec 22, 1999, 5:57:10 PM (26 years ago)

Author:

ansari

Message:

add loadobj method

Location:

trunk/SophyaExt/FitsIOServer

Files:

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

trunk/SophyaExt/FitsIOServer/fitsioserver.cc

@@ -12 +12 @@
 #include "machdefs.h" 
 #include <iostream.h>
+
+#ifdef __USE_STD_IOSTREAM
+#include <sstream>
+#endif
+
 #include <list>
 #include <string>
@@ -19 +24 @@
 #include "strutil.h"
 
+AnyDataObj* FitsIoServer::loadobj(char flnm[],int hdunum)
+{
+  // pointer to the FITS file, defined in fitsio.h
+  fitsfile *fptr;
+
+  // initialize status before calling fitsio routines
+  int status= 0;
+  fits_open_file(&fptr,flnm,READONLY,&status);
+  if( status ) printerror( status );
+
+  // move to the specified HDU number
+  int hdutype;
+  fits_movabs_hdu(fptr,hdunum,&hdutype,&status);
+  if( status ) printerror( status );
+
+  if(hdutype == BINARY_TBL)
+    {
+      cout << " Reading a FITS binary table in HDU : " << hdunum << endl;
+
+      // get number of keywords
+      int nkeys = 0;
+      int keypos= 0;
+      fits_get_hdrpos(fptr,&nkeys,&keypos,&status);
+      if( status ) printerror( status );
+
+      // get number of fields (columns) in the table
+      int tfields= 0;
+      fits_get_num_cols(fptr,&tfields,&status);
+      if( status ) printerror( status );
+
+      // only a table with ONE column is allowed
+      if(tfields != 1)
+        throw IOExc("FitsIOServer::dataobj  ERROR more than one column !");
+
+      // get the datatype of the values
+      int DTYPE;
+      long repeat,width;
+      fits_get_coltype(fptr,1,&DTYPE,&repeat,&width,&status);
+      if( status ) printerror( status );
+
+      // check if the keyword ORDERING exists
+      bool ordering= check_keyword(fptr,nkeys,"ORDERING");
+
+      // check if the keyword NSIDE exists
+      int nside= 2;
+      if(check_keyword(fptr,nkeys,"NSIDE"))
+        {
+          fits_read_key(fptr,TINT,"NSIDE",&nside,NULL,&status);
+          if( status ) printerror( status );
+        }
+
+      // close the file
+      fits_close_file(fptr,&status);
+      if( status )  printerror( status );
+
+      if(ordering)
+        {
+          if(DTYPE == TDOUBLE)
+            {
+              SphereGorski<double> *sph= new SphereGorski<double>(nside);
+              load(*sph,flnm,hdunum);
+              return sph;
+            }
+      
+          else if(DTYPE == TFLOAT)
+            {
+              SphereGorski<float> *sph= new SphereGorski<float>(nside);
+              load(*sph,flnm,hdunum);
+              return sph;
+            }
+
+          else
+            {
+              cout << " FitsIOServer::dataobj:: DTYPE= " << DTYPE << endl;
+              throw IOExc("datatype code not yet programmed");
+            }
+        }
+      else
+        {
+          NTuple *ntpl= new NTuple;
+          load(*ntpl,flnm,hdunum);
+          return ntpl;
+        }
+    }
+  else if(hdutype == IMAGE_HDU)
+    {
+      cout << " Reading a FITS image in HDU : " << hdunum << endl;
+
+      // bits per pixels
+      int bitpix;
+      fits_read_key(fptr,TINT,"BITPIX",&bitpix,NULL,&status);
+      if( status )  printerror( status );
+
+      // number of dimensions in the FITS array
+      int naxis= 0;
+      fits_read_key(fptr,TINT,"NAXIS",&naxis,NULL,&status);
+      if( status ) printerror( status );
+
+      // read the NAXIS1 and NAXIS2 keyword to get image size
+      long naxes[3]= {0,0,0};
+      int nfound;
+      fits_read_keys_lng(fptr,"NAXIS",1,naxis,naxes,&nfound,&status);
+      if( status ) printerror( status );
+
+      int nrows= (int)naxes[0];
+      int ncols= 0;
+
+      if(naxis == 1) ncols= 1;
+      if(naxis == 2) ncols= (int)naxes[1];
+      if(naxis == 3 && naxes[2] < 2)
+        {
+          naxis= 2;
+          ncols= (int)naxes[1];
+        }
+
+      // close the file
+      fits_close_file(fptr,&status);
+      if( status )  printerror( status );
+      
+      if(bitpix ==  DOUBLE_IMG)
+        {
+          TMatrix<double> *mat=new TMatrix<double>(nrows,ncols);
+          load(*mat,flnm);
+          if(naxis == 1)
+            {
+              TVector<double> *vect=new TVector<double>(*mat);
+              delete mat;
+              return vect;
+            }
+          else if(naxis == 2)
+            return mat;
+        }
+      else if(bitpix ==  FLOAT_IMG)
+        {
+          TMatrix<float> *mat=new TMatrix<float>(nrows,ncols);
+          load(*mat,flnm);
+          if(naxis == 1)
+            {
+              TVector<float> *vect=new TVector<float>(*mat);
+              delete mat;
+              return vect;
+            }
+          else if(naxis == 2)
+            return mat;
+        }
+      else if(bitpix ==  LONG_IMG)
+        {
+          TMatrix<int> *mat=new TMatrix<int>(nrows,ncols);
+          load(*mat,flnm);
+          if(naxis == 1)
+            {
+              TVector<int> *vect=new TVector<int>(*mat);
+              delete mat;
+              return vect;
+            }
+          else if(naxis == 2)
+            return mat;
+        }
+    }
+  else
+    throw IOExc("Error:: this HDU is not a FITS image or binary table");
+
+  return NULL;
+}
+
 void FitsIoServer::load(TMatrix<double>& mat,char flnm[]) 
 {
-  int nbrows=0;
-  int nbcols=0;
-  FITS_tab_typ_ = TDOUBLE;
-  long naxis;
-  int n1, n2, n3;
+  int nbrows= 0;
+  int nbcols= 0;
+  FITS_tab_typ_= TDOUBLE;
+  int naxis;
+  int n1,n2,n3;
   DVList dvl;
-  planck_read_img(flnm, naxis, n1, n2, n3, dvl);
-
-  nbrows=n1;
-  nbcols=n2;
-  if (naxis == 1) nbcols=1;
-  if (naxis > 2 && n3>1) 
-    {
-      cout<<"  FitsIOServer : le fichier fits n'est pas une matrice, naxis= " <<naxis<< endl;
+  planck_read_img(flnm,naxis,n1,n2,n3,dvl);
+
+  nbrows= n1;
+  nbcols= n2;
+  if(naxis == 1) nbcols= 1;
+  if(naxis > 2 && n3 > 1) 
+    {
+      cout << " naxis = " << naxis << endl;
+      throw IOExc("FitsIOServer::load() this Fits file is not a matrix");
     }
 
   // number of components
-  if (mat.NRows() != nbrows || mat.NCols() != nbcols )
-    {
-      //    cout << " found " << nbrows << " rows ";
-      //   cout << "  expected " << mat.NRows()  << endl;
-      //  cout << " found " << nbcols << " columns " ;
-      // cout << "  expected " <<  mat.NCols()  << endl;
+  if(mat.NRows() != nbrows || mat.NCols() != nbcols )
+    {
       mat.ReSize(nbrows,nbcols);
-      cout << " resize the vector to nbrows=  " << nbrows  << " nbcols= " << nbcols << endl;
-    }
-  int ij=0;
-  for (int j=0; j< nbcols; j++) 
-    for (int i = 0; i < nbrows; i++)  mat(i,j) = (double)r_8tab_[ij++];
+      cout<<" FitsIOServer::load resize the matrix";
+      cout<<" nbrows= "<<nbrows<<" nbcols= "<<nbcols<<endl;
+    }
+
+  int ij= 0;
+  for(int j = 0; j < nbcols; j++) 
+    for (int i = 0; i < nbrows; i++)  mat(i,j)= (double)r_8tab_[ij++];
 }
 
 void FitsIoServer::load(TMatrix<float>& mat,char flnm[]) 
 {
-  int nbrows=0;
-  int nbcols=0;
-  FITS_tab_typ_ = TFLOAT;
-  long naxis;
-  int n1, n2, n3;
+  int nbrows= 0;
+  int nbcols= 0;
+  FITS_tab_typ_= TFLOAT;
+  int naxis;
+  int n1,n2,n3;
   DVList dvl;
-  planck_read_img(flnm, naxis, n1, n2, n3, dvl);
-
-  nbrows=n1;
-  nbcols=n2;
-  if (naxis == 1) nbcols=1;
-  if (naxis > 2 && n3>1) 
-    {
-      cout<<"  FitsIOServer : le fichier fits n'est pas une matrice, naxis= " <<naxis<< endl;
+  planck_read_img(flnm,naxis,n1,n2,n3,dvl);
+
+  nbrows= n1;
+  nbcols= n2;
+  if(naxis == 1) nbcols= 1;
+  if(naxis > 2 && n3 > 1) 
+    {
+      cout << " naxis = " << naxis << endl;
+      throw IOExc("FitsIOServer::load() this Fits file is not a matrix");
     }
 
   // number of components
-  if (mat.NRows() != nbrows || mat.NCols() != nbcols )
-    {
-      //      cout << " found " << nbrows << " rows ";
-      //    cout << "  expected " << mat.NRows()  << endl;
-      //  cout << " found " << nbcols << " columns " ;
-      //  cout << "  expected " <<  mat.NCols()  << endl;
+  if(mat.NRows() != nbrows || mat.NCols() != nbcols )
+    {
       mat.ReSize(nbrows,nbcols);
-      cout << " resize the vector to nbrows=  " << nbrows  << " nbcols= " << nbcols << endl;
-    }
-  int ij=0;
-  for (int j=0; j< nbcols; j++) 
-    for (int i = 0; i < nbrows; i++)  mat(i,j) = (float)r_4tab_[ij++];
+      cout<<" FitsIOServer::load resize the matrix";
+      cout<<" nbrows= "<<nbrows<<" nbcols= "<<nbcols<<endl;
+    }
+
+  int ij= 0;
+  for(int j = 0; j < nbcols; j++) 
+    for(int i = 0; i < nbrows; i++)  mat(i,j)= (float)r_4tab_[ij++];
 }
 
 void FitsIoServer::load(TMatrix<int_4>& mat,char flnm[]) 
 {
-  int nbrows=0;
-  int nbcols=0;
-  FITS_tab_typ_ = TINT;
-  long naxis;
-  int n1, n2, n3;
+  int nbrows= 0;
+  int nbcols= 0;
+  FITS_tab_typ_= TINT;
+  int naxis;
+  int n1,n2,n3;
   DVList dvl;
-  planck_read_img(flnm, naxis, n1, n2, n3, dvl);
-
-  nbrows=n1;
-  nbcols=n2;
-  if (naxis == 1) nbcols=1;
-  if (naxis > 2 && n3>1) 
-    {
-      cout<<"  FitsIOServer : le fichier fits n'est pas une matrice, naxis= " <<naxis<< endl;
+  planck_read_img(flnm,naxis,n1,n2,n3,dvl);
+
+  nbrows= n1;
+  nbcols= n2;
+  if(naxis == 1) nbcols= 1;
+  if(naxis > 2 && n3 > 1) 
+    {
+      cout << " naxis = " << naxis << endl;
+      throw IOExc("FitsIOServer::load() this Fits file is not a matrix");
     }
 
   // number of components
-  if (mat.NRows() != nbrows || mat.NCols() != nbcols )
-    {
-      //     cout << " found " << nbrows << " rows ";
-      //   cout << "  expected " << mat.NRows()  << endl;
-      //  cout << " found " << nbcols << " columns " ;
-      //  cout << "  expected " <<  mat.NCols()  << endl;
+  if(mat.NRows() != nbrows || mat.NCols() != nbcols )
+    {
       mat.ReSize(nbrows,nbcols);
-      cout << " resize the vector to nbrows=  " << nbrows  << " nbcols= " << nbcols << endl;
-    }
-  int ij=0;
-  for (int j=0; j< nbcols; j++) 
-    for (int i = 0; i < nbrows; i++)  mat(i,j) = (int_4)i_4tab_[ij++];
-}
-
-void FitsIoServer::load(NTuple& ntpl,char flnm[],int hdunum) 
-
-  //********************************************************/
-  //* move to the HDU which has the specified number hdunum*/
-  //* in the FITS file, perform read operations and write  */
-  //* the elements in an NTuple.                           */
-  //********************************************************/
-{
-
-  // pointer to the FITS file, defined in fitsio.h 
+      cout<<" FitsIOServer::load resize the matrix";
+      cout<<" nbrows= "<<nbrows<<" nbcols= "<<nbcols<<endl;
+    }
+
+  int ij= 0;
+  for(int j = 0; j < nbcols; j++) 
+    for(int i = 0; i < nbrows; i++)  mat(i,j)= (int_4)i_4tab_[ij++];
+}
+
+void FitsIoServer::load(NTuple& ntpl,char flnm[],int hdunum)
+
+  // *****************************************************
+  // move to the HDU which has the specified number hdunum
+  // in the FITS file, read data values form an ASCII or  
+  // binary table and write the elements in an NTuple.
+  // Only TFLOAT or TDOUBLE datatype are allowed
+  // *****************************************************
+
+{
+  // pointer to the FITS file, defined in fitsio.h
   fitsfile *fptr;
-  int status = 0;
-  if( fits_open_file(&fptr,flnm,READONLY,&status) ) 
-    printerror( status );
+  int status= 0;
+  fits_open_file(&fptr,flnm,READONLY,&status);
+  if( status ) printerror( status );
 
   // move to the HDU
-  int hdutype;
-  if( fits_movabs_hdu(fptr,hdunum,&hdutype,&status) ) 
-    printerror( status );
+  int hdutype= 0;
+  fits_movabs_hdu(fptr,hdunum,&hdutype,&status);
+  if( status ) printerror( status );
+
+  if(hdutype == ASCII_TBL)
+    printf("\nReading ASCII table in HDU %d:\n",hdunum);
+  else if(hdutype == BINARY_TBL)
+    printf("\nReading binary table in HDU %d:\n",hdunum);
+  else
+    {
+      printf("Error:: this HDU is not an ASCII or binary table\n");
+      throw IOExc("FitsIoServer::load(NTuple& ," + (string)flnm + ") Error");
+    }
+
+  // get the number of columns
+  int ncols= 0;
+  fits_get_num_cols(fptr,&ncols,&status);
+  if( status ) printerror( status );
+
+  // get the number of rows
+  long naxis2= 0;
+  fits_get_num_rows(fptr,&naxis2,&status);
+  if( status ) printerror( status );
+  int nrows= (int)naxis2;
+
+  // get the datatype of the values and the min repeat count
+  list<int> tfields;
+  long repeat;
+  int ntest= 0;
+  for(int ii = 0; ii < ncols; ii++)
+    {
+      int DTYPE;
+      long rept,width;
+      fits_get_coltype(fptr,ii+1,&DTYPE,&rept,&width,&status);
+      if(DTYPE == TFLOAT || DTYPE == TDOUBLE)
+        {
+          tfields.push_back(ii+1);
+          if(ntest == 0) repeat= rept;
+          if(rept < repeat) repeat= rept;
+          ntest++;
+        }
+    }
+  
+  if(tfields.empty())
+    {
+      cout << " No data values of type TFLOAT or TDOUBLE" << endl;
+      throw IOExc("FitsIoServer::load(NTuple&) Exit");
+    }
+  else
+    cout << " nbre tfields= " << tfields.size() << " rep= " << repeat << endl; 
+
+  // get input elements to create the NTuple
+  char **clname;
+  clname= new char*[tfields.size()];
+  for(int ii = 0; ii < tfields.size(); ii++)
+    {
+      ostringstream oss; oss<<"Col_"<<ii+1;
+      string ss= oss.str();
+      clname[ii]= new char[FLEN_VALUE];
+      strcpy(clname[ii],ss.data());
+    }
+  
+  // create a NTuple
+  NTuple nt0(tfields.size(),clname);
+  for(int ii = 0; ii < tfields.size(); ii++)
+    delete [] clname[ii];
+  delete [] clname;
+
+  // value to represent undefined array elements
+  int anull= 0;
+  float fnull = FLOATNULLVALUE;
+  float value[1];
+
+  // read elements from columns and fill NTuple
+  for (int k = 0; k < nrows; k++)
+    {
+      for(int i = 0; i < repeat; i++)
+        {
+          int j= 0;
+          float *xnt= new float[tfields.size()];
+          list<int>::iterator jtr;
+          for(jtr= tfields.begin(); jtr != tfields.end(); jtr++)
+            {
+              fits_read_col(fptr,TFLOAT,*jtr,k+1,i+1,1,&fnull,value,
+                            &anull,&status);
+              if( status ) printerror(status,"fits_read_col");
+              xnt[j++]= value[0];
+            }
+          nt0.Fill(xnt);
+          delete[] xnt;
+        }
+    } 
 
   // get number of keywords
   int nkeys,keypos;
-  if( fits_get_hdrpos(fptr,&nkeys,&keypos,&status) )
-    printerror( status );
-
-
-  if (hdutype == BINARY_TBL)
-    printf("\nReading binary table in HDU %d:\n",hdunum);
-  else
-    {
-      printf("Error:: this HDU is not a binary table\n");
-      throw IOExc("FitsIoServer::load(NTuple& ," + (string)flnm + ") Error Not a bin table !");
-//      exit ( status );
-    }
-
-  // number of columns
-  int tfields;
-  if( fits_get_num_cols(fptr,&tfields,&status) )
-    printerror( status );
-
-  // to get table size
-  long naxis[2];
-  int nfound;
-   if( fits_read_keys_lng(fptr, "NAXIS", 1, 2, naxis, &nfound, &status) ) 
-     printerror( status );
-   int nrows= naxis[1];
-
-   //Information about each column
-   char **ttype, **tform;
-   ttype= new char*[tfields];
-   tform= new char*[tfields];
-   int ii;
-   for(ii = 0; ii < tfields; ii++)
-     {
-       ttype[ii]= new char[FLEN_VALUE];
-       tform[ii]= new char[FLEN_VALUE];
-     }
-
-   // number of TTYPEn,TFORMn and TUNITn keywords found
-   int num= 0;
-
-   // read the column names from the TTYPEn keywords
-   fits_read_keys_str(fptr,"TTYPE",1,tfields,ttype,&nfound,&status);
-   num += nfound;
-   // read the column names from the TFORMn keywords
-   fits_read_keys_str(fptr,"TFORM",1,tfields,tform,&nfound,&status);
-   num += nfound;
-
-   //for(int ii = 0; ii < tfields; ii++)
-   //  printf("\nColumn name & Format %8s %8s", ttype[ii],tform[ii]);
+  fits_get_hdrpos(fptr,&nkeys,&keypos,&status);
+  if( status ) printerror( status );
+
+  // put other reserved keywords in a DVList object
+  char keyname[LEN_KEYWORD]= "";
+  char strval[FLEN_VALUE]= "";
+  char dtype;
+  char card[FLEN_CARD];  
+  char *comkey = "COMMENT";
+
+  // shift with the number of mandatory keywords
+  int num= 8;
+
+  for(int j = num+1; j <= nkeys; j++)
+    {
+      fits_read_keyn(fptr,j,card,strval,NULL,&status);
+      if(status) printerror(status);
+
+      strncpy(keyname,card,LEN_KEYWORD-1);
+      
+      if(strncmp(keyname,comkey,LEN_KEYWORD-1) != 0 && strlen(keyname) != 0 
+         && strlen(strval) != 0)
+        {
+          fits_get_keytype(strval,&dtype,&status);
+          if(status)  printerror(status);
+          
+          strip(keyname, 'B',' '); 
+          strip(strval, 'B',' ');
+          strip(strval, 'B','\'');
+          switch( dtype )  
+            { 
+            case 'C':
+              nt0.Info()[keyname]= strval;        
+              break;
+            case 'I':
+              int ival;
+              fits_read_key(fptr,TINT,keyname,&ival,NULL,&status);
+              nt0.Info()[keyname]= ival;
+              break;
+            case 'L':
+              int ilog;
+              if(strncmp(strval,"T",1) == 0) ilog= 1;
+              else ilog= 0;
+              nt0.Info()[keyname]= ilog;
+              break;
+            case 'F':
+              double dval;
+              fits_read_key(fptr,TDOUBLE,keyname,&dval,NULL,&status);
+              nt0.Info()[keyname]= dval;
+              break;
+            }
+        }
+    }
   
-   // select the columns with float data values 
-   int typecode;
-   long repeat,width;
-   list<int> column;
-   for(ii = 0; ii < tfields; ii++)
-     {
-       fits_binary_tform(tform[ii],&typecode, &repeat, &width, &status);
-       //printf("\n%4s %3d %2ld %2ld", tform[ii], typecode, repeat, width);
-       if(typecode == TFLOAT) 
-         column.push_back(ii+1);
-     }
-
-   // get input elements to create the NTuple
-   char **clname;
-   clname= new char*[column.size()];
-   for(ii = 0; ii < column.size(); ii++)
-     clname[ii]= new char[FLEN_VALUE];
-
-   list<int>::iterator itr;
-   int index= 0;
-   for(itr= column.begin(); itr != column.end(); itr++)
-     strcpy(clname[index++],ttype[*itr-1]);
-
-   for(ii = 0; ii < tfields; ii++)
-     {
-       delete [] ttype[ii];
-       delete [] tform[ii];
-     }
-   delete [] ttype;
-   delete [] tform;
-   
-
-   // check if the specified keyword BLK exists
-   int blk= 512;
-   if(check_keyword(fptr,nkeys,"BLK"))
-     {
-       if( fits_read_key(fptr,TINT,"BLK",&blk,NULL,&status) )
-         printerror( status );
-     }
-   // create a NTuple
-   NTuple nt0(column.size(),clname,blk);
-
-   for(ii = 0; ii < column.size(); ii++)
-     delete [] clname[ii];
-   delete [] clname;
-
-   float value[1];
-   long felem = 1;
-   char strnull[10];
-   strcpy(strnull, " ");
-   int anynull= 0;
-
-   // value to represent undefined array elements
-   float floatnull = FLOATNULLVALUE;
-
-   // read elements from columns and fill NTuple
-   for (int k = 0; k < nrows; k++)
-     {
-       int j= 0;
-       float *xnt= new float[column.size()];
-       list<int>::iterator itr;
-       for(itr= column.begin(); itr != column.end(); itr++)
-         {
-           fits_read_col(fptr,TFLOAT,*itr,k+1,felem,1,&floatnull,value,
-                         &anynull, &status);
-           xnt[j++]= value[0];
-         }
-       nt0.Fill(xnt);
-       delete[] xnt;
-     }
-
-   // the TUNITn keywords are optional, if they exist they are put 
-   // in the DVLIst object
-   char keyname[LEN_KEYWORD]= "";
-   char strval[FLEN_VALUE];
-   for(ii = 0; ii < tfields; ii++)
-     {
-       fits_make_keyn("TUNIT",ii+1,keyname,&status);
-       if(check_keyword(fptr,nkeys,keyname))
-         {
-           num++;
-           if( fits_read_key_str(fptr,keyname,strval,NULL,&status) )
-             printerror(status);
-           strip (keyname, 'B',' '); 
-           strip(strval, 'B',' ');
-           nt0.Info()[keyname]= strval;
-         }
-     }
-
-   // add the number of mandatory keywords of a binary table
-   num += 8;
-
-   // put names and values of other reserved keywords in a DVList object
-   char comment[FLEN_COMMENT];
-   char dtype;
-   for(int j = num+1; j <= nkeys; j++)
-     {
-       fits_read_keyn(fptr,j,keyname,strval,comment,&status);
-       fits_get_keytype(strval,&dtype,&status);
-       strip (keyname, 'B',' '); 
-       strip(strval, 'B',' ');
-       strip(strval, 'B','\'');
-       
-       switch( dtype ) 
-         {
-         case 'C':
-           nt0.Info()[keyname]= strval;   
-           break;
-         case 'I':
-           int ival;
-           ctoi(strval,&ival);
-           nt0.Info()[keyname]= (int_4)ival;
-           break;
-         case 'F':
-           double dval;
-           ctof(strval,&dval);
-           nt0.Info()[keyname]= dval;
-           break;
-         }
-     }
-
-   // copy in the input NTuple
-   ntpl= nt0;
-
-   if( fits_close_file(fptr, &status) )
-     printerror(status);
-
-   printf("\n");
-   return;
-}
-
+  // copy in the input NTuple
+  ntpl= nt0;
+  
+  // close the file
+  fits_close_file(fptr,&status);
+  if(status) printerror(status);
+}
 
 void FitsIoServer::load(SphericalMap<double>& sph, char flnm[]) 
@@ -325 +463 @@
   int npixels=0;
   int nside=0;
-  long naxis;
+  int naxis;
   int n1, n2, n3;
 
@@ -365 +503 @@
 
 
-void FitsIoServer::load(SphereGorski<float>& sph, char flnm[], int hdunum) 
-{
-  int npixels=0;
-  int nside=0;
-
-  FITS_tab_typ_ = TFLOAT;
-
+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();
+
+  planck_read_bntbl(flnm,hdunum,npixels,dvl);
   nside= dvl.GetI("NSIDE");
+
   const char* ordering= dvl.GetS("ORDERING").c_str();
+
   char* ring = "RING";
-  if (  strncmp(ordering, ring,4) != 0) 
-
-    //  if (ordering!="RING    ")
-    {
-      cout << " numerotation non RING" << endl;
-    }
+  if(strncmp(ordering,ring,4) != 0) 
+    cout << " numerotation non RING" << endl;
+
+
   // number of pixels in the sphere
-  if (sph.NbPixels() != npixels)
-    {
-      //DBG      cout << " found " << npixels << " pixels" << endl;
-      //DBG      cout << " expected " << sph.NbPixels() << endl;
+  if(sph.NbPixels() != npixels)
+    {
       if (nside <= 0 )  
         {
           cout<<" FITSIOSERVER: no resolution parameter on fits file "<<endl;
-          throw IOExc("FitsIoServer::load(SphericalMap<float>& ," + (string)flnm + ", ) Error No resolution parameter !");
-//        exit(0);
+          throw IOExc("FitsIoServer::load(SphereGorski<float>& ," + (string)flnm + ", ) Error No resolution parameter !");
         }
       if (nside != sph.SizeIndex())
         { 
           sph.Resize(nside);
-          cout << "FitsIoServer::load(SphereGorski<float> ...) ReSizing to NSide= " << nside  << endl;
-        }
-      //      else
-      //        {
-      //          cout << " FITSIOSERVER : same resolution, surprising ! " << endl;
-//        exit(0);   $CHECK$ - Il ne faut pas sortir, me semble-t-il , Reza 20/11/99
-      //        }
-    }
-  for (int j = 0; j < sph.NbPixels(); j++) 
-    {
-      sph(j)= (float)r_4tab_[j];
-    }
-}
-void FitsIoServer::load(SphereGorski<double>& sph, char flnm[], int hdunum) 
+          cout << " FitsIoServer::load(SphereGorski<float> ...)";
+          cout << "  ReSizing to NSide= " << nside << endl;
+        }
+    }
+  
+  for(int j = 0; j < sph.NbPixels(); j++) sph(j)= (float)r_4tab_[j];
+}
+
+void FitsIoServer::load(SphereGorski<double>& sph,char flnm[],int hdunum) 
+{
+  int npixels= 0;
+  int nside  = 0;
+
+  FITS_tab_typ_ = TDOUBLE;
+  DVList dvl;
+  planck_read_bntbl(flnm,hdunum,npixels,dvl);
+
+  nside= dvl.GetI("NSIDE");
+
+  const char* ordering= dvl.GetS("ORDERING").c_str();
+
+  char* ring = "RING";
+  if(strncmp(ordering,ring,4) != 0) 
+    cout << " numerotation non RING" << endl;
+
+  // number of pixels in the sphere
+  if(sph.NbPixels() != npixels)
+    {
+      if(nside <= 0) 
+        throw IOExc("FitsIoServer::load(SphereGorski<double>& ," + (string)flnm + ", )  No resol parameter !");
+
+      if(nside != sph.SizeIndex())
+        { 
+          sph.Resize(nside);
+          cout << " FitsIoServer::load(SphereGorski<double> ...)";
+          cout << "  ReSizing to NSide= " << nside << endl;
+        }
+    }
+  
+  for(int j = 0; j < sph.NbPixels(); j++) sph(j)= (double)r_8tab_[j];
+}
+
+void FitsIoServer::load(SphericalMap<float>& sph, char flnm[]) 
 {
   int npixels=0;
   int nside=0;
-
-  FITS_tab_typ_ = TDOUBLE;
-
+  int naxis;
+  int n1, n2, n3;
   DVList dvl;
-  planck_read_bntbl(flnm, hdunum, npixels, dvl);
-  //dvl.Print();
+
+  FITS_tab_typ_ = TFLOAT;
+
+  planck_read_img(flnm, naxis, n1, n2, n3, dvl);
+  if (naxis != 1) 
+    {
+      cout << " le fichier fits n'est pas une sphere, naxis= " << naxis << endl;
+    }
+  npixels=n1;
   nside= dvl.GetI("NSIDE");
-  const char* ordering= dvl.GetS("ORDERING").c_str();
-  char* ring = "RING";
-  if (  strncmp(ordering, ring,4) != 0) 
-
-    //  if (ordering!="RING    ")
-    {
-      cout << " numerotation non RING" << endl;
-    }
+
   // number of pixels in the sphere
   if (sph.NbPixels() != npixels)
@@ -438 +600 @@
         {
           cout<<" FITSIOSERVER: no resolution parameter on fits file "<<endl;
-          throw IOExc("FitsIoServer::load(SphereGorski<double>& ," + (string)flnm + ", )  No resol parameter !");
+          throw IOExc("FitsIoServer::load(SphericalMap<float>& ," + (string)flnm + ") No resolution param !");
 //        exit(0);
         }
@@ -449 +611 @@
         {
           cout << " FITSIOSERVER : same resolution, surprising ! " << endl;
-//        exit(0);   $CHECK$ - ne pas sortir , Reza 20/11/99
-        }
-    }
-  for (int j = 0; j < sph.NbPixels(); j++) sph(j)= (double)r_8tab_[j];
-}
-
-void FitsIoServer::load(SphericalMap<float>& sph, char flnm[]) 
-{
-  int npixels=0;
-  int nside=0;
-  long naxis;
-  int n1, n2, n3;
-  DVList dvl;
-
-  FITS_tab_typ_ = TFLOAT;
-
-  planck_read_img(flnm, naxis, n1, n2, n3, dvl);
-  if (naxis != 1) 
-    {
-      cout << " le fichier fits n'est pas une sphere, naxis= " << naxis << endl;
-    }
-  npixels=n1;
-  nside= dvl.GetI("NSIDE");
-
-  // 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;
-          throw IOExc("FitsIoServer::load(SphericalMap<float>& ," + (string)flnm + ") No resolution param !");
-//        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);  $CHECK$ - Ne pas sortir , Reza 20/11/99
         }
@@ -503 +622 @@
   int nbcols=0;
   FITS_tab_typ_ = TDOUBLE;
-  long naxis;
+  int naxis;
   int n1, n2, n3;
   DVList dvl;
@@ -545 +664 @@
 {
   FITS_tab_typ_ = TFLOAT;
-  long naxis;
+  int naxis;
   int siz_x;
   int siz_y;
@@ -585 +704 @@
 {
   FITS_tab_typ_ = TINT;
-  long naxis;
+  int naxis;
   int siz_x;
   int siz_y;
@@ -1235 +1354 @@
 
   // primary header
-  int bitpix=LONG_IMG;
-  int naxis=0;
-  fits_create_img(fptr, bitpix, naxis, NULL, &status);
+  //  int bitpix=LONG_IMG;
+  //  int naxis=0;
+  // fits_create_img(fptr, bitpix, naxis, NULL, &status);
   // write the current date
-  fits_write_date(fptr, &status);
+  //  fits_write_date(fptr, &status);
   //DBG  cerr << " DBG - Apres write_date status = " << status << endl;
-  if( status )  printerror( status );
+  // if( status )  printerror( status );
 
 
@@ -1367 +1486 @@
 }
 
-void FitsIoServer::planck_read_img(char flnm[], long& naxis,int& n1, int& n2, int& n3, DVList& dvl)
-{
-  int status=0;
-  long bitpix;
-  long naxes[3]={0,0,0};
-  char* comment=NULL;
+void FitsIoServer::planck_read_img(char flnm[],int &naxis,int &n1,int &n2,int &n3,DVList &dvl)
+{
+  int status= 0;
 
   // 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, "BITPIX", &bitpix, comment, &status);
-  if( status )  printerror( status );
-  fits_read_key_lng(fptr, "NAXIS", &naxis, comment, &status);
-  if( status )  printerror( status );
+  fits_open_file(&fptr,flnm,READONLY,&status);
+  if( status )  printerror( status );
+
+  // bits per pixels
+  int bitpix;
+  fits_read_key(fptr,TINT,"BITPIX",&bitpix,NULL,&status);
+  if( status )  printerror( status );
+
+  // number of dimensions in the FITS array
+  fits_read_key(fptr,TINT,"NAXIS",&naxis,NULL,&status);
+  if( status )  printerror( status );
+
+  // read the NAXIS1,NAXIS2 and NAXIS3 keyword to get image size
+  long naxes[3]= {0,0,0};
   int nfound;
-  int nkeys=(int)naxis;
-  fits_read_keys_lng(fptr, "NAXIS", 1, nkeys, naxes, &nfound, &status); 
-  if( status )  printerror( status );
-
-  n1 = naxes[0] ;
-  n2 = naxes[1] ;
-  n3 = naxes[2] ;
-
+  fits_read_keys_lng(fptr,"NAXIS",1,naxis,naxes,&nfound,&status); 
+  if( status )  printerror( status );
+
+  n1 = (int)naxes[0];
+  n2 = (int)naxes[1];
+  n3 = (int)naxes[2];
 
   long nelements= naxes[0];
-  if (naxis >=2) nelements*=naxes[1];
-  if (naxis == 3)  nelements*=naxes[2];
-  int anynull=0;
-  r_8 dnullval=DOUBLENULLVALUE;
-  r_4 fnullval=FLOATNULLVALUE;
-  int_4 inullval=0;
+  if(naxis >= 2) nelements*= naxes[1];
+  if(naxis == 3) nelements*= naxes[2];
+
+  int anull= 0;
+  r_8 dnull= DOUBLENULLVALUE;
+  r_4 fnull= FLOATNULLVALUE;
+  int_4 inull= 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_) 
+  switch (FITS_tab_typ_) 
     {
     case TDOUBLE :
-      if (bitpix !=  DOUBLE_IMG) 
-        {
-          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_ = NULL; }
+      if(bitpix != DOUBLE_IMG) 
+        {
+          cout << " FitsIOServer:: the data type on fits file " << flnm;
+          cout << " is not double, " << "conversion to double will";
+          cout << " be achieved by cfitsio lib" << endl; 
+        }
+      if(r_8tab_ != NULL) { delete [] r_8tab_; r_8tab_= NULL; }
       r_8tab_=new r_8[nelements];
-      fits_read_img(fptr, TDOUBLE, 1, nelements,  &dnullval, r_8tab_, 
-                    &anynull, &status);
+      fits_read_img(fptr,TDOUBLE,1,nelements,&dnull,r_8tab_,&anull,&status);
       if( status )  printerror( status );
       break;
+
     case TFLOAT :
-      if (bitpix !=  FLOAT_IMG) 
-        {
-          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_ = NULL; } 
+      if(bitpix != FLOAT_IMG) 
+        {
+          cout << " FitsIOServer:: the data type on fits file " << flnm;
+          cout << " is not float, " << "conversion to float will";
+          cout << " be achieved by cfitsio lib" << endl; 
+        }
+      if(r_4tab_ != NULL) { delete [] r_4tab_; r_4tab_= NULL; } 
       r_4tab_=new r_4[nelements];
-      fits_read_img(fptr, TFLOAT, 1, nelements,  &fnullval, r_4tab_, 
-                    &anynull, &status);
+      fits_read_img(fptr,TFLOAT,1,nelements,&fnull,r_4tab_,&anull,&status);
+      if( status )  printerror( status );
       //SHV: remove useless print...
-      //      for (int kk=0; kk<nelements; kk++) cout << r_4tab_[kk] << endl;
-      cout << " anynull= " << anynull << endl;
+      //     for (int kk=0; kk<nelements; kk++) cout << r_4tab_[kk] << endl;
+      break;
+
+    case TINT :
+      if(bitpix != LONG_IMG) 
+        {
+          cout << " FitsIOServer:: the data type on fits file " << flnm;
+          cout << " is not long, " << "conversion to long will";
+          cout << " be achieved by cfitsio lib" << endl; 
+        }
+      if (i_4tab_ != NULL) { delete [] i_4tab_; i_4tab_= NULL; }
+      i_4tab_=new int_4[nelements];
+      fits_read_img(fptr,TINT,1,nelements,&inull,i_4tab_,&anull, &status);
       if( status )  printerror( status );
       break;
 
-
-    case TINT :
-      if (bitpix !=  LONG_IMG) 
-        {
-          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_ = NULL; }
-      i_4tab_=new int_4[nelements];
-      fits_read_img(fptr, TINT, 1, nelements,  &inullval, i_4tab_, 
-                    &anynull, &status);
-      if( status )  printerror( status );
+    default :
+      cout << " FitsIOServer:: datatype code " << FITS_tab_typ_;
+      cout << " FitsIOServer::planck_read_img  not yet programmed" << endl;
       break;
-
-
-    default :
-      cout << " FitsIOServer::read_img : type non traite: " << FITS_tab_typ_  << endl;
-      break;
-    }
+    }
+
   status = 0;
+
+  int nkeys;
+  int keypos;
+  fits_get_hdrpos(fptr,&nkeys,&keypos,&status);
+  if( status )  printerror( status );
+
   char card[FLEN_CARD];
   char keyname[LEN_KEYWORD]= "";
   char strval[FLEN_VALUE];
-  char comment1[FLEN_COMMENT];
   char dtype;
-  char * comkey = "COMMENT";
-  int keypos;
-   fits_get_hdrpos(fptr,&nkeys,&keypos,&status);
-   if( status )  printerror( status );
-   //   cout << " planck_read_img : nombre de mots-cles : " << nkeys << endl;
-   for(int j = 1; j <= nkeys; j++)
-     {
-       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);
-            //          cout<<" keyname= "<< keyname <<" dtype= "<<dtype <<endl; 
-            strip (keyname, 'B',' '); 
-            strip(strval, '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 'L':
-                int ilog;
-                if (strncmp(strval,"T" ,1) ==0) ilog=1;
-                else ilog=0;
-                dvl[keyname]= (int_4)ilog;
-                break;
-              case 'F':
-                double dval;
-                ctof(strval,&dval);
-                dvl[keyname]= dval;
-                break;
-              default :
-                cout << " FitsIOServer : type de mot-cle bizarre: " << keyname  <<" dtype= "<<dtype << endl;
-                break;
-              }
-          }
-     }
+  char *comkey = "COMMENT";
+
+  for(int j = 1; j <= nkeys; j++)
+    {
+      fits_read_keyn(fptr,j,card,strval,NULL,&status);
+      strncpy(keyname,card,LEN_KEYWORD-1);
+
+      if(strncmp(keyname,comkey,LEN_KEYWORD-1) != 0 && strlen(keyname) != 0) 
+        {
+          fits_get_keytype(strval,&dtype,&status);
+          strip(keyname, 'B',' '); 
+          strip(strval, '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 'L':
+              int ilog;
+              if (strncmp(strval,"T" ,1) ==0) ilog=1;
+              else ilog=0;
+              dvl[keyname]= (int_4)ilog;
+              break;
+            case 'F':
+              double dval;
+              ctof(strval,&dval);
+              dvl[keyname]= dval;
+              break;
+            default :
+              cout << " FitsIOServer:: type de mot-cle bizarre";
+              cout << " keyname= " << keyname  <<" dtype= "<< dtype << endl;
+              break;
+            }
+        }
+    }
 
   // close the file
-  status=0;
+  status= 0;
   fits_close_file(fptr, &status);
   if( status )  printerror( status );
 }
-
 
 void FitsIoServer::planck_read_bntbl(char flnm[], int hdunum, int& npixels, DVList& dvl)
@@ -1657 +1781 @@
   //char bidon[LEN_KEYWORD];
   char * comkey = "COMMENT";
+  char blank[LEN_KEYWORD]= "";
 
    for(int j = 1; j <= nkeys; j++)
      {
-       //       fits_read_record(fptr, j, card, &status);
-       // strncpy(keyname,card,LEN_KEYWORD-1);
-       // cout << " bidon= " << keyname << endl;
+       //fits_read_record(fptr, j, card, &status);
+       //strncpy(keyname,card,LEN_KEYWORD-1);
+       //cout << " cle= " << keyname << endl;
        // if ( strncmp(keyname,comkey ,LEN_KEYWORD-1) != 0)
-       fits_read_keyn(fptr,j,card,strval,comment1,&status);
+
+       fits_read_keyn(fptr,j,card,strval,NULL,&status);
        strncpy(keyname,card,LEN_KEYWORD-1);
 
-        if ( strncmp(keyname,comkey ,LEN_KEYWORD-1) != 0)
-          {
-            fits_get_keytype(strval,&dtype,&status);
-            //  cout<<" keyname= "<< keyname <<" dtype= "<<dtype <<endl; 
-            strip (keyname, 'B',' '); 
-            strip(strval, '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 'L':
-                int ilog;
-                if (strncmp(strval,"T" ,1) ==0) ilog=1;
-                else ilog=0;
-                dvl[keyname]= (int_4)ilog;
-                break;
-              case 'F':
-                double dval;
-                ctof(strval,&dval);
-                dvl[keyname]= dval;
-                break;
-              default :
-                cout << " FitsIOServer : type de mot-cle bizarre: "  << keyname <<" dtype= "<<dtype << endl;
-                break;
-              }
-          }
+       if(strncmp(keyname,comkey,LEN_KEYWORD-1) != 0 && strlen(keyname) != 0)  
+         {
+           fits_get_keytype(strval,&dtype,&status);
+           //   cout<<" keyname= "<< keyname <<" dtype= "<<dtype <<endl; 
+           strip (keyname, 'B',' '); 
+           strip(strval, '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 'L':
+               int ilog;
+               if (strncmp(strval,"T" ,1) ==0) ilog=1;
+               else ilog=0;
+               dvl[keyname]= (int_4)ilog;
+               break;
+             case 'F':
+               double dval;
+               ctof(strval,&dval);
+               dvl[keyname]= dval;
+               break;
+             default :
+               cout << " FitsIOServer : type de mot-cle bizarre: "  << keyname <<" dtype= "<<dtype << endl;
+               break;
+             }
+         }
      }
-
+   
    // close the file
    status=0;

trunk/SophyaExt/FitsIOServer/fitsioserver.h

@@ -11 +11 @@
 #include "cimage.h"
 #include "dvlist.h"
-
+#include "anydataobj.h"
 #include <fitsio.h>   /* entre <> pour que mkmflib fonctionne correctement  */
 
@@ -21 +21 @@
 public :
 
-  FitsIoServer() : FITS_tab_typ_(0),i_4tab_(NULL),r_4tab_(NULL),r_8tab_(NULL) {;}
+  FitsIoServer() : FITS_tab_typ_(0),i_4tab_(NULL),r_4tab_(NULL),r_8tab_(NULL) 
+  {;}
  ~FitsIoServer()
    {
@@ -56 +57 @@
  void Mollweide_picture_projection(SphericalMap<double>& sph, char flnm[]);
  void picture(LocalMap<double>& lcm, char flnm[]);
-
-
-
  void readheader(char flnm[]);
  
+ virtual AnyDataObj* loadobj(char flnm[],int hdunum = 1);
+
  private :
-  
 
  // creer, ecrire une imageformat FITS, a partir des tableaux de donnees 
  // dtab_, ftab_ prealablement remplis 
- void planck_write_img(char flnm[],int naxis,int n1,int n2,int n3,DVList& dvl);
+ void planck_write_img(char flnm[],int naxis,int n1,int n2,int n3,DVList &dvl);
 
- void planck_read_img(char flnm[],long& naxis,int& n1,int& n2,int& n3,DVList& dvl);
+ void planck_read_img(char flnm[],int &naxis,int &n1,int &n2,int &n3,DVList &dvl);
  void planck_read_bntbl(char flnm[], int hdunum, int& npixels, DVList& dvl);
  void planck_write_bntbl(char flnm[], int npixels, char typeOfContent[], char extname[], char comment1[], DVList& dvl);