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

Last change on this file since 378 was 378, checked in by MatthewWhiting, 17 years ago

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