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

Last change on this file since 621 was 578, checked in by MatthewWhiting, 15 years ago

Fixing cfitsio problems in writing reconstructed cubes. Solve ticket #57.
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          this->objectList->at(o).getPixelSet();
155        for(unsigned int p=0;p<voxlist.size();p++){
156          int pixelpos = voxlist[p].getX() + this->axisDim[0]*voxlist[p].getY() +
157            this->axisDim[0]*this->axisDim[1]*voxlist[p].getZ();
158          if(this->par.getFlagMaskWithObjectNum()) mask[pixelpos] = this->objectList->at(o).getID();
159          else mask[pixelpos] = 1;
160        }
161      }
162      status=0;
163      fits_write_pix(fptrNew, TSHORT, fpixel, this->numPixels, mask, &status);
164      if(status){
165        duchampError("saveMask","Fits Error 8:");
166        fits_report_error(stderr,status);
167      }
168      status = 0;
169      fits_close_file(fptrNew, &status);
170      if (status){
171        duchampError("saveMask","Fits Error 9:");
172        fits_report_error(stderr, status);
173      }
174
175      delete [] mask;
176
177    }
178
179    delete [] fpixel;
180
181  }
182 
183  //---------------------------------------------------------------------------
184
185  void Cube::saveSmoothedCube()
186  {
187    /// @brief
188    ///   A function to save the smoothed arrays to a FITS file.
189    ///   Additional header keywords are written as well, indicating the
190    ///   width of the Hanning filter or the dimensions of the Gaussian
191    ///   kernel.
192    ///   The file is always written -- if the filename (as calculated
193    ///    based on the parameters) exists, then it is overwritten.
194 
195    float blankval = this->par.getBlankPixVal();
196
197    int status = 0;  /* MUST initialize status */
198    fitsfile *fptrOld, *fptrNew;         
199    fits_open_file(&fptrOld,this->par.getFullImageFile().c_str(),READONLY,&status);
200    if (status) fits_report_error(stderr, status);
201 
202    if(this->par.getFlagOutputSmooth()){
203      std::string fileout = "!" + this->par.outputSmoothFile();
204      // the ! is there so that it writes over an existing file.
205
206      status = 0;
207      fits_create_file(&fptrNew,fileout.c_str(),&status);
208      if (status)
209        fits_report_error(stderr, status);
210      else {
211        status = 0;
212        fits_movabs_hdu(fptrOld, 1, NULL, &status);
213        if (status) fits_report_error(stderr, status);
214        status = 0;
215        fits_copy_header(fptrOld, fptrNew, &status);
216        if (status) fits_report_error(stderr, status);
217
218        writeSmoothHeaderInfo(fptrNew, this->par);
219
220        if(this->par.getFlagBlankPix())
221          fits_write_imgnull(fptrNew, TFLOAT, 1, this->numPixels, this->recon, &blankval, &status);
222        else
223          fits_write_img(fptrNew, TFLOAT, 1, this->numPixels, this->recon, &status);
224
225        status = 0;
226        fits_close_file(fptrNew, &status);
227        if (status) fits_report_error(stderr, status);
228
229      }
230    }
231
232  }
233
234  //---------------------------------------------------------------------------
235
236  void Cube::saveReconstructedCube()
237  {
238    /// @details
239    ///  A function to save the reconstructed and/or residual arrays.
240    ///   A number of header keywords are written as well, indicating the
241    ///    nature of the reconstruction that has been done.
242    ///   The file is always written -- if the filename (as calculated
243    ///    based on the recon parameters) exists, then it is overwritten.
244 
245    float blankval = this->par.getBlankPixVal();
246
247    int status = 0;  /* MUST initialize status */
248    fitsfile *fptrOld, *fptrNew;         
249    fits_open_file(&fptrOld,this->par.getFullImageFile().c_str(),READONLY,&status);
250    if (status) fits_report_error(stderr, status);
251 
252    if(this->par.getFlagOutputRecon()){
253      std::string fileout = "!" + this->par.outputReconFile();
254      // the ! is there so that it writes over an existing file.
255
256      status = 0;
257      fits_create_file(&fptrNew,fileout.c_str(),&status);
258      if (status)
259        fits_report_error(stderr, status);
260      else
261        {
262          status = 0;
263          fits_movabs_hdu(fptrOld, 1, NULL, &status);
264          if (status) fits_report_error(stderr, status);
265          status = 0;
266          fits_copy_header(fptrOld, fptrNew, &status);
267          if (status) fits_report_error(stderr, status);
268
269          writeReconHeaderInfo(fptrNew, this->par, "recon");
270
271          if(this->par.getFlagBlankPix())
272            fits_write_imgnull(fptrNew, TFLOAT, 1, this->numPixels, this->recon, &blankval, &status);
273          else 
274            fits_write_img(fptrNew, TFLOAT, 1, this->numPixels, this->recon, &status);
275
276          status = 0;
277          fits_close_file(fptrNew, &status);
278          if (status) fits_report_error(stderr, status);
279        }
280    }
281
282
283    if(this->par.getFlagOutputResid()){
284      float *resid = new float[this->numPixels];
285      for(int i=0;i<this->numPixels;i++)
286        resid[i] = this->array[i] - this->recon[i];
287
288      std::string fileout = "!" + this->par.outputResidFile();
289      // the ! is there so that it writes over an existing file.
290      status = 0;
291      fits_create_file(&fptrNew,fileout.c_str(),&status);
292      if (status)
293        fits_report_error(stderr, status);
294      else
295        {
296          status = 0;
297          fits_movabs_hdu(fptrOld, 1, NULL, &status);
298          if (status) fits_report_error(stderr, status);
299          status = 0;
300          fits_copy_header(fptrOld, fptrNew, &status);
301          if (status) fits_report_error(stderr, status);
302
303          writeReconHeaderInfo(fptrNew, this->par, "resid");
304
305          if(this->par.getFlagBlankPix())
306            fits_write_imgnull(fptrNew, TFLOAT, 1, this->numPixels, resid, &blankval, &status);
307          else 
308            fits_write_img(fptrNew, TFLOAT, 1, this->numPixels, resid, &status);
309
310          fits_close_file(fptrNew, &status);
311        }
312      delete [] resid;
313    }
314
315
316  }
317
318  //---------------------------------------------------------------------------
319
320  void writeReconHeaderInfo(fitsfile *fptr, Param &par, std::string nature)
321  {
322    /// @details
323    ///   A simple function that writes all the necessary keywords and comments
324    ///    to the FITS header pointed to by fptr.
325    ///   The keyword names and comments are taken from duchamp.hh
326    ///   The parameter "nature" indicates what type of file is being written:
327    ///    should be either "recon" or "resid".
328
329    int status = 0;
330    std::string explanation = "",ReconResid="";
331
332    fits_write_history(fptr, (char *)header_reconHistory1.c_str(), &status);
333                                   
334    fits_write_history(fptr, (char *)header_reconHistory2.c_str(), &status);
335
336    fits_write_history(fptr, (char *)header_reconHistory_input.c_str(), &status);
337
338    fits_write_history(fptr, (char *)par.getImageFile().c_str(), &status);
339
340    if(par.getFlagSubsection()){
341      fits_write_comment(fptr,(char *)header_reconSubsection_comment.c_str(),
342                         &status);
343      fits_write_key(fptr, TSTRING, (char *)keyword_subsection.c_str(),
344                     (char *)par.getSubsection().c_str(),
345                     (char *)comment_subsection.c_str(), &status);
346    }
347   
348    fits_write_comment(fptr, (char *)header_atrous_comment.c_str(), &status);
349
350    float valf = par.getAtrousCut();
351    fits_write_key(fptr, TFLOAT, (char *)keyword_snrRecon.c_str(), &valf,
352                   (char *)comment_snrRecon.c_str(), &status);
353
354    int vali = par.getReconDim();
355    fits_write_key(fptr, TINT, (char *)keyword_reconDim.c_str(), &vali,
356                   (char *)comment_reconDim.c_str(), &status);
357
358    vali = par.getMinScale();
359    fits_write_key(fptr, TINT, (char *)keyword_scaleMin.c_str(), &vali,
360                   (char *)comment_scaleMin.c_str(), &status);
361
362    vali = par.getFilterCode();
363    fits_write_key(fptr, TINT, (char *)keyword_filterCode.c_str(), &vali,
364                   (char *)comment_filterCode.c_str(), &status);
365
366    if(nature == "recon"){
367      explanation = "Duchamp: This is the RECONSTRUCTED cube";
368      ReconResid = "RECON";
369    }
370    else if(nature == "resid"){
371      explanation = "Duchamp: This is the RESIDUAL cube";
372      ReconResid = "RESID";
373    }
374    else duchampWarning("write_header_info","explanation not present!\n");
375    fits_write_comment(fptr, (char *)explanation.c_str(), &status);
376    fits_write_key(fptr, TSTRING, (char *)keyword_ReconResid.c_str(),
377                   (char *)ReconResid.c_str(),
378                   (char *)comment_ReconResid.c_str(), &status);
379
380  }
381
382  //---------------------------------------------------------------------------
383
384  void writeSmoothHeaderInfo(fitsfile *fptr, Param &par)
385  {
386    /// @details
387    ///   A simple function that writes all the necessary keywords and comments
388    ///    to the FITS header pointed to by fptr.
389    ///   The keyword names and comments are taken from duchamp.hh
390
391    int status = 0;
392
393    fits_write_history(fptr, (char *)header_smoothHistory.c_str(), &status);
394    status = 0;
395    fits_write_history(fptr, (char *)header_smoothHistory_input.c_str(),&status);
396    status = 0;
397    fits_write_history(fptr, (char *)par.getImageFile().c_str(), &status);
398
399    if(par.getFlagSubsection()){
400      status = 0;
401      fits_write_comment(fptr,(char *)header_smoothSubsection_comment.c_str(),
402                         &status);
403      status = 0;
404      fits_write_key(fptr, TSTRING, (char *)keyword_subsection.c_str(),
405                     (char *)par.getSubsection().c_str(),
406                     (char *)comment_subsection.c_str(), &status);
407    }
408   
409    if(par.getSmoothType()=="spatial"){
410      // if kernMin is negative (not defined), make it equal to kernMaj
411      if(par.getKernMin() < 0) par.setKernMin(par.getKernMaj());
412
413      fits_write_key(fptr, TSTRING, (char *)keyword_smoothtype.c_str(),
414                     (char *)header_smoothSpatial.c_str(),
415                     (char *)comment_smoothtype.c_str(), &status);
416      float valf = par.getKernMaj();
417      fits_write_key(fptr, TFLOAT, (char *)keyword_kernmaj.c_str(), &valf,
418                     (char *)comment_kernmaj.c_str(), &status);
419      valf = par.getKernMin();
420      fits_write_key(fptr, TFLOAT, (char *)keyword_kernmin.c_str(), &valf,
421                     (char *)comment_kernmin.c_str(), &status);
422      valf = par.getKernPA();
423      fits_write_key(fptr, TFLOAT, (char *)keyword_kernpa.c_str(), &valf,
424                     (char *)comment_kernpa.c_str(), &status);
425    }
426    else if(par.getSmoothType()=="spectral"){
427      fits_write_key(fptr, TSTRING, (char *)keyword_smoothtype.c_str(),
428                     (char *)header_smoothSpectral.c_str(),
429                     (char *)comment_smoothtype.c_str(), &status);
430      int vali = par.getHanningWidth();
431      fits_write_key(fptr, TINT, (char *)keyword_hanningwidth.c_str(), &vali,
432                     (char *)comment_hanningwidth.c_str(), &status);
433    }
434  }
435
436  //---------------------------------------------------------------------------
437
438  void writeMaskHeaderInfo(fitsfile *fptr, Param &par)
439  {
440    /// @details
441    ///   A simple function that writes all the necessary keywords and comments
442    ///    to the FITS header pointed to by fptr.
443    ///   The keyword names and comments are taken from duchamp.hh
444
445    int status = 0;
446
447    fits_write_history(fptr, (char *)header_maskHistory.c_str(), &status);
448    status = 0;
449    fits_write_history(fptr, (char *)header_maskHistory_input.c_str(),&status);
450    status = 0;
451    fits_write_history(fptr, (char *)par.getImageFile().c_str(), &status);
452
453    if(par.getFlagSubsection()){
454      status = 0;
455      fits_write_comment(fptr,(char *)header_maskSubsection_comment.c_str(),
456                         &status);
457      status = 0;
458      fits_write_key(fptr, TSTRING, (char *)keyword_subsection.c_str(),
459                     (char *)par.getSubsection().c_str(),
460                     (char *)comment_subsection.c_str(), &status);
461    }
462  }
463
464}
Note: See TracBrowser for help on using the repository browser.