source: Sophya/trunk/SophyaExt/FitsIOServer/fitsioserver.cc@ 531

Last change on this file since 531 was 531, checked in by ansari, 26 years ago

dvl[NSIDE] = (int_4) au lieu de dvl[NSIDE] = (long) 28-OCT-99 GLM
File size: 36.5 KB
Line 
1//************************************************************************
2// Class for loadind and saving from FITS-formatted file to DPC objects
3// (G. Le Meur ; Francois Touze) OCT. 99
4//
5// methods 'load(X& x, char f[])' get from FITS file "f" a DPC object x
6// from DPC class X.
7// methods 'save(X& x, char f[])' save a DPC object x from DPC // class X
8// onto a FITS file "f" .
9
10//************************************************************************
11
12#include <iostream.h>
13#include <list>
14
15#include "fitsioserver.h"
16#include "strutil.h"
17
18void FitsIoServer::load(TMatrix<double>& mat,char flnm[])
19{
20 int nbrows=0;
21 int nbcols=0;
22 FITS_tab_typ_ = TDOUBLE;
23 long naxis;
24 int n1, n2, n3;
25 DVList dvl;
26 planck_read_img(flnm, naxis, n1, n2, n3, dvl);
27
28 nbrows=n1;
29 nbcols=n2;
30 if (naxis == 1) nbcols=1;
31 if (naxis > 2)
32 {
33 cout<<" FitsIOServer : le fichier fits n'est pas une matrice, naxis= " <<naxis<< endl;
34 }
35
36 // number of components
37 if (mat.NRows() != nbrows || mat.NCols() != nbcols )
38 {
39 cout << " found " << nbrows << " rows ";
40 cout << " expected " << mat.NRows() << endl;
41 cout << " found " << nbcols << " columns " ;
42 cout << " expected " << mat.NCols() << endl;
43 mat.ReSize(nbrows,nbcols);
44 cout << " resize the vector to nbrows= " << nbrows << " nbcols= " << nbcols << endl;
45 }
46 int ij=0;
47 for (int j=0; j< nbcols; j++)
48 for (int i = 0; i < nbrows; i++) mat(i,j) = (double)r_8tab_[ij++];
49}
50
51void FitsIoServer::load(NTuple& ntpl,char flnm[],int hdunum)
52
53 //********************************************************/
54 //* move to the HDU which has the specified number hdunum*/
55 //* in the FITS file, perform read operations and write */
56 //* the elements in an NTuple. */
57 //********************************************************/
58{
59
60 // pointer to the FITS file, defined in fitsio.h
61 fitsfile *fptr;
62 int status = 0;
63 if( fits_open_file(&fptr,flnm,READONLY,&status) )
64 printerror( status );
65
66 // move to the HDU
67 int hdutype;
68 if( fits_movabs_hdu(fptr,hdunum,&hdutype,&status) )
69 printerror( status );
70
71 // get number of keywords
72 int nkeys,keypos;
73 if( fits_get_hdrpos(fptr,&nkeys,&keypos,&status) )
74 printerror( status );
75
76
77 if (hdutype == BINARY_TBL)
78 printf("\nReading binary table in HDU %d:\n",hdunum);
79 else
80 {
81 printf("Error:: this HDU is not a binary table\n");
82 exit ( status );
83 }
84
85 // number of columns
86 int tfields;
87 if( fits_get_num_cols(fptr,&tfields,&status) )
88 printerror( status );
89
90 // to get table size
91 long naxis[2];
92 int nfound;
93 if( fits_read_keys_lng(fptr, "NAXIS", 1, 2, naxis, &nfound, &status) )
94 printerror( status );
95 int nrows= naxis[1];
96
97 //Information about each column
98 char **ttype, **tform;
99 ttype= new char*[tfields];
100 tform= new char*[tfields];
101 for( int ii = 0; ii < tfields; ii++)
102 {
103 ttype[ii]= new char[FLEN_VALUE];
104 tform[ii]= new char[FLEN_VALUE];
105 }
106
107 // number of TTYPEn,TFORMn and TUNITn keywords found
108 int num= 0;
109
110 // read the column names from the TTYPEn keywords
111 fits_read_keys_str(fptr,"TTYPE",1,tfields,ttype,&nfound,&status);
112 num += nfound;
113 // read the column names from the TFORMn keywords
114 fits_read_keys_str(fptr,"TFORM",1,tfields,tform,&nfound,&status);
115 num += nfound;
116
117 //for(int ii = 0; ii < tfields; ii++)
118 // printf("\nColumn name & Format %8s %8s", ttype[ii],tform[ii]);
119
120 // select the columns with float data values
121 int typecode;
122 long repeat,width;
123 list<int> column;
124 for( int ii = 0; ii < tfields; ii++)
125 {
126 fits_binary_tform(tform[ii],&typecode, &repeat, &width, &status);
127 //printf("\n%4s %3d %2ld %2ld", tform[ii], typecode, repeat, width);
128 if(typecode == TFLOAT)
129 column.push_back(ii+1);
130 }
131
132 // get input elements to create the NTuple
133 char **clname;
134 clname= new char*[column.size()];
135 for(int ii = 0; ii < column.size(); ii++)
136 clname[ii]= new char[FLEN_VALUE];
137
138 list<int>::iterator itr;
139 int index= 0;
140 for(itr= column.begin(); itr != column.end(); itr++)
141 strcpy(clname[index++],ttype[*itr-1]);
142
143 for( int ii = 0; ii < tfields; ii++)
144 {
145 delete [] ttype[ii];
146 delete [] tform[ii];
147 }
148 delete [] ttype;
149 delete [] tform;
150
151
152 // check if the specified keyword BLK exists
153 int blk= 512;
154 if(check_keyword(fptr,nkeys,"BLK"))
155 {
156 if( fits_read_key(fptr,TINT,"BLK",&blk,NULL,&status) )
157 printerror( status );
158 }
159 // create a NTuple
160 NTuple nt0(column.size(),clname,blk);
161
162 for(int ii = 0; ii < column.size(); ii++)
163 delete [] clname[ii];
164 delete [] clname;
165
166 float value[1];
167 long felem = 1;
168 char strnull[10];
169 strcpy(strnull, " ");
170 int anynull= 0;
171 long longnull = 0;
172
173 // value to represent undefined array elements
174 float floatnull = FLOATNULLVALUE;
175
176 // read elements from columns and fill NTuple
177 for (int k = 0; k < nrows; k++)
178 {
179 int j= 0;
180 float *xnt= new float[column.size()];
181 list<int>::iterator itr;
182 for(itr= column.begin(); itr != column.end(); itr++)
183 {
184 fits_read_col(fptr,TFLOAT,*itr,k+1,felem,1,&floatnull,value,
185 &anynull, &status);
186 xnt[j++]= value[0];
187 }
188 nt0.Fill(xnt);
189 delete[] xnt;
190 }
191
192 // the TUNITn keywords are optional, if they exist they are put
193 // in the DVLIst object
194 char keyname[LEN_KEYWORD]= "";
195 char strval[FLEN_VALUE];
196 for(int ii = 0; ii < tfields; ii++)
197 {
198 fits_make_keyn("TUNIT",ii+1,keyname,&status);
199 if(check_keyword(fptr,nkeys,keyname))
200 {
201 num++;
202 if( fits_read_key_str(fptr,keyname,strval,NULL,&status) )
203 printerror(status);
204 strip (keyname, 'B',' ');
205 strip(strval, 'B',' ');
206 nt0.Info()[keyname]= strval;
207 }
208 }
209
210 // add the number of mandatory keywords of a binary table
211 num += 8;
212
213 // put names and values of other reserved keywords in a DVList object
214 char comment[FLEN_COMMENT];
215 char dtype;
216 for(int j = num+1; j <= nkeys; j++)
217 {
218 fits_read_keyn(fptr,j,keyname,strval,comment,&status);
219 fits_get_keytype(strval,&dtype,&status);
220 strip (keyname, 'B',' ');
221 strip(strval, 'B',' ');
222
223 switch( dtype )
224 {
225 case 'C':
226 nt0.Info()[keyname]= strval;
227 break;
228 case 'I':
229 int ival;
230 ctoi(strval,&ival);
231 nt0.Info()[keyname]= (int_4)ival;
232 break;
233 case 'F':
234 double dval;
235 ctof(strval,&dval);
236 nt0.Info()[keyname]= dval;
237 break;
238 }
239 }
240
241 // copy in the input NTuple
242 ntpl= nt0;
243
244 if( fits_close_file(fptr, &status) )
245 printerror(status);
246
247 printf("\n");
248 return;
249}
250
251
252void FitsIoServer::load(SphericalMap<double>& sph, char flnm[])
253{
254 int npixels=0;
255 int nside=0;
256 long naxis;
257 int n1, n2, n3;
258
259 FITS_tab_typ_ = TDOUBLE;
260
261 DVList dvl;
262 planck_read_img(flnm, naxis, n1, n2, n3, dvl);
263 if (naxis != 1)
264 {
265 cout << " le fichier fits n'est pas une sphere, naxis= " << naxis << endl;
266 }
267 npixels=n1;
268 nside= dvl.GetI("NSIDE");
269
270 // number of pixels in the sphere
271 if (sph.NbPixels() != npixels)
272 {
273 cout << " found " << npixels << " pixels" << endl;
274 cout << " expected " << sph.NbPixels() << endl;
275 if (nside <= 0 )
276 {
277 cout<<" FITSIOSERVER: no resolution parameter on fits file "<<endl;
278 exit(0);
279 }
280 if (nside != sph.SizeIndex())
281 {
282 sph.Resize(nside);
283 cout << " resizing the sphere to nside= " << nside << endl;
284 }
285 else
286 {
287 cout << " FITSIOSERVER : same resolution, surprising ! " << endl;
288 exit(0);
289 }
290 }
291 for (int j = 0; j < sph.NbPixels(); j++) sph(j)= (double)r_8tab_[j];
292}
293void FitsIoServer::load(SphericalMap<float>& sph, char flnm[])
294{
295 int npixels=0;
296 int nside=0;
297 long naxis;
298 int n1, n2, n3;
299 DVList dvl;
300
301 FITS_tab_typ_ = TFLOAT;
302
303 planck_read_img(flnm, naxis, n1, n2, n3, dvl);
304 if (naxis != 1)
305 {
306 cout << " le fichier fits n'est pas une sphere, naxis= " << naxis << endl;
307 }
308 npixels=n1;
309 nside= dvl.GetI("NSIDE");
310
311 // number of pixels in the sphere
312 if (sph.NbPixels() != npixels)
313 {
314 cout << " found " << npixels << " pixels" << endl;
315 cout << " expected " << sph.NbPixels() << endl;
316 if (nside <= 0 )
317 {
318 cout<<" FITSIOSERVER: no resolution parameter on fits file "<<endl;
319 exit(0);
320 }
321 if (nside != sph.SizeIndex())
322 {
323 sph.Resize(nside);
324 cout << " resizing the sphere to nside= " << nside << endl;
325 }
326 else
327 {
328 cout << " FITSIOSERVER : same resolution, surprising ! " << endl;
329 exit(0);
330 }
331 }
332 for (int j = 0; j < sph.NbPixels(); j++) sph(j)= (float)r_4tab_[j];
333}
334
335void FitsIoServer::load(LocalMap<double>& lcm, char flnm[])
336{
337 int nbrows=0;
338 int nbcols=0;
339 FITS_tab_typ_ = TDOUBLE;
340 long naxis;
341 int n1, n2, n3;
342 DVList dvl;
343 planck_read_img(flnm, naxis, n1, n2, n3, dvl);
344
345 nbrows=n1;
346 nbcols=n2;
347 if (naxis != 2)
348 {
349 cout<<" FitsIOServer : le fichier fits n'est pas une localmap, naxis= "<<naxis<<endl;
350 }
351 float theta0 = dvl.GetD("THETA0");
352 float phi0 = dvl.GetD("PHI0");
353 int x0 = dvl.GetI("X0");
354 int y0 = dvl.GetI("Y0");
355 int xo = dvl.GetI("XO");
356 int yo = dvl.GetI("YO");
357 float anglex=dvl.GetD("ANGLEX");
358 float angley=dvl.GetD("ANGLEY");
359 float angle = dvl.GetD("ANGLE");
360
361 // number of components
362
363 if (lcm.Size_x() != nbrows || lcm.Size_y() != nbcols )
364 {
365 cout << " found " << nbrows << " x-pixels ";
366 cout << " expected " << lcm.Size_x() << endl;
367 cout << " found " << nbcols << " y-pixels " ;
368 cout << " expected " << lcm.Size_y() << endl;
369 lcm.ReSize(nbrows,nbcols);
370 cout << " resizing the map to x-size= " << nbrows << " y-size= " << nbcols << endl;
371 }
372 lcm.SetSize(anglex, angley);
373 lcm.SetOrigin(theta0, phi0, x0, y0, angle);
374 int ij=0;
375 for (int j=0; j< nbcols; j++)
376 for (int i = 0; i < nbrows; i++) lcm(i,j) = (double)r_8tab_[ij++];
377}
378
379void FitsIoServer::load(ImageR4& DpcImg,char flnm[])
380{
381 FITS_tab_typ_ = TFLOAT;
382 long naxis;
383 int siz_x;
384 int siz_y;
385 int n3;
386 DVList dvl;
387
388 planck_read_img(flnm, naxis, siz_x, siz_y, n3, dvl);
389
390
391 if (naxis != 2)
392 {
393 cout << " FitsIOServer : naxis= " << naxis << " not equal to 2 ?" << endl;
394 }
395
396 DpcImg.Allocate(siz_x, siz_y, 0., 0);
397 float* pixelPtr=DpcImg.ImagePtr();
398 // verifications de type
399 PBaseDataTypes dataT=DataType((r_4)0);
400 int TypeDonnees=dvl.GetI("DATATYPE");
401 if (int(dataT) != TypeDonnees)
402 {
403 cout << " FitsIOServer : parameter DATATYPE on file " << flnm << " is not float " << endl;
404 cout << " eventual conversion to float is achieved by cfitsio lib " << endl;
405 }
406
407 memcpy(pixelPtr, r_4tab_, siz_x*siz_y*DataSize(dataT));
408 const char* nom=dvl.GetS("NAME").c_str();
409 DpcImg.SetNameId(dvl.GetI("IDENT"), nom);
410 DpcImg.SetOrg(dvl.GetI("ORG_X"), dvl.GetI("ORG_Y"));
411 DpcImg.SetPxSize(dvl.GetD("PIXSZ_X"), dvl.GetD("PIXSZ_Y"));
412 DpcImg.SetAtt(dvl.GetI("NBNUL"), dvl.GetI("NBSAT"),
413 dvl.GetD("MINPIX"), dvl.GetD("MAXPIX"), dvl.GetD("MOYPIX"),
414 dvl.GetD("SIGPIX"),dvl.GetD("FOND"), dvl.GetD("SIGFON"));
415
416}
417
418
419void FitsIoServer::load(ImageI4& DpcImg,char flnm[])
420{
421 FITS_tab_typ_ = TINT;
422 long naxis;
423 int siz_x;
424 int siz_y;
425 int n3;
426 DVList dvl;
427
428 // pointer to the FITS file, defined in fitsio.h
429 //fitsfile *fptr;
430 // initialize status before calling fitsio routines
431 //int status = 0;
432
433 //fits_open_file(&fptr, flnm, READONLY, &status);
434 //if( status ) printerror( status );
435 planck_read_img(flnm, naxis, siz_x, siz_y, n3, dvl);
436 // close the file
437 //fits_close_file(fptr, &status);
438 //if( status ) printerror( status );
439
440 if (naxis != 2)
441 {
442 cout << " FitsIOServer : naxis= " << naxis << " not equal to 2 ?" << endl;
443 }
444
445 DpcImg.Allocate(siz_x, siz_y, 0., 0);
446 int_4* pixelPtr=DpcImg.ImagePtr();
447
448
449 // verifications de type
450 PBaseDataTypes dataT=DataType((int_4)0);
451 int TypeDonnees=dvl.GetI("DATATYPE");
452 if (int(dataT) != TypeDonnees)
453 {
454 cout << " FitsIOServer : parameter DATATYPE on file " << flnm << " is not int_4 " << endl;
455 cout << " eventual conversion to int_4 is achieved by cfitsio lib " << endl;
456 }
457 memcpy(pixelPtr, i_4tab_, siz_x*siz_y*DataSize(dataT));
458 const char* nom=dvl.GetS("NAME").c_str();
459 DpcImg.SetNameId(dvl.GetI("IDENT"), nom);
460 DpcImg.SetOrg(dvl.GetI("ORG_X"), dvl.GetI("ORG_Y"));
461 DpcImg.SetPxSize(dvl.GetD("PIXSZ_X"), dvl.GetD("PIXSZ_Y"));
462 DpcImg.SetAtt(dvl.GetI("NBNUL"), dvl.GetI("NBSAT"),
463 dvl.GetD("MINPIX"), dvl.GetD("MAXPIX"), dvl.GetD("MOYPIX"),
464 dvl.GetD("SIGPIX"),dvl.GetD("FOND"), dvl.GetD("SIGFON"));
465}
466
467
468void FitsIoServer::save( TMatrix<double>& mat, char filename[])
469{
470 int nbrows = mat.NRows();
471 int nbcols = mat.NCols();
472 long naxis = nbcols > 1 ? 2 : 1;
473 cout << " nombre de lignes : " << nbrows << " colonnes " << nbcols << endl;
474 FITS_tab_typ_ = TDOUBLE;
475 if (r_8tab_ != NULL ) delete[] r_8tab_;
476 r_8tab_=new r_8[nbrows*nbcols];
477
478
479 int ij=0;
480 for (int j=0; j< nbcols; j++)
481 for (int i = 0; i < nbrows; i++) r_8tab_[ij++]= (r_8)mat(i,j);
482
483 DVList dvl;
484
485 planck_write_img(filename, naxis, nbrows, nbcols, 0, dvl);
486 delete[] r_8tab_;
487}
488
489void FitsIoServer::save(NTuple& ntpl,char flnm[])
490
491 //****************************************************/
492 //* read the elements of the NTuple ntpl, and create */
493 //* a FITS file with a binary table extension */
494 //****************************************************/
495{
496
497 // delete old file if it already exists
498 remove(flnm);
499
500 // pointer to the FITS file, defined in fitsio.h
501 fitsfile *fptr;
502 int status = 0;
503 if( fits_create_file(&fptr, flnm, &status) )
504 printerror( status );
505
506 // table will have tfields columns
507 int tfields= ntpl.NVar();
508
509 // table will have nrows rows
510 int nrows= ntpl.NEntry();
511
512 // extension name
513 char extname[] = "NTuple_Binary_tbl";
514
515 // define the name, and the datatype for the tfields columns
516 char **ttype, **tform;
517 ttype= new char*[tfields];
518 tform= new char*[tfields];
519 for(int i = 0; i < tfields; i++)
520 {
521 ttype[i]= new char[FLEN_VALUE];
522 strcpy(ttype[i], ntpl.NomIndex(i));
523 tform[i]= new char[FLEN_VALUE];
524 strcpy(tform[i], "1E");
525 }
526
527 // create a new empty binary table onto the FITS file
528 // physical units if they exist, are defined in the DVList object
529 // so the null pointer is given for the tunit parameters.
530 if ( fits_create_tbl(fptr,BINARY_TBL,nrows,tfields,ttype,tform,
531 NULL,extname,&status) )
532 printerror( status );
533
534 for( int ii = 0; ii < tfields; ii++)
535 {
536 delete [] ttype[ii];
537 delete [] tform[ii];
538 }
539 delete [] ttype;
540 delete [] tform;
541
542
543 // first row in table to write
544 int firstrow = 1;
545
546 // first element in row (ignored in ASCII tables)
547 int firstelem = 1;
548
549 for(int i = 0; i < tfields; i++)
550 {
551 float *dens= new float[nrows];
552 for(int j = 0; j < nrows; j++)
553 {
554 dens[j]= ntpl.GetVal(j,i);
555 }
556 fits_write_col(fptr,TFLOAT,i+1,firstrow,firstelem,nrows,dens,&status);
557 delete []dens;
558 }
559
560 // number of blocks per event
561 int blk= ntpl.BLock();
562 fits_write_key(fptr,TINT,"BLK",&blk,"number of blocks per evt",&status);
563
564 // get names and values from the join DVList object
565 DVList dvl= ntpl.Info();
566 dvl.Print();
567 DVList::ValList::const_iterator it;
568 for(it = dvl.Begin(); it != dvl.End(); it++)
569 {
570 char keytype= (*it).second.typ;
571 char keyname[9]="";
572 strncpy(keyname,(*it).first.substr(0,64).c_str(),8);
573 char comment[FLEN_COMMENT]="";
574
575 switch (keytype)
576 {
577 case 'I' :
578 {
579 int ival=(*it).second.mtv.iv;
580 strcpy(comment,"I entier");
581 fits_write_key(fptr,TINT,keyname,&ival,comment,&status);
582 break;
583 }
584 case 'D' :
585 {
586 double dval=(*it).second.mtv.dv;
587 strcpy(comment,"D double");
588 fits_write_key(fptr,TDOUBLE,keyname,&dval,comment,&status);
589 break;
590 }
591 case 'S' :
592 {
593 char strval[]="";
594 strcpy(strval,(*it).second.mtv.strv);
595 strcpy(comment,"S character string");
596 fits_write_key(fptr,TSTRING,keyname,&strval,comment,&status);
597 break;
598 }
599 }
600 }
601
602 //close the FITS file
603 if ( fits_close_file(fptr, &status) )
604 printerror( status );
605 return;
606}
607
608
609
610
611//void FitsIoServer::save(SphericalMap<double>& sph, char filename[])
612void FitsIoServer::save(SphericalMap<double>& sph, char filename[])
613{
614 int npixels = sph.NbPixels();
615 FITS_tab_typ_ = TDOUBLE;
616 if (r_8tab_ != NULL ) delete[] r_8tab_;
617 r_8tab_=new r_8[npixels];
618
619
620 for (int j = 0; j < npixels; j++) r_8tab_[j]= (r_8)sph(j);
621 DVList dvl;
622 dvl["NSIDE"] = (int_4) sph.SizeIndex();
623 dvl["ORDERING"]=sph.TypeOfMap();
624
625 planck_write_img(filename, 1, npixels, 0, 0, dvl);
626
627 // decider ulterieurement ce qu'on fait de ce qui suit, specifique
628 // pour l'instant, aux spheres gorski.
629
630 /*
631 // write supplementary keywords
632 fits_write_comment(fptr, " ", &status);
633 if( status ) printerror( status );
634
635
636 strcpy(comment,"HEALPIX Pixelisation");
637 strcpy(svalue, "HEALPIX");
638 fits_write_key(fptr, TSTRING, "PIXTYPE", svalue, comment, &status);
639 if( status ) printerror( status );
640
641
642 strcpy(comment,"pixel ordering scheme, either RING or NESTED");
643 strcpy(svalue, "RING");
644 fits_write_key(fptr, TSTRING, "ORDERING", svalue, comment, &status);
645 if( status ) printerror( status );
646
647
648 strcpy(comment,"Random generator seed");
649 int iseed= sph.iseed();
650 fits_write_key(fptr, TINT, "RANDSEED", &iseed, comment, &status);
651 if( status ) printerror( status );
652
653
654 strcpy(comment,"--------------------------------------------");
655 fits_write_comment(fptr, comment, &status);
656 if( status ) printerror( status );
657
658 strcpy(comment," Above keywords are still likely to change !");
659 fits_write_comment(fptr, comment, &status);
660 if( status ) printerror( status );
661
662 strcpy(comment,"--------------------------------------------");
663 fits_write_comment(fptr, comment, &status);
664 if( status ) printerror( status );
665
666 */
667
668 delete[] r_8tab_;
669}
670
671void FitsIoServer::save(SphericalMap<float>& sph, char filename[])
672{
673 int npixels = sph.NbPixels();
674 FITS_tab_typ_ = TFLOAT;
675 if (r_4tab_ != NULL ) delete[] r_4tab_;
676 r_4tab_=new r_4[npixels];
677
678
679 for (int j = 0; j < npixels; j++) r_4tab_[j]= (r_4)sph(j);
680 DVList dvl;
681 dvl["NSIDE"] = (int_4)sph.SizeIndex();
682 dvl["ORDERING"]=sph.TypeOfMap();
683
684 planck_write_img(filename, 1, npixels, 0, 0, dvl);
685 delete[] r_4tab_;
686
687}
688
689void FitsIoServer::save(LocalMap<double>& locm, char filename[])
690{
691 int nbrows = locm.Size_x();
692 int nbcols = locm.Size_y();
693 long naxis = 2;
694 cout << " nombre de pts en x : " << nbrows << " en y " << nbcols << endl;
695 FITS_tab_typ_ = TDOUBLE;
696 if (r_8tab_ != NULL ) delete[] r_8tab_;
697 r_8tab_=new r_8[nbrows*nbcols];
698
699 int ij=0;
700 for (int j=0; j< nbcols; j++)
701 for (int i = 0; i < nbrows; i++) r_8tab_[ij++]= (r_8)locm(i,j);
702
703 DVList dvl;
704 dvl["NSIDE"] = (int_4) locm.SizeIndex();
705 dvl["ORDERING"]=locm.TypeOfMap();
706 double theta0;
707 double phi0;
708 int x0;
709 int y0;
710 double angle;
711 locm.Origin(theta0,phi0,x0,y0,angle);
712 double anglex;
713 double angley;
714 locm.Aperture(anglex,angley);
715 dvl["THETA0"] = theta0;
716 dvl["PHI0"] = phi0;
717 dvl["X0"] = (int_4)x0;
718 dvl["Y0"] = (int_4)y0;
719 dvl["ANGLE"] = angle;
720 dvl["ANGLEX"] = anglex;
721 dvl["ANGLEY"] = angley;
722 planck_write_img(filename, naxis, nbrows, nbcols, 0, dvl);
723 delete[] r_8tab_;}
724
725
726void FitsIoServer::save(const ImageR4& DpcImg,char flnm[])
727{
728 long naxis=2;
729 int siz_x = DpcImg.XSize();
730 int siz_y = DpcImg.YSize();
731 int nbpixels = siz_x*siz_y;
732 FITS_tab_typ_ = TFLOAT;
733
734
735
736 // write FITS image
737 DVList dvl;
738
739 dvl["DATATYPE"] = (int_4)DpcImg.PixelType();
740 dvl["ORG_X"] = DpcImg.XOrg();
741 dvl["ORG_Y"] = DpcImg.YOrg();
742 dvl["PIXSZ_X"] = DpcImg.XPxSize();
743 dvl["PIXSZ_Y"] = DpcImg.YPxSize();
744 dvl["IDENT"] = DpcImg.Ident();
745
746 //dvl["NAME"] = DpcImg.Nom();
747
748 // j utilise la methode SetS parce que ses parametres sont const et
749 // que l'argument DpcImg est const dans la presente method.
750 // (dans dvlist, l'operateur surcharge [] renvoie MuTyV&, ce qui
751 // est non const et provoque un warning sur mac (CodeWarrior)
752 dvl.SetS("NAME",DpcImg.Nom());
753
754 dvl["NBSAT"] = DpcImg.nbSat;
755 dvl["NBNUL"] = DpcImg.nbNul;
756 dvl["MINPIX"] = DpcImg.minPix;
757 dvl["MAXPIX"] = DpcImg.maxPix;
758 dvl["MOYPIX"] = DpcImg.moyPix;
759 dvl["SIGPIX"] = DpcImg.sigPix;
760 dvl["FOND"] = DpcImg.fond;
761 dvl["SIGFON"] = DpcImg.sigmaFond;
762
763 // get the values of the DpcImage
764 if (r_4tab_ != NULL) delete [] r_4tab_;
765 r_4tab_=new r_4[siz_x*siz_y];
766 PBaseDataTypes dataT=DataType((r_4)0);
767 memcpy( r_4tab_, DpcImg.ImagePtr(), siz_x*siz_y*DataSize(dataT));
768 planck_write_img(flnm, naxis, siz_x, siz_y, 0, dvl);
769
770
771 delete [] r_4tab_;
772}
773
774
775void FitsIoServer::save(const ImageI4& DpcImg,char flnm[])
776{
777 long naxis=2;
778 int siz_x = DpcImg.XSize();
779 int siz_y = DpcImg.YSize();
780 int nbpixels = siz_x*siz_y;
781 FITS_tab_typ_ = TINT;
782
783 // pointer to the FITS file, defined in fitsio.h
784 //fitsfile *fptr;
785
786 // initialize status before calling fitsio routines
787 //int status = 0;
788
789 // delete old file if it already exists
790 //remove(flnm);
791
792 // create new FITS file
793 //fits_create_file(&fptr, flnm, &status);
794 //if( status ) printerror( status );
795
796
797 // write FITS image
798 DVList dvl;
799
800 dvl["DATATYPE"] = (int_4)DpcImg.PixelType();
801 dvl["ORG_X"] = DpcImg.XOrg();
802 dvl["ORG_Y"] = DpcImg.YOrg();
803 dvl["PIXSZ_X"] = DpcImg.XPxSize();
804 dvl["PIXSZ_Y"] = DpcImg.YPxSize();
805 dvl["IDENT"] = DpcImg.Ident();
806
807 //dvl["NAME"] = DpcImg.Nom();
808 // j utilise la methode SetS parce que ses parametres sont const et
809 // que l'argument DpcImg est const dans la presente method.
810 // (dans dvlist, l'operateur surcharge [] renvoie MuTyV&, ce qui
811 // est non const et provoque un warning sur mac (CodeWarrior)
812 dvl.SetS("NAME",DpcImg.Nom());
813
814 dvl["NBSAT"] = DpcImg.nbSat;
815 dvl["NBNUL"] = DpcImg.nbNul;
816 dvl["MINPIX"] = DpcImg.minPix;
817 dvl["MAXPIX"] = DpcImg.maxPix;
818 dvl["MOYPIX"] = DpcImg.moyPix;
819 dvl["SIGPIX"] = DpcImg.sigPix;
820 dvl["FOND"] = DpcImg.fond;
821 dvl["SIGFON"] = DpcImg.sigmaFond;
822
823 // get the values of the DpcImage
824 if (i_4tab_ != NULL) delete [] i_4tab_;
825 i_4tab_=new int_4[siz_x*siz_y];
826 PBaseDataTypes dataT=DataType((int_4)0);
827 memcpy( i_4tab_, DpcImg.ImagePtr(), siz_x*siz_y*DataSize(dataT));
828
829
830 planck_write_img(flnm, naxis, siz_x, siz_y, 0, dvl);
831
832 // close the file
833 //fits_close_file(fptr, &status);
834 //if( status ) printerror( status );
835
836 delete [] i_4tab_;
837}
838
839
840
841
842void FitsIoServer::planck_write_img(char flnm[], int naxis,int n1, int n2, int n3, DVList& dvl)
843{
844
845
846 // pointer to the FITS file, defined in fitsio.h
847 fitsfile *fptr;
848
849 // initialize status before calling fitsio routines
850 int status = 0;
851
852 // delete old file if it already exists
853 remove(flnm);
854
855 // create new FITS file
856 fits_create_file(&fptr, flnm, &status);
857 if( status ) printerror( status );
858
859 int bitpix=0;
860 if ( FITS_tab_typ_ == TDOUBLE) bitpix = DOUBLE_IMG;
861 if ( FITS_tab_typ_ == TFLOAT) bitpix = FLOAT_IMG;
862 if ( FITS_tab_typ_ == TINT) bitpix = LONG_IMG;
863 long naxes[3];
864 naxes[0] = n1;
865 naxes[1] = n2;
866 naxes[2] = n3;
867 if (n2 > 0 && naxis < 2) cout << " FitsIoServer:: n2 = " << n2 << " seems to be incompatible with naxis = " << naxis << endl;
868 if (n3 > 0 && naxis < 3) cout << " FitsIoServer:: n3 = " << n3 << " seems to be incompatible with naxis = " << naxis << endl;
869 fits_create_img(fptr, bitpix, naxis, naxes, &status);
870 if( status ) printerror( status );
871
872
873
874 long nelements= naxes[0];
875 if (naxis >=2) nelements*=naxes[1];
876 if (naxis == 3) nelements*=naxes[2];
877 switch (FITS_tab_typ_)
878 {
879 case TDOUBLE :
880 fits_write_img(fptr, TDOUBLE, 1, nelements, r_8tab_, &status);
881 if( status ) printerror( status );
882 break;
883 case TFLOAT :
884 fits_write_img(fptr, TFLOAT, 1, nelements, r_4tab_, &status);
885 if( status ) printerror( status );
886 break;
887 case TINT :
888 fits_write_img(fptr, TINT, 1, nelements, i_4tab_, &status);
889 if( status ) printerror( status );
890 break;
891 default :
892 cout << " FitsIOServer : type de tableau non traite en ecriture " << endl;
893 break;
894 }
895
896 // write the current date
897 fits_write_date(fptr, &status);
898 if( status ) printerror( status );
899 // on convient d ecrire apres la date, les mots cles utilisateur.
900 // la date servira de repere pour relire ces mots cles.
901
902 //dvl.Print();
903 DVList::ValList::const_iterator it;
904 for(it = dvl.Begin(); it != dvl.End(); it++)
905 {
906 int datatype= key_type_PL2FITS( (*it).second.typ);
907 char keyname[]="";
908 int flen_keyword = 9;
909 // FLEN_KEYWORD est la longueur max d'un mot-cle. Il doit y avoir une
910 // erreur dans la librairie fits qui donne FLEN_KEYWORD=72
911 // contrairement a la notice qui donne FLEN_KEYWORD=9 (ch. 5, p.39)
912 strncpy(keyname, (*it).first.substr(0,64).c_str(),flen_keyword-1);
913 char comment[FLEN_COMMENT]="";
914 int ival=0;
915 double dval=0.;
916 char strval[]="";
917 switch (datatype)
918 {
919 case TINT :
920 ival=(*it).second.mtv.iv;
921 strcpy(comment,"I entier");
922 fits_write_key(fptr, datatype, keyname, &ival, comment, &status);
923 break;
924 case TDOUBLE :
925 dval=(*it).second.mtv.dv;
926 strcpy(comment,"D double");
927 fits_write_key(fptr, datatype, keyname, &dval, comment, &status);
928 break;
929 case TSTRING :
930 strcpy(strval, (*it).second.mtv.strv);
931 strcpy(comment,"S character string");
932 fits_write_key(fptr, datatype, keyname, &strval, comment, &status);
933 break;
934 default :
935 cout << " FitsIOServer : probleme dans type mot cle optionnel" << endl;
936 break;
937 }
938 if( status ) printerror( status );
939 }
940 // close the file
941 fits_close_file(fptr, &status);
942 if( status ) printerror( status );
943}
944
945void FitsIoServer::planck_read_img(char flnm[], long& naxis,int& n1, int& n2, int& n3, DVList& dvl)
946{
947 int status=0;
948 long bitpix;
949 long naxes[3]={0,0,0};
950 char* comment=NULL;
951
952 // pointer to the FITS file, defined in fitsio.h
953 fitsfile *fptr;
954 // initialize status before calling fitsio routines
955 fits_open_file(&fptr, flnm, READONLY, &status);
956 if( status ) printerror( status );
957
958
959 fits_read_key_lng(fptr, "BITPIX", &bitpix, comment, &status);
960 if( status ) printerror( status );
961 fits_read_key_lng(fptr, "NAXIS", &naxis, comment, &status);
962 if( status ) printerror( status );
963 int nfound;
964 int nkeys=(int)naxis;
965 fits_read_keys_lng(fptr, "NAXIS", 1, nkeys, naxes, &nfound, &status);
966 if( status ) printerror( status );
967
968 n1 = naxes[0] ;
969 n2 = naxes[1] ;
970 n3 = naxes[2] ;
971
972
973 long nelements= naxes[0];
974 if (naxis >=2) nelements*=naxes[1];
975 if (naxis == 3) nelements*=naxes[2];
976 int anynull;
977 r_8 dnullval=0.;
978 r_4 fnullval=0.;
979 int_4 inullval=0;
980 // on laisse a fits le soin de convertir le type du tableau lu vers
981 // le type suppose par l'utilisateur de fitsioserver
982 //
983 switch ( FITS_tab_typ_)
984 {
985 case TDOUBLE :
986 if (bitpix != DOUBLE_IMG)
987 {
988 cout << " FitsIOServer : the data type on fits file " << flnm << " is not double, "
989 << " conversion to double will be achieved by cfitsio lib " << endl;
990 }
991 if (r_8tab_ != NULL) delete [] r_8tab_;
992 r_8tab_=new r_8[nelements];
993 fits_read_img(fptr, TDOUBLE, 1, nelements, &dnullval, r_8tab_,
994 &anynull, &status);
995 if( status ) printerror( status );
996 break;
997 case TFLOAT :
998 if (bitpix != FLOAT_IMG)
999 {
1000 cout << " FitsIOServer : the data type on fits file " << flnm << " is not float, "
1001 << " conversion to float will be achieved by cfitsio lib " << endl;
1002 }
1003 if (r_4tab_ != NULL) delete [] r_4tab_;
1004 r_4tab_=new r_4[nelements];
1005 fits_read_img(fptr, TFLOAT, 1, nelements, &fnullval, r_4tab_,
1006 &anynull, &status);
1007 if( status ) printerror( status );
1008 break;
1009
1010
1011 case TINT :
1012 if (bitpix != LONG_IMG)
1013 {
1014 cout << " FitsIOServer : the data type on fits file " << flnm << " is not long, "
1015 << " conversion to long will be achieved by cfitsio lib " << endl;
1016 }
1017 if (i_4tab_ != NULL) delete [] i_4tab_;
1018 i_4tab_=new int_4[nelements];
1019 fits_read_img(fptr, TINT, 1, nelements, &inullval, i_4tab_,
1020 &anynull, &status);
1021 if( status ) printerror( status );
1022 break;
1023
1024
1025 default :
1026 cout << " FitsIOServer::read_img : type non traite: " << FITS_tab_typ_ << endl;
1027 break;
1028 }
1029 status = 0;
1030 char card[FLEN_CARD];
1031 int num = 0;
1032 char comment2[FLEN_COMMENT] = "x";
1033 char keyname[]= "";
1034 char datekey[]= "DATE";
1035 char endkey[] = "END";
1036 char typ='x';
1037 int ival;
1038 double dval;
1039 char strval[]="";
1040 // on a convenu que les mots cles utilisateur sont apres le mot cle DATE
1041 // on va jusqu'au mot cle DATE
1042 int flen_keyword = 9;
1043 // FLEN_KEYWORD est la longueur max d'un mot-cle. Il doit y avoir une
1044 // erreur dans la librairie fits qui donne FLEN_KEYWORD=72
1045 // contrairement a la notice qui donne FLEN_KEYWORD=9 (ch. 5, p.39)
1046 while (status == 0 && strncmp(keyname, datekey,4) != 0 )
1047 {
1048 num++;
1049 fits_read_record(fptr, num, card, &status);
1050 strncpy(keyname,card,flen_keyword-1);
1051 }
1052 if (status != 0 )
1053 {
1054 cout << " fitsio::planck_read_img : erreur, mot cle DATE absent " << endl;
1055 }
1056 // on recupere la liste des mots-cles utilisateurs
1057 while (status == 0)
1058 {
1059 num++;
1060 // on lit un record pour recuperer le nom du mot-cle
1061 fits_read_record(fptr, num, card, &status);
1062 strncpy(keyname,card,flen_keyword-1);
1063 char value[FLEN_VALUE];
1064 // on recupere le premier caractere du commentaire, qui contient
1065 // le code du type de la valeur
1066 // (tant que ce n est pas le mot cle END)
1067 fits_read_keyword(fptr, keyname, value, comment2, &status);
1068 if ( strncmp(keyname, endkey,flen_keyword-1) != 0)
1069 {
1070 typ = comment2[0];
1071 // quand le type est connu, on lit la valeur
1072 switch (typ)
1073 {
1074 case 'I' :
1075 fits_read_key(fptr, TINT, keyname, &ival, comment2, &status);
1076 if( status ) printerror( status );
1077 strip (keyname, 'B',' ');
1078 dvl[keyname] = (int_4)ival;
1079 break;
1080 case 'D' :
1081 fits_read_key(fptr, TDOUBLE, keyname, &dval, comment2, &status);
1082 if( status ) printerror( status );
1083 strip (keyname, 'B',' ');
1084 dvl[keyname] = dval;
1085 break;
1086 case 'S' :
1087 fits_read_key(fptr, TSTRING, keyname, strval, comment2, &status);
1088 if( status ) printerror( status );
1089 strip (keyname, 'B',' ');
1090 strip(strval, 'B',' ');
1091 dvl[keyname]=strval;
1092 break;
1093 default :
1094 cout << " FITSIOSERVER::planck_read_img : type de donnee non prevu " << endl;
1095 break;
1096 }
1097 }
1098 }
1099
1100
1101 // close the file
1102 status=0;
1103 fits_close_file(fptr, &status);
1104 if( status ) printerror( status );
1105}
1106
1107
1108
1109// projects a SphericalMap<double>, according sinus-method, and saves onto
1110// a FITS-file
1111void FitsIoServer::sinus_picture_projection(SphericalMap<double>& sph, char filename[])
1112{
1113
1114 long naxes[2]={600, 300};
1115 float* map =new float[ 600*300 ];
1116 int npixels= naxes[0]*naxes[1];
1117
1118 cout << " image FITS en projection SINUS" << endl;
1119 // table will have npixels rows
1120 for(int j=0; j < npixels; j++) map[j]=0.;
1121 for(int j=0; j<naxes[1]; j++)
1122 {
1123 double yd = (j+0.5)/naxes[1]-0.5;
1124 double theta = (0.5-yd)*Pi;
1125 double facteur=1./sin(theta);
1126 for(int i=0; i<naxes[0]; i++)
1127 {
1128 double xa = (i+0.5)/naxes[0]-0.5;
1129 double phi = 2.*Pi*xa*facteur+Pi;
1130 float th=float(theta);
1131 float ff=float(phi);
1132 if (phi<2*Pi && phi>= 0)
1133 {
1134 map[j*naxes[0]+i] = sph.PixValSph(th, ff);
1135 }
1136 }
1137 }
1138
1139 write_picture(naxes, map, filename);
1140 delete [] map;
1141}
1142
1143// projects a SphericalMap<double>, according sinus-method, and saves onto
1144// a FITS-file
1145void FitsIoServer::sinus_picture_projection(SphericalMap<float>& sph, char filename[])
1146{
1147 // le code de cete methode duplique celui de la precedente, la seule
1148 //difference etant le type de sphere en entree. Ces methodes de projection
1149 // sont provisoires, et ne servent que pour les tests. C est pourquoi je
1150 // ne me suis pas casse la tete, pour l instant
1151
1152 long naxes[2]={600, 300};
1153 float* map = new float[ 600*300 ];
1154 int npixels= naxes[0]*naxes[1];
1155
1156 cout << " image FITS en projection SINUS" << endl;
1157 // table will have npixels rows
1158 for(int j=0; j < npixels; j++) map[j]=0.;
1159 for(int j=0; j<naxes[1]; j++)
1160 {
1161 double yd = (j+0.5)/naxes[1]-0.5;
1162 double theta = (0.5-yd)*Pi;
1163 double facteur=1./sin(theta);
1164 for(int i=0; i<naxes[0]; i++)
1165 {
1166 double xa = (i+0.5)/naxes[0]-0.5;
1167 double phi = 2.*Pi*xa*facteur+Pi;
1168 float th=float(theta);
1169 float ff=float(phi);
1170 if (phi<2*Pi && phi>= 0)
1171 {
1172 map[j*naxes[0]+i] = sph.PixValSph(th, ff);
1173 }
1174 }
1175 }
1176
1177 write_picture(naxes, map, filename);
1178 delete [] map;
1179
1180}
1181
1182// saves a (LocalMap<double> on a FITS-file in order to be visualized
1183// (for tests)
1184void FitsIoServer::picture(LocalMap<double>& lcm, char filename[])
1185{
1186
1187 long naxes[2];
1188 naxes[0] = lcm.Size_x();
1189 naxes[1] = lcm.Size_y();
1190 int npixels= naxes[0]*naxes[1];
1191 float* map = new float[npixels];
1192
1193 // table will have npixels rows
1194 for(int j=0; j < npixels; j++) map[j]=0.;
1195 for(int j=0; j<naxes[1]; j++)
1196 {
1197 for(int i=0; i<naxes[0]; i++)
1198 {
1199 map[j*naxes[0]+i] = lcm(i, j);
1200 }
1201 }
1202
1203 write_picture(naxes, map, filename);
1204 delete [] map;
1205}
1206
1207
1208
1209void FitsIoServer::write_picture(long* naxes, float* map, char* filename) const
1210{
1211
1212 int bitpix = FLOAT_IMG;
1213 long naxis = 2;
1214
1215 //pointer to the FITS file, defined in fitsio.h
1216 fitsfile *fptr;
1217 // delete old file if it already exists
1218 remove(filename);
1219 // initialize status before calling fitsio routines
1220 int status = 0;
1221
1222 // create new FITS file
1223 fits_create_file(&fptr, filename, &status);
1224 if( status ) printerror( status );
1225
1226 // write the required header keywords
1227 fits_create_img(fptr, bitpix, naxis, naxes, &status);
1228 if( status ) printerror( status );
1229
1230 // write the current date
1231 fits_write_date(fptr, &status);
1232 if( status ) printerror( status );
1233
1234
1235 // first row in table to write
1236 long firstrow = 1;
1237 // first element in row
1238 long firstelem = 1;
1239 int colnum = 1;
1240 int nelements=naxes[0]*naxes[1];
1241 fits_write_img(fptr, TFLOAT, firstelem, nelements, map, &status);
1242 if( status ) printerror( status );
1243
1244 // close the file
1245 fits_close_file(fptr, &status);
1246 if( status ) printerror( status );
1247}
1248
1249
1250bool FitsIoServer::check_keyword(fitsfile *fptr,int nkeys,char keyword[])
1251
1252 //*****************************************************/
1253 //* check if the specified keyword exits in the CHU */
1254 //*****************************************************/
1255{
1256
1257 bool KEY_EXIST = false;
1258 int status = 0;
1259 char strbide[FLEN_VALUE];
1260 char keybide[LEN_KEYWORD]= "";
1261 for(int jj = 1; jj <= nkeys; jj++)
1262 {
1263 if( fits_read_keyn(fptr,jj,keybide,strbide,NULL,&status) )
1264 printerror( status );
1265 if( !strcmp(keybide,keyword) )
1266 {
1267 KEY_EXIST= true;
1268 break;
1269 }
1270 }
1271 return(KEY_EXIST);
1272}
1273
1274void FitsIoServer::readheader ( char filename[] )
1275
1276 //**********************************************************************/
1277 //* Print out all the header keywords in all extensions of a FITS file */
1278 //**********************************************************************/
1279{
1280
1281 // standard string lengths defined in fitsioc.h
1282 char card[FLEN_CARD];
1283
1284 // pointer to the FITS file, defined in fitsio.h
1285 fitsfile *fptr;
1286
1287 int status = 0;
1288 if ( fits_open_file(&fptr, filename, READONLY, &status) )
1289 printerror( status );
1290
1291 // attempt to move to next HDU, until we get an EOF error
1292 int hdutype;
1293 for (int ii = 1; !(fits_movabs_hdu(fptr,ii,&hdutype,&status));ii++)
1294 {
1295 if (hdutype == ASCII_TBL)
1296 printf("\nReading ASCII table in HDU %d:\n", ii);
1297 else if (hdutype == BINARY_TBL)
1298 printf("\nReading binary table in HDU %d:\n", ii);
1299 else if (hdutype == IMAGE_HDU)
1300 printf("\nReading FITS image in HDU %d:\n", ii);
1301 else
1302 {
1303 printf("Error: unknown type of this HDU \n");
1304 printerror( status );
1305 }
1306
1307 // get the number of keywords
1308 int nkeys, keypos;
1309 if ( fits_get_hdrpos(fptr, &nkeys, &keypos, &status) )
1310 printerror( status );
1311
1312 printf("Header listing for HDU #%d:\n", ii);
1313 for (int jj = 1; jj <= nkeys; jj++)
1314 {
1315 if ( fits_read_record(fptr, jj, card, &status) )
1316 printerror( status );
1317
1318 // print the keyword card
1319 printf("%s\n", card);
1320 }
1321 printf("END\n\n");
1322 }
1323
1324 // got the expected EOF error; reset = 0
1325 if (status == END_OF_FILE)
1326 status = 0;
1327 else
1328 printerror( status );
1329
1330 if ( fits_close_file(fptr, &status) )
1331 printerror( status );
1332
1333 return;
1334}
1335
1336void FitsIoServer::printerror(int status) const
1337
1338 //*****************************************************/
1339 //* Print out cfitsio error messages and exit program */
1340 //*****************************************************/
1341{
1342
1343 // print out cfitsio error messages and exit program
1344 if( status )
1345 {
1346 // print error report
1347 fits_report_error(stderr, status);
1348 // terminate the program, returning error status
1349 exit( status );
1350 }
1351 return;
1352}
1353
1354
1355
Note: See TracBrowser for help on using the repository browser.