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

Last change on this file since 221 was 220, checked in by Matthew Whiting, 17 years ago
  • Two new files: plots.cc and feedback.cc. Introduced to separate the declarations and definitions for various classes.
  • Mostly just improving the documentation for use with Doxygen.
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
14void writeReconHeaderInfo(fitsfile *fptr, Param &par, string nature);
15void writeSmoothHeaderInfo(fitsfile *fptr, Param &par);
16
17void Cube::saveSmoothedCube()
18{
19  /**
20   *   A function to save the Hanning-smoothed arrays to a FITS file.
21   *   Additional header keywords are written as well, indicating the
22   *    width of the Hanning filter.
23   *   The file is always written -- if the filename (as calculated
24   *    based on the parameters) exists, then it is overwritten.
25   */
26 
27  int newbitpix = FLOAT_IMG;
28
29  float blankval = this->par.getBlankPixVal();
30
31  long *fpixel = new long[this->numDim];
32  for(int i=0;i<this->numDim;i++) fpixel[i]=1;
33  int status = 0;  /* MUST initialize status */
34  fitsfile *fptrOld, *fptrNew;         
35  fits_open_file(&fptrOld,this->par.getImageFile().c_str(),READONLY,&status);
36  if (status) fits_report_error(stderr, status);
37 
38  if(this->par.getFlagOutputSmooth()){
39    string fileout = "!" + this->par.outputSmoothFile();
40    // the ! is there so that it writes over an existing file.
41
42    status = 0;
43    fits_create_file(&fptrNew,fileout.c_str(),&status);
44    if (status)
45      fits_report_error(stderr, status);
46    else {
47      status = 0;
48      fits_movabs_hdu(fptrOld, 1, NULL, &status);
49      if (status) fits_report_error(stderr, status);
50      status = 0;
51      fits_copy_header(fptrOld, fptrNew, &status);
52      if (status) fits_report_error(stderr, status);
53      char *comment = new char[80];
54      long dud;
55      status = 0;
56      fits_update_key(fptrNew, TINT, "BITPIX", &newbitpix,
57                      "number of bits per data pixel", &status);
58      if (status) fits_report_error(stderr, status);
59      status = 0;
60      float bscale=1., bzero=0.;
61      fits_update_key(fptrNew, TFLOAT, "BSCALE", &bscale,
62                      "PHYSICAL = PIXEL*BSCALE + BZERO", &status);
63      fits_update_key(fptrNew, TFLOAT, "BZERO", &bzero, "", &status);
64      fits_set_bscale(fptrNew, 1., 0., &status);
65      if (status) fits_report_error(stderr, status);
66      // Need to correct the dimensions, if we have subsectioned the image
67      if(this->par.getFlagSubsection()){
68        fits_read_key(fptrOld, TLONG, "NAXIS1", &dud, comment, &status);
69        fits_update_key(fptrNew, TLONG, "NAXIS1",
70                        &(this->axisDim[0]), comment, &status);
71        fits_read_key(fptrOld, TLONG, "NAXIS2", &dud, comment, &status);
72        fits_update_key(fptrNew, TLONG, "NAXIS2",
73                        &(this->axisDim[1]), comment, &status);
74        fits_read_key(fptrOld, TLONG, "NAXIS3", &dud, comment, &status);
75        fits_update_key(fptrNew, TLONG, "NAXIS3",
76                        &(this->axisDim[2]), comment, &status);
77      }
78
79      writeSmoothHeaderInfo(fptrNew, this->par);
80
81      if(this->par.getFlagBlankPix())
82        fits_write_pixnull(fptrNew, TFLOAT, fpixel, this->numPixels,
83                           this->recon, &blankval, &status);
84      else fits_write_pix(fptrNew, TFLOAT, fpixel, this->numPixels,
85                          this->recon, &status);
86
87      status = 0;
88      fits_close_file(fptrNew, &status);
89      if (status) fits_report_error(stderr, status);
90    }
91  }
92
93}
94
95
96void Cube::saveReconstructedCube()
97{
98  /**
99   *  A function to save the reconstructed and/or residual arrays.
100   *   A number of header keywords are written as well, indicating the
101   *    nature of the reconstruction that has been done.
102   *   The file is always written -- if the filename (as calculated
103   *    based on the recon parameters) exists, then it is overwritten.
104   */
105 
106  int newbitpix = FLOAT_IMG;
107
108  float *resid = new float[this->numPixels];
109  for(int i=0;i<this->numPixels;i++)
110    resid[i] = this->array[i] - this->recon[i];
111  float blankval = this->par.getBlankPixVal();
112
113  long *fpixel = new long[this->numDim];
114  for(int i=0;i<this->numDim;i++) fpixel[i]=1;
115  int status = 0;  /* MUST initialize status */
116  fitsfile *fptrOld, *fptrNew;         
117  fits_open_file(&fptrOld,this->par.getImageFile().c_str(),READONLY,&status);
118  if (status) fits_report_error(stderr, status);
119 
120  if(this->par.getFlagOutputRecon()){
121    string fileout = "!" + this->par.outputReconFile();
122    // the ! is there so that it writes over an existing file.
123
124    status = 0;
125    fits_create_file(&fptrNew,fileout.c_str(),&status);
126    if (status)
127      fits_report_error(stderr, status);
128    else
129      {
130        status = 0;
131        fits_movabs_hdu(fptrOld, 1, NULL, &status);
132        if (status) fits_report_error(stderr, status);
133        status = 0;
134        fits_copy_header(fptrOld, fptrNew, &status);
135        if (status) fits_report_error(stderr, status);
136        char *comment = new char[80];
137        long dud;
138        status = 0;
139        fits_update_key(fptrNew, TINT, "BITPIX", &newbitpix,
140                        "number of bits per data pixel", &status);
141        if (status) fits_report_error(stderr, status);
142        status = 0;
143        float bscale=1., bzero=0.;
144        fits_update_key(fptrNew, TFLOAT, "BSCALE", &bscale,
145                        "PHYSICAL = PIXEL*BSCALE + BZERO", &status);
146        fits_update_key(fptrNew, TFLOAT, "BZERO", &bzero, "", &status);
147        fits_set_bscale(fptrNew, 1., 0., &status);
148        if (status) fits_report_error(stderr, status);
149        // Need to correct the dimensions, if we have subsectioned the image
150        if(this->par.getFlagSubsection()){
151          fits_read_key(fptrOld, TLONG, "NAXIS1", &dud, comment, &status);
152          fits_update_key(fptrNew, TLONG, "NAXIS1",
153                          &(this->axisDim[0]), comment, &status);
154          fits_read_key(fptrOld, TLONG, "NAXIS2", &dud, comment, &status);
155          fits_update_key(fptrNew, TLONG, "NAXIS2",
156                          &(this->axisDim[1]), comment, &status);
157          fits_read_key(fptrOld, TLONG, "NAXIS3", &dud, comment, &status);
158          fits_update_key(fptrNew, TLONG, "NAXIS3",
159                          &(this->axisDim[2]), comment, &status);
160        }
161
162        writeReconHeaderInfo(fptrNew, this->par, "recon");
163
164        if(this->par.getFlagBlankPix())
165          fits_write_pixnull(fptrNew, TFLOAT, fpixel, this->numPixels,
166                             this->recon, &blankval, &status);
167        else  fits_write_pix(fptrNew, TFLOAT, fpixel, this->numPixels,
168                             this->recon, &status);
169
170        status = 0;
171        fits_close_file(fptrNew, &status);
172        if (status) fits_report_error(stderr, status);
173      }
174  }
175
176
177  if(this->par.getFlagOutputResid()){
178    string fileout = "!" + this->par.outputResidFile();
179    // the ! is there so that it writes over an existing file.
180    status = 0;
181    fits_create_file(&fptrNew,fileout.c_str(),&status);
182    if (status)
183      fits_report_error(stderr, status);
184    else
185      {
186        status = 0;
187        fits_movabs_hdu(fptrOld, 1, NULL, &status);
188        if (status) fits_report_error(stderr, status);
189        status = 0;
190        fits_copy_header(fptrOld, fptrNew, &status);
191        if (status) fits_report_error(stderr, status);
192
193        status = 0;
194        fits_update_key(fptrNew, TINT, "BITPIX", &newbitpix,
195                        "number of bits per data pixel", &status);
196        if (status) fits_report_error(stderr, status);
197        status = 0;
198        float bscale=1., bzero=0.;
199        fits_update_key(fptrNew, TFLOAT, "BSCALE", &bscale,
200                        "PHYSICAL = PIXEL*BSCALE + BZERO", &status);
201        fits_update_key(fptrNew, TFLOAT, "BZERO", &bzero, "", &status);
202        fits_set_bscale(fptrNew, 1., 0., &status);
203        if (status) fits_report_error(stderr, status);
204
205        // Need to correct the dimensions, if we have subsectioned the image...
206        char *comment = new char[80];
207        long dud;
208        if(this->pars().getFlagSubsection()){
209          fits_read_key(fptrOld, TLONG, "NAXIS1", &dud, comment, &status);
210          fits_update_key(fptrNew, TLONG, "NAXIS1",
211                          &(this->axisDim[0]), comment, &status);
212          fits_read_key(fptrOld, TLONG, "NAXIS2", &dud, comment, &status);
213          fits_update_key(fptrNew, TLONG, "NAXIS2",
214                          &(this->axisDim[1]), comment, &status);
215          fits_read_key(fptrOld, TLONG, "NAXIS3", &dud, comment, &status);
216          fits_update_key(fptrNew, TLONG, "NAXIS3",
217                          &(this->axisDim[2]), comment, &status);
218        }
219
220        writeReconHeaderInfo(fptrNew, this->par, "resid");
221
222        if(this->par.getFlagBlankPix())
223          fits_write_pixnull(fptrNew, TFLOAT, fpixel, this->numPixels,
224                             resid, &blankval, &status);
225        else  fits_write_pix(fptrNew, TFLOAT, fpixel, this->numPixels,
226                             resid, &status);
227
228        fits_close_file(fptrNew, &status);
229      }
230  }
231
232  delete [] resid;
233  delete [] fpixel;
234
235}
236
237
238void writeReconHeaderInfo(fitsfile *fptr, Param &par, string nature)
239{
240  /**
241   *  writeReconHeaderInfo(fptr, par, nature)
242   *
243   *   A simple function that writes all the necessary keywords and comments
244   *    to the FITS header pointed to by fptr.
245   *   The keyword names and comments are taken from duchamp.hh
246   *   The parameter "nature" indicates what type of file is being written:
247   *    should be either "recon" or "resid".
248   */
249
250
251  int status = 0;
252  char *comment, *keyword;
253  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   *  writeSmoothHeaderInfo(fptr, par)
309   *
310   *   A simple function that writes all the necessary keywords and comments
311   *    to the FITS header pointed to by fptr.
312   *   The keyword names and comments are taken from duchamp.hh
313   */
314
315
316  int status = 0;
317  char *comment, *keyword;
318
319  fits_write_history(fptr, (char *)header_smoothHistory.c_str(), &status);
320  status = 0;
321  fits_write_history(fptr, (char *)header_smoothHistory_input.c_str(),&status);
322  status = 0;
323  fits_write_history(fptr, (char *)par.getImageFile().c_str(), &status);
324
325  if(par.getFlagSubsection()){
326    status = 0;
327    fits_write_comment(fptr,(char *)header_smoothSubsection_comment.c_str(),
328                       &status);
329    status = 0;
330    fits_write_key(fptr, TSTRING, (char *)keyword_subsection.c_str(),
331                   (char *)par.getSubsection().c_str(),
332                   (char *)comment_subsection.c_str(), &status);
333  }
334   
335  int val = par.getHanningWidth();
336  fits_write_key(fptr, TINT, (char *)keyword_hanningwidth.c_str(), &val,
337                 (char *)comment_hanningwidth.c_str(), &status);
338
339}
Note: See TracBrowser for help on using the repository browser.