source: tags/release-0.9/Cubes/saveImage.cc @ 813

Last change on this file since 813 was 65, checked in by Matthew Whiting, 18 years ago

Changing the references in mainDuchamp.cc, saveImage.cc and getImage.cc (in
the v0.9 release) to fitsio.h so that they compile with the new version of
CFITSIO:
*remove reference to fitsio.h in mainDuchamp.cc
*move fitsio.h reference to after the wcslib include statements in other two
files.
File size: 3.0 KB
Line 
1#include <iostream>
2#include <sstream>
3#include <string>
4#include <wcs.h>
5#include <wcshdr.h>
6#include <fitsio.h>
7#include <Cubes/cubes.hh>
8
9void Cube::saveReconstructedCube()
10{
11 
12  float *resid = new float[this->numPixels];
13  for(int i=0;i<this->numPixels;i++) resid[i] = this->array[i] - this->recon[i];
14  float blankval = this->par.getBlankPixVal();
15
16  short int maxdim=3;
17  long *fpixel = new long[maxdim];
18  for(int i=0;i<maxdim;i++) fpixel[i]=1;
19  long *dimAxes = new long[maxdim];
20  for(int i=0;i<maxdim;i++) dimAxes[i]=1;
21  long nelements;
22  int bitpix,numAxes,anynul;
23  int status = 0,  nkeys, ct=1;  /* MUST initialize status */
24  fitsfile *fptrOld, *fptrNew;         
25  fits_open_file(&fptrOld,this->par.getImageFile().c_str(),READONLY,&status);
26  if (status) fits_report_error(stderr, status);
27 
28  if(this->par.getFlagOutputRecon()){
29    string fileout1 = this->par.getImageFile();
30    fileout1 = fileout1.substr(0,fileout1.size()-5); // remove the ".fits" on the end.
31    std::stringstream ss1;
32    ss1 << ".RECON"<<this->par.getAtrousCut()<<".fits";
33    fileout1 += ss1.str();
34    string syscall = "rm -f " + fileout1;
35    system(syscall.c_str());
36
37    status = 0;
38    fits_create_file(&fptrNew,fileout1.c_str(),&status);
39    if (status) fits_report_error(stderr, status);
40    status = 0;
41    fits_movabs_hdu(fptrOld, 1, NULL, &status);
42    if (status) fits_report_error(stderr, status);
43    status = 0;
44    fits_copy_hdu(fptrOld, fptrNew, 0, &status);
45    if (status) fits_report_error(stderr, status);
46    //////////////////////////////
47    ////// NEED TO ADD NEW HEADER KEYS
48    ////// INDICATING WHAT THE FILE IS...
49    //////////////////////////////
50
51    if(this->par.getFlagBlankPix())
52      fits_write_pixnull(fptrNew, TFLOAT, fpixel, this->numPixels, recon, &blankval, &status);
53    else  fits_write_pix(fptrNew, TFLOAT, fpixel, this->numPixels, recon, &status);
54
55    fits_close_file(fptrNew, &status);
56  }
57
58  if(this->par.getFlagOutputResid()){
59    string fileout2 = this->par.getImageFile();
60    fileout2 = fileout2.substr(0,fileout2.size()-5);
61    std::stringstream ss2;
62    ss2 << ".RESID"<<this->par.getAtrousCut()<<".fits";
63    fileout2 += ss2.str();
64    string syscall = "rm -f " + fileout2;
65    system(syscall.c_str()); 
66
67    status = 0;
68    fits_create_file(&fptrNew,fileout2.c_str(),&status);
69    if (status) fits_report_error(stderr, status);
70    status = 0;
71    fits_movabs_hdu(fptrOld, 1, NULL, &status);
72    if (status) fits_report_error(stderr, status);
73    status = 0;
74    fits_copy_hdu(fptrOld, fptrNew, 0, &status);
75    if (status) fits_report_error(stderr, status);
76    //////////////////////////////
77    ////// NEED TO ADD NEW HEADER KEYS
78    ////// INDICATING WHAT THE FILE IS...
79    //////////////////////////////
80
81    if(this->par.getFlagBlankPix())
82      fits_write_pixnull(fptrNew, TFLOAT, fpixel, this->numPixels, resid, &blankval, &status);
83    else  fits_write_pix(fptrNew, TFLOAT, fpixel, this->numPixels, resid, &status);
84
85    fits_close_file(fptrNew, &status);
86  }
87
88  delete [] recon;
89  delete [] resid;
90  delete [] fpixel;
91  delete [] dimAxes;
92
93}
Note: See TracBrowser for help on using the repository browser.