Ignore:
Timestamp:
Oct 19, 1999, 12:11:06 PM (26 years ago)
Author:
ansari
Message:

load et save d'images 19-OCT-99 GLM

File:
1 edited

Legend:

Unmodified
Added
Removed

  • trunk/SophyaExt/FitsIOServer/fitsioserver.cc

    r471 r477  
    11
    2 #include "nbmath.h"
    3 #include <math.h>
    42#include <iostream.h>
     3#include <list>
     4
    55#include "fitsioserver.h"
    6 #include "dvlist.h"
    76#include "strutil.h"
    87
    9   typedef map<string, MuTyV, less<string> >  ValList;
    10 void FitsIoServer::load(TMatrix<double>& mat, char flnm[])
     8void FitsIoServer::load(TMatrix<double>& mat,char flnm[])
    119{
    1210  int nbrows=0;
    … …  
    2927  nbcols=n2;
    3028  if (naxis == 1) nbcols=1;
    31   // cout << " lecture  du dvlist par load " << endl;
    32   dvl.Print();
    33   ValList::const_iterator it;
     29  //dvl.Print();
     30  DVList::ValList::const_iterator it;
    3431  for(it = dvl.Begin(); it != dvl.End(); it++) 
    3532    {
    … …  
    5249          break;
    5350        default :
    54           cout << " probleme dans type mot cle optionnel" << endl;
     51          cout << " FitsIOServer : probleme dans type mot cle optionnel" << endl;
    5552          break;
    5653        }
    … …  
    5956  if (naxis > 2)
    6057    {
    61       cout << " le fichier fits n'est pas une matrice, naxis= " << naxis << endl;
    62     }
     58      cout<<"  FitsIOServer : le fichier fits n'est pas une matrice, naxis= " <<naxis<< endl;
     59    }
     60
    6361  // number of components
    64 
    6562  if (mat.NRows() != nbrows || mat.NCols() != nbcols )
    6663    {
    … …  
    7572  for (int j=0; j< nbcols; j++)
    7673    for (int i = 0; i < nbrows; i++)  mat(i,j) = dtab_[ij++];
    77   // cout << " chargement termine " << endl;
    7874}
    7975
    … …  
    8480  FITS_tab_typ_ = TDOUBLE;
    8581  read_sphere(flnm, npixels, nside);
     82
    8683  // number of pixels in the sphere
    8784  if (sph.NbPixels() != npixels)
    … …  
    9188      if (nside <= 0 )
    9289        {
    93           cout << " FITSIOSERVER : no resolution parameter on fits file " << endl;
     90          cout<<" FITSIOSERVER: no resolution parameter on fits file "<<endl;
    9491          exit(0);
    9592        }
    … …  
    9794        {
    9895          sph.Resize(nside);
    99           cout << " resize the sphere to nside=  " << nside  << endl;
     96          cout << " resizing the sphere to nside=  " << nside  << endl;
    10097        }
    10198      else
    … …  
    106103    }
    107104  for (int j = 0; j < sph.NbPixels(); j++) sph(j)= dtab_[j];
    108   // cout << " chargement termine " << endl;
    109105}
    110106void FitsIoServer::load(SphericalMap<float>& sph, char flnm[])
    … …  
    121117      if (nside <= 0 )
    122118        {
    123           cout << " FITSIOSERVER : no resolution parameter on fits file " << endl;
     119          cout<<" FITSIOSERVER: no resolution parameter on fits file "<<endl;
    124120          exit(0);
    125121        }
    … …  
    127123        {
    128124          sph.Resize(nside);
    129           cout << " resize the sphere to nside=  " << nside  << endl;
     125          cout << " resizing the sphere to nside=  " << nside  << endl;
    130126        }
    131127      else
    … …  
    136132    }
    137133  for (int j = 0; j < sph.NbPixels(); j++) sph(j)= ftab_[j];
    138   //  cout << " chargement termine " << endl;
    139134}
    140135
    … …  
    162157  if (naxis != 2)
    163158    {
    164       cout << " le fichier fits n'est pas une localmap, naxis= " << naxis << endl;
    165     }
    166   // cout << " lecture  du dvlist par load " << endl;
    167   dvl.Print();
    168   ValList::const_iterator it;
     159      cout<<" FitsIOServer : le fichier fits n'est pas une localmap, naxis= "<<naxis<<endl;
     160    }
     161  //dvl.Print();
     162  DVList::ValList::const_iterator it;
    169163  for(it = dvl.Begin(); it != dvl.End(); it++) 
    170164    {
    … …  
    187181          break;
    188182        default :
    189           cout << " probleme dans type mot cle optionnel" << endl;
     183          cout << " FitsIOServer : probleme dans type mot cle optionnel" << endl;
    190184          break;
    191185        }
    192186      if (!strcmp(keyname, "NSIDE") )
    193187        {
    194           // cout << " j'ai trouve NSIDE " << endl;
    195188          nside = ival;
    196189        }
    197190      //if (!strcmp(keyname, "ORDERING") )
    198191      //        {
    199       //          cout << " j'ai trouve ORDERING " << endl;
     192      //  cout << " j'ai trouve ORDERING " << endl;
    200193      //        }
    201194    }
    … …  
    209202  float angley=dvl.GetD("ANGLEY");
    210203  float angle = dvl.GetD("ANGLE");
    211   //cout << " theta0= " << theta0 << endl;
    212   //cout << " phi0= " << phi0 << endl;
    213   //cout << " x0= " << x0 << endl;
    214   //cout << " y0= " << y0 << endl;
    215   //cout << " angle= " << angle << endl;
    216   //cout << " anglex= " << anglex << endl;
    217   //cout << " angley= " << angley << endl;
    218204
    219205  // number of components
    … …  
    226212      cout << "  expected " <<   lcm.Size_y()  << endl;
    227213      lcm.ReSize(nbrows,nbcols);
    228       cout << " resize the map to x-size=  " << nbrows  << " y-size= " << nbcols << endl;
     214      cout << " resizing the map to x-size=  " << nbrows  << " y-size= " << nbcols << endl;
    229215    }
    230216  lcm.SetSize(anglex, angley);
    … …  
    233219  for (int j=0; j< nbcols; j++)
    234220    for (int i = 0; i < nbrows; i++)  lcm(i,j) = dtab_[ij++];
    235   //cout << " chargement termine " << endl;
    236 }
    237 
    238 void FitsIoServer::printerror(int status) const
    239 {
    240   //Print out cfitsio error messages and exit program
    241   if( status )
    242     {
    243       //print error report 
    244       fits_report_error(stderr, status);
    245       //terminate the program, returning error status
    246       exit( status );
    247     }
    248   return;
    249 }
     221}
     222
     223
    250224void FitsIoServer::save( TMatrix<double>& mat, char filename[])   
    251225{
    … …  
    253227  int nbcols = mat.NCols();
    254228  long naxis = nbcols > 1 ? 2 : 1;
    255   //cout << " nombre de lignes : " << nbrows << " colonnes " << nbcols << endl;
     229  cout << " nombre de lignes : " << nbrows << " colonnes " << nbcols << endl;
    256230  FITS_tab_typ_ = TDOUBLE;
    257231  if (dtab_ != NULL ) delete[] dtab_;
    … …  
    278252  DVList dvl;
    279253
    280   //  cout << " save  : debut ecriture " << endl;
    281254  planck_write_img(fptr, naxis, nbrows, nbcols, 0, dvl);
    282   //  cout << " save  : fin ecriture " << endl;
    283255
    284256
    … …  
    326298  if( status ) printerror( status );
    327299
    328   char comment[] = "Generated by synfast";
    329   char svalue [] = "SYNFAST ";
    330   strcat(svalue, version);
    331   fits_write_key(fptr, TSTRING, "ORIG-MAP", svalue, comment, &status);
    332   if( status )  printerror( status );
    333300
    334301  strcpy(comment,"HEALPIX Pixelisation");
    … …  
    343310  if( status )  printerror( status );
    344311
    345   strcpy(comment,"Maximum multipole l used in map synthesis");
    346   int lmax= sph.nlmax();
    347   fits_write_key(fptr, TINT, "MAX-LPOL", &lmax, comment, &status);
    348   if( status )  printerror( status );
    349312
    350313  strcpy(comment,"Random generator seed");
    … …  
    353316  if( status )  printerror( status );
    354317
    355   strcpy(comment,"FWHM in DEGREES");
    356   float fwhm_deg= sph.fwhm()/60.0;
    357   fits_write_key(fptr, TFLOAT, "FWHM", &fwhm_deg, comment, &status);
    358   if( status )  printerror( status );
    359318
    360319  strcpy(comment,"--------------------------------------------");
    … …  
    370329  if( status ) printerror( status );
    371330 
    372   //   write information about the power spectrum
    373   char card[]="";
    374   sph.powfile(card);
    375   strcpy(comment,"Input power spectrum in : ");
    376   strcat(comment, card);
    377   fits_write_comment(fptr, comment, &status);
    378   if( status ) printerror( status );
    379 
    380   strcpy(comment,"Quadrupole :");
    381   fits_write_comment(fptr, comment, &status);
    382   if( status ) printerror( status );
    383 
    384   strcpy(comment,"C(2)= ");
    385   sprintf(card,"%g",sph.quadrupole());
    386   strcat(comment, card);
    387   fits_write_comment(fptr, comment, &status);
    388   if( status ) printerror( status );
    389 
    390   strcpy(comment,"keywords relative to input power spectrum ");
    391   fits_write_comment(fptr, comment, &status);
    392   if( status ) printerror( status );
    393   strcpy(comment,"have not yet been defined for PLANCK Surveyor Project");
    394   fits_write_comment(fptr, comment, &status);
    395   if( status ) printerror( status );
    396331  */
    397332
    … …  
    467402  dvl["NSIDE"] = locm.SizeIndex();
    468403  dvl["ORDERING"]=locm.TypeOfMap();
    469   float theta0;
    470   float phi0;
     404  double theta0;
     405  double phi0;
    471406  int x0;
    472407  int y0;
    473   float angle;
    474   locm.Origin(theta0, phi0,x0, y0, angle);
    475   float anglex;
    476   float angley;
    477   locm.Aperture(anglex, angley);
     408  double angle;
     409  locm.Origin(theta0,phi0,x0,y0,angle);
     410  double anglex;
     411  double angley;
     412  locm.Aperture(anglex,angley);
    478413  dvl["THETA0"] = theta0;
    479414  dvl["PHI0"] = phi0;
    … …  
    491426
    492427}
    493 
    494 
    495428void FitsIoServer::planck_write_img( fitsfile *fptr, int naxis,int n1, int n2, int n3, DVList& dvl)
    496429{
    … …  
    500433  if ( FITS_tab_typ_ == TDOUBLE) bitpix = DOUBLE_IMG;
    501434  if ( FITS_tab_typ_ == TFLOAT)  bitpix = FLOAT_IMG;
     435  if ( FITS_tab_typ_ == TINT)   bitpix = LONG_IMG;
    502436  long naxes[3];
    503437  naxes[0] = n1;
    … …  
    524458      if( status ) printerror( status );
    525459      break;
     460    case TINT :
     461      fits_write_img(fptr, TINT, 1, nelements, i_4tab_,  &status);
     462      if( status ) printerror( status );
     463      break;
    526464    default :
    527       cout << " fitsioserver : type de tableau non traite en ecriture " << endl;
     465      cout << " FitsIOServer : type de tableau non traite en ecriture " << endl;
    528466      break;
    529467    }
    … …  
    535473  // la date servira de repere pour relire ces mots cles.
    536474
    537   dvl.Print();
    538   ValList::const_iterator it;
     475  //dvl.Print();
     476  DVList::ValList::const_iterator it;
    539477  for(it = dvl.Begin(); it != dvl.End(); it++) 
    540478    {
    … …  
    568506          break;
    569507        default :
    570           cout << " probleme dans type mot cle optionnel" << endl;
     508          cout << " FitsIOServer : probleme dans type mot cle optionnel" << endl;
    571509          break;
    572510        }
    … …  
    574512    }
    575513}
     514
    576515void FitsIoServer::planck_read_img( fitsfile *fptr, long& naxis,int& n1, int& n2, int& n3, DVList& dvl)
    577516{
    … …  
    607546    {
    608547    case TDOUBLE :
     548      if (bitpix !=  DOUBLE_IMG)
     549        {
     550          cout << " FitsIOServer : the data type on fits file is not double, "
     551               << " conversion to double will be achieved by cfitsio lib " << endl;
     552        }
    609553      if (dtab_ != NULL) delete [] dtab_;
    610554      dtab_=new double[nelements];
    … …  
    614558      break;
    615559    case TFLOAT :
     560      if (bitpix !=  FLOAT_IMG)
     561        {
     562          cout << " FitsIOServer : the data type on fits file is not float, "
     563               << " conversion to float will be achieved by cfitsio lib " << endl;
     564        }
    616565      if (ftab_ != NULL) delete [] ftab_;
    617566      ftab_=new float[nelements];
    … …  
    620569      if( status )  printerror( status );
    621570      break;
     571
     572
     573    case TINT :
     574      if (bitpix !=  LONG_IMG)
     575        {
     576          cout << " FitsIOServer : the data type on fits file is not long, "
     577               << " conversion to long will be achieved by cfitsio lib " << endl;
     578        }
     579      if (i_4tab_ != NULL) delete [] i_4tab_;
     580      i_4tab_=new int_4[nelements];
     581      fits_read_img(fptr, TINT, 1, nelements,  &fnullval, i_4tab_,
     582                    &anynull, &status);
     583      if( status )  printerror( status );
     584      break;
     585
     586
    622587    default :
    623       cout << " fitsio::read_img : type non traite: " << FITS_tab_typ_  << endl;
     588      cout << " FitsIOServer::read_img : type non traite: " << FITS_tab_typ_  << endl;
    624589      break;
    625590    }
    … …  
    806771
    807772  npixels=n1;
    808   dvl.Print();
    809   ValList::const_iterator it;
     773  //dvl.Print();
     774  DVList::ValList::const_iterator it;
    810775  for(it = dvl.Begin(); it != dvl.End(); it++) 
    811776    {
    … …  
    828793          break;
    829794        default :
    830           cout << " probleme dans type mot cle optionnel" << endl;
     795          cout << " FitsIOServer : probleme dans type mot cle optionnel" << endl;
    831796          break;
    832797        }
    833798      if (!strcmp(keyname, "NSIDE") )
    834799        {
    835           //  cout << " j'ai trouve NSIDE " << endl;
    836800          nside = ival;
    837801        }
    838       // if (!strcmp(keyname, "ORDERING") )
    839       ///       {
     802      //if (!strcmp(keyname, "ORDERING") )
     803      //        {
    840804      //          cout << " j'ai trouve ORDERING " << endl;
    841805      //        }
    … …  
    854818  fitsfile *fptr; 
    855819  int bitpix = FLOAT_IMG;
    856   long naxis = 2;   
     820  long naxis = 2;
    857821
    858822  // delete old file if it already exists
    … …  
    889853}
    890854
    891 
     855void FitsIoServer::save(NTuple& ntpl,char flnm[])
     856
     857  //****************************************************/
     858  //* read the elements of the NTuple ntpl, and create */
     859  //* a FITS file with a binary table extension        */
     860  //****************************************************/
     861{
     862
     863  // delete old file if it already exists
     864  remove(flnm);
     865
     866  // pointer to the FITS file, defined in fitsio.h
     867  fitsfile *fptr;
     868  int status = 0;
     869  if( fits_create_file(&fptr, flnm, &status) )
     870    printerror( status );
     871
     872  // table will have tfields columns
     873  int tfields= ntpl.NVar();
     874
     875  // table will have nrows rows
     876  int nrows= ntpl.NEntry();
     877
     878  // extension name
     879  char extname[] = "NTuple_Binary_tbl";
     880
     881  // define the name, and the datatype for the tfields columns
     882  char **ttype, **tform;
     883  ttype= new char*[tfields];
     884  tform= new char*[tfields];
     885  for(int i = 0; i < tfields; i++)
     886    {
     887      ttype[i]= new char[FLEN_VALUE];
     888      strcpy(ttype[i], ntpl.NomIndex(i));
     889      tform[i]= new char[FLEN_VALUE];
     890      strcpy(tform[i], "1E");
     891    }
     892
     893  // create a new empty binary table onto the FITS file
     894  // physical units if they exist, are defined in the DVList object
     895  // so the null pointer is given for the tunit parameters.
     896  if ( fits_create_tbl(fptr,BINARY_TBL,nrows,tfields,ttype,tform,
     897                       NULL,extname,&status) )
     898    printerror( status );
     899
     900  // first row in table to write
     901  int firstrow  = 1;
     902
     903  // first element in row  (ignored in ASCII tables)
     904  int firstelem = 1;
     905
     906  for(int i = 0; i < tfields; i++)
     907    {
     908      float *dens= new float[nrows];
     909      for(int j = 0; j < nrows; j++)
     910        {
     911          dens[j]= ntpl.GetVal(j,i);
     912        }
     913      fits_write_col(fptr,TFLOAT,i+1,firstrow,firstelem,nrows,dens,&status);
     914      delete []dens;
     915    }
     916
     917  // number of blocks per event
     918  int blk= ntpl.BLock();
     919  fits_write_key(fptr,TINT,"BLK",&blk,"number of blocks per evt",&status);
     920
     921  // get names and values from the join DVList object
     922  DVList dvl= ntpl.Info();
     923  dvl.Print();
     924  DVList::ValList::const_iterator it;
     925  for(it = dvl.Begin(); it != dvl.End(); it++)
     926    {
     927      char keytype= (*it).second.typ;
     928      char keyname[9]="";
     929      strncpy(keyname,(*it).first.substr(0,64).c_str(),8);
     930      char comment[FLEN_COMMENT]="";
     931
     932      switch (keytype)
     933        {
     934        case 'I' :
     935          {
     936            int ival=(*it).second.mtv.iv;
     937            strcpy(comment,"I entier");
     938            fits_write_key(fptr,TINT,keyname,&ival,comment,&status);
     939            break;
     940          }
     941        case 'D' :
     942          {
     943            double dval=(*it).second.mtv.dv;
     944            strcpy(comment,"D double");
     945            fits_write_key(fptr,TDOUBLE,keyname,&dval,comment,&status);
     946            break;
     947          }
     948        case 'S' :
     949          {
     950            char strval[]="";
     951            strcpy(strval,(*it).second.mtv.strv);
     952            strcpy(comment,"S character string");
     953            fits_write_key(fptr,TSTRING,keyname,&strval,comment,&status);
     954            break;
     955          }
     956        }
     957    }
     958
     959  //close the FITS file
     960  if ( fits_close_file(fptr, &status) )
     961    printerror( status );
     962  return;
     963} 
     964
     965void FitsIoServer::load(NTuple& ntpl,char flnm[],int hdunum)
     966
     967  //********************************************************/
     968  //* move to the HDU which has the specified number hdunum*/
     969  //* in the FITS file, perform read operations and write  */
     970  //* the elements in an NTuple.                           */
     971  //********************************************************/
     972{
     973
     974  // pointer to the FITS file, defined in fitsio.h
     975  fitsfile *fptr;
     976  int status = 0;
     977  if( fits_open_file(&fptr,flnm,READONLY,&status) )
     978    printerror( status );
     979
     980  // move to the HDU
     981  int hdutype;
     982  if( fits_movabs_hdu(fptr,hdunum,&hdutype,&status) )
     983    printerror( status );
     984
     985  // get number of keywords
     986  int nkeys,keypos;
     987  if( fits_get_hdrpos(fptr,&nkeys,&keypos,&status) )
     988    printerror( status );
     989
     990  /*
     991  printf("Move to extensions by name and version number: (ffmnhd)\n");
     992  char* binname= "Test-BINTABLE";
     993  int extvers = 1;
     994  fits_movnam_hdu(fptr,ANY_HDU, binname,extvers,&status);
     995  fits_get_hdu_num(fptr,&hdunum);
     996  printf(" %s, %ld = hdu %d, %d\n", binname, extvers, hdunum, status);
     997  */
     998
     999  if (hdutype == BINARY_TBL)
     1000    printf("\nReading binary table in HDU %d:\n",hdunum);
     1001  else
     1002    {
     1003      printf("Error:: this HDU is not a binary table\n");
     1004      exit ( status );
     1005    }
     1006
     1007  // number of columns
     1008  int tfields;
     1009  if( fits_get_num_cols(fptr,&tfields,&status) )
     1010    printerror( status );
     1011
     1012  // to get table size
     1013  long naxis[2];
     1014  int nfound;
     1015   if( fits_read_keys_lng(fptr, "NAXIS", 1, 2, naxis, &nfound, &status) )
     1016     printerror( status );
     1017   int nrows= naxis[1];
     1018
     1019   printf("\nInformation about each column:\n");
     1020
     1021   //Information about each column
     1022   char **ttype, **tform;
     1023   ttype= new char*[tfields];
     1024   tform= new char*[tfields];
     1025   for( int ii = 0; ii < tfields; ii++)
     1026     {
     1027       ttype[ii]= new char[FLEN_VALUE];
     1028       tform[ii]= new char[FLEN_VALUE];
     1029     }
     1030
     1031   // number of TTYPEn,TFORMn and TUNITn keywords found
     1032   int num= 0;
     1033
     1034   // read the column names from the TTYPEn keywords
     1035   fits_read_keys_str(fptr,"TTYPE",1,tfields,ttype,&nfound,&status);
     1036   num += nfound;
     1037   // read the column names from the TFORMn keywords
     1038   fits_read_keys_str(fptr,"TFORM",1,tfields,tform,&nfound,&status);
     1039   num += nfound;
     1040
     1041   for(int ii = 0; ii < tfields; ii++)
     1042     printf("\nColumn name & Format %8s %8s", ttype[ii],tform[ii]);
     1043 
     1044   // select the columns with float data values (typecode = 42)
     1045   int typecode;
     1046   long repeat,width;
     1047   list<int> column;
     1048   for( int ii = 0; ii < tfields; ii++)
     1049     {
     1050       fits_binary_tform(tform[ii],&typecode, &repeat, &width, &status);
     1051       //printf("\n%4s %3d %2ld %2ld", tform[ii], typecode, repeat, width);
     1052       if(typecode == 42)
     1053         column.push_back(ii+1);
     1054     }
     1055
     1056   // get input elements to create the NTuple
     1057   char **clname;
     1058   clname= new char*[column.size()];
     1059   for(int ii = 0; ii < column.size(); ii++)
     1060     clname[ii]= new char[FLEN_VALUE];
     1061
     1062   list<int>::iterator itr;
     1063   int index= 0;
     1064   for(itr= column.begin(); itr != column.end(); itr++)
     1065     strcpy(clname[index++],ttype[*itr-1]);
     1066
     1067   // check if the specified keyword BLK exists
     1068   int blk= 512;
     1069   if(check_keyword(fptr,nkeys,"BLK"))
     1070     {
     1071       if( fits_read_key(fptr,TINT,"BLK",&blk,NULL,&status) )
     1072         printerror( status );
     1073     }
     1074 
     1075   // create a NTuple
     1076   NTuple nt0(column.size(),clname,blk);
     1077
     1078   float value[1];
     1079   long felem = 1;
     1080   char strnull[10];
     1081   strcpy(strnull, " ");
     1082   int anynull= 0;
     1083   long longnull  = 0;
     1084
     1085   // value to represent undefined array elements
     1086   float floatnull = FLOATNULLVALUE;
     1087
     1088   // read elements from columns and fill NTuple
     1089   for (int k = 0; k < nrows; k++)
     1090     {
     1091       int j= 0;
     1092       float *xnt= new float[column.size()];
     1093       list<int>::iterator itr;
     1094       for(itr= column.begin(); itr != column.end(); itr++)
     1095         {
     1096           fits_read_col(fptr,TFLOAT,*itr,k+1,felem,1,&floatnull,value,
     1097                         &anynull, &status);
     1098           xnt[j++]= value[0];
     1099         }
     1100       nt0.Fill(xnt);
     1101       delete[] xnt;
     1102     }
     1103
     1104   // the TUNITn keywords are optional, if they exist they are put
     1105   // in the DVLIst object
     1106   char keyname[LEN_KEYWORD]= "";
     1107   char strval[FLEN_VALUE];
     1108   for(int ii = 0; ii < tfields; ii++)
     1109     {
     1110       fits_make_keyn("TUNIT",ii+1,keyname,&status);
     1111       if(check_keyword(fptr,nkeys,keyname))
     1112         {
     1113           num++;
     1114           if( fits_read_key_str(fptr,keyname,strval,NULL,&status) )
     1115             printerror(status);
     1116           strip (keyname, 'B',' ');
     1117           strip(strval, 'B',' ');
     1118           nt0.Info()[keyname]= strval;
     1119         }
     1120     }
     1121
     1122   // add the number of mandatory keywords of a binary table
     1123   num += 8;
     1124
     1125   // put names and values of other reserved keywords in a DVList object
     1126   char comment[FLEN_COMMENT];
     1127   char dtype;
     1128   for(int j = num+1; j <= nkeys; j++)
     1129     {
     1130       fits_read_keyn(fptr,j,keyname,strval,comment,&status);
     1131       fits_get_keytype(strval,&dtype,&status);
     1132       strip (keyname, 'B',' ');
     1133       strip(strval, 'B',' ');
     1134       
     1135       switch( dtype )
     1136         {
     1137         case 'C':
     1138           nt0.Info()[keyname]= strval;   
     1139           break;
     1140         case 'I':
     1141           int ival;
     1142           ctoi(strval,&ival);
     1143           nt0.Info()[keyname]= ival;
     1144           break;
     1145         case 'F':
     1146           double dval;
     1147           ctof(strval,&dval);
     1148           nt0.Info()[keyname]= dval;
     1149           break;
     1150         }
     1151     }
     1152
     1153   // copy in the input NTuple
     1154   ntpl= nt0;
     1155
     1156   if( fits_close_file(fptr, &status) )
     1157     printerror(status);
     1158
     1159   printf("\n");
     1160   return;
     1161}
     1162
     1163bool FitsIoServer::check_keyword(fitsfile *fptr,int nkeys,char keyword[])
     1164
     1165  //*****************************************************/
     1166  //* check if the specified keyword exits in the CHU   */
     1167  //*****************************************************/
     1168{
     1169
     1170  bool KEY_EXIST = false;
     1171  int status = 0;
     1172  char strbide[FLEN_VALUE];
     1173  char keybide[LEN_KEYWORD]= "";
     1174  for(int jj = 1; jj <= nkeys; jj++)
     1175    {
     1176      if( fits_read_keyn(fptr,jj,keybide,strbide,NULL,&status) )
     1177        printerror( status );
     1178      if( !strcmp(keybide,keyword) )
     1179        {
     1180          KEY_EXIST= true;
     1181          break;
     1182        }
     1183    }
     1184  return(KEY_EXIST);
     1185}
     1186
     1187void FitsIoServer::readheader ( char filename[] )
     1188
     1189  //**********************************************************************/
     1190  //* Print out all the header keywords in all extensions of a FITS file */
     1191  //**********************************************************************/
     1192{
     1193
     1194  // standard string lengths defined in fitsioc.h
     1195  char card[FLEN_CARD];
     1196
     1197  // pointer to the FITS file, defined in fitsio.h
     1198  fitsfile *fptr;
     1199
     1200  int status = 0;
     1201  if ( fits_open_file(&fptr, filename, READONLY, &status) )
     1202    printerror( status );
     1203 
     1204  // attempt to move to next HDU, until we get an EOF error
     1205  int hdutype;
     1206  for (int ii = 1; !(fits_movabs_hdu(fptr,ii,&hdutype,&status));ii++)
     1207    {
     1208      if (hdutype == ASCII_TBL)
     1209        printf("\nReading ASCII table in HDU %d:\n", ii);
     1210      else if (hdutype == BINARY_TBL)
     1211        printf("\nReading binary table in HDU %d:\n", ii);
     1212      else if (hdutype == IMAGE_HDU)
     1213        printf("\nReading FITS image in HDU %d:\n", ii);
     1214      else
     1215        {
     1216          printf("Error: unknown type of this HDU \n");
     1217          printerror( status );
     1218        }
     1219
     1220      // get the number of keywords
     1221      int nkeys, keypos;
     1222      if ( fits_get_hdrpos(fptr, &nkeys, &keypos, &status) )
     1223        printerror( status );
     1224     
     1225      printf("Header listing for HDU #%d:\n", ii);
     1226      for (int jj = 1; jj <= nkeys; jj++) 
     1227        {
     1228          if ( fits_read_record(fptr, jj, card, &status) )
     1229            printerror( status );
     1230
     1231          // print the keyword card
     1232          printf("%s\n", card);
     1233        }
     1234      printf("END\n\n");
     1235    }
     1236
     1237  // got the expected EOF error; reset = 0
     1238  if (status == END_OF_FILE)   
     1239    status = 0;             
     1240  else
     1241    printerror( status );
     1242
     1243  if ( fits_close_file(fptr, &status) )
     1244    printerror( status );
     1245 
     1246  return;
     1247}
     1248
     1249void FitsIoServer::printerror(int status) const
     1250
     1251  //*****************************************************/
     1252  //* Print out cfitsio error messages and exit program */
     1253  //*****************************************************/
     1254{
     1255
     1256  // print out cfitsio error messages and exit program
     1257  if( status )
     1258    {
     1259      // print error report 
     1260      fits_report_error(stderr, status);
     1261      // terminate the program, returning error status
     1262      exit( status );
     1263    }
     1264  return;
     1265}
     1266
     1267void FitsIoServer::save(ImageR4& DpcImg,char flnm[])
     1268{
     1269  long naxis=2;
     1270  int siz_x = DpcImg.XSize();
     1271  int siz_y = DpcImg.YSize();
     1272  int nbpixels = siz_x*siz_y;
     1273  FITS_tab_typ_ = TFLOAT;
     1274
     1275  // pointer to the FITS file, defined in fitsio.h
     1276  fitsfile *fptr;
     1277
     1278  // initialize status before calling fitsio routines
     1279  int status = 0;         
     1280
     1281  // delete old file if it already exists
     1282  remove(flnm);
     1283
     1284  // create new FITS file 
     1285  fits_create_file(&fptr, flnm, &status);
     1286  if( status )  printerror( status );
     1287
     1288  // get the values of the DpcImage
     1289  if (ftab_ != NULL) delete [] ftab_;
     1290  ftab_=DpcImg.ImagePtr();
     1291
     1292  // write FITS image
     1293  DVList dvl;
     1294
     1295  dvl["DATATYPE"] = int(DpcImg.PixelType());
     1296  dvl["ORG_X"] = DpcImg.XOrg();
     1297  dvl["ORG_Y"] = DpcImg.YOrg();
     1298  dvl["PIXSZ_X"] = DpcImg.XPxSize();
     1299  dvl["PIXSZ_Y"] = DpcImg.YPxSize();
     1300  dvl["IDENT"] = DpcImg.Ident();
     1301  dvl["NAME"] = DpcImg.Nom();
     1302  dvl["NBSAT"] = DpcImg.nbSat;
     1303  dvl["NBNUL"] = DpcImg.nbNul;
     1304  dvl["MINPIX"] = DpcImg.minPix;
     1305  dvl["MAXPIX"] = DpcImg.maxPix;
     1306  dvl["MOYPIX"] = DpcImg.moyPix;
     1307  dvl["SIGPIX"] = DpcImg.sigPix;
     1308  dvl["FOND"] = DpcImg.fond;
     1309  dvl["SIGFON"] = DpcImg.sigmaFond;
     1310
     1311
     1312  planck_write_img(fptr, naxis, siz_x, siz_y, 0, dvl);
     1313
     1314  // close the file
     1315  fits_close_file(fptr, &status);
     1316  if( status )  printerror( status );
     1317
     1318  ftab_=NULL;
     1319}
     1320
     1321void FitsIoServer::load(ImageR4& DpcImg,char flnm[])
     1322{
     1323  FITS_tab_typ_ = TFLOAT;
     1324  long naxis;
     1325  int siz_x;
     1326  int siz_y;
     1327  int n3;
     1328  DVList dvl;
     1329
     1330  // pointer to the FITS file, defined in fitsio.h
     1331  fitsfile *fptr;
     1332  // initialize status before calling fitsio routines
     1333  int status = 0;
     1334       
     1335  fits_open_file(&fptr, flnm, READONLY, &status);
     1336  if( status )  printerror( status );
     1337
     1338  planck_read_img(fptr, naxis, siz_x, siz_y, n3, dvl);
     1339
     1340  // close the file
     1341  fits_close_file(fptr, &status);
     1342  if( status )  printerror( status );
     1343 
     1344  if (naxis != 2)
     1345    {
     1346      cout << " FitsIOServer : naxis= " << naxis << " not equal to 2 ?" << endl;
     1347    }
     1348
     1349  DpcImg.Allocate(siz_x, siz_y, 0., 0);
     1350  float* pixelPtr=DpcImg.ImagePtr();
     1351  //for (int i=0; i < siz_x*siz_y; i++) pixelPtr[i]=ftab_[i];
     1352  // verifications de type
     1353  PBaseDataTypes  dataT=DataType((r_4)0);
     1354  int TypeDonnees=dvl.GetI("DATATYPE");
     1355  if (int(dataT) !=  TypeDonnees)
     1356    {
     1357      cout << " FitsIOServer : parameter DATATYPE on file is not float " << endl;
     1358      cout << "  eventual conversion to float is achieved by cfitsio lib " << endl;
     1359    }
     1360
     1361  memcpy(pixelPtr, ftab_, siz_x*siz_y*DataSize(dataT));
     1362  const char* nom=dvl.GetS("NAME").c_str();
     1363  DpcImg.SetNameId(dvl.GetI("IDENT"), nom);
     1364  DpcImg.SetOrg(dvl.GetI("ORG_X"), dvl.GetI("ORG_Y"));
     1365  DpcImg.SetPxSize(dvl.GetD("PIXSZ_X"), dvl.GetD("PIXSZ_Y"));   
     1366  DpcImg.SetAtt(dvl.GetI("NBNUL"), dvl.GetI("NBSAT"),
     1367       dvl.GetD("MINPIX"), dvl.GetD("MAXPIX"), dvl.GetD("MOYPIX"),
     1368       dvl.GetD("SIGPIX"),dvl.GetD("FOND"), dvl.GetD("SIGFON"));
     1369
     1370}
     1371
     1372void FitsIoServer::save(ImageI4& DpcImg,char flnm[])
     1373{
     1374  long naxis=2;
     1375  int siz_x = DpcImg.XSize();
     1376  int siz_y = DpcImg.YSize();
     1377  int nbpixels = siz_x*siz_y;
     1378  FITS_tab_typ_ = TINT;
     1379
     1380  // pointer to the FITS file, defined in fitsio.h
     1381  fitsfile *fptr;
     1382
     1383  // initialize status before calling fitsio routines
     1384  int status = 0;         
     1385
     1386  // delete old file if it already exists
     1387  remove(flnm);
     1388
     1389  // create new FITS file 
     1390  fits_create_file(&fptr, flnm, &status);
     1391  if( status )  printerror( status );
     1392
     1393  // get the values of the DpcImage
     1394  if (i_4tab_ != NULL) delete [] i_4tab_;
     1395  i_4tab_=DpcImg.ImagePtr();
     1396
     1397  // write FITS image
     1398  DVList dvl;
     1399
     1400  dvl["DATATYPE"] = int(DpcImg.PixelType());
     1401  dvl["ORG_X"] = DpcImg.XOrg();
     1402  dvl["ORG_Y"] = DpcImg.YOrg();
     1403  dvl["PIXSZ_X"] = DpcImg.XPxSize();
     1404  dvl["PIXSZ_Y"] = DpcImg.YPxSize();
     1405  dvl["IDENT"] = DpcImg.Ident();
     1406  dvl["NAME"] = DpcImg.Nom();
     1407  dvl["NBSAT"] = DpcImg.nbSat;
     1408  dvl["NBNUL"] = DpcImg.nbNul;
     1409  dvl["MINPIX"] = DpcImg.minPix;
     1410  dvl["MAXPIX"] = DpcImg.maxPix;
     1411  dvl["MOYPIX"] = DpcImg.moyPix;
     1412  dvl["SIGPIX"] = DpcImg.sigPix;
     1413  dvl["FOND"] = DpcImg.fond;
     1414  dvl["SIGFON"] = DpcImg.sigmaFond;
     1415
     1416
     1417  planck_write_img(fptr, naxis, siz_x, siz_y, 0, dvl);
     1418
     1419  // close the file
     1420  fits_close_file(fptr, &status);
     1421  if( status )  printerror( status );
     1422
     1423  i_4tab_=NULL;
     1424}
     1425
     1426void FitsIoServer::load(ImageI4& DpcImg,char flnm[])
     1427{
     1428  FITS_tab_typ_ = TINT;
     1429  long naxis;
     1430  int siz_x;
     1431  int siz_y;
     1432  int n3;
     1433  DVList dvl;
     1434
     1435  // pointer to the FITS file, defined in fitsio.h
     1436  fitsfile *fptr;
     1437  // initialize status before calling fitsio routines
     1438  int status = 0;
     1439       
     1440  fits_open_file(&fptr, flnm, READONLY, &status);
     1441  if( status )  printerror( status );
     1442  planck_read_img(fptr, naxis, siz_x, siz_y, n3, dvl);
     1443
     1444  // close the file
     1445  fits_close_file(fptr, &status);
     1446  if( status )  printerror( status );
     1447 
     1448  if (naxis != 2)
     1449    {
     1450      cout << " FitsIOServer : naxis= " << naxis << " not equal to 2 ?" << endl;
     1451    }
     1452
     1453  DpcImg.Allocate(siz_x, siz_y, 0., 0);
     1454  int_4* pixelPtr=DpcImg.ImagePtr();
     1455
     1456
     1457  // verifications de type
     1458  PBaseDataTypes  dataT=DataType((int_4)0);
     1459  int TypeDonnees=dvl.GetI("DATATYPE");
     1460  if (int(dataT) !=  TypeDonnees)
     1461    {
     1462      cout << " FitsIOServer : parameter DATATYPE on file is not int_4 " << endl;
     1463      cout << "  eventual conversion to float is achieved by cfitsio lib " << endl;
     1464    }
     1465  memcpy(pixelPtr, i_4tab_, siz_x*siz_y*DataSize(dataT));
     1466  const char* nom=dvl.GetS("NAME").c_str();
     1467  DpcImg.SetNameId(dvl.GetI("IDENT"), nom);
     1468  DpcImg.SetOrg(dvl.GetI("ORG_X"), dvl.GetI("ORG_Y"));
     1469  DpcImg.SetPxSize(dvl.GetD("PIXSZ_X"), dvl.GetD("PIXSZ_Y"));   
     1470  DpcImg.SetAtt(dvl.GetI("NBNUL"), dvl.GetI("NBSAT"),
     1471       dvl.GetD("MINPIX"), dvl.GetD("MAXPIX"), dvl.GetD("MOYPIX"),
     1472       dvl.GetD("SIGPIX"),dvl.GetD("FOND"), dvl.GetD("SIGFON"));
     1473}
     1474
     1475
Note: See TracChangeset for help on using the changeset viewer.