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

Last change on this file since 208 was 208, checked in by Matthew Whiting, 18 years ago
  • Enabled saving and reading in of a smoothed array, in manner directly analogous to that for the recon array.
    • New file : src/Cubes/readSmooth.cc
    • The other new functions go in existing files eg. saveImage.cc
    • Renamed some functions (like writeHeader...) to be more obvious what they do.
    • The reading in is taken care of by new function Cube::readSavedArrays() -- handles both smoothed and recon'd arrays.
    • Associated parameters in Param class
    • Clarified names of FITS header strings in duchamp.hh.
  • Updated the documentation to describe the ability to smooth a cube.
  • Added description of feedback mechanisms in the Install appendix.
  • Also, Hanning class improved to guard against memory leaks.


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