source: tags/release-1.1.8/src/Cubes/saveImage.cc

Last change on this file was 594, checked in by MatthewWhiting, 15 years ago

Saving changes to FITSIO stuff, from [573], [578] and [579]
File size: 16.5 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 <wcslib/wcs.h>
33#include <wcslib/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/duchamp.hh>
40#include <duchamp/Cubes/cubes.hh>
41#include <duchamp/PixelMap/Voxel.hh>
42#include <duchamp/PixelMap/Object3D.hh>
43
44namespace duchamp
45{
46
47  /// @brief Write FITS headers in correct format for reconstructed array output
48  void writeReconHeaderInfo(fitsfile *fptr, Param &par, std::string nature);
49
50  /// @brief Write FITS headers in correct format for smoothed array output
51  void writeSmoothHeaderInfo(fitsfile *fptr, Param &par);
52
53  /// @brief Write FITS headers in correct format for mask array output
54  void writeMaskHeaderInfo(fitsfile *fptr, Param &par);
55  //---------------------------------------------------------------------------
56
57  void Cube::saveMaskCube()
58  {
59    /// @details
60    ///  A function to save a mask to a FITS file, indicating where the
61    ///  detections where made. The value of the detected pixels is
62    ///  determined by the flagMaskWithObjectNum parameter: if true,
63    ///  the value of the pixels is given by the corresponding object
64    ///  ID number; if false, they take the value 1 for all
65    ///  objects. Pixels not in a detected object have the value 0.
66
67    int newbitpix = SHORT_IMG;
68    long *fpixel = new long[this->numDim];
69    for(int i=0;i<this->header().WCS().naxis;i++) fpixel[i]=1;
70    int status = 0;  /* MUST initialize status */
71    fitsfile *fptrOld, *fptrNew;         
72    fits_open_file(&fptrOld,this->par.getImageFile().c_str(),READONLY,&status);
73    if (status){
74      duchampError("saveMask","Fits Error 1:");
75      fits_report_error(stderr, status);
76    }
77
78    std::string fileout = "!" + this->par.outputMaskFile();
79    // the ! is there so that it writes over an existing file.
80
81    status = 0;
82    fits_create_file(&fptrNew,fileout.c_str(),&status);
83    if (status){
84      duchampError("saveMask","Fits Error 2:");
85      fits_report_error(stderr, status);
86    }
87    else {
88      status = 0;
89      fits_movabs_hdu(fptrOld, 1, NULL, &status);
90      if (status){
91        duchampError("saveMask","Fits Error 3:");
92        fits_report_error(stderr, status);
93      }
94      status = 0;
95      fits_copy_header(fptrOld, fptrNew, &status);
96      if (status){
97        duchampError("saveMask","Fits Error 4:");
98        fits_report_error(stderr, status);
99      }
100      status = 0;
101      fits_update_key(fptrNew, TINT, "BITPIX", &newbitpix,
102                      "number of bits per data pixel", &status);
103      if (status){
104        duchampError("saveMask","Fits Error 5:");
105        fits_report_error(stderr, status);
106      }
107      float bscale=1., bzero=0.;
108      fits_update_key(fptrNew, TFLOAT, "BSCALE", &bscale, "", &status);
109      fits_update_key(fptrNew, TFLOAT, "BZERO", &bzero, "", &status);
110      fits_set_bscale(fptrNew, 1, 0, &status);
111      if (status){
112        duchampError("saveMask","Fits Error 6:");
113        fits_report_error(stderr, status);
114      }
115      std::string newunits;
116      if(this->par.getFlagMaskWithObjectNum())
117        newunits = "Object ID";
118      else
119        newunits = "Detection flag";
120      fits_update_key(fptrNew, TSTRING, "BUNIT", (char *)newunits.c_str(), "", &status);
121      if (status){
122        duchampError("saveMask","Fits Error 7:");
123        fits_report_error(stderr, status);
124      }
125      char *comment = new char[80];
126      long dud;
127      // Need to correct the dimensions, if we have subsectioned the image
128      if(this->par.getFlagSubsection()){
129        std::stringstream naxis;
130        fits_read_key(fptrOld, TLONG, (char *)naxis.str().c_str(), &dud, comment, &status);
131        fits_update_key(fptrNew, TLONG, (char *)naxis.str().c_str(),
132                        &(this->axisDim[0]), comment, &status);
133        naxis.str("");
134        naxis << "NAXIS" << this->head.WCS().lat;
135        fits_read_key(fptrOld, TLONG, (char *)naxis.str().c_str(), &dud, comment, &status);
136        fits_update_key(fptrNew, TLONG, (char *)naxis.str().c_str(),
137                        &(this->axisDim[1]), comment, &status);
138        naxis.str("");
139        naxis << "NAXIS" << this->head.WCS().spec;
140        fits_read_key(fptrOld, TLONG, (char *)naxis.str().c_str(), &dud, comment, &status);
141        fits_update_key(fptrNew, TLONG, (char *)naxis.str().c_str(),
142                        &(this->axisDim[2]), comment, &status);
143      }
144
145      delete [] comment;
146
147      writeMaskHeaderInfo(fptrNew, this->par);
148       
149      short *mask = new short[this->numPixels];
150      for(int i=0;i<this->numPixels;i++) mask[i]=0;
151      for(unsigned int o=0;o<this->objectList->size();o++){
152        std::vector<PixelInfo::Voxel> voxlist =
153          this->objectList->at(o).pixels().getPixelSet();
154        for(unsigned int p=0;p<voxlist.size();p++){
155          int pixelpos = voxlist[p].getX() + this->axisDim[0]*voxlist[p].getY() +
156            this->axisDim[0]*this->axisDim[1]*voxlist[p].getZ();
157          if(this->par.getFlagMaskWithObjectNum()) mask[pixelpos] = this->objectList->at(o).getID();
158          else mask[pixelpos] = 1;
159        }
160      }
161      status=0;
162      fits_write_pix(fptrNew, TSHORT, fpixel, this->numPixels, mask, &status);
163      if(status){
164        duchampError("saveMask","Fits Error 8:");
165        fits_report_error(stderr,status);
166      }
167      status = 0;
168      fits_close_file(fptrNew, &status);
169      if (status){
170        duchampError("saveMask","Fits Error 9:");
171        fits_report_error(stderr, status);
172      }
173
174      delete [] mask;
175
176    }
177
178    delete [] fpixel;
179
180  }
181 
182  //---------------------------------------------------------------------------
183
184  void Cube::saveSmoothedCube()
185  {
186    /// @brief
187    ///   A function to save the smoothed arrays to a FITS file.
188    ///   Additional header keywords are written as well, indicating the
189    ///   width of the Hanning filter or the dimensions of the Gaussian
190    ///   kernel.
191    ///   The file is always written -- if the filename (as calculated
192    ///    based on the parameters) exists, then it is overwritten.
193 
194    float blankval = this->par.getBlankPixVal();
195
196    int status = 0;  /* MUST initialize status */
197    fitsfile *fptrOld, *fptrNew;         
198    fits_open_file(&fptrOld,this->par.getFullImageFile().c_str(),READONLY,&status);
199    if (status) fits_report_error(stderr, status);
200 
201    if(this->par.getFlagOutputSmooth()){
202      std::string fileout = "!" + this->par.outputSmoothFile();
203      // the ! is there so that it writes over an existing file.
204
205      status = 0;
206      fits_create_file(&fptrNew,fileout.c_str(),&status);
207      if (status)
208        fits_report_error(stderr, status);
209      else {
210        status = 0;
211        fits_movabs_hdu(fptrOld, 1, NULL, &status);
212        if (status) fits_report_error(stderr, status);
213        status = 0;
214        fits_copy_header(fptrOld, fptrNew, &status);
215        if (status) fits_report_error(stderr, status);
216
217        writeSmoothHeaderInfo(fptrNew, this->par);
218
219        if(this->par.getFlagBlankPix())
220          fits_write_imgnull(fptrNew, TFLOAT, 1, this->numPixels, this->recon, &blankval, &status);
221        else
222          fits_write_img(fptrNew, TFLOAT, 1, this->numPixels, this->recon, &status);
223
224        status = 0;
225        fits_close_file(fptrNew, &status);
226        if (status) fits_report_error(stderr, status);
227
228      }
229    }
230
231  }
232
233  //---------------------------------------------------------------------------
234
235  void Cube::saveReconstructedCube()
236  {
237    /// @details
238    ///  A function to save the reconstructed and/or residual arrays.
239    ///   A number of header keywords are written as well, indicating the
240    ///    nature of the reconstruction that has been done.
241    ///   The file is always written -- if the filename (as calculated
242    ///    based on the recon parameters) exists, then it is overwritten.
243 
244    float blankval = this->par.getBlankPixVal();
245
246    int status = 0;  /* MUST initialize status */
247    fitsfile *fptrOld, *fptrNew;         
248    fits_open_file(&fptrOld,this->par.getFullImageFile().c_str(),READONLY,&status);
249    if (status) fits_report_error(stderr, status);
250 
251    if(this->par.getFlagOutputRecon()){
252      std::string fileout = "!" + this->par.outputReconFile();
253      // the ! is there so that it writes over an existing file.
254
255      status = 0;
256      fits_create_file(&fptrNew,fileout.c_str(),&status);
257      if (status)
258        fits_report_error(stderr, status);
259      else
260        {
261          status = 0;
262          fits_movabs_hdu(fptrOld, 1, NULL, &status);
263          if (status) fits_report_error(stderr, status);
264          status = 0;
265          fits_copy_header(fptrOld, fptrNew, &status);
266          if (status) fits_report_error(stderr, status);
267
268          writeReconHeaderInfo(fptrNew, this->par, "recon");
269
270          if(this->par.getFlagBlankPix())
271            fits_write_imgnull(fptrNew, TFLOAT, 1, this->numPixels, this->recon, &blankval, &status);
272          else 
273            fits_write_img(fptrNew, TFLOAT, 1, this->numPixels, this->recon, &status);
274
275          status = 0;
276          fits_close_file(fptrNew, &status);
277          if (status) fits_report_error(stderr, status);
278        }
279    }
280
281
282    if(this->par.getFlagOutputResid()){
283      float *resid = new float[this->numPixels];
284      for(int i=0;i<this->numPixels;i++)
285        resid[i] = this->array[i] - this->recon[i];
286
287      std::string fileout = "!" + this->par.outputResidFile();
288      // the ! is there so that it writes over an existing file.
289      status = 0;
290      fits_create_file(&fptrNew,fileout.c_str(),&status);
291      if (status)
292        fits_report_error(stderr, status);
293      else
294        {
295          status = 0;
296          fits_movabs_hdu(fptrOld, 1, NULL, &status);
297          if (status) fits_report_error(stderr, status);
298          status = 0;
299          fits_copy_header(fptrOld, fptrNew, &status);
300          if (status) fits_report_error(stderr, status);
301
302          writeReconHeaderInfo(fptrNew, this->par, "resid");
303
304          if(this->par.getFlagBlankPix())
305            fits_write_imgnull(fptrNew, TFLOAT, 1, this->numPixels, resid, &blankval, &status);
306          else 
307            fits_write_img(fptrNew, TFLOAT, 1, this->numPixels, resid, &status);
308
309          fits_close_file(fptrNew, &status);
310        }
311      delete [] resid;
312    }
313
314
315  }
316
317  //---------------------------------------------------------------------------
318
319  void writeReconHeaderInfo(fitsfile *fptr, Param &par, std::string nature)
320  {
321    /// @details
322    ///   A simple function that writes all the necessary keywords and comments
323    ///    to the FITS header pointed to by fptr.
324    ///   The keyword names and comments are taken from duchamp.hh
325    ///   The parameter "nature" indicates what type of file is being written:
326    ///    should be either "recon" or "resid".
327
328    int status = 0;
329    std::string explanation = "",ReconResid="";
330
331    fits_write_history(fptr, (char *)header_reconHistory1.c_str(), &status);
332                                   
333    fits_write_history(fptr, (char *)header_reconHistory2.c_str(), &status);
334
335    fits_write_history(fptr, (char *)header_reconHistory_input.c_str(), &status);
336
337    fits_write_history(fptr, (char *)par.getImageFile().c_str(), &status);
338
339    if(par.getFlagSubsection()){
340      fits_write_comment(fptr,(char *)header_reconSubsection_comment.c_str(),
341                         &status);
342      fits_write_key(fptr, TSTRING, (char *)keyword_subsection.c_str(),
343                     (char *)par.getSubsection().c_str(),
344                     (char *)comment_subsection.c_str(), &status);
345    }
346   
347    fits_write_comment(fptr, (char *)header_atrous_comment.c_str(), &status);
348
349    float valf = par.getAtrousCut();
350    fits_write_key(fptr, TFLOAT, (char *)keyword_snrRecon.c_str(), &valf,
351                   (char *)comment_snrRecon.c_str(), &status);
352
353    int vali = par.getReconDim();
354    fits_write_key(fptr, TINT, (char *)keyword_reconDim.c_str(), &vali,
355                   (char *)comment_reconDim.c_str(), &status);
356
357    vali = par.getMinScale();
358    fits_write_key(fptr, TINT, (char *)keyword_scaleMin.c_str(), &vali,
359                   (char *)comment_scaleMin.c_str(), &status);
360
361    vali = par.getFilterCode();
362    fits_write_key(fptr, TINT, (char *)keyword_filterCode.c_str(), &vali,
363                   (char *)comment_filterCode.c_str(), &status);
364
365    if(nature == "recon"){
366      explanation = "Duchamp: This is the RECONSTRUCTED cube";
367      ReconResid = "RECON";
368    }
369    else if(nature == "resid"){
370      explanation = "Duchamp: This is the RESIDUAL cube";
371      ReconResid = "RESID";
372    }
373    else duchampWarning("write_header_info","explanation not present!\n");
374    fits_write_comment(fptr, (char *)explanation.c_str(), &status);
375    fits_write_key(fptr, TSTRING, (char *)keyword_ReconResid.c_str(),
376                   (char *)ReconResid.c_str(),
377                   (char *)comment_ReconResid.c_str(), &status);
378
379  }
380
381  //---------------------------------------------------------------------------
382
383  void writeSmoothHeaderInfo(fitsfile *fptr, Param &par)
384  {
385    /// @details
386    ///   A simple function that writes all the necessary keywords and comments
387    ///    to the FITS header pointed to by fptr.
388    ///   The keyword names and comments are taken from duchamp.hh
389
390    int status = 0;
391
392    fits_write_history(fptr, (char *)header_smoothHistory.c_str(), &status);
393    status = 0;
394    fits_write_history(fptr, (char *)header_smoothHistory_input.c_str(),&status);
395    status = 0;
396    fits_write_history(fptr, (char *)par.getImageFile().c_str(), &status);
397
398    if(par.getFlagSubsection()){
399      status = 0;
400      fits_write_comment(fptr,(char *)header_smoothSubsection_comment.c_str(),
401                         &status);
402      status = 0;
403      fits_write_key(fptr, TSTRING, (char *)keyword_subsection.c_str(),
404                     (char *)par.getSubsection().c_str(),
405                     (char *)comment_subsection.c_str(), &status);
406    }
407   
408    if(par.getSmoothType()=="spatial"){
409      // if kernMin is negative (not defined), make it equal to kernMaj
410      if(par.getKernMin() < 0) par.setKernMin(par.getKernMaj());
411
412      fits_write_key(fptr, TSTRING, (char *)keyword_smoothtype.c_str(),
413                     (char *)header_smoothSpatial.c_str(),
414                     (char *)comment_smoothtype.c_str(), &status);
415      float valf = par.getKernMaj();
416      fits_write_key(fptr, TFLOAT, (char *)keyword_kernmaj.c_str(), &valf,
417                     (char *)comment_kernmaj.c_str(), &status);
418      valf = par.getKernMin();
419      fits_write_key(fptr, TFLOAT, (char *)keyword_kernmin.c_str(), &valf,
420                     (char *)comment_kernmin.c_str(), &status);
421      valf = par.getKernPA();
422      fits_write_key(fptr, TFLOAT, (char *)keyword_kernpa.c_str(), &valf,
423                     (char *)comment_kernpa.c_str(), &status);
424    }
425    else if(par.getSmoothType()=="spectral"){
426      fits_write_key(fptr, TSTRING, (char *)keyword_smoothtype.c_str(),
427                     (char *)header_smoothSpectral.c_str(),
428                     (char *)comment_smoothtype.c_str(), &status);
429      int vali = par.getHanningWidth();
430      fits_write_key(fptr, TINT, (char *)keyword_hanningwidth.c_str(), &vali,
431                     (char *)comment_hanningwidth.c_str(), &status);
432    }
433  }
434
435  //---------------------------------------------------------------------------
436
437  void writeMaskHeaderInfo(fitsfile *fptr, Param &par)
438  {
439    /// @details
440    ///   A simple function that writes all the necessary keywords and comments
441    ///    to the FITS header pointed to by fptr.
442    ///   The keyword names and comments are taken from duchamp.hh
443
444    int status = 0;
445
446    fits_write_history(fptr, (char *)header_maskHistory.c_str(), &status);
447    status = 0;
448    fits_write_history(fptr, (char *)header_maskHistory_input.c_str(),&status);
449    status = 0;
450    fits_write_history(fptr, (char *)par.getImageFile().c_str(), &status);
451
452    if(par.getFlagSubsection()){
453      status = 0;
454      fits_write_comment(fptr,(char *)header_maskSubsection_comment.c_str(),
455                         &status);
456      status = 0;
457      fits_write_key(fptr, TSTRING, (char *)keyword_subsection.c_str(),
458                     (char *)par.getSubsection().c_str(),
459                     (char *)comment_subsection.c_str(), &status);
460    }
461  }
462
463}
Note: See TracBrowser for help on using the repository browser.