source: trunk/Cubes/saveImage.cc @ 84

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

Added ability to save and read back in the MW parameters. If a reconstructed
cube has different MW range to that being considered, it would cause problems.
ReadRecon? now just returns a FAILURE and will reconstruct the cube anew.
File size: 6.7 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 <duchamp.hh>
8#include <Cubes/cubes.hh>
9
10void write_header_info(fitsfile *fptr, Param &par, string nature);
11void Cube::saveReconstructedCube()
12{
13 
14  int newbitpix = FLOAT_IMG;
15
16  float *resid = new float[this->numPixels];
17  for(int i=0;i<this->numPixels;i++) resid[i] = this->array[i] - this->recon[i];
18  float blankval = this->par.getBlankPixVal();
19
20  long *fpixel = new long[this->numDim];
21  for(int i=0;i<this->numDim;i++) fpixel[i]=1;
22  int status = 0;  /* MUST initialize status */
23  fitsfile *fptrOld, *fptrNew;         
24  fits_open_file(&fptrOld,this->par.getImageFile().c_str(),READONLY,&status);
25  if (status) fits_report_error(stderr, status);
26 
27  if(this->par.getFlagOutputRecon()){
28    string fileout = "!" + this->par.getImageFile(); // ! so that it writes over an existing file.
29    fileout = fileout.substr(0,fileout.size()-5); // remove the ".fits" on the end.
30    std::stringstream ss;
31    ss << ".RECON"<<this->par.getAtrousCut()<<".fits";
32    fileout += ss.str();
33
34    status = 0;
35    fits_create_file(&fptrNew,fileout.c_str(),&status);
36    if (status) fits_report_error(stderr, status);
37    status = 0;
38    fits_movabs_hdu(fptrOld, 1, NULL, &status);
39    if (status) fits_report_error(stderr, status);
40    status = 0;
41    fits_copy_header(fptrOld, fptrNew, &status);
42    if (status) fits_report_error(stderr, status);
43    char *comment = new char[80];
44    long dud;
45    status = 0;
46    fits_update_key(fptrNew, TINT, "BITPIX", &newbitpix, "number of bits per data pixel", &status);
47    if (status) fits_report_error(stderr, status);
48    status = 0;
49    float bscale=1., bzero=0.;
50    fits_update_key(fptrNew, TFLOAT, "BSCALE", &bscale, "PHYSICAL = PIXEL*BSCALE + BZERO", &status);
51    fits_update_key(fptrNew, TFLOAT, "BZERO", &bzero, "", &status);
52    fits_set_bscale(fptrNew, 1., 0., &status);
53    if (status) fits_report_error(stderr, status);
54    // Need to correct the dimensions, if we have subsectioned the image...   
55    if(this->par.getFlagSubsection()){
56      fits_read_key(fptrOld, TLONG, "NAXIS1", &dud, comment, &status);
57      fits_update_key(fptrNew, TLONG, "NAXIS1", &(this->axisDim[0]), comment, &status);
58      fits_read_key(fptrOld, TLONG, "NAXIS2", &dud, comment, &status);
59      fits_update_key(fptrNew, TLONG, "NAXIS2", &(this->axisDim[1]), comment, &status);
60      fits_read_key(fptrOld, TLONG, "NAXIS3", &dud, comment, &status);
61      fits_update_key(fptrNew, TLONG, "NAXIS3", &(this->axisDim[2]), comment, &status);
62    }
63
64    write_header_info(fptrNew, this->par, "recon");
65
66    if(this->par.getFlagBlankPix())
67      fits_write_pixnull(fptrNew, TFLOAT, fpixel, this->numPixels, this->recon, &blankval, &status);
68    else  fits_write_pix(fptrNew, TFLOAT, fpixel, this->numPixels, this->recon, &status);
69
70    status = 0;
71    fits_close_file(fptrNew, &status);
72    if (status) fits_report_error(stderr, status);
73  }
74
75  if(this->par.getFlagOutputResid()){
76    string fileout = "!" + this->par.getImageFile(); // ! so that it writes over an existing file.
77    fileout = fileout.substr(0,fileout.size()-5);
78    std::stringstream ss;
79    ss << ".RESID"<<this->par.getAtrousCut()<<".fits";
80    fileout += ss.str();
81
82    status = 0;
83    fits_create_file(&fptrNew,fileout.c_str(),&status);
84    if (status) fits_report_error(stderr, status);
85    status = 0;
86    fits_movabs_hdu(fptrOld, 1, NULL, &status);
87    if (status) fits_report_error(stderr, status);
88    status = 0;
89    fits_copy_header(fptrOld, fptrNew, &status);
90    if (status) fits_report_error(stderr, status);
91
92    status = 0;
93    fits_update_key(fptrNew, TINT, "BITPIX", &newbitpix, "number of bits per data pixel", &status);
94    if (status) fits_report_error(stderr, status);
95    status = 0;
96    float bscale=1., bzero=0.;
97    fits_update_key(fptrNew, TFLOAT, "BSCALE", &bscale, "PHYSICAL = PIXEL*BSCALE + BZERO", &status);
98    fits_update_key(fptrNew, TFLOAT, "BZERO", &bzero, "", &status);
99    fits_set_bscale(fptrNew, 1., 0., &status);
100    if (status) fits_report_error(stderr, status);
101
102    // Need to correct the dimensions, if we have subsectioned the image...
103    char *comment = new char[80];
104    long dud;
105    if(this->pars().getFlagSubsection()){
106      fits_read_key(fptrOld, TLONG, "NAXIS1", &dud, comment, &status);
107      fits_update_key(fptrNew, TLONG, "NAXIS1", &(this->axisDim[0]), comment, &status);
108      fits_read_key(fptrOld, TLONG, "NAXIS2", &dud, comment, &status);
109      fits_update_key(fptrNew, TLONG, "NAXIS2", &(this->axisDim[1]), comment, &status);
110      fits_read_key(fptrOld, TLONG, "NAXIS3", &dud, comment, &status);
111      fits_update_key(fptrNew, TLONG, "NAXIS3", &(this->axisDim[2]), comment, &status);
112    }
113
114    write_header_info(fptrNew, this->par, "resid");
115
116    if(this->par.getFlagBlankPix())
117      fits_write_pixnull(fptrNew, TFLOAT, fpixel, this->numPixels, resid, &blankval, &status);
118    else  fits_write_pix(fptrNew, TFLOAT, fpixel, this->numPixels, resid, &status);
119
120    fits_close_file(fptrNew, &status);
121  }
122
123  delete [] recon;
124  delete [] resid;
125  delete [] fpixel;
126
127}
128
129
130void write_header_info(fitsfile *fptr, Param &par, string nature)
131{
132  int status = 0;
133  char *comment, *keyword;
134  string explanation = "",ReconResid="";
135
136  fits_write_history(fptr, (char *)header_history1.c_str(), &status);
137                                   
138  fits_write_history(fptr, (char *)header_history2.c_str(), &status);
139
140  fits_write_comment(fptr, (char *)header_comment.c_str(), &status);
141
142  float valf = par.getAtrousCut();
143  fits_write_key(fptr, TFLOAT, (char *)keyword_snrRecon.c_str(), &valf,
144                 (char *)comment_snrRecon.c_str(), &status);
145
146  int vali = par.getMinScale();
147  fits_write_key(fptr, TINT, (char *)keyword_scaleMin.c_str(), &vali,
148                 (char *)comment_scaleMin.c_str(), &status);
149
150  vali = par.getFilterCode();
151  fits_write_key(fptr, TINT, (char *)keyword_filterCode.c_str(), &vali,
152                 (char *)comment_filterCode.c_str(), &status);
153
154  if(par.getFlagMW()){
155    fits_write_comment(fptr, (char *)flagMW_comment.c_str(), &status);
156    vali = par.getMinMW();
157    fits_write_key(fptr, TINT, (char *)keyword_minMW.c_str(), &vali,
158                   (char *)comment_minMW.c_str(), &status);
159    vali = par.getMaxMW();
160    fits_write_key(fptr, TINT, (char *)keyword_maxMW.c_str(), &vali,
161                   (char *)comment_maxMW.c_str(), &status);
162  }
163
164  if(nature == "recon"){
165    explanation = "Duchamp: This is the RECONSTRUCTED cube";
166    ReconResid = "RECON";
167  }
168  else if(nature == "resid"){
169    explanation = "Duchamp: This is the RESIDUAL cube";
170    ReconResid = "RESID";
171  }
172  else std::cerr << "WARNING : write_header_info : explanation not present!\n";
173  fits_write_comment(fptr, (char *)explanation.c_str(), &status);
174  fits_write_key(fptr, TSTRING, (char *)keyword_ReconResid.c_str(),
175                 (char *)ReconResid.c_str(), (char *)comment_ReconResid.c_str(), &status);
176
177}
Note: See TracBrowser for help on using the repository browser.