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

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

load et save d'images 19-OCT-99 GLM
File size: 39.1 KB
Line 
1
2#include <iostream.h>
3#include <list>
4
5#include "fitsioserver.h"
6#include "strutil.h"
7
8void FitsIoServer::load(TMatrix<double>& mat,char flnm[])
9{
10 int nbrows=0;
11 int nbcols=0;
12 FITS_tab_typ_ = TDOUBLE;
13 int status=0;
14 long naxis;
15 int n1, n2, n3;
16 DVList dvl;
17 // pointer to the FITS file, defined in fitsio.h
18 fitsfile *fptr;
19 fits_open_file(&fptr, flnm, READONLY, &status);
20 if( status ) printerror( status );
21 planck_read_img(fptr, naxis, n1, n2, n3, dvl);
22 // close the file
23 fits_close_file(fptr, &status);
24 if( status ) printerror( status );
25
26 nbrows=n1;
27 nbcols=n2;
28 if (naxis == 1) nbcols=1;
29 //dvl.Print();
30 DVList::ValList::const_iterator it;
31 for(it = dvl.Begin(); it != dvl.End(); it++)
32 {
33 char datatype= (*it).second.typ;
34 char keyname[]="";
35 strcpy(keyname, (*it).first.substr(0,64).c_str());
36 int ival=0;
37 double dval=0.;
38 char strval[]="";
39 switch (datatype)
40 {
41 case 'I' :
42 ival=(*it).second.mtv.iv;
43 break;
44 case 'D' :
45 dval=(*it).second.mtv.dv;
46 break;
47 case 'S' :
48 strcpy(strval, (*it).second.mtv.strv);
49 break;
50 default :
51 cout << " FitsIOServer : probleme dans type mot cle optionnel" << endl;
52 break;
53 }
54 }
55
56 if (naxis > 2)
57 {
58 cout<<" FitsIOServer : le fichier fits n'est pas une matrice, naxis= " <<naxis<< endl;
59 }
60
61 // number of components
62 if (mat.NRows() != nbrows || mat.NCols() != nbcols )
63 {
64 cout << " found " << nbrows << " rows ";
65 cout << " expected " << mat.NRows() << endl;
66 cout << " found " << nbcols << " columns " ;
67 cout << " expected " << mat.NCols() << endl;
68 mat.ReSize(nbrows,nbcols);
69 cout << " resize the vector to nbrows= " << nbrows << " nbcols= " << nbcols << endl;
70 }
71 int ij=0;
72 for (int j=0; j< nbcols; j++)
73 for (int i = 0; i < nbrows; i++) mat(i,j) = dtab_[ij++];
74}
75
76void FitsIoServer::load(SphericalMap<double>& sph, char flnm[])
77{
78 int npixels=0;
79 int nside=0;
80 FITS_tab_typ_ = TDOUBLE;
81 read_sphere(flnm, npixels, nside);
82
83 // number of pixels in the sphere
84 if (sph.NbPixels() != npixels)
85 {
86 cout << " found " << npixels << " pixels" << endl;
87 cout << " expected " << sph.NbPixels() << endl;
88 if (nside <= 0 )
89 {
90 cout<<" FITSIOSERVER: no resolution parameter on fits file "<<endl;
91 exit(0);
92 }
93 if (nside != sph.SizeIndex())
94 {
95 sph.Resize(nside);
96 cout << " resizing the sphere to nside= " << nside << endl;
97 }
98 else
99 {
100 cout << " FITSIOSERVER : same resolution, surprising ! " << endl;
101 exit(0);
102 }
103 }
104 for (int j = 0; j < sph.NbPixels(); j++) sph(j)= dtab_[j];
105}
106void FitsIoServer::load(SphericalMap<float>& sph, char flnm[])
107{
108 int npixels=0;
109 int nside=0;
110 FITS_tab_typ_ = TFLOAT;
111 read_sphere(flnm, npixels, nside);
112 // number of pixels in the sphere
113 if (sph.NbPixels() != npixels)
114 {
115 cout << " found " << npixels << " pixels" << endl;
116 cout << " expected " << sph.NbPixels() << endl;
117 if (nside <= 0 )
118 {
119 cout<<" FITSIOSERVER: no resolution parameter on fits file "<<endl;
120 exit(0);
121 }
122 if (nside != sph.SizeIndex())
123 {
124 sph.Resize(nside);
125 cout << " resizing the sphere to nside= " << nside << endl;
126 }
127 else
128 {
129 cout << " FITSIOSERVER : same resolution, surprising ! " << endl;
130 exit(0);
131 }
132 }
133 for (int j = 0; j < sph.NbPixels(); j++) sph(j)= ftab_[j];
134}
135
136void FitsIoServer::load(LocalMap<double>& lcm, char flnm[])
137{
138 int nside;
139 int nbrows=0;
140 int nbcols=0;
141 FITS_tab_typ_ = TDOUBLE;
142 int status=0;
143 long naxis;
144 int n1, n2, n3;
145 DVList dvl;
146 // pointer to the FITS file, defined in fitsio.h
147 fitsfile *fptr;
148 fits_open_file(&fptr, flnm, READONLY, &status);
149 if( status ) printerror( status );
150 planck_read_img(fptr, naxis, n1, n2, n3, dvl);
151 // close the file
152 fits_close_file(fptr, &status);
153 if( status ) printerror( status );
154
155 nbrows=n1;
156 nbcols=n2;
157 if (naxis != 2)
158 {
159 cout<<" FitsIOServer : le fichier fits n'est pas une localmap, naxis= "<<naxis<<endl;
160 }
161 //dvl.Print();
162 DVList::ValList::const_iterator it;
163 for(it = dvl.Begin(); it != dvl.End(); it++)
164 {
165 char datatype= (*it).second.typ;
166 char keyname[]="";
167 strcpy(keyname, (*it).first.substr(0,64).c_str());
168 int ival=0;
169 double dval=0.;
170 char strval[]="";
171 switch (datatype)
172 {
173 case 'I' :
174 ival=(*it).second.mtv.iv;
175 break;
176 case 'D' :
177 dval=(*it).second.mtv.dv;
178 break;
179 case 'S' :
180 strcpy(strval, (*it).second.mtv.strv);
181 break;
182 default :
183 cout << " FitsIOServer : probleme dans type mot cle optionnel" << endl;
184 break;
185 }
186 if (!strcmp(keyname, "NSIDE") )
187 {
188 nside = ival;
189 }
190 //if (!strcmp(keyname, "ORDERING") )
191 // {
192 // cout << " j'ai trouve ORDERING " << endl;
193 // }
194 }
195 float theta0 = dvl.GetD("THETA0");
196 float phi0 = dvl.GetD("PHI0");
197 int x0 = dvl.GetI("X0");
198 int y0 = dvl.GetI("Y0");
199 int xo = dvl.GetI("XO");
200 int yo = dvl.GetI("YO");
201 float anglex=dvl.GetD("ANGLEX");
202 float angley=dvl.GetD("ANGLEY");
203 float angle = dvl.GetD("ANGLE");
204
205 // number of components
206
207 if (lcm.Size_x() != nbrows || lcm.Size_y() != nbcols )
208 {
209 cout << " found " << nbrows << " x-pixels ";
210 cout << " expected " << lcm.Size_x() << endl;
211 cout << " found " << nbcols << " y-pixels " ;
212 cout << " expected " << lcm.Size_y() << endl;
213 lcm.ReSize(nbrows,nbcols);
214 cout << " resizing the map to x-size= " << nbrows << " y-size= " << nbcols << endl;
215 }
216 lcm.SetSize(anglex, angley);
217 lcm.SetOrigin(theta0, phi0, x0, y0, angle);
218 int ij=0;
219 for (int j=0; j< nbcols; j++)
220 for (int i = 0; i < nbrows; i++) lcm(i,j) = dtab_[ij++];
221}
222
223
224void FitsIoServer::save( TMatrix<double>& mat, char filename[])
225{
226 int nbrows = mat.NRows();
227 int nbcols = mat.NCols();
228 long naxis = nbcols > 1 ? 2 : 1;
229 cout << " nombre de lignes : " << nbrows << " colonnes " << nbcols << endl;
230 FITS_tab_typ_ = TDOUBLE;
231 if (dtab_ != NULL ) delete[] dtab_;
232 dtab_=new double[nbrows*nbcols];
233
234
235 //pointer to the FITS file, defined in fitsio.h
236 fitsfile *fptr;
237
238
239 // initialize status before calling fitsio routines
240 int status = 0;
241
242 // delete old file if it already exists
243 remove(filename);
244
245 // create new FITS file
246 fits_create_file(&fptr, filename, &status);
247 if( status ) printerror( status );
248 int ij=0;
249 for (int j=0; j< nbcols; j++)
250 for (int i = 0; i < nbrows; i++) dtab_[ij++]= mat(i,j);
251
252 DVList dvl;
253
254 planck_write_img(fptr, naxis, nbrows, nbcols, 0, dvl);
255
256
257 // close the file
258 fits_close_file(fptr, &status);
259 if( status ) printerror( status );
260
261}
262
263//void FitsIoServer::save(SphericalMap<double>& sph, char filename[])
264void FitsIoServer::save(SphericalMap<double>& sph, char filename[])
265{
266 int npixels = sph.NbPixels();
267 FITS_tab_typ_ = TDOUBLE;
268 if (dtab_ != NULL ) delete[] dtab_;
269 dtab_=new double[npixels];
270
271 //pointer to the FITS file, defined in fitsio.h
272 fitsfile *fptr;
273
274
275 // initialize status before calling fitsio routines
276 int status = 0;
277
278 // delete old file if it already exists
279 remove(filename);
280
281 // create new FITS file
282 fits_create_file(&fptr, filename, &status);
283 if( status ) printerror( status );
284
285 for (int j = 0; j < npixels; j++) dtab_[j]= sph(j);
286 DVList dvl;
287 dvl["NSIDE"] = sph.SizeIndex();
288 dvl["ORDERING"]=sph.TypeOfMap();
289
290 planck_write_img(fptr, 1, npixels, 0, 0, dvl);
291
292 // decider ulterieurement ce qu'on fait de ce qui suit, specifique
293 // pour l'instant, aux spheres gorski. Actuellement, on ne sauve pas
294 // les informations d'analyse harmonique
295 /*
296 // write supplementary keywords
297 fits_write_comment(fptr, " ", &status);
298 if( status ) printerror( status );
299
300
301 strcpy(comment,"HEALPIX Pixelisation");
302 strcpy(svalue, "HEALPIX");
303 fits_write_key(fptr, TSTRING, "PIXTYPE", svalue, comment, &status);
304 if( status ) printerror( status );
305
306
307 strcpy(comment,"pixel ordering scheme, either RING or NESTED");
308 strcpy(svalue, "RING");
309 fits_write_key(fptr, TSTRING, "ORDERING", svalue, comment, &status);
310 if( status ) printerror( status );
311
312
313 strcpy(comment,"Random generator seed");
314 int iseed= sph.iseed();
315 fits_write_key(fptr, TINT, "RANDSEED", &iseed, comment, &status);
316 if( status ) printerror( status );
317
318
319 strcpy(comment,"--------------------------------------------");
320 fits_write_comment(fptr, comment, &status);
321 if( status ) printerror( status );
322
323 strcpy(comment," Above keywords are still likely to change !");
324 fits_write_comment(fptr, comment, &status);
325 if( status ) printerror( status );
326
327 strcpy(comment,"--------------------------------------------");
328 fits_write_comment(fptr, comment, &status);
329 if( status ) printerror( status );
330
331 */
332
333 // close the file
334 fits_close_file(fptr, &status);
335 if( status ) printerror( status );
336
337}
338
339void FitsIoServer::save(SphericalMap<float>& sph, char filename[])
340{
341 int npixels = sph.NbPixels();
342 FITS_tab_typ_ = TFLOAT;
343 if (ftab_ != NULL ) delete[] ftab_;
344 ftab_=new float[npixels];
345
346
347 //pointer to the FITS file, defined in fitsio.h
348 fitsfile *fptr;
349
350
351 // initialize status before calling fitsio routines
352 int status = 0;
353
354 // delete old file if it already exists
355 remove(filename);
356
357 // create new FITS file
358 fits_create_file(&fptr, filename, &status);
359 if( status ) printerror( status );
360 for (int j = 0; j < npixels; j++) ftab_[j]= sph(j);
361 DVList dvl;
362 dvl["NSIDE"] = sph.SizeIndex();
363 dvl["ORDERING"]=sph.TypeOfMap();
364
365 planck_write_img(fptr, 1, npixels, 0, 0, dvl);
366
367
368 // close the file
369 fits_close_file(fptr, &status);
370 if( status ) printerror( status );
371
372}
373
374void FitsIoServer::save(LocalMap<double>& locm, char filename[])
375{
376 int nbrows = locm.Size_x();
377 int nbcols = locm.Size_y();
378 long naxis = 2;
379 cout << " nombre de pts en x : " << nbrows << " en y " << nbcols << endl;
380 FITS_tab_typ_ = TDOUBLE;
381 if (dtab_ != NULL ) delete[] dtab_;
382 dtab_=new double[nbrows*nbcols];
383
384 //pointer to the FITS file, defined in fitsio.h
385 fitsfile *fptr;
386
387
388 // initialize status before calling fitsio routines
389 int status = 0;
390
391 // delete old file if it already exists
392 remove(filename);
393
394 // create new FITS file
395 fits_create_file(&fptr, filename, &status);
396 if( status ) printerror( status );
397 int ij=0;
398 for (int j=0; j< nbcols; j++)
399 for (int i = 0; i < nbrows; i++) dtab_[ij++]= locm(i,j);
400
401 DVList dvl;
402 dvl["NSIDE"] = locm.SizeIndex();
403 dvl["ORDERING"]=locm.TypeOfMap();
404 double theta0;
405 double phi0;
406 int x0;
407 int y0;
408 double angle;
409 locm.Origin(theta0,phi0,x0,y0,angle);
410 double anglex;
411 double angley;
412 locm.Aperture(anglex,angley);
413 dvl["THETA0"] = theta0;
414 dvl["PHI0"] = phi0;
415 dvl["X0"] = x0;
416 dvl["Y0"] = y0;
417 dvl["ANGLE"] = angle;
418 dvl["ANGLEX"] = anglex;
419 dvl["ANGLEY"] = angley;
420 planck_write_img(fptr, naxis, nbrows, nbcols, 0, dvl);
421
422
423 // close the file
424 fits_close_file(fptr, &status);
425 if( status ) printerror( status );
426
427}
428void FitsIoServer::planck_write_img( fitsfile *fptr, int naxis,int n1, int n2, int n3, DVList& dvl)
429{
430
431 int status=0;
432 int bitpix=0;
433 if ( FITS_tab_typ_ == TDOUBLE) bitpix = DOUBLE_IMG;
434 if ( FITS_tab_typ_ == TFLOAT) bitpix = FLOAT_IMG;
435 if ( FITS_tab_typ_ == TINT) bitpix = LONG_IMG;
436 long naxes[3];
437 naxes[0] = n1;
438 naxes[1] = n2;
439 naxes[2] = n3;
440 if (n2 > 0 && naxis < 2) cout << " FitsIoServer:: n2 = " << n2 << " seems to be incompatible with naxis = " << naxis << endl;
441 if (n3 > 0 && naxis < 3) cout << " FitsIoServer:: n3 = " << n3 << " seems to be incompatible with naxis = " << naxis << endl;
442 fits_create_img(fptr, bitpix, naxis, naxes, &status);
443 if( status ) printerror( status );
444
445
446
447 long nelements= naxes[0];
448 if (naxis >=2) nelements*=naxes[1];
449 if (naxis == 3) nelements*=naxes[2];
450 switch (FITS_tab_typ_)
451 {
452 case TDOUBLE :
453 fits_write_img(fptr, TDOUBLE, 1, nelements, dtab_, &status);
454 if( status ) printerror( status );
455 break;
456 case TFLOAT :
457 fits_write_img(fptr, TFLOAT, 1, nelements, ftab_, &status);
458 if( status ) printerror( status );
459 break;
460 case TINT :
461 fits_write_img(fptr, TINT, 1, nelements, i_4tab_, &status);
462 if( status ) printerror( status );
463 break;
464 default :
465 cout << " FitsIOServer : type de tableau non traite en ecriture " << endl;
466 break;
467 }
468
469 // write the current date
470 fits_write_date(fptr, &status);
471 if( status ) printerror( status );
472 // on convient d ecrire apres la date, les mots cles utilisateur.
473 // la date servira de repere pour relire ces mots cles.
474
475 //dvl.Print();
476 DVList::ValList::const_iterator it;
477 for(it = dvl.Begin(); it != dvl.End(); it++)
478 {
479 int datatype= key_type_PL2FITS( (*it).second.typ);
480 char keyname[]="";
481 int flen_keyword = 9;
482 // FLEN_KEYWORD est la longueur max d'un mot-cle. Il doit y avoir une
483 // erreur dans la librairie fits qui donne FLEN_KEYWORD=72
484 // contrairement a la notice qui donne FLEN_KEYWORD=9 (ch. 5, p.39)
485 strncpy(keyname, (*it).first.substr(0,64).c_str(),flen_keyword-1);
486 char comment[FLEN_COMMENT]="";
487 int ival=0;
488 double dval=0.;
489 char strval[]="";
490 switch (datatype)
491 {
492 case TINT :
493 ival=(*it).second.mtv.iv;
494 strcpy(comment,"I entier");
495 fits_write_key(fptr, datatype, keyname, &ival, comment, &status);
496 break;
497 case TDOUBLE :
498 dval=(*it).second.mtv.dv;
499 strcpy(comment,"D double");
500 fits_write_key(fptr, datatype, keyname, &dval, comment, &status);
501 break;
502 case TSTRING :
503 strcpy(strval, (*it).second.mtv.strv);
504 strcpy(comment,"S character string");
505 fits_write_key(fptr, datatype, keyname, &strval, comment, &status);
506 break;
507 default :
508 cout << " FitsIOServer : probleme dans type mot cle optionnel" << endl;
509 break;
510 }
511 if( status ) printerror( status );
512 }
513}
514
515void FitsIoServer::planck_read_img( fitsfile *fptr, long& naxis,int& n1, int& n2, int& n3, DVList& dvl)
516{
517 int status=0;
518 long bitpix;
519 long naxes[3]={0,0,0};
520 char* comment=NULL;
521
522 fits_read_key_lng(fptr, "BITPIX", &bitpix, comment, &status);
523 if( status ) printerror( status );
524 fits_read_key_lng(fptr, "NAXIS", &naxis, comment, &status);
525 if( status ) printerror( status );
526 int nfound;
527 int nkeys=(int)naxis;
528 fits_read_keys_lng(fptr, "NAXIS", 1, nkeys, naxes, &nfound, &status);
529 if( status ) printerror( status );
530
531 n1 = naxes[0] ;
532 n2 = naxes[1] ;
533 n3 = naxes[2] ;
534
535
536 long nelements= naxes[0];
537 if (naxis >=2) nelements*=naxes[1];
538 if (naxis == 3) nelements*=naxes[2];
539 int anynull;
540 double dnullval=0.;
541 float fnullval=0.;
542 // on laisse a fits le soin de convertir le type du tableau lu vers
543 // le type suppose par l'utilisateur de fitsioserver
544 //
545 switch ( FITS_tab_typ_)
546 {
547 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 }
553 if (dtab_ != NULL) delete [] dtab_;
554 dtab_=new double[nelements];
555 fits_read_img(fptr, TDOUBLE, 1, nelements, &dnullval, dtab_,
556 &anynull, &status);
557 if( status ) printerror( status );
558 break;
559 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 }
565 if (ftab_ != NULL) delete [] ftab_;
566 ftab_=new float[nelements];
567 fits_read_img(fptr, TFLOAT, 1, nelements, &fnullval, ftab_,
568 &anynull, &status);
569 if( status ) printerror( status );
570 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
587 default :
588 cout << " FitsIOServer::read_img : type non traite: " << FITS_tab_typ_ << endl;
589 break;
590 }
591 status = 0;
592 char card[FLEN_CARD];
593 int num = 0;
594 char comment2[FLEN_COMMENT] = "x";
595 char keyname[]= "";
596 char datekey[]= "DATE";
597 char endkey[] = "END";
598 char typ='x';
599 int ival;
600 double dval;
601 char strval[]="";
602 // on a convenu que les mots cles utilisateur sont apres le mot cle DATE
603 // on va jusqu'au mot cle DATE
604 int flen_keyword = 9;
605 // FLEN_KEYWORD est la longueur max d'un mot-cle. Il doit y avoir une
606 // erreur dans la librairie fits qui donne FLEN_KEYWORD=72
607 // contrairement a la notice qui donne FLEN_KEYWORD=9 (ch. 5, p.39)
608 while (status == 0 && strncmp(keyname, datekey,4) != 0 )
609 {
610 num++;
611 fits_read_record(fptr, num, card, &status);
612 strncpy(keyname,card,flen_keyword-1);
613 }
614 if (status != 0 )
615 {
616 cout << " fitsio::planck_read_img : erreur, mot cle DATE absent " << endl;
617 }
618 // on recupere la liste des mots-cles utilisateurs
619 while (status == 0)
620 {
621 num++;
622 // on lit un record pour recuperer le nom du mot-cle
623 fits_read_record(fptr, num, card, &status);
624 strncpy(keyname,card,flen_keyword-1);
625 char value[FLEN_VALUE];
626 // on recupere le premier caractere du commentaire, qui contient
627 // le code du type de la valeur
628 // (tant que ce n est pas le mot cle END)
629 fits_read_keyword(fptr, keyname, value, comment2, &status);
630 if ( strncmp(keyname, endkey,flen_keyword-1) != 0)
631 {
632 typ = comment2[0];
633 // quand le type est connu, on lit la valeur
634 switch (typ)
635 {
636 case 'I' :
637 fits_read_key(fptr, TINT, keyname, &ival, comment2, &status);
638 if( status ) printerror( status );
639 strip (keyname, 'B',' ');
640 dvl[keyname] = ival;
641 break;
642 case 'D' :
643 fits_read_key(fptr, TDOUBLE, keyname, &dval, comment2, &status);
644 if( status ) printerror( status );
645 strip (keyname, 'B',' ');
646 dvl[keyname] = dval;
647 break;
648 case 'S' :
649 fits_read_key(fptr, TSTRING, keyname, strval, comment2, &status);
650 if( status ) printerror( status );
651 strip (keyname, 'B',' ');
652 strip(strval, 'B',' ');
653 dvl[keyname]=strval;
654 break;
655 default :
656 cout << " FITSIOSERVER::planck_read_img : type de donnee non prevu " << endl;
657 break;
658 }
659 }
660 }
661
662}
663
664
665void FitsIoServer::sinus_picture_projection(SphericalMap<double>& sph, char filename[])
666{
667
668 long naxes[2]={600, 300};
669 float map [ 600*300 ];
670 int npixels= naxes[0]*naxes[1];
671
672 cout << " image FITS en projection SINUS" << endl;
673 // table will have npixels rows
674 for(int j=0; j < npixels; j++) map[j]=0.;
675 for(int j=0; j<naxes[1]; j++)
676 {
677 double yd = (j+0.5)/naxes[1]-0.5;
678 double theta = (0.5-yd)*Pi;
679 double facteur=1./sin(theta);
680 for(int i=0; i<naxes[0]; i++)
681 {
682 double xa = (i+0.5)/naxes[0]-0.5;
683 double phi = 2.*Pi*xa*facteur+Pi;
684 float th=float(theta);
685 float ff=float(phi);
686 if (phi<2*Pi && phi>= 0)
687 {
688 map[j*naxes[0]+i] = sph.PixValSph(th, ff);
689 }
690 }
691 }
692
693 write_picture(naxes, map, filename);
694
695}
696void FitsIoServer::sinus_picture_projection(SphericalMap<float>& sph, char filename[])
697{
698 // le code de cete methode duplique celui de la precedente, la seule
699 //difference etant le type de sphere en entree. Ces methodes de projection
700 // sont provisoires, et ne servent que pour les tests. C est pourquoi je
701 // ne me suis pas casse la tete, pour l instant
702
703 long naxes[2]={600, 300};
704 float map [ 600*300 ];
705 int npixels= naxes[0]*naxes[1];
706
707 cout << " image FITS en projection SINUS" << endl;
708 // table will have npixels rows
709 for(int j=0; j < npixels; j++) map[j]=0.;
710 for(int j=0; j<naxes[1]; j++)
711 {
712 double yd = (j+0.5)/naxes[1]-0.5;
713 double theta = (0.5-yd)*Pi;
714 double facteur=1./sin(theta);
715 for(int i=0; i<naxes[0]; i++)
716 {
717 double xa = (i+0.5)/naxes[0]-0.5;
718 double phi = 2.*Pi*xa*facteur+Pi;
719 float th=float(theta);
720 float ff=float(phi);
721 if (phi<2*Pi && phi>= 0)
722 {
723 map[j*naxes[0]+i] = sph.PixValSph(th, ff);
724 }
725 }
726 }
727
728 write_picture(naxes, map, filename);
729
730
731}
732
733void FitsIoServer::picture(LocalMap<double>& lcm, char filename[])
734{
735
736 long naxes[2];
737 naxes[0] = lcm.Size_x();
738 naxes[1] = lcm.Size_y();
739 int npixels= naxes[0]*naxes[1];
740 float* map = new float[npixels];
741
742 // table will have npixels rows
743 for(int j=0; j < npixels; j++) map[j]=0.;
744 for(int j=0; j<naxes[1]; j++)
745 {
746 for(int i=0; i<naxes[0]; i++)
747 {
748 map[j*naxes[0]+i] = lcm(i, j);
749 }
750 }
751
752 write_picture(naxes, map, filename);
753 delete [] map;
754}
755
756void FitsIoServer::read_sphere(char flnm[], int& npixels, int& nside)
757{
758 int status=0;
759 long naxis;
760 nside=0;
761 int n1, n2, n3;
762 DVList dvl;
763 // pointer to the FITS file, defined in fitsio.h
764 fitsfile *fptr;
765 fits_open_file(&fptr, flnm, READONLY, &status);
766 if( status ) printerror( status );
767 planck_read_img(fptr, naxis, n1, n2, n3, dvl);
768 // close the file
769 fits_close_file(fptr, &status);
770 if( status ) printerror( status );
771
772 npixels=n1;
773 //dvl.Print();
774 DVList::ValList::const_iterator it;
775 for(it = dvl.Begin(); it != dvl.End(); it++)
776 {
777 char datatype= (*it).second.typ;
778 char keyname[]="";
779 strcpy(keyname, (*it).first.substr(0,64).c_str());
780 int ival=0;
781 double dval=0.;
782 char strval[]="";
783 switch (datatype)
784 {
785 case 'I' :
786 ival=(*it).second.mtv.iv;
787 break;
788 case 'D' :
789 dval=(*it).second.mtv.dv;
790 break;
791 case 'S' :
792 strcpy(strval, (*it).second.mtv.strv);
793 break;
794 default :
795 cout << " FitsIOServer : probleme dans type mot cle optionnel" << endl;
796 break;
797 }
798 if (!strcmp(keyname, "NSIDE") )
799 {
800 nside = ival;
801 }
802 //if (!strcmp(keyname, "ORDERING") )
803 // {
804 // cout << " j'ai trouve ORDERING " << endl;
805 // }
806 }
807
808 if (naxis != 1)
809 {
810 cout << " le fichier fits n'est pas une sphere, naxis= " << naxis << endl;
811 }
812}
813
814
815void FitsIoServer::write_picture(long* naxes, float* map, char* filename) const
816{
817 //pointer to the FITS file, defined in fitsio.h
818 fitsfile *fptr;
819 int bitpix = FLOAT_IMG;
820 long naxis = 2;
821
822 // delete old file if it already exists
823 remove(filename);
824
825 // initialize status before calling fitsio routines
826 int status = 0;
827
828 // create new FITS file
829 fits_create_file(&fptr, filename, &status);
830 if( status ) printerror( status );
831
832 // write the required header keywords
833 fits_create_img(fptr, bitpix, naxis, naxes, &status);
834 if( status ) printerror( status );
835
836 // write the current date
837 fits_write_date(fptr, &status);
838 if( status ) printerror( status );
839
840
841 // first row in table to write
842 long firstrow = 1;
843 // first element in row
844 long firstelem = 1;
845 int colnum = 1;
846 int nelements=naxes[0]*naxes[1];
847 fits_write_img(fptr, TFLOAT, firstelem, nelements, map, &status);
848 if( status ) printerror( status );
849
850 // close the file
851 fits_close_file(fptr, &status);
852 if( status ) printerror( status );
853}
854
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 TracBrowser for help on using the repository browser.