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

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

ajout methode load pour spheregorski
File size: 44.0 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}
293
294
295void FitsIoServer::load(SphereGorski<float>& sph, char flnm[], int hdunum)
296{
297 int npixels=0;
298 int nside=0;
299
300 FITS_tab_typ_ = TFLOAT;
301
302 DVList dvl;
303 planck_read_bntbl(flnm, hdunum, npixels, dvl);
304 //dvl.Print();
305 nside= dvl.GetI("NSIDE");
306 const char* ordering= dvl.GetS("ORDERING").c_str();
307 char* ring = "'RING";
308 if ( strncmp(ordering, ring,3) != 0)
309
310 // if (ordering!="RING ")
311 {
312 cout << " numerotation non RING" << endl;
313 }
314 // number of pixels in the sphere
315 if (sph.NbPixels() != npixels)
316 {
317 cout << " found " << npixels << " pixels" << endl;
318 cout << " expected " << sph.NbPixels() << endl;
319 if (nside <= 0 )
320 {
321 cout<<" FITSIOSERVER: no resolution parameter on fits file "<<endl;
322 exit(0);
323 }
324 if (nside != sph.SizeIndex())
325 {
326 sph.Resize(nside);
327 cout << " resizing the sphere to nside= " << nside << endl;
328 }
329 else
330 {
331 cout << " FITSIOSERVER : same resolution, surprising ! " << endl;
332 exit(0);
333 }
334 }
335 for (int j = 0; j < sph.NbPixels(); j++) sph(j)= (float)r_4tab_[j];
336}
337
338void FitsIoServer::load(SphericalMap<float>& sph, char flnm[])
339{
340 int npixels=0;
341 int nside=0;
342 long naxis;
343 int n1, n2, n3;
344 DVList dvl;
345
346 FITS_tab_typ_ = TFLOAT;
347
348 planck_read_img(flnm, naxis, n1, n2, n3, dvl);
349 if (naxis != 1)
350 {
351 cout << " le fichier fits n'est pas une sphere, naxis= " << naxis << endl;
352 }
353 npixels=n1;
354 nside= dvl.GetI("NSIDE");
355
356 // number of pixels in the sphere
357 if (sph.NbPixels() != npixels)
358 {
359 cout << " found " << npixels << " pixels" << endl;
360 cout << " expected " << sph.NbPixels() << endl;
361 if (nside <= 0 )
362 {
363 cout<<" FITSIOSERVER: no resolution parameter on fits file "<<endl;
364 exit(0);
365 }
366 if (nside != sph.SizeIndex())
367 {
368 sph.Resize(nside);
369 cout << " resizing the sphere to nside= " << nside << endl;
370 }
371 else
372 {
373 cout << " FITSIOSERVER : same resolution, surprising ! " << endl;
374 exit(0);
375 }
376 }
377 for (int j = 0; j < sph.NbPixels(); j++) sph(j)= (float)r_4tab_[j];
378}
379
380void FitsIoServer::load(LocalMap<double>& lcm, char flnm[])
381{
382 int nbrows=0;
383 int nbcols=0;
384 FITS_tab_typ_ = TDOUBLE;
385 long naxis;
386 int n1, n2, n3;
387 DVList dvl;
388 planck_read_img(flnm, naxis, n1, n2, n3, dvl);
389
390 nbrows=n1;
391 nbcols=n2;
392 if (naxis != 2)
393 {
394 cout<<" FitsIOServer : le fichier fits n'est pas une localmap, naxis= "<<naxis<<endl;
395 }
396 float theta0 = dvl.GetD("THETA0");
397 float phi0 = dvl.GetD("PHI0");
398 int x0 = dvl.GetI("X0");
399 int y0 = dvl.GetI("Y0");
400 int xo = dvl.GetI("XO");
401 int yo = dvl.GetI("YO");
402 float anglex=dvl.GetD("ANGLEX");
403 float angley=dvl.GetD("ANGLEY");
404 float angle = dvl.GetD("ANGLE");
405
406 // number of components
407
408 if (lcm.Size_x() != nbrows || lcm.Size_y() != nbcols )
409 {
410 cout << " found " << nbrows << " x-pixels ";
411 cout << " expected " << lcm.Size_x() << endl;
412 cout << " found " << nbcols << " y-pixels " ;
413 cout << " expected " << lcm.Size_y() << endl;
414 lcm.ReSize(nbrows,nbcols);
415 cout << " resizing the map to x-size= " << nbrows << " y-size= " << nbcols << endl;
416 }
417 lcm.SetSize(anglex, angley);
418 lcm.SetOrigin(theta0, phi0, x0, y0, angle);
419 int ij=0;
420 for (int j=0; j< nbcols; j++)
421 for (int i = 0; i < nbrows; i++) lcm(i,j) = (double)r_8tab_[ij++];
422}
423
424void FitsIoServer::load(ImageR4& DpcImg,char flnm[])
425{
426 FITS_tab_typ_ = TFLOAT;
427 long naxis;
428 int siz_x;
429 int siz_y;
430 int n3;
431 DVList dvl;
432
433 planck_read_img(flnm, naxis, siz_x, siz_y, n3, dvl);
434
435
436 if (naxis != 2)
437 {
438 cout << " FitsIOServer : naxis= " << naxis << " not equal to 2 ?" << endl;
439 }
440
441 DpcImg.Allocate(siz_x, siz_y, 0., 0);
442 float* pixelPtr=DpcImg.ImagePtr();
443 // verifications de type
444 PBaseDataTypes dataT=DataType((r_4)0);
445 int TypeDonnees=dvl.GetI("DATATYPE");
446 if (int(dataT) != TypeDonnees)
447 {
448 cout << " FitsIOServer : parameter DATATYPE on file " << flnm << " is not float " << endl;
449 cout << " eventual conversion to float is achieved by cfitsio lib " << endl;
450 }
451
452 memcpy(pixelPtr, r_4tab_, siz_x*siz_y*DataSize(dataT));
453 const char* nom=dvl.GetS("NAME").c_str();
454 DpcImg.SetNameId(dvl.GetI("IDENT"), nom);
455 DpcImg.SetOrg(dvl.GetI("ORG_X"), dvl.GetI("ORG_Y"));
456 DpcImg.SetPxSize(dvl.GetD("PIXSZ_X"), dvl.GetD("PIXSZ_Y"));
457 DpcImg.SetAtt(dvl.GetI("NBNUL"), dvl.GetI("NBSAT"),
458 dvl.GetD("MINPIX"), dvl.GetD("MAXPIX"), dvl.GetD("MOYPIX"),
459 dvl.GetD("SIGPIX"),dvl.GetD("FOND"), dvl.GetD("SIGFON"));
460
461}
462
463
464void FitsIoServer::load(ImageI4& DpcImg,char flnm[])
465{
466 FITS_tab_typ_ = TINT;
467 long naxis;
468 int siz_x;
469 int siz_y;
470 int n3;
471 DVList dvl;
472
473 // pointer to the FITS file, defined in fitsio.h
474 //fitsfile *fptr;
475 // initialize status before calling fitsio routines
476 //int status = 0;
477
478 //fits_open_file(&fptr, flnm, READONLY, &status);
479 //if( status ) printerror( status );
480 planck_read_img(flnm, naxis, siz_x, siz_y, n3, dvl);
481 // close the file
482 //fits_close_file(fptr, &status);
483 //if( status ) printerror( status );
484
485 if (naxis != 2)
486 {
487 cout << " FitsIOServer : naxis= " << naxis << " not equal to 2 ?" << endl;
488 }
489
490 DpcImg.Allocate(siz_x, siz_y, 0., 0);
491 int_4* pixelPtr=DpcImg.ImagePtr();
492
493
494 // verifications de type
495 PBaseDataTypes dataT=DataType((int_4)0);
496 int TypeDonnees=dvl.GetI("DATATYPE");
497 if (int(dataT) != TypeDonnees)
498 {
499 cout << " FitsIOServer : parameter DATATYPE on file " << flnm << " is not int_4 " << endl;
500 cout << " eventual conversion to int_4 is achieved by cfitsio lib " << endl;
501 }
502 memcpy(pixelPtr, i_4tab_, siz_x*siz_y*DataSize(dataT));
503 const char* nom=dvl.GetS("NAME").c_str();
504 DpcImg.SetNameId(dvl.GetI("IDENT"), nom);
505 DpcImg.SetOrg(dvl.GetI("ORG_X"), dvl.GetI("ORG_Y"));
506 DpcImg.SetPxSize(dvl.GetD("PIXSZ_X"), dvl.GetD("PIXSZ_Y"));
507 DpcImg.SetAtt(dvl.GetI("NBNUL"), dvl.GetI("NBSAT"),
508 dvl.GetD("MINPIX"), dvl.GetD("MAXPIX"), dvl.GetD("MOYPIX"),
509 dvl.GetD("SIGPIX"),dvl.GetD("FOND"), dvl.GetD("SIGFON"));
510}
511
512
513void FitsIoServer::save( TMatrix<double>& mat, char filename[])
514{
515 int nbrows = mat.NRows();
516 int nbcols = mat.NCols();
517 long naxis = nbcols > 1 ? 2 : 1;
518 cout << " nombre de lignes : " << nbrows << " colonnes " << nbcols << endl;
519 FITS_tab_typ_ = TDOUBLE;
520 if (r_8tab_ != NULL ) delete[] r_8tab_;
521 r_8tab_=new r_8[nbrows*nbcols];
522
523
524 int ij=0;
525 for (int j=0; j< nbcols; j++)
526 for (int i = 0; i < nbrows; i++) r_8tab_[ij++]= (r_8)mat(i,j);
527
528 DVList dvl;
529
530 planck_write_img(filename, naxis, nbrows, nbcols, 0, dvl);
531 delete[] r_8tab_;
532}
533
534void FitsIoServer::save(NTuple& ntpl,char flnm[])
535
536 //****************************************************/
537 //* read the elements of the NTuple ntpl, and create */
538 //* a FITS file with a binary table extension */
539 //****************************************************/
540{
541
542 // delete old file if it already exists
543 remove(flnm);
544
545 // pointer to the FITS file, defined in fitsio.h
546 fitsfile *fptr;
547 int status = 0;
548 if( fits_create_file(&fptr, flnm, &status) )
549 printerror( status );
550
551 // table will have tfields columns
552 int tfields= ntpl.NVar();
553
554 // table will have nrows rows
555 int nrows= ntpl.NEntry();
556
557 // extension name
558 char extname[] = "NTuple_Binary_tbl";
559
560 // define the name, and the datatype for the tfields columns
561 char **ttype, **tform;
562 ttype= new char*[tfields];
563 tform= new char*[tfields];
564 for(int i = 0; i < tfields; i++)
565 {
566 ttype[i]= new char[FLEN_VALUE];
567 strcpy(ttype[i], ntpl.NomIndex(i));
568 tform[i]= new char[FLEN_VALUE];
569 strcpy(tform[i], "1E");
570 }
571
572 // create a new empty binary table onto the FITS file
573 // physical units if they exist, are defined in the DVList object
574 // so the null pointer is given for the tunit parameters.
575 if ( fits_create_tbl(fptr,BINARY_TBL,nrows,tfields,ttype,tform,
576 NULL,extname,&status) )
577 printerror( status );
578
579 for( int ii = 0; ii < tfields; ii++)
580 {
581 delete [] ttype[ii];
582 delete [] tform[ii];
583 }
584 delete [] ttype;
585 delete [] tform;
586
587
588 // first row in table to write
589 int firstrow = 1;
590
591 // first element in row (ignored in ASCII tables)
592 int firstelem = 1;
593
594 for(int i = 0; i < tfields; i++)
595 {
596 float *dens= new float[nrows];
597 for(int j = 0; j < nrows; j++)
598 {
599 dens[j]= ntpl.GetVal(j,i);
600 }
601 fits_write_col(fptr,TFLOAT,i+1,firstrow,firstelem,nrows,dens,&status);
602 delete []dens;
603 }
604
605 // number of blocks per event
606 int blk= ntpl.BLock();
607 fits_write_key(fptr,TINT,"BLK",&blk,"number of blocks per evt",&status);
608
609 // get names and values from the join DVList object
610 DVList dvl= ntpl.Info();
611 dvl.Print();
612 DVList::ValList::const_iterator it;
613 for(it = dvl.Begin(); it != dvl.End(); it++)
614 {
615 char keytype= (*it).second.typ;
616 char keyname[9]="";
617 strncpy(keyname,(*it).first.substr(0,64).c_str(),8);
618 char comment[FLEN_COMMENT]="";
619
620 switch (keytype)
621 {
622 case 'I' :
623 {
624 int ival=(*it).second.mtv.iv;
625 strcpy(comment,"I entier");
626 fits_write_key(fptr,TINT,keyname,&ival,comment,&status);
627 break;
628 }
629 case 'D' :
630 {
631 double dval=(*it).second.mtv.dv;
632 strcpy(comment,"D double");
633 fits_write_key(fptr,TDOUBLE,keyname,&dval,comment,&status);
634 break;
635 }
636 case 'S' :
637 {
638 char strval[]="";
639 strcpy(strval,(*it).second.mtv.strv);
640 strcpy(comment,"S character string");
641 fits_write_key(fptr,TSTRING,keyname,&strval,comment,&status);
642 break;
643 }
644 }
645 }
646
647 //close the FITS file
648 if ( fits_close_file(fptr, &status) )
649 printerror( status );
650 return;
651}
652
653
654
655
656//void FitsIoServer::save(SphericalMap<double>& sph, char filename[])
657void FitsIoServer::save(SphericalMap<double>& sph, char filename[])
658{
659 int npixels = sph.NbPixels();
660 FITS_tab_typ_ = TDOUBLE;
661 if (r_8tab_ != NULL ) delete[] r_8tab_;
662 r_8tab_=new r_8[npixels];
663
664
665 for (int j = 0; j < npixels; j++) r_8tab_[j]= (r_8)sph(j);
666 DVList dvl;
667 dvl["NSIDE"] = (int_4) sph.SizeIndex();
668 dvl["ORDERING"]=sph.TypeOfMap();
669
670 planck_write_img(filename, 1, npixels, 0, 0, dvl);
671
672 // decider ulterieurement ce qu'on fait de ce qui suit, specifique
673 // pour l'instant, aux spheres gorski.
674
675 /*
676 // write supplementary keywords
677 fits_write_comment(fptr, " ", &status);
678 if( status ) printerror( status );
679
680
681 strcpy(comment,"HEALPIX Pixelisation");
682 strcpy(svalue, "HEALPIX");
683 fits_write_key(fptr, TSTRING, "PIXTYPE", svalue, comment, &status);
684 if( status ) printerror( status );
685
686
687 strcpy(comment,"pixel ordering scheme, either RING or NESTED");
688 strcpy(svalue, "RING");
689 fits_write_key(fptr, TSTRING, "ORDERING", svalue, comment, &status);
690 if( status ) printerror( status );
691
692
693 strcpy(comment,"Random generator seed");
694 int iseed= sph.iseed();
695 fits_write_key(fptr, TINT, "RANDSEED", &iseed, comment, &status);
696 if( status ) printerror( status );
697
698
699 strcpy(comment,"--------------------------------------------");
700 fits_write_comment(fptr, comment, &status);
701 if( status ) printerror( status );
702
703 strcpy(comment," Above keywords are still likely to change !");
704 fits_write_comment(fptr, comment, &status);
705 if( status ) printerror( status );
706
707 strcpy(comment,"--------------------------------------------");
708 fits_write_comment(fptr, comment, &status);
709 if( status ) printerror( status );
710
711 */
712
713 delete[] r_8tab_;
714}
715
716void FitsIoServer::save(SphericalMap<float>& sph, char filename[])
717{
718 int npixels = sph.NbPixels();
719 FITS_tab_typ_ = TFLOAT;
720 if (r_4tab_ != NULL ) delete[] r_4tab_;
721 r_4tab_=new r_4[npixels];
722
723
724 for (int j = 0; j < npixels; j++) r_4tab_[j]= (r_4)sph(j);
725 DVList dvl;
726 dvl["NSIDE"] = (int_4)sph.SizeIndex();
727 dvl["ORDERING"]=sph.TypeOfMap();
728
729 planck_write_img(filename, 1, npixels, 0, 0, dvl);
730 delete[] r_4tab_;
731
732}
733
734void FitsIoServer::save(LocalMap<double>& locm, char filename[])
735{
736 int nbrows = locm.Size_x();
737 int nbcols = locm.Size_y();
738 long naxis = 2;
739 cout << " nombre de pts en x : " << nbrows << " en y " << nbcols << endl;
740 FITS_tab_typ_ = TDOUBLE;
741 if (r_8tab_ != NULL ) delete[] r_8tab_;
742 r_8tab_=new r_8[nbrows*nbcols];
743
744 int ij=0;
745 for (int j=0; j< nbcols; j++)
746 for (int i = 0; i < nbrows; i++) r_8tab_[ij++]= (r_8)locm(i,j);
747
748 DVList dvl;
749 dvl["NSIDE"] = (int_4) locm.SizeIndex();
750 dvl["ORDERING"]=locm.TypeOfMap();
751 double theta0;
752 double phi0;
753 int x0;
754 int y0;
755 double angle;
756 locm.Origin(theta0,phi0,x0,y0,angle);
757 double anglex;
758 double angley;
759 locm.Aperture(anglex,angley);
760 dvl["THETA0"] = theta0;
761 dvl["PHI0"] = phi0;
762 dvl["X0"] = (int_4)x0;
763 dvl["Y0"] = (int_4)y0;
764 dvl["ANGLE"] = angle;
765 dvl["ANGLEX"] = anglex;
766 dvl["ANGLEY"] = angley;
767 planck_write_img(filename, naxis, nbrows, nbcols, 0, dvl);
768 delete[] r_8tab_;}
769
770
771void FitsIoServer::save(const ImageR4& DpcImg,char flnm[])
772{
773 long naxis=2;
774 int siz_x = DpcImg.XSize();
775 int siz_y = DpcImg.YSize();
776 int nbpixels = siz_x*siz_y;
777 FITS_tab_typ_ = TFLOAT;
778
779
780
781 // write FITS image
782 DVList dvl;
783
784 dvl["DATATYPE"] = (int_4)DpcImg.PixelType();
785 dvl["ORG_X"] = DpcImg.XOrg();
786 dvl["ORG_Y"] = DpcImg.YOrg();
787 dvl["PIXSZ_X"] = DpcImg.XPxSize();
788 dvl["PIXSZ_Y"] = DpcImg.YPxSize();
789 dvl["IDENT"] = DpcImg.Ident();
790
791 //dvl["NAME"] = DpcImg.Nom();
792
793 // j utilise la methode SetS parce que ses parametres sont const et
794 // que l'argument DpcImg est const dans la presente method.
795 // (dans dvlist, l'operateur surcharge [] renvoie MuTyV&, ce qui
796 // est non const et provoque un warning sur mac (CodeWarrior)
797 dvl.SetS("NAME",DpcImg.Nom());
798
799 dvl["NBSAT"] = DpcImg.nbSat;
800 dvl["NBNUL"] = DpcImg.nbNul;
801 dvl["MINPIX"] = DpcImg.minPix;
802 dvl["MAXPIX"] = DpcImg.maxPix;
803 dvl["MOYPIX"] = DpcImg.moyPix;
804 dvl["SIGPIX"] = DpcImg.sigPix;
805 dvl["FOND"] = DpcImg.fond;
806 dvl["SIGFON"] = DpcImg.sigmaFond;
807
808 // get the values of the DpcImage
809 if (r_4tab_ != NULL) delete [] r_4tab_;
810 r_4tab_=new r_4[siz_x*siz_y];
811 PBaseDataTypes dataT=DataType((r_4)0);
812 memcpy( r_4tab_, DpcImg.ImagePtr(), siz_x*siz_y*DataSize(dataT));
813 planck_write_img(flnm, naxis, siz_x, siz_y, 0, dvl);
814
815
816 delete [] r_4tab_;
817}
818
819
820void FitsIoServer::save(const ImageI4& DpcImg,char flnm[])
821{
822 long naxis=2;
823 int siz_x = DpcImg.XSize();
824 int siz_y = DpcImg.YSize();
825 int nbpixels = siz_x*siz_y;
826 FITS_tab_typ_ = TINT;
827
828 // pointer to the FITS file, defined in fitsio.h
829 //fitsfile *fptr;
830
831 // initialize status before calling fitsio routines
832 //int status = 0;
833
834 // delete old file if it already exists
835 //remove(flnm);
836
837 // create new FITS file
838 //fits_create_file(&fptr, flnm, &status);
839 //if( status ) printerror( status );
840
841
842 // write FITS image
843 DVList dvl;
844
845 dvl["DATATYPE"] = (int_4)DpcImg.PixelType();
846 dvl["ORG_X"] = DpcImg.XOrg();
847 dvl["ORG_Y"] = DpcImg.YOrg();
848 dvl["PIXSZ_X"] = DpcImg.XPxSize();
849 dvl["PIXSZ_Y"] = DpcImg.YPxSize();
850 dvl["IDENT"] = DpcImg.Ident();
851
852 //dvl["NAME"] = DpcImg.Nom();
853 // j utilise la methode SetS parce que ses parametres sont const et
854 // que l'argument DpcImg est const dans la presente method.
855 // (dans dvlist, l'operateur surcharge [] renvoie MuTyV&, ce qui
856 // est non const et provoque un warning sur mac (CodeWarrior)
857 dvl.SetS("NAME",DpcImg.Nom());
858
859 dvl["NBSAT"] = DpcImg.nbSat;
860 dvl["NBNUL"] = DpcImg.nbNul;
861 dvl["MINPIX"] = DpcImg.minPix;
862 dvl["MAXPIX"] = DpcImg.maxPix;
863 dvl["MOYPIX"] = DpcImg.moyPix;
864 dvl["SIGPIX"] = DpcImg.sigPix;
865 dvl["FOND"] = DpcImg.fond;
866 dvl["SIGFON"] = DpcImg.sigmaFond;
867
868 // get the values of the DpcImage
869 if (i_4tab_ != NULL) delete [] i_4tab_;
870 i_4tab_=new int_4[siz_x*siz_y];
871 PBaseDataTypes dataT=DataType((int_4)0);
872 memcpy( i_4tab_, DpcImg.ImagePtr(), siz_x*siz_y*DataSize(dataT));
873
874
875 planck_write_img(flnm, naxis, siz_x, siz_y, 0, dvl);
876
877 // close the file
878 //fits_close_file(fptr, &status);
879 //if( status ) printerror( status );
880
881 delete [] i_4tab_;
882}
883
884
885
886
887void FitsIoServer::planck_write_img(char flnm[], int naxis,int n1, int n2, int n3, DVList& dvl)
888{
889
890
891 // pointer to the FITS file, defined in fitsio.h
892 fitsfile *fptr;
893
894 // initialize status before calling fitsio routines
895 int status = 0;
896
897 // delete old file if it already exists
898 remove(flnm);
899
900 // create new FITS file
901 fits_create_file(&fptr, flnm, &status);
902 if( status ) printerror( status );
903
904 int bitpix=0;
905 if ( FITS_tab_typ_ == TDOUBLE) bitpix = DOUBLE_IMG;
906 if ( FITS_tab_typ_ == TFLOAT) bitpix = FLOAT_IMG;
907 if ( FITS_tab_typ_ == TINT) bitpix = LONG_IMG;
908 long naxes[3];
909 naxes[0] = n1;
910 naxes[1] = n2;
911 naxes[2] = n3;
912 if (n2 > 0 && naxis < 2) cout << " FitsIoServer:: n2 = " << n2 << " seems to be incompatible with naxis = " << naxis << endl;
913 if (n3 > 0 && naxis < 3) cout << " FitsIoServer:: n3 = " << n3 << " seems to be incompatible with naxis = " << naxis << endl;
914 fits_create_img(fptr, bitpix, naxis, naxes, &status);
915 if( status ) printerror( status );
916
917
918
919 long nelements= naxes[0];
920 if (naxis >=2) nelements*=naxes[1];
921 if (naxis == 3) nelements*=naxes[2];
922 switch (FITS_tab_typ_)
923 {
924 case TDOUBLE :
925 fits_write_img(fptr, TDOUBLE, 1, nelements, r_8tab_, &status);
926 if( status ) printerror( status );
927 break;
928 case TFLOAT :
929 fits_write_img(fptr, TFLOAT, 1, nelements, r_4tab_, &status);
930 if( status ) printerror( status );
931 break;
932 case TINT :
933 fits_write_img(fptr, TINT, 1, nelements, i_4tab_, &status);
934 if( status ) printerror( status );
935 break;
936 default :
937 cout << " FitsIOServer : type de tableau non traite en ecriture " << endl;
938 break;
939 }
940
941 // write the current date
942 fits_write_date(fptr, &status);
943 if( status ) printerror( status );
944 // on convient d ecrire apres la date, les mots cles utilisateur.
945 // la date servira de repere pour relire ces mots cles.
946
947 //dvl.Print();
948 DVList::ValList::const_iterator it;
949 for(it = dvl.Begin(); it != dvl.End(); it++)
950 {
951 int datatype= key_type_PL2FITS( (*it).second.typ);
952 char keyname[]="";
953 int flen_keyword = 9;
954 // FLEN_KEYWORD est la longueur max d'un mot-cle. Il doit y avoir une
955 // erreur dans la librairie fits qui donne FLEN_KEYWORD=72
956 // contrairement a la notice qui donne FLEN_KEYWORD=9 (ch. 5, p.39)
957 strncpy(keyname, (*it).first.substr(0,64).c_str(),flen_keyword-1);
958 char comment[FLEN_COMMENT]="";
959 int ival=0;
960 double dval=0.;
961 char strval[]="";
962 switch (datatype)
963 {
964 case TINT :
965 ival=(*it).second.mtv.iv;
966 strcpy(comment,"I entier");
967 fits_write_key(fptr, datatype, keyname, &ival, comment, &status);
968 break;
969 case TDOUBLE :
970 dval=(*it).second.mtv.dv;
971 strcpy(comment,"D double");
972 fits_write_key(fptr, datatype, keyname, &dval, comment, &status);
973 break;
974 case TSTRING :
975 strcpy(strval, (*it).second.mtv.strv);
976 strcpy(comment,"S character string");
977 fits_write_key(fptr, datatype, keyname, &strval, comment, &status);
978 break;
979 default :
980 cout << " FitsIOServer : probleme dans type mot cle optionnel" << endl;
981 break;
982 }
983 if( status ) printerror( status );
984 }
985 // close the file
986 fits_close_file(fptr, &status);
987 if( status ) printerror( status );
988}
989
990void FitsIoServer::planck_read_img(char flnm[], long& naxis,int& n1, int& n2, int& n3, DVList& dvl)
991{
992 int status=0;
993 long bitpix;
994 long naxes[3]={0,0,0};
995 char* comment=NULL;
996
997 // pointer to the FITS file, defined in fitsio.h
998 fitsfile *fptr;
999 // initialize status before calling fitsio routines
1000 fits_open_file(&fptr, flnm, READONLY, &status);
1001 if( status ) printerror( status );
1002
1003
1004 fits_read_key_lng(fptr, "BITPIX", &bitpix, comment, &status);
1005 if( status ) printerror( status );
1006 fits_read_key_lng(fptr, "NAXIS", &naxis, comment, &status);
1007 if( status ) printerror( status );
1008 int nfound;
1009 int nkeys=(int)naxis;
1010 fits_read_keys_lng(fptr, "NAXIS", 1, nkeys, naxes, &nfound, &status);
1011 if( status ) printerror( status );
1012
1013 n1 = naxes[0] ;
1014 n2 = naxes[1] ;
1015 n3 = naxes[2] ;
1016
1017
1018 long nelements= naxes[0];
1019 if (naxis >=2) nelements*=naxes[1];
1020 if (naxis == 3) nelements*=naxes[2];
1021 int anynull;
1022 r_8 dnullval=0.;
1023 r_4 fnullval=0.;
1024 int_4 inullval=0;
1025 // on laisse a fits le soin de convertir le type du tableau lu vers
1026 // le type suppose par l'utilisateur de fitsioserver
1027 //
1028 switch ( FITS_tab_typ_)
1029 {
1030 case TDOUBLE :
1031 if (bitpix != DOUBLE_IMG)
1032 {
1033 cout << " FitsIOServer : the data type on fits file " << flnm << " is not double, "
1034 << " conversion to double will be achieved by cfitsio lib " << endl;
1035 }
1036 if (r_8tab_ != NULL) delete [] r_8tab_;
1037 r_8tab_=new r_8[nelements];
1038 fits_read_img(fptr, TDOUBLE, 1, nelements, &dnullval, r_8tab_,
1039 &anynull, &status);
1040 if( status ) printerror( status );
1041 break;
1042 case TFLOAT :
1043 if (bitpix != FLOAT_IMG)
1044 {
1045 cout << " FitsIOServer : the data type on fits file " << flnm << " is not float, "
1046 << " conversion to float will be achieved by cfitsio lib " << endl;
1047 }
1048 if (r_4tab_ != NULL) delete [] r_4tab_;
1049 r_4tab_=new r_4[nelements];
1050 fits_read_img(fptr, TFLOAT, 1, nelements, &fnullval, r_4tab_,
1051 &anynull, &status);
1052 if( status ) printerror( status );
1053 break;
1054
1055
1056 case TINT :
1057 if (bitpix != LONG_IMG)
1058 {
1059 cout << " FitsIOServer : the data type on fits file " << flnm << " is not long, "
1060 << " conversion to long will be achieved by cfitsio lib " << endl;
1061 }
1062 if (i_4tab_ != NULL) delete [] i_4tab_;
1063 i_4tab_=new int_4[nelements];
1064 fits_read_img(fptr, TINT, 1, nelements, &inullval, i_4tab_,
1065 &anynull, &status);
1066 if( status ) printerror( status );
1067 break;
1068
1069
1070 default :
1071 cout << " FitsIOServer::read_img : type non traite: " << FITS_tab_typ_ << endl;
1072 break;
1073 }
1074 status = 0;
1075 char card[FLEN_CARD];
1076 int num = 0;
1077 char comment2[FLEN_COMMENT] = "x";
1078 char keyname[]= "";
1079 char datekey[]= "DATE";
1080 char endkey[] = "END";
1081 char typ='x';
1082 int ival;
1083 double dval;
1084 char strval[]="";
1085 // on a convenu que les mots cles utilisateur sont apres le mot cle DATE
1086 // on va jusqu'au mot cle DATE
1087 int flen_keyword = 9;
1088 // FLEN_KEYWORD est la longueur max d'un mot-cle. Il doit y avoir une
1089 // erreur dans la librairie fits qui donne FLEN_KEYWORD=72
1090 // contrairement a la notice qui donne FLEN_KEYWORD=9 (ch. 5, p.39)
1091 while (status == 0 && strncmp(keyname, datekey,4) != 0 )
1092 {
1093 num++;
1094 fits_read_record(fptr, num, card, &status);
1095 strncpy(keyname,card,flen_keyword-1);
1096 }
1097 if (status != 0 )
1098 {
1099 cout << " fitsio::planck_read_img : erreur, mot cle DATE absent " << endl;
1100 }
1101 // on recupere la liste des mots-cles utilisateurs
1102 while (status == 0)
1103 {
1104 num++;
1105 // on lit un record pour recuperer le nom du mot-cle
1106 fits_read_record(fptr, num, card, &status);
1107 strncpy(keyname,card,flen_keyword-1);
1108 char value[FLEN_VALUE];
1109 // on recupere le premier caractere du commentaire, qui contient
1110 // le code du type de la valeur
1111 // (tant que ce n est pas le mot cle END)
1112 fits_read_keyword(fptr, keyname, value, comment2, &status);
1113 if ( strncmp(keyname, endkey,flen_keyword-1) != 0)
1114 {
1115 typ = comment2[0];
1116 // quand le type est connu, on lit la valeur
1117 switch (typ)
1118 {
1119 case 'I' :
1120 fits_read_key(fptr, TINT, keyname, &ival, comment2, &status);
1121 if( status ) printerror( status );
1122 strip (keyname, 'B',' ');
1123 dvl[keyname] = (int_4)ival;
1124 break;
1125 case 'D' :
1126 fits_read_key(fptr, TDOUBLE, keyname, &dval, comment2, &status);
1127 if( status ) printerror( status );
1128 strip (keyname, 'B',' ');
1129 dvl[keyname] = dval;
1130 break;
1131 case 'S' :
1132 fits_read_key(fptr, TSTRING, keyname, strval, comment2, &status);
1133 if( status ) printerror( status );
1134 strip (keyname, 'B',' ');
1135 strip(strval, 'B',' ');
1136 dvl[keyname]=strval;
1137 break;
1138 default :
1139 cout << " FITSIOSERVER::planck_read_img : type de donnee non prevu " << endl;
1140 break;
1141 }
1142 }
1143 }
1144
1145
1146 // close the file
1147 status=0;
1148 fits_close_file(fptr, &status);
1149 if( status ) printerror( status );
1150}
1151
1152
1153void FitsIoServer::planck_read_bntbl(char flnm[], int hdunum, int& npixels, DVList& dvl)
1154{
1155 int status=0;
1156 int nkeys,keypos;
1157 int hdutype;
1158 int tfields;
1159 int datype;
1160 long lastpix;
1161 long repeat, width;
1162 long nrows;
1163 long extend;
1164 char* comment=NULL;
1165
1166 // pointer to the FITS file, defined in fitsio.h
1167 fitsfile *fptr;
1168 // initialize status before calling fitsio routines
1169 fits_open_file(&fptr, flnm, READONLY, &status);
1170 if( status ) printerror( status );
1171 fits_read_key_lng(fptr, "EXTEND", &extend, comment, &status);
1172 if( status ) printerror( status );
1173 if (extend!=1)
1174 {
1175 cout << "FitsIoServer:: le fichier fits ne contient pas d'extension binary table" << endl;
1176 return;
1177 }
1178 fits_movabs_hdu(fptr, hdunum,&hdutype,&status);
1179 if( status ) printerror( status );
1180 if (hdutype!=BINARY_TBL)
1181 {
1182 cout << "FitsIoServer:: this HDU is not a binary table " << endl;
1183 exit(status);
1184 }
1185 char xtension[FLEN_VALUE];
1186 fits_read_key_str(fptr,"XTENSION",xtension,NULL,&status);
1187 if( status ) printerror( status );
1188 cout << "xtension= " << xtension << endl;
1189
1190 char binta[] = "BINTABLE";
1191 if ( strncmp(xtension, binta,8) != 0)
1192 // if (xtension !="BINTABLE")
1193 {
1194 cout << "FitsIoServer:: not a binary table " << endl;
1195 exit(status);
1196 }
1197 fits_get_hdrpos(fptr,&nkeys,&keypos,&status);
1198 if( status ) printerror( status );
1199 //cout << " nombre de mots-cles : " << nkeys << endl;
1200 fits_get_num_cols(fptr, &tfields, &status);
1201 cout << "tfields= " << tfields << endl;
1202 if (tfields != 1)
1203 {
1204 cout << "FitsIoServer:: il y a plus d'une colonne" << endl;
1205 return;
1206 }
1207 fits_get_num_rows(fptr, &nrows, &status);
1208 //cout << "nblignes= " << nrows << endl;
1209 fits_get_coltype(fptr, 1, &datype, &repeat, &width, &status);
1210 if( status ) printerror( status );
1211 //cout << " type de donnees : " << datype << endl;
1212 //cout << " repeat : " << repeat << endl;
1213 //cout << " width : " << width << endl;
1214 fits_read_key_lng(fptr, "LASTPIX", &lastpix, comment, &status);
1215 if( status ) printerror( status );
1216
1217 long nelements= nrows*repeat;
1218 if (nelements!=lastpix+1)
1219 {
1220 cout << " erreur sur longueur du vecteur " << endl;
1221 cout << " nelements= " << nelements << " lastpix+1=" << lastpix+1 << endl;
1222 }
1223 npixels=nelements;
1224 int anynull;
1225 r_8 dnullval=0.;
1226 r_4 fnullval=0.;
1227 int_4 inullval=0;
1228 // on laisse a fits le soin de convertir le type du tableau lu vers
1229 // le type suppose par l'utilisateur de fitsioserver
1230 //
1231 switch ( FITS_tab_typ_)
1232 {
1233 case TDOUBLE :
1234 if (datype != TDOUBLE)
1235 {
1236 cout << " FitsIOServer : the data type on fits file " << flnm << " is not double, "
1237 << " conversion to double will be achieved by cfitsio lib " << endl;
1238 }
1239 if (r_8tab_ != NULL) delete [] r_8tab_;
1240 r_8tab_=new r_8[ npixels];
1241 fits_read_col(fptr, TDOUBLE, 1, 1, 1, nelements, &dnullval,
1242 r_8tab_,
1243 &anynull, &status);
1244 if( status ) printerror( status );
1245 break;
1246 case TFLOAT :
1247 if (datype != TFLOAT)
1248 {
1249
1250 cout << " FitsIOServer : the data type on fits file " << flnm << " is not float, "
1251 << " conversion to float will be achieved by cfitsio lib " << endl;
1252 }
1253 if (r_4tab_ != NULL) delete [] r_4tab_;
1254 r_4tab_=new r_4[nelements];
1255 fits_read_col(fptr, TFLOAT, 1, 1, 1, nelements, &fnullval,
1256 r_4tab_, &anynull, &status);
1257 if( status ) printerror( status );
1258 break;
1259
1260
1261 case TINT :
1262 if (datype != TLONG)
1263 {
1264 cout << " FitsIOServer : the data type on fits file " << flnm << " is not long, "
1265 << " conversion to long will be achieved by cfitsio lib " << endl;
1266 }
1267 if (i_4tab_ != NULL) delete [] i_4tab_;
1268 i_4tab_=new int_4[nelements];
1269 fits_read_col(fptr, TLONG, 1, 1, 1, nelements, &inullval,
1270 i_4tab_, &anynull, &status);
1271 //fits_read_img(fptr, TINT, 1, nelements, &inullval, i_4tab_,
1272 // &anynull, &status);
1273 //if( status ) printerror( status );
1274 break;
1275
1276
1277 default :
1278 cout << " FitsIOServer::read_img : type non traite: " << FITS_tab_typ_ << endl;
1279 break;
1280 }
1281 char card[FLEN_CARD];
1282 char keyname[LEN_KEYWORD]= "";
1283 char strval[FLEN_VALUE];
1284 char comment1[FLEN_COMMENT];
1285 char dtype;
1286 //char bidon[LEN_KEYWORD];
1287 char comkey[] = "COMMENT";
1288
1289 for(int j = 1; j <= nkeys; j++)
1290 {
1291 // fits_read_record(fptr, j, card, &status);
1292 // strncpy(keyname,card,LEN_KEYWORD-1);
1293 // cout << " bidon= " << keyname << endl;
1294 // if ( strncmp(keyname,comkey ,LEN_KEYWORD-1) != 0)
1295 fits_read_keyn(fptr,j,card,strval,comment1,&status);
1296 strncpy(keyname,card,LEN_KEYWORD-1);
1297
1298 if ( strncmp(keyname,comkey ,LEN_KEYWORD-1) != 0)
1299 {
1300 fits_get_keytype(strval,&dtype,&status);
1301 strip (keyname, 'B',' ');
1302 strip(strval, 'B',' ');
1303 switch( dtype )
1304 {
1305 case 'C':
1306 dvl[keyname]= strval;
1307 break;
1308 case 'I':
1309 int ival;
1310 ctoi(strval,&ival);
1311 dvl[keyname]= (int_4)ival;
1312 break;
1313 case 'F':
1314 double dval;
1315 ctof(strval,&dval);
1316 dvl[keyname]= dval;
1317 break;
1318 default :
1319 cout << " FitsIOServer : mot-cle bizerre " << endl;
1320 break;
1321 }
1322 }
1323 }
1324
1325 // close the file
1326 status=0;
1327 fits_close_file(fptr, &status);
1328 if( status ) printerror( status );
1329}
1330
1331
1332// projects a SphericalMap<double>, according sinus-method, and saves onto
1333// a FITS-file
1334void FitsIoServer::sinus_picture_projection(SphericalMap<double>& sph, char filename[])
1335{
1336
1337 long naxes[2]={600, 300};
1338 float* map =new float[ 600*300 ];
1339 int npixels= naxes[0]*naxes[1];
1340
1341 cout << " image FITS en projection SINUS" << endl;
1342 // table will have npixels rows
1343 for(int j=0; j < npixels; j++) map[j]=0.;
1344 for(int j=0; j<naxes[1]; j++)
1345 {
1346 double yd = (j+0.5)/naxes[1]-0.5;
1347 double theta = (0.5-yd)*Pi;
1348 double facteur=1./sin(theta);
1349 for(int i=0; i<naxes[0]; i++)
1350 {
1351 double xa = (i+0.5)/naxes[0]-0.5;
1352 double phi = 2.*Pi*xa*facteur+Pi;
1353 float th=float(theta);
1354 float ff=float(phi);
1355 if (phi<2*Pi && phi>= 0)
1356 {
1357 map[j*naxes[0]+i] = sph.PixValSph(th, ff);
1358 }
1359 }
1360 }
1361
1362 write_picture(naxes, map, filename);
1363 delete [] map;
1364}
1365
1366// projects a SphericalMap<double>, according sinus-method, and saves onto
1367// a FITS-file
1368void FitsIoServer::sinus_picture_projection(SphericalMap<float>& sph, char filename[])
1369{
1370 // le code de cete methode duplique celui de la precedente, la seule
1371 //difference etant le type de sphere en entree. Ces methodes de projection
1372 // sont provisoires, et ne servent que pour les tests. C est pourquoi je
1373 // ne me suis pas casse la tete, pour l instant
1374
1375 long naxes[2]={600, 300};
1376 float* map = new float[ 600*300 ];
1377 int npixels= naxes[0]*naxes[1];
1378
1379 cout << " image FITS en projection SINUS" << endl;
1380 // table will have npixels rows
1381 for(int j=0; j < npixels; j++) map[j]=0.;
1382 for(int j=0; j<naxes[1]; j++)
1383 {
1384 double yd = (j+0.5)/naxes[1]-0.5;
1385 double theta = (0.5-yd)*Pi;
1386 double facteur=1./sin(theta);
1387 for(int i=0; i<naxes[0]; i++)
1388 {
1389 double xa = (i+0.5)/naxes[0]-0.5;
1390 double phi = 2.*Pi*xa*facteur+Pi;
1391 float th=float(theta);
1392 float ff=float(phi);
1393 if (phi<2*Pi && phi>= 0)
1394 {
1395 map[j*naxes[0]+i] = sph.PixValSph(th, ff);
1396 }
1397 }
1398 }
1399
1400 write_picture(naxes, map, filename);
1401 delete [] map;
1402
1403}
1404
1405// projects a SphericalMap<double>, according Mollweide-method, and saves onto
1406// a FITS-file
1407void FitsIoServer::Mollweide_picture_projection(SphericalMap<float>& sph, char filename[])
1408{
1409 // le code de cete methode duplique celui de la precedente, la seule
1410 //difference etant le type de sphere en entree. Ces methodes de projection
1411 // sont provisoires, et ne servent que pour les tests. C est pourquoi je
1412 // ne me suis pas casse la tete, pour l instant
1413
1414 long naxes[2]={600, 300};
1415 float* map = new float[ 600*300 ];
1416 int npixels= naxes[0]*naxes[1];
1417
1418 cout << " image FITS en projection MOLLWEIDE" << endl;
1419 // table will have npixels rows
1420 for(int j=0; j < npixels; j++) map[j]=0.;
1421 for(int j=0; j<naxes[1]; j++)
1422 {
1423 double yd = (j+0.5)/naxes[1]-0.5;
1424 double facteur=2.*Pi/sin(acos(yd*2));
1425 double theta = (0.5-yd)*Pi;
1426 for(int i=0; i<naxes[0]; i++)
1427 {
1428 double xa = (i+0.5)/naxes[0]-0.5;
1429 double phi = xa*facteur+Pi;
1430 float th=float(theta);
1431 float ff=float(phi);
1432 if (phi<2*Pi && phi>= 0)
1433 {
1434 map[j*naxes[0]+i] = sph.PixValSph(th, ff);
1435 }
1436 }
1437 }
1438
1439 write_picture(naxes, map, filename);
1440 delete [] map;
1441
1442}
1443
1444
1445
1446// saves a (LocalMap<double> on a FITS-file in order to be visualized
1447// (for tests)
1448void FitsIoServer::picture(LocalMap<double>& lcm, char filename[])
1449{
1450
1451 long naxes[2];
1452 naxes[0] = lcm.Size_x();
1453 naxes[1] = lcm.Size_y();
1454 int npixels= naxes[0]*naxes[1];
1455 float* map = new float[npixels];
1456
1457 // table will have npixels rows
1458 for(int j=0; j < npixels; j++) map[j]=0.;
1459 for(int j=0; j<naxes[1]; j++)
1460 {
1461 for(int i=0; i<naxes[0]; i++)
1462 {
1463 map[j*naxes[0]+i] = lcm(i, j);
1464 }
1465 }
1466
1467 write_picture(naxes, map, filename);
1468 delete [] map;
1469}
1470
1471
1472
1473void FitsIoServer::write_picture(long* naxes, float* map, char* filename) const
1474{
1475
1476 int bitpix = FLOAT_IMG;
1477 long naxis = 2;
1478
1479 //pointer to the FITS file, defined in fitsio.h
1480 fitsfile *fptr;
1481 // delete old file if it already exists
1482 remove(filename);
1483 // initialize status before calling fitsio routines
1484 int status = 0;
1485
1486 // create new FITS file
1487 fits_create_file(&fptr, filename, &status);
1488 if( status ) printerror( status );
1489
1490 // write the required header keywords
1491 fits_create_img(fptr, bitpix, naxis, naxes, &status);
1492 if( status ) printerror( status );
1493
1494 // write the current date
1495 fits_write_date(fptr, &status);
1496 if( status ) printerror( status );
1497
1498
1499 // first row in table to write
1500 long firstrow = 1;
1501 // first element in row
1502 long firstelem = 1;
1503 int colnum = 1;
1504 int nelements=naxes[0]*naxes[1];
1505 fits_write_img(fptr, TFLOAT, firstelem, nelements, map, &status);
1506 if( status ) printerror( status );
1507
1508 // close the file
1509 fits_close_file(fptr, &status);
1510 if( status ) printerror( status );
1511}
1512
1513
1514bool FitsIoServer::check_keyword(fitsfile *fptr,int nkeys,char keyword[])
1515
1516 //*****************************************************/
1517 //* check if the specified keyword exits in the CHU */
1518 //*****************************************************/
1519{
1520
1521 bool KEY_EXIST = false;
1522 int status = 0;
1523 char strbide[FLEN_VALUE];
1524 char keybide[LEN_KEYWORD]= "";
1525 for(int jj = 1; jj <= nkeys; jj++)
1526 {
1527 if( fits_read_keyn(fptr,jj,keybide,strbide,NULL,&status) )
1528 printerror( status );
1529 if( !strcmp(keybide,keyword) )
1530 {
1531 KEY_EXIST= true;
1532 break;
1533 }
1534 }
1535 return(KEY_EXIST);
1536}
1537
1538void FitsIoServer::readheader ( char filename[] )
1539
1540 //**********************************************************************/
1541 //* Print out all the header keywords in all extensions of a FITS file */
1542 //**********************************************************************/
1543{
1544
1545 // standard string lengths defined in fitsioc.h
1546 char card[FLEN_CARD];
1547
1548 // pointer to the FITS file, defined in fitsio.h
1549 fitsfile *fptr;
1550
1551 int status = 0;
1552 if ( fits_open_file(&fptr, filename, READONLY, &status) )
1553 printerror( status );
1554
1555 // attempt to move to next HDU, until we get an EOF error
1556 int hdutype;
1557 for (int ii = 1; !(fits_movabs_hdu(fptr,ii,&hdutype,&status));ii++)
1558 {
1559 if (hdutype == ASCII_TBL)
1560 printf("\nReading ASCII table in HDU %d:\n", ii);
1561 else if (hdutype == BINARY_TBL)
1562 printf("\nReading binary table in HDU %d:\n", ii);
1563 else if (hdutype == IMAGE_HDU)
1564 printf("\nReading FITS image in HDU %d:\n", ii);
1565 else
1566 {
1567 printf("Error: unknown type of this HDU \n");
1568 printerror( status );
1569 }
1570
1571 // get the number of keywords
1572 int nkeys, keypos;
1573 if ( fits_get_hdrpos(fptr, &nkeys, &keypos, &status) )
1574 printerror( status );
1575
1576 printf("Header listing for HDU #%d:\n", ii);
1577 for (int jj = 1; jj <= nkeys; jj++)
1578 {
1579 if ( fits_read_record(fptr, jj, card, &status) )
1580 printerror( status );
1581
1582 // print the keyword card
1583 printf("%s\n", card);
1584 }
1585 printf("END\n\n");
1586 }
1587
1588 // got the expected EOF error; reset = 0
1589 if (status == END_OF_FILE)
1590 status = 0;
1591 else
1592 printerror( status );
1593
1594 if ( fits_close_file(fptr, &status) )
1595 printerror( status );
1596
1597 return;
1598}
1599
1600void FitsIoServer::printerror(int status) const
1601
1602 //*****************************************************/
1603 //* Print out cfitsio error messages and exit program */
1604 //*****************************************************/
1605{
1606
1607 // print out cfitsio error messages and exit program
1608 if( status )
1609 {
1610 // print error report
1611 fits_report_error(stderr, status);
1612 // terminate the program, returning error status
1613 exit( status );
1614 }
1615 return;
1616}
1617
1618
1619
Note: See TracBrowser for help on using the repository browser.