source: trunk/src/Cubes/saveImage.cc @ 272

Last change on this file since 272 was 270, checked in by Matthew Whiting, 17 years ago

A large set of changes, each of which small ones from compiling with the -Wall flag (plus the -Wno-sign-compare flag -- as we don't care about warnings about comparing ints and unsigned ints).
File size: 11.4 KB
Line 
1#include <iostream>
2#include <sstream>
3#include <string>
4#include <wcs.h>
5#include <wcshdr.h>
6#define WCSLIB_GETWCSTAB
7// define this so that we don't try and redefine wtbarr
8// (this is a problem when using cfitsio v.3 and g++ v.4)
9
10#include <fitsio.h>
11#include <duchamp.hh>
12#include <Cubes/cubes.hh>
13
14/** Write FITS headers in correct format for reconstructed array output */
15void writeReconHeaderInfo(fitsfile *fptr, Param &par, std::string nature);
16
17/** Write FITS headers in correct format for smoothed array output */
18void writeSmoothHeaderInfo(fitsfile *fptr, Param &par);
19
20void Cube::saveSmoothedCube()
21{
22  /**
23   *   A function to save the Hanning-smoothed arrays to a FITS file.
24   *   Additional header keywords are written as well, indicating the
25   *    width of the Hanning filter.
26   *   The file is always written -- if the filename (as calculated
27   *    based on the parameters) exists, then it is overwritten.
28   */
29 
30  int newbitpix = FLOAT_IMG;
31
32  float blankval = this->par.getBlankPixVal();
33
34  long *fpixel = new long[this->numDim];
35  for(int i=0;i<this->numDim;i++) fpixel[i]=1;
36  int status = 0;  /* MUST initialize status */
37  fitsfile *fptrOld, *fptrNew;         
38  fits_open_file(&fptrOld,this->par.getImageFile().c_str(),READONLY,&status);
39  if (status) fits_report_error(stderr, status);
40 
41  if(this->par.getFlagOutputSmooth()){
42    std::string fileout = "!" + this->par.outputSmoothFile();
43    // the ! is there so that it writes over an existing file.
44
45    status = 0;
46    fits_create_file(&fptrNew,fileout.c_str(),&status);
47    if (status)
48      fits_report_error(stderr, status);
49    else {
50      status = 0;
51      fits_movabs_hdu(fptrOld, 1, NULL, &status);
52      if (status) fits_report_error(stderr, status);
53      status = 0;
54      fits_copy_header(fptrOld, fptrNew, &status);
55      if (status) fits_report_error(stderr, status);
56      char *comment = new char[80];
57      long dud;
58      status = 0;
59      fits_update_key(fptrNew, TINT, "BITPIX", &newbitpix,
60                      "number of bits per data pixel", &status);
61      if (status) fits_report_error(stderr, status);
62      status = 0;
63      float bscale=1., bzero=0.;
64      fits_update_key(fptrNew, TFLOAT, "BSCALE", &bscale,
65                      "PHYSICAL = PIXEL*BSCALE + BZERO", &status);
66      fits_update_key(fptrNew, TFLOAT, "BZERO", &bzero, "", &status);
67      fits_set_bscale(fptrNew, 1., 0., &status);
68      if (status) fits_report_error(stderr, status);
69      // Need to correct the dimensions, if we have subsectioned the image
70      if(this->par.getFlagSubsection()){
71        fits_read_key(fptrOld, TLONG, "NAXIS1", &dud, comment, &status);
72        fits_update_key(fptrNew, TLONG, "NAXIS1",
73                        &(this->axisDim[0]), comment, &status);
74        fits_read_key(fptrOld, TLONG, "NAXIS2", &dud, comment, &status);
75        fits_update_key(fptrNew, TLONG, "NAXIS2",
76                        &(this->axisDim[1]), comment, &status);
77        fits_read_key(fptrOld, TLONG, "NAXIS3", &dud, comment, &status);
78        fits_update_key(fptrNew, TLONG, "NAXIS3",
79                        &(this->axisDim[2]), comment, &status);
80      }
81
82      writeSmoothHeaderInfo(fptrNew, this->par);
83
84      if(this->par.getFlagBlankPix())
85        fits_write_pixnull(fptrNew, TFLOAT, fpixel, this->numPixels,
86                           this->recon, &blankval, &status);
87      else fits_write_pix(fptrNew, TFLOAT, fpixel, this->numPixels,
88                          this->recon, &status);
89
90      status = 0;
91      fits_close_file(fptrNew, &status);
92      if (status) fits_report_error(stderr, status);
93    }
94  }
95
96}
97
98
99void Cube::saveReconstructedCube()
100{
101  /**
102   *  A function to save the reconstructed and/or residual arrays.
103   *   A number of header keywords are written as well, indicating the
104   *    nature of the reconstruction that has been done.
105   *   The file is always written -- if the filename (as calculated
106   *    based on the recon parameters) exists, then it is overwritten.
107   */
108 
109  int newbitpix = FLOAT_IMG;
110
111  float *resid = new float[this->numPixels];
112  for(int i=0;i<this->numPixels;i++)
113    resid[i] = this->array[i] - this->recon[i];
114  float blankval = this->par.getBlankPixVal();
115
116  long *fpixel = new long[this->numDim];
117  for(int i=0;i<this->numDim;i++) fpixel[i]=1;
118  int status = 0;  /* MUST initialize status */
119  fitsfile *fptrOld, *fptrNew;         
120  fits_open_file(&fptrOld,this->par.getImageFile().c_str(),READONLY,&status);
121  if (status) fits_report_error(stderr, status);
122 
123  if(this->par.getFlagOutputRecon()){
124    std::string fileout = "!" + this->par.outputReconFile();
125    // the ! is there so that it writes over an existing file.
126
127    status = 0;
128    fits_create_file(&fptrNew,fileout.c_str(),&status);
129    if (status)
130      fits_report_error(stderr, status);
131    else
132      {
133        status = 0;
134        fits_movabs_hdu(fptrOld, 1, NULL, &status);
135        if (status) fits_report_error(stderr, status);
136        status = 0;
137        fits_copy_header(fptrOld, fptrNew, &status);
138        if (status) fits_report_error(stderr, status);
139        char *comment = new char[80];
140        long dud;
141        status = 0;
142        fits_update_key(fptrNew, TINT, "BITPIX", &newbitpix,
143                        "number of bits per data pixel", &status);
144        if (status) fits_report_error(stderr, status);
145        status = 0;
146        float bscale=1., bzero=0.;
147        fits_update_key(fptrNew, TFLOAT, "BSCALE", &bscale,
148                        "PHYSICAL = PIXEL*BSCALE + BZERO", &status);
149        fits_update_key(fptrNew, TFLOAT, "BZERO", &bzero, "", &status);
150        fits_set_bscale(fptrNew, 1., 0., &status);
151        if (status) fits_report_error(stderr, status);
152        // Need to correct the dimensions, if we have subsectioned the image
153        if(this->par.getFlagSubsection()){
154          fits_read_key(fptrOld, TLONG, "NAXIS1", &dud, comment, &status);
155          fits_update_key(fptrNew, TLONG, "NAXIS1",
156                          &(this->axisDim[0]), comment, &status);
157          fits_read_key(fptrOld, TLONG, "NAXIS2", &dud, comment, &status);
158          fits_update_key(fptrNew, TLONG, "NAXIS2",
159                          &(this->axisDim[1]), comment, &status);
160          fits_read_key(fptrOld, TLONG, "NAXIS3", &dud, comment, &status);
161          fits_update_key(fptrNew, TLONG, "NAXIS3",
162                          &(this->axisDim[2]), comment, &status);
163        }
164
165        writeReconHeaderInfo(fptrNew, this->par, "recon");
166
167        if(this->par.getFlagBlankPix())
168          fits_write_pixnull(fptrNew, TFLOAT, fpixel, this->numPixels,
169                             this->recon, &blankval, &status);
170        else  fits_write_pix(fptrNew, TFLOAT, fpixel, this->numPixels,
171                             this->recon, &status);
172
173        status = 0;
174        fits_close_file(fptrNew, &status);
175        if (status) fits_report_error(stderr, status);
176      }
177  }
178
179
180  if(this->par.getFlagOutputResid()){
181    std::string fileout = "!" + this->par.outputResidFile();
182    // the ! is there so that it writes over an existing file.
183    status = 0;
184    fits_create_file(&fptrNew,fileout.c_str(),&status);
185    if (status)
186      fits_report_error(stderr, status);
187    else
188      {
189        status = 0;
190        fits_movabs_hdu(fptrOld, 1, NULL, &status);
191        if (status) fits_report_error(stderr, status);
192        status = 0;
193        fits_copy_header(fptrOld, fptrNew, &status);
194        if (status) fits_report_error(stderr, status);
195
196        status = 0;
197        fits_update_key(fptrNew, TINT, "BITPIX", &newbitpix,
198                        "number of bits per data pixel", &status);
199        if (status) fits_report_error(stderr, status);
200        status = 0;
201        float bscale=1., bzero=0.;
202        fits_update_key(fptrNew, TFLOAT, "BSCALE", &bscale,
203                        "PHYSICAL = PIXEL*BSCALE + BZERO", &status);
204        fits_update_key(fptrNew, TFLOAT, "BZERO", &bzero, "", &status);
205        fits_set_bscale(fptrNew, 1., 0., &status);
206        if (status) fits_report_error(stderr, status);
207
208        // Need to correct the dimensions, if we have subsectioned the image...
209        char *comment = new char[80];
210        long dud;
211        if(this->pars().getFlagSubsection()){
212          fits_read_key(fptrOld, TLONG, "NAXIS1", &dud, comment, &status);
213          fits_update_key(fptrNew, TLONG, "NAXIS1",
214                          &(this->axisDim[0]), comment, &status);
215          fits_read_key(fptrOld, TLONG, "NAXIS2", &dud, comment, &status);
216          fits_update_key(fptrNew, TLONG, "NAXIS2",
217                          &(this->axisDim[1]), comment, &status);
218          fits_read_key(fptrOld, TLONG, "NAXIS3", &dud, comment, &status);
219          fits_update_key(fptrNew, TLONG, "NAXIS3",
220                          &(this->axisDim[2]), comment, &status);
221        }
222
223        writeReconHeaderInfo(fptrNew, this->par, "resid");
224
225        if(this->par.getFlagBlankPix())
226          fits_write_pixnull(fptrNew, TFLOAT, fpixel, this->numPixels,
227                             resid, &blankval, &status);
228        else  fits_write_pix(fptrNew, TFLOAT, fpixel, this->numPixels,
229                             resid, &status);
230
231        fits_close_file(fptrNew, &status);
232      }
233  }
234
235  delete [] resid;
236  delete [] fpixel;
237
238}
239
240
241void writeReconHeaderInfo(fitsfile *fptr, Param &par, std::string nature)
242{
243  /**
244   *   A simple function that writes all the necessary keywords and comments
245   *    to the FITS header pointed to by fptr.
246   *   The keyword names and comments are taken from duchamp.hh
247   *   The parameter "nature" indicates what type of file is being written:
248   *    should be either "recon" or "resid".
249   */
250
251
252  int status = 0;
253  std::string explanation = "",ReconResid="";
254
255  fits_write_history(fptr, (char *)header_reconHistory1.c_str(), &status);
256                                   
257  fits_write_history(fptr, (char *)header_reconHistory2.c_str(), &status);
258
259  fits_write_history(fptr, (char *)header_reconHistory_input.c_str(), &status);
260
261  fits_write_history(fptr, (char *)par.getImageFile().c_str(), &status);
262
263  if(par.getFlagSubsection()){
264    fits_write_comment(fptr,(char *)header_reconSubsection_comment.c_str(),
265                       &status);
266    fits_write_key(fptr, TSTRING, (char *)keyword_subsection.c_str(),
267                   (char *)par.getSubsection().c_str(),
268                   (char *)comment_subsection.c_str(), &status);
269  }
270   
271  fits_write_comment(fptr, (char *)header_atrous_comment.c_str(), &status);
272
273  float valf = par.getAtrousCut();
274  fits_write_key(fptr, TFLOAT, (char *)keyword_snrRecon.c_str(), &valf,
275                 (char *)comment_snrRecon.c_str(), &status);
276
277  int vali = par.getReconDim();
278  fits_write_key(fptr, TINT, (char *)keyword_reconDim.c_str(), &vali,
279                 (char *)comment_reconDim.c_str(), &status);
280
281  vali = par.getMinScale();
282  fits_write_key(fptr, TINT, (char *)keyword_scaleMin.c_str(), &vali,
283                 (char *)comment_scaleMin.c_str(), &status);
284
285  vali = par.getFilterCode();
286  fits_write_key(fptr, TINT, (char *)keyword_filterCode.c_str(), &vali,
287                 (char *)comment_filterCode.c_str(), &status);
288
289  if(nature == "recon"){
290    explanation = "Duchamp: This is the RECONSTRUCTED cube";
291    ReconResid = "RECON";
292  }
293  else if(nature == "resid"){
294    explanation = "Duchamp: This is the RESIDUAL cube";
295    ReconResid = "RESID";
296  }
297  else duchampWarning("write_header_info","explanation not present!\n");
298  fits_write_comment(fptr, (char *)explanation.c_str(), &status);
299  fits_write_key(fptr, TSTRING, (char *)keyword_ReconResid.c_str(),
300                 (char *)ReconResid.c_str(),
301                 (char *)comment_ReconResid.c_str(), &status);
302
303}
304
305void writeSmoothHeaderInfo(fitsfile *fptr, Param &par)
306{
307  /**
308   *   A simple function that writes all the necessary keywords and comments
309   *    to the FITS header pointed to by fptr.
310   *   The keyword names and comments are taken from duchamp.hh
311   */
312
313
314  int status = 0;
315
316  fits_write_history(fptr, (char *)header_smoothHistory.c_str(), &status);
317  status = 0;
318  fits_write_history(fptr, (char *)header_smoothHistory_input.c_str(),&status);
319  status = 0;
320  fits_write_history(fptr, (char *)par.getImageFile().c_str(), &status);
321
322  if(par.getFlagSubsection()){
323    status = 0;
324    fits_write_comment(fptr,(char *)header_smoothSubsection_comment.c_str(),
325                       &status);
326    status = 0;
327    fits_write_key(fptr, TSTRING, (char *)keyword_subsection.c_str(),
328                   (char *)par.getSubsection().c_str(),
329                   (char *)comment_subsection.c_str(), &status);
330  }
331   
332  int val = par.getHanningWidth();
333  fits_write_key(fptr, TINT, (char *)keyword_hanningwidth.c_str(), &val,
334                 (char *)comment_hanningwidth.c_str(), &status);
335
336}
Note: See TracBrowser for help on using the repository browser.