source: tags/release-1.1/src/Cubes/saveImage.cc @ 1391

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

Adding distribution text at the start of each file...
File size: 13.8 KB
Line 
1// -----------------------------------------------------------------------
2// saveImage.cc: Write a wavelet-reconstructed or smoothed array to a
3//               FITS file.
4// -----------------------------------------------------------------------
5// Copyright (C) 2006, Matthew Whiting, ATNF
6//
7// This program is free software; you can redistribute it and/or modify it
8// under the terms of the GNU General Public License as published by the
9// Free Software Foundation; either version 2 of the License, or (at your
10// option) any later version.
11//
12// Duchamp is distributed in the hope that it will be useful, but WITHOUT
13// ANY WARRANTY; without even the implied warranty of MERCHANTABILITY or
14// FITNESS FOR A PARTICULAR PURPOSE.  See the GNU General Public License
15// for more details.
16//
17// You should have received a copy of the GNU General Public License
18// along with Duchamp; if not, write to the Free Software Foundation,
19// Inc., 59 Temple Place, Suite 330, Boston, MA 02111-1307, USA
20//
21// Correspondence concerning Duchamp may be directed to:
22//    Internet email: Matthew.Whiting [at] atnf.csiro.au
23//    Postal address: Dr. Matthew Whiting
24//                    Australia Telescope National Facility, CSIRO
25//                    PO Box 76
26//                    Epping NSW 1710
27//                    AUSTRALIA
28// -----------------------------------------------------------------------
29#include <iostream>
30#include <sstream>
31#include <string>
32#include <wcs.h>
33#include <wcshdr.h>
34#define WCSLIB_GETWCSTAB
35// define this so that we don't try and redefine wtbarr
36// (this is a problem when using cfitsio v.3 and g++ v.4)
37
38#include <fitsio.h>
39#include <duchamp.hh>
40#include <Cubes/cubes.hh>
41
42/** Write FITS headers in correct format for reconstructed array output */
43void writeReconHeaderInfo(fitsfile *fptr, Param &par, std::string nature);
44
45/** Write FITS headers in correct format for smoothed array output */
46void writeSmoothHeaderInfo(fitsfile *fptr, Param &par);
47
48void Cube::saveSmoothedCube()
49{
50  /**
51   *   A function to save the smoothed arrays to a FITS file.
52   *   Additional header keywords are written as well, indicating the
53   *   width of the Hanning filter or the dimensions of the Gaussian
54   *   kernel.
55   *   The file is always written -- if the filename (as calculated
56   *    based on the parameters) exists, then it is overwritten.
57   */
58 
59  int newbitpix = FLOAT_IMG;
60
61  float blankval = this->par.getBlankPixVal();
62
63  long *fpixel = new long[this->numDim];
64  for(int i=0;i<this->numDim;i++) fpixel[i]=1;
65  int status = 0;  /* MUST initialize status */
66  fitsfile *fptrOld, *fptrNew;         
67  fits_open_file(&fptrOld,this->par.getImageFile().c_str(),READONLY,&status);
68  if (status) fits_report_error(stderr, status);
69 
70  if(this->par.getFlagOutputSmooth()){
71    std::string fileout = "!" + this->par.outputSmoothFile();
72    // the ! is there so that it writes over an existing file.
73
74    status = 0;
75    fits_create_file(&fptrNew,fileout.c_str(),&status);
76    if (status)
77      fits_report_error(stderr, status);
78    else {
79      status = 0;
80      fits_movabs_hdu(fptrOld, 1, NULL, &status);
81      if (status) fits_report_error(stderr, status);
82      status = 0;
83      fits_copy_header(fptrOld, fptrNew, &status);
84      if (status) fits_report_error(stderr, status);
85      char *comment = new char[80];
86      long dud;
87      status = 0;
88      fits_update_key(fptrNew, TINT, "BITPIX", &newbitpix,
89                      "number of bits per data pixel", &status);
90      if (status) fits_report_error(stderr, status);
91      status = 0;
92      float bscale=1., bzero=0.;
93      fits_update_key(fptrNew, TFLOAT, "BSCALE", &bscale,
94                      "PHYSICAL = PIXEL*BSCALE + BZERO", &status);
95      fits_update_key(fptrNew, TFLOAT, "BZERO", &bzero, "", &status);
96      fits_set_bscale(fptrNew, 1., 0., &status);
97      if (status) fits_report_error(stderr, status);
98      // Need to correct the dimensions, if we have subsectioned the image
99      if(this->par.getFlagSubsection()){
100        fits_read_key(fptrOld, TLONG, "NAXIS1", &dud, comment, &status);
101        fits_update_key(fptrNew, TLONG, "NAXIS1",
102                        &(this->axisDim[0]), comment, &status);
103        fits_read_key(fptrOld, TLONG, "NAXIS2", &dud, comment, &status);
104        fits_update_key(fptrNew, TLONG, "NAXIS2",
105                        &(this->axisDim[1]), comment, &status);
106        fits_read_key(fptrOld, TLONG, "NAXIS3", &dud, comment, &status);
107        fits_update_key(fptrNew, TLONG, "NAXIS3",
108                        &(this->axisDim[2]), comment, &status);
109      }
110
111      writeSmoothHeaderInfo(fptrNew, this->par);
112
113      if(this->par.getFlagBlankPix())
114        fits_write_pixnull(fptrNew, TFLOAT, fpixel, this->numPixels,
115                           this->recon, &blankval, &status);
116      else fits_write_pix(fptrNew, TFLOAT, fpixel, this->numPixels,
117                          this->recon, &status);
118
119      status = 0;
120      fits_close_file(fptrNew, &status);
121      if (status) fits_report_error(stderr, status);
122    }
123  }
124
125}
126
127
128void Cube::saveReconstructedCube()
129{
130  /**
131   *  A function to save the reconstructed and/or residual arrays.
132   *   A number of header keywords are written as well, indicating the
133   *    nature of the reconstruction that has been done.
134   *   The file is always written -- if the filename (as calculated
135   *    based on the recon parameters) exists, then it is overwritten.
136   */
137 
138  int newbitpix = FLOAT_IMG;
139
140  float *resid = new float[this->numPixels];
141  for(int i=0;i<this->numPixels;i++)
142    resid[i] = this->array[i] - this->recon[i];
143  float blankval = this->par.getBlankPixVal();
144
145  long *fpixel = new long[this->numDim];
146  for(int i=0;i<this->numDim;i++) fpixel[i]=1;
147  int status = 0;  /* MUST initialize status */
148  fitsfile *fptrOld, *fptrNew;         
149  fits_open_file(&fptrOld,this->par.getImageFile().c_str(),READONLY,&status);
150  if (status) fits_report_error(stderr, status);
151 
152  if(this->par.getFlagOutputRecon()){
153    std::string fileout = "!" + this->par.outputReconFile();
154    // the ! is there so that it writes over an existing file.
155
156    status = 0;
157    fits_create_file(&fptrNew,fileout.c_str(),&status);
158    if (status)
159      fits_report_error(stderr, status);
160    else
161      {
162        status = 0;
163        fits_movabs_hdu(fptrOld, 1, NULL, &status);
164        if (status) fits_report_error(stderr, status);
165        status = 0;
166        fits_copy_header(fptrOld, fptrNew, &status);
167        if (status) fits_report_error(stderr, status);
168        char *comment = new char[80];
169        long dud;
170        status = 0;
171        fits_update_key(fptrNew, TINT, "BITPIX", &newbitpix,
172                        "number of bits per data pixel", &status);
173        if (status) fits_report_error(stderr, status);
174        status = 0;
175        float bscale=1., bzero=0.;
176        fits_update_key(fptrNew, TFLOAT, "BSCALE", &bscale,
177                        "PHYSICAL = PIXEL*BSCALE + BZERO", &status);
178        fits_update_key(fptrNew, TFLOAT, "BZERO", &bzero, "", &status);
179        fits_set_bscale(fptrNew, 1., 0., &status);
180        if (status) fits_report_error(stderr, status);
181        // Need to correct the dimensions, if we have subsectioned the image
182        if(this->par.getFlagSubsection()){
183          fits_read_key(fptrOld, TLONG, "NAXIS1", &dud, comment, &status);
184          fits_update_key(fptrNew, TLONG, "NAXIS1",
185                          &(this->axisDim[0]), comment, &status);
186          fits_read_key(fptrOld, TLONG, "NAXIS2", &dud, comment, &status);
187          fits_update_key(fptrNew, TLONG, "NAXIS2",
188                          &(this->axisDim[1]), comment, &status);
189          fits_read_key(fptrOld, TLONG, "NAXIS3", &dud, comment, &status);
190          fits_update_key(fptrNew, TLONG, "NAXIS3",
191                          &(this->axisDim[2]), comment, &status);
192        }
193
194        writeReconHeaderInfo(fptrNew, this->par, "recon");
195
196        if(this->par.getFlagBlankPix())
197          fits_write_pixnull(fptrNew, TFLOAT, fpixel, this->numPixels,
198                             this->recon, &blankval, &status);
199        else  fits_write_pix(fptrNew, TFLOAT, fpixel, this->numPixels,
200                             this->recon, &status);
201
202        status = 0;
203        fits_close_file(fptrNew, &status);
204        if (status) fits_report_error(stderr, status);
205      }
206  }
207
208
209  if(this->par.getFlagOutputResid()){
210    std::string fileout = "!" + this->par.outputResidFile();
211    // the ! is there so that it writes over an existing file.
212    status = 0;
213    fits_create_file(&fptrNew,fileout.c_str(),&status);
214    if (status)
215      fits_report_error(stderr, status);
216    else
217      {
218        status = 0;
219        fits_movabs_hdu(fptrOld, 1, NULL, &status);
220        if (status) fits_report_error(stderr, status);
221        status = 0;
222        fits_copy_header(fptrOld, fptrNew, &status);
223        if (status) fits_report_error(stderr, status);
224
225        status = 0;
226        fits_update_key(fptrNew, TINT, "BITPIX", &newbitpix,
227                        "number of bits per data pixel", &status);
228        if (status) fits_report_error(stderr, status);
229        status = 0;
230        float bscale=1., bzero=0.;
231        fits_update_key(fptrNew, TFLOAT, "BSCALE", &bscale,
232                        "PHYSICAL = PIXEL*BSCALE + BZERO", &status);
233        fits_update_key(fptrNew, TFLOAT, "BZERO", &bzero, "", &status);
234        fits_set_bscale(fptrNew, 1., 0., &status);
235        if (status) fits_report_error(stderr, status);
236
237        // Need to correct the dimensions, if we have subsectioned the image...
238        char *comment = new char[80];
239        long dud;
240        if(this->pars().getFlagSubsection()){
241          fits_read_key(fptrOld, TLONG, "NAXIS1", &dud, comment, &status);
242          fits_update_key(fptrNew, TLONG, "NAXIS1",
243                          &(this->axisDim[0]), comment, &status);
244          fits_read_key(fptrOld, TLONG, "NAXIS2", &dud, comment, &status);
245          fits_update_key(fptrNew, TLONG, "NAXIS2",
246                          &(this->axisDim[1]), comment, &status);
247          fits_read_key(fptrOld, TLONG, "NAXIS3", &dud, comment, &status);
248          fits_update_key(fptrNew, TLONG, "NAXIS3",
249                          &(this->axisDim[2]), comment, &status);
250        }
251
252        writeReconHeaderInfo(fptrNew, this->par, "resid");
253
254        if(this->par.getFlagBlankPix())
255          fits_write_pixnull(fptrNew, TFLOAT, fpixel, this->numPixels,
256                             resid, &blankval, &status);
257        else  fits_write_pix(fptrNew, TFLOAT, fpixel, this->numPixels,
258                             resid, &status);
259
260        fits_close_file(fptrNew, &status);
261      }
262  }
263
264  delete [] resid;
265  delete [] fpixel;
266
267}
268
269
270void writeReconHeaderInfo(fitsfile *fptr, Param &par, std::string nature)
271{
272  /**
273   *   A simple function that writes all the necessary keywords and comments
274   *    to the FITS header pointed to by fptr.
275   *   The keyword names and comments are taken from duchamp.hh
276   *   The parameter "nature" indicates what type of file is being written:
277   *    should be either "recon" or "resid".
278   */
279
280
281  int status = 0;
282  std::string explanation = "",ReconResid="";
283
284  fits_write_history(fptr, (char *)header_reconHistory1.c_str(), &status);
285                                   
286  fits_write_history(fptr, (char *)header_reconHistory2.c_str(), &status);
287
288  fits_write_history(fptr, (char *)header_reconHistory_input.c_str(), &status);
289
290  fits_write_history(fptr, (char *)par.getImageFile().c_str(), &status);
291
292  if(par.getFlagSubsection()){
293    fits_write_comment(fptr,(char *)header_reconSubsection_comment.c_str(),
294                       &status);
295    fits_write_key(fptr, TSTRING, (char *)keyword_subsection.c_str(),
296                   (char *)par.getSubsection().c_str(),
297                   (char *)comment_subsection.c_str(), &status);
298  }
299   
300  fits_write_comment(fptr, (char *)header_atrous_comment.c_str(), &status);
301
302  float valf = par.getAtrousCut();
303  fits_write_key(fptr, TFLOAT, (char *)keyword_snrRecon.c_str(), &valf,
304                 (char *)comment_snrRecon.c_str(), &status);
305
306  int vali = par.getReconDim();
307  fits_write_key(fptr, TINT, (char *)keyword_reconDim.c_str(), &vali,
308                 (char *)comment_reconDim.c_str(), &status);
309
310  vali = par.getMinScale();
311  fits_write_key(fptr, TINT, (char *)keyword_scaleMin.c_str(), &vali,
312                 (char *)comment_scaleMin.c_str(), &status);
313
314  vali = par.getFilterCode();
315  fits_write_key(fptr, TINT, (char *)keyword_filterCode.c_str(), &vali,
316                 (char *)comment_filterCode.c_str(), &status);
317
318  if(nature == "recon"){
319    explanation = "Duchamp: This is the RECONSTRUCTED cube";
320    ReconResid = "RECON";
321  }
322  else if(nature == "resid"){
323    explanation = "Duchamp: This is the RESIDUAL cube";
324    ReconResid = "RESID";
325  }
326  else duchampWarning("write_header_info","explanation not present!\n");
327  fits_write_comment(fptr, (char *)explanation.c_str(), &status);
328  fits_write_key(fptr, TSTRING, (char *)keyword_ReconResid.c_str(),
329                 (char *)ReconResid.c_str(),
330                 (char *)comment_ReconResid.c_str(), &status);
331
332}
333
334void writeSmoothHeaderInfo(fitsfile *fptr, Param &par)
335{
336  /**
337   *   A simple function that writes all the necessary keywords and comments
338   *    to the FITS header pointed to by fptr.
339   *   The keyword names and comments are taken from duchamp.hh
340   */
341
342
343  int status = 0;
344
345  fits_write_history(fptr, (char *)header_smoothHistory.c_str(), &status);
346  status = 0;
347  fits_write_history(fptr, (char *)header_smoothHistory_input.c_str(),&status);
348  status = 0;
349  fits_write_history(fptr, (char *)par.getImageFile().c_str(), &status);
350
351  if(par.getFlagSubsection()){
352    status = 0;
353    fits_write_comment(fptr,(char *)header_smoothSubsection_comment.c_str(),
354                       &status);
355    status = 0;
356    fits_write_key(fptr, TSTRING, (char *)keyword_subsection.c_str(),
357                   (char *)par.getSubsection().c_str(),
358                   (char *)comment_subsection.c_str(), &status);
359  }
360   
361  if(par.getSmoothType()=="spatial"){
362    // if kernMin is negative (not defined), make it equal to kernMaj
363    if(par.getKernMin() < 0) par.setKernMin(par.getKernMaj());
364
365    fits_write_key(fptr, TSTRING, (char *)keyword_smoothtype.c_str(),
366                   (char *)header_smoothSpatial.c_str(),
367                   (char *)comment_smoothtype.c_str(), &status);
368    float valf = par.getKernMaj();
369    fits_write_key(fptr, TFLOAT, (char *)keyword_kernmaj.c_str(), &valf,
370                   (char *)comment_kernmaj.c_str(), &status);
371    valf = par.getKernMin();
372    fits_write_key(fptr, TFLOAT, (char *)keyword_kernmin.c_str(), &valf,
373                   (char *)comment_kernmin.c_str(), &status);
374    valf = par.getKernPA();
375    fits_write_key(fptr, TFLOAT, (char *)keyword_kernpa.c_str(), &valf,
376                   (char *)comment_kernpa.c_str(), &status);
377  }
378  else if(par.getSmoothType()=="spectral"){
379    fits_write_key(fptr, TSTRING, (char *)keyword_smoothtype.c_str(),
380                   (char *)header_smoothSpectral.c_str(),
381                   (char *)comment_smoothtype.c_str(), &status);
382    int vali = par.getHanningWidth();
383    fits_write_key(fptr, TINT, (char *)keyword_hanningwidth.c_str(), &vali,
384                   (char *)comment_hanningwidth.c_str(), &status);
385  }
386}
Note: See TracBrowser for help on using the repository browser.