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

Last change on this file since 221 was 220, checked in by Matthew Whiting, 18 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
RevLine 
[3]1#include <iostream>
2#include <sstream>
3#include <string>
4#include <wcs.h>
5#include <wcshdr.h>
[146]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
[66]10#include <fitsio.h>
[71]11#include <duchamp.hh>
[3]12#include <Cubes/cubes.hh>
13
[208]14void writeReconHeaderInfo(fitsfile *fptr, Param &par, string nature);
15void writeSmoothHeaderInfo(fitsfile *fptr, Param &par);
[86]16
[208]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
[3]96void Cube::saveReconstructedCube()
97{
[86]98  /**
[220]99   *  A function to save the reconstructed and/or residual arrays.
[146]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.
[86]104   */
[3]105 
[71]106  int newbitpix = FLOAT_IMG;
107
[3]108  float *resid = new float[this->numPixels];
[146]109  for(int i=0;i<this->numPixels;i++)
110    resid[i] = this->array[i] - this->recon[i];
[3]111  float blankval = this->par.getBlankPixVal();
112
[71]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 */
[3]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()){
[146]121    string fileout = "!" + this->par.outputReconFile();
122    // the ! is there so that it writes over an existing file.
[3]123
124    status = 0;
[71]125    fits_create_file(&fptrNew,fileout.c_str(),&status);
[103]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;
[146]139        fits_update_key(fptrNew, TINT, "BITPIX", &newbitpix,
140                        "number of bits per data pixel", &status);
[103]141        if (status) fits_report_error(stderr, status);
142        status = 0;
143        float bscale=1., bzero=0.;
[146]144        fits_update_key(fptrNew, TFLOAT, "BSCALE", &bscale,
145                        "PHYSICAL = PIXEL*BSCALE + BZERO", &status);
[103]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);
[146]149        // Need to correct the dimensions, if we have subsectioned the image
[103]150        if(this->par.getFlagSubsection()){
151          fits_read_key(fptrOld, TLONG, "NAXIS1", &dud, comment, &status);
[146]152          fits_update_key(fptrNew, TLONG, "NAXIS1",
153                          &(this->axisDim[0]), comment, &status);
[103]154          fits_read_key(fptrOld, TLONG, "NAXIS2", &dud, comment, &status);
[146]155          fits_update_key(fptrNew, TLONG, "NAXIS2",
156                          &(this->axisDim[1]), comment, &status);
[103]157          fits_read_key(fptrOld, TLONG, "NAXIS3", &dud, comment, &status);
[146]158          fits_update_key(fptrNew, TLONG, "NAXIS3",
159                          &(this->axisDim[2]), comment, &status);
[103]160        }
[46]161
[208]162        writeReconHeaderInfo(fptrNew, this->par, "recon");
[3]163
[103]164        if(this->par.getFlagBlankPix())
[146]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);
[3]169
[103]170        status = 0;
171        fits_close_file(fptrNew, &status);
172        if (status) fits_report_error(stderr, status);
173      }
[3]174  }
175
[103]176
[3]177  if(this->par.getFlagOutputResid()){
[146]178    string fileout = "!" + this->par.outputResidFile();
179    // the ! is there so that it writes over an existing file.
[3]180    status = 0;
[71]181    fits_create_file(&fptrNew,fileout.c_str(),&status);
[103]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);
[46]192
[103]193        status = 0;
[146]194        fits_update_key(fptrNew, TINT, "BITPIX", &newbitpix,
195                        "number of bits per data pixel", &status);
[103]196        if (status) fits_report_error(stderr, status);
197        status = 0;
198        float bscale=1., bzero=0.;
[146]199        fits_update_key(fptrNew, TFLOAT, "BSCALE", &bscale,
200                        "PHYSICAL = PIXEL*BSCALE + BZERO", &status);
[103]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);
[71]204
[103]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);
[146]210          fits_update_key(fptrNew, TLONG, "NAXIS1",
211                          &(this->axisDim[0]), comment, &status);
[103]212          fits_read_key(fptrOld, TLONG, "NAXIS2", &dud, comment, &status);
[146]213          fits_update_key(fptrNew, TLONG, "NAXIS2",
214                          &(this->axisDim[1]), comment, &status);
[103]215          fits_read_key(fptrOld, TLONG, "NAXIS3", &dud, comment, &status);
[146]216          fits_update_key(fptrNew, TLONG, "NAXIS3",
217                          &(this->axisDim[2]), comment, &status);
[103]218        }
[3]219
[208]220        writeReconHeaderInfo(fptrNew, this->par, "resid");
[71]221
[103]222        if(this->par.getFlagBlankPix())
[146]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);
[3]227
[103]228        fits_close_file(fptrNew, &status);
229      }
[3]230  }
231
232  delete [] resid;
233  delete [] fpixel;
234
235}
[71]236
237
[208]238void writeReconHeaderInfo(fitsfile *fptr, Param &par, string nature)
[71]239{
[86]240  /**
[208]241   *  writeReconHeaderInfo(fptr, par, nature)
[86]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
[71]251  int status = 0;
252  char *comment, *keyword;
253  string explanation = "",ReconResid="";
254
[208]255  fits_write_history(fptr, (char *)header_reconHistory1.c_str(), &status);
[71]256                                   
[208]257  fits_write_history(fptr, (char *)header_reconHistory2.c_str(), &status);
[71]258
[208]259  fits_write_history(fptr, (char *)header_reconHistory_input.c_str(), &status);
[71]260
[105]261  fits_write_history(fptr, (char *)par.getImageFile().c_str(), &status);
262
263  if(par.getFlagSubsection()){
[208]264    fits_write_comment(fptr,(char *)header_reconSubsection_comment.c_str(),
265                       &status);
[105]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
[71]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
[103]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();
[71]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  }
[120]297  else duchampWarning("write_header_info","explanation not present!\n");
[71]298  fits_write_comment(fptr, (char *)explanation.c_str(), &status);
299  fits_write_key(fptr, TSTRING, (char *)keyword_ReconResid.c_str(),
[146]300                 (char *)ReconResid.c_str(),
301                 (char *)comment_ReconResid.c_str(), &status);
[71]302
303}
[208]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.