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

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

Copying changesets [539]-[543] onto release-1.1.8
File size: 20.2 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    int newbitpix = FLOAT_IMG;
195
196    float blankval = this->par.getBlankPixVal();
197
198    long *fpixel = new long[this->numDim];
199    for(int i=0;i<this->numDim;i++) fpixel[i]=1;
200    int status = 0;  /* MUST initialize status */
201    fitsfile *fptrOld, *fptrNew;         
202    fits_open_file(&fptrOld,this->par.getImageFile().c_str(),READONLY,&status);
203    if (status) fits_report_error(stderr, status);
204 
205    if(this->par.getFlagOutputSmooth()){
206      std::string fileout = "!" + this->par.outputSmoothFile();
207      // the ! is there so that it writes over an existing file.
208
209      status = 0;
210      fits_create_file(&fptrNew,fileout.c_str(),&status);
211      if (status)
212        fits_report_error(stderr, status);
213      else {
214        status = 0;
215        fits_movabs_hdu(fptrOld, 1, NULL, &status);
216        if (status) fits_report_error(stderr, status);
217        status = 0;
218        fits_copy_header(fptrOld, fptrNew, &status);
219        if (status) fits_report_error(stderr, status);
220        char *comment = new char[80];
221        long dud;
222        status = 0;
223        fits_update_key(fptrNew, TINT, "BITPIX", &newbitpix,
224                        "number of bits per data pixel", &status);
225        if (status) fits_report_error(stderr, status);
226        status = 0;
227        float bscale=1., bzero=0.;
228        fits_update_key(fptrNew, TFLOAT, "BSCALE", &bscale,
229                        "PHYSICAL = PIXEL*BSCALE + BZERO", &status);
230        fits_update_key(fptrNew, TFLOAT, "BZERO", &bzero, "", &status);
231        fits_set_bscale(fptrNew, 1., 0., &status);
232        if (status) fits_report_error(stderr, status);
233        // Need to correct the dimensions, if we have subsectioned the image
234        if(this->par.getFlagSubsection()){
235          fits_read_key(fptrOld, TLONG, "NAXIS1", &dud, comment, &status);
236          fits_update_key(fptrNew, TLONG, "NAXIS1",
237                          &(this->axisDim[0]), comment, &status);
238          fits_read_key(fptrOld, TLONG, "NAXIS2", &dud, comment, &status);
239          fits_update_key(fptrNew, TLONG, "NAXIS2",
240                          &(this->axisDim[1]), comment, &status);
241          fits_read_key(fptrOld, TLONG, "NAXIS3", &dud, comment, &status);
242          fits_update_key(fptrNew, TLONG, "NAXIS3",
243                          &(this->axisDim[2]), comment, &status);
244        }
245
246        delete [] comment;
247
248        writeSmoothHeaderInfo(fptrNew, this->par);
249
250        if(this->par.getFlagBlankPix())
251          fits_write_pixnull(fptrNew, TFLOAT, fpixel, this->numPixels,
252                             this->recon, &blankval, &status);
253        else fits_write_pix(fptrNew, TFLOAT, fpixel, this->numPixels,
254                            this->recon, &status);
255
256        status = 0;
257        fits_close_file(fptrNew, &status);
258        if (status) fits_report_error(stderr, status);
259
260      }
261    }
262
263    delete [] fpixel;
264
265  }
266
267  //---------------------------------------------------------------------------
268
269  void Cube::saveReconstructedCube()
270  {
271    /// @details
272    ///  A function to save the reconstructed and/or residual arrays.
273    ///   A number of header keywords are written as well, indicating the
274    ///    nature of the reconstruction that has been done.
275    ///   The file is always written -- if the filename (as calculated
276    ///    based on the recon parameters) exists, then it is overwritten.
277 
278    int newbitpix = FLOAT_IMG;
279
280    float *resid = new float[this->numPixels];
281    for(int i=0;i<this->numPixels;i++)
282      resid[i] = this->array[i] - this->recon[i];
283    float blankval = this->par.getBlankPixVal();
284
285    long *fpixel = new long[this->numDim];
286    for(int i=0;i<this->numDim;i++) fpixel[i]=1;
287    int status = 0;  /* MUST initialize status */
288    fitsfile *fptrOld, *fptrNew;         
289    fits_open_file(&fptrOld,this->par.getImageFile().c_str(),READONLY,&status);
290    if (status) fits_report_error(stderr, status);
291 
292    if(this->par.getFlagOutputRecon()){
293      std::string fileout = "!" + this->par.outputReconFile();
294      // the ! is there so that it writes over an existing file.
295
296      status = 0;
297      fits_create_file(&fptrNew,fileout.c_str(),&status);
298      if (status)
299        fits_report_error(stderr, status);
300      else
301        {
302          status = 0;
303          fits_movabs_hdu(fptrOld, 1, NULL, &status);
304          if (status) fits_report_error(stderr, status);
305          status = 0;
306          fits_copy_header(fptrOld, fptrNew, &status);
307          if (status) fits_report_error(stderr, status);
308          char *comment = new char[80];
309          long dud;
310          status = 0;
311          fits_update_key(fptrNew, TINT, "BITPIX", &newbitpix,
312                          "number of bits per data pixel", &status);
313          if (status) fits_report_error(stderr, status);
314          status = 0;
315          float bscale=1., bzero=0.;
316          fits_update_key(fptrNew, TFLOAT, "BSCALE", &bscale,
317                          "PHYSICAL = PIXEL*BSCALE + BZERO", &status);
318          fits_update_key(fptrNew, TFLOAT, "BZERO", &bzero, "", &status);
319          fits_set_bscale(fptrNew, 1., 0., &status);
320          if (status) fits_report_error(stderr, status);
321          // Need to correct the dimensions, if we have subsectioned the image
322          if(this->par.getFlagSubsection()){
323            fits_read_key(fptrOld, TLONG, "NAXIS1", &dud, comment, &status);
324            fits_update_key(fptrNew, TLONG, "NAXIS1",
325                            &(this->axisDim[0]), comment, &status);
326            fits_read_key(fptrOld, TLONG, "NAXIS2", &dud, comment, &status);
327            fits_update_key(fptrNew, TLONG, "NAXIS2",
328                            &(this->axisDim[1]), comment, &status);
329            fits_read_key(fptrOld, TLONG, "NAXIS3", &dud, comment, &status);
330            fits_update_key(fptrNew, TLONG, "NAXIS3",
331                            &(this->axisDim[2]), comment, &status);
332          }
333
334          delete [] comment;
335
336          writeReconHeaderInfo(fptrNew, this->par, "recon");
337
338          if(this->par.getFlagBlankPix())
339            fits_write_pixnull(fptrNew, TFLOAT, fpixel, this->numPixels,
340                               this->recon, &blankval, &status);
341          else  fits_write_pix(fptrNew, TFLOAT, fpixel, this->numPixels,
342                               this->recon, &status);
343
344          status = 0;
345          fits_close_file(fptrNew, &status);
346          if (status) fits_report_error(stderr, status);
347        }
348    }
349
350
351    if(this->par.getFlagOutputResid()){
352      std::string fileout = "!" + this->par.outputResidFile();
353      // the ! is there so that it writes over an existing file.
354      status = 0;
355      fits_create_file(&fptrNew,fileout.c_str(),&status);
356      if (status)
357        fits_report_error(stderr, status);
358      else
359        {
360          status = 0;
361          fits_movabs_hdu(fptrOld, 1, NULL, &status);
362          if (status) fits_report_error(stderr, status);
363          status = 0;
364          fits_copy_header(fptrOld, fptrNew, &status);
365          if (status) fits_report_error(stderr, status);
366
367          status = 0;
368          fits_update_key(fptrNew, TINT, "BITPIX", &newbitpix,
369                          "number of bits per data pixel", &status);
370          if (status) fits_report_error(stderr, status);
371          status = 0;
372          float bscale=1., bzero=0.;
373          fits_update_key(fptrNew, TFLOAT, "BSCALE", &bscale,
374                          "PHYSICAL = PIXEL*BSCALE + BZERO", &status);
375          fits_update_key(fptrNew, TFLOAT, "BZERO", &bzero, "", &status);
376          fits_set_bscale(fptrNew, 1., 0., &status);
377          if (status) fits_report_error(stderr, status);
378
379          // Need to correct the dimensions, if we have subsectioned the image...
380          char *comment = new char[80];
381          long dud;
382          if(this->pars().getFlagSubsection()){
383            fits_read_key(fptrOld, TLONG, "NAXIS1", &dud, comment, &status);
384            fits_update_key(fptrNew, TLONG, "NAXIS1",
385                            &(this->axisDim[0]), comment, &status);
386            fits_read_key(fptrOld, TLONG, "NAXIS2", &dud, comment, &status);
387            fits_update_key(fptrNew, TLONG, "NAXIS2",
388                            &(this->axisDim[1]), comment, &status);
389            fits_read_key(fptrOld, TLONG, "NAXIS3", &dud, comment, &status);
390            fits_update_key(fptrNew, TLONG, "NAXIS3",
391                            &(this->axisDim[2]), comment, &status);
392          }
393
394          delete [] comment;
395
396          writeReconHeaderInfo(fptrNew, this->par, "resid");
397
398          if(this->par.getFlagBlankPix())
399            fits_write_pixnull(fptrNew, TFLOAT, fpixel, this->numPixels,
400                               resid, &blankval, &status);
401          else  fits_write_pix(fptrNew, TFLOAT, fpixel, this->numPixels,
402                               resid, &status);
403
404          fits_close_file(fptrNew, &status);
405        }
406    }
407
408    delete [] resid;
409    delete [] fpixel;
410
411  }
412
413  //---------------------------------------------------------------------------
414
415  void writeReconHeaderInfo(fitsfile *fptr, Param &par, std::string nature)
416  {
417    /// @details
418    ///   A simple function that writes all the necessary keywords and comments
419    ///    to the FITS header pointed to by fptr.
420    ///   The keyword names and comments are taken from duchamp.hh
421    ///   The parameter "nature" indicates what type of file is being written:
422    ///    should be either "recon" or "resid".
423
424    int status = 0;
425    std::string explanation = "",ReconResid="";
426
427    fits_write_history(fptr, (char *)header_reconHistory1.c_str(), &status);
428                                   
429    fits_write_history(fptr, (char *)header_reconHistory2.c_str(), &status);
430
431    fits_write_history(fptr, (char *)header_reconHistory_input.c_str(), &status);
432
433    fits_write_history(fptr, (char *)par.getImageFile().c_str(), &status);
434
435    if(par.getFlagSubsection()){
436      fits_write_comment(fptr,(char *)header_reconSubsection_comment.c_str(),
437                         &status);
438      fits_write_key(fptr, TSTRING, (char *)keyword_subsection.c_str(),
439                     (char *)par.getSubsection().c_str(),
440                     (char *)comment_subsection.c_str(), &status);
441    }
442   
443    fits_write_comment(fptr, (char *)header_atrous_comment.c_str(), &status);
444
445    float valf = par.getAtrousCut();
446    fits_write_key(fptr, TFLOAT, (char *)keyword_snrRecon.c_str(), &valf,
447                   (char *)comment_snrRecon.c_str(), &status);
448
449    int vali = par.getReconDim();
450    fits_write_key(fptr, TINT, (char *)keyword_reconDim.c_str(), &vali,
451                   (char *)comment_reconDim.c_str(), &status);
452
453    vali = par.getMinScale();
454    fits_write_key(fptr, TINT, (char *)keyword_scaleMin.c_str(), &vali,
455                   (char *)comment_scaleMin.c_str(), &status);
456
457    vali = par.getFilterCode();
458    fits_write_key(fptr, TINT, (char *)keyword_filterCode.c_str(), &vali,
459                   (char *)comment_filterCode.c_str(), &status);
460
461    if(nature == "recon"){
462      explanation = "Duchamp: This is the RECONSTRUCTED cube";
463      ReconResid = "RECON";
464    }
465    else if(nature == "resid"){
466      explanation = "Duchamp: This is the RESIDUAL cube";
467      ReconResid = "RESID";
468    }
469    else duchampWarning("write_header_info","explanation not present!\n");
470    fits_write_comment(fptr, (char *)explanation.c_str(), &status);
471    fits_write_key(fptr, TSTRING, (char *)keyword_ReconResid.c_str(),
472                   (char *)ReconResid.c_str(),
473                   (char *)comment_ReconResid.c_str(), &status);
474
475  }
476
477  //---------------------------------------------------------------------------
478
479  void writeSmoothHeaderInfo(fitsfile *fptr, Param &par)
480  {
481    /// @details
482    ///   A simple function that writes all the necessary keywords and comments
483    ///    to the FITS header pointed to by fptr.
484    ///   The keyword names and comments are taken from duchamp.hh
485
486    int status = 0;
487
488    fits_write_history(fptr, (char *)header_smoothHistory.c_str(), &status);
489    status = 0;
490    fits_write_history(fptr, (char *)header_smoothHistory_input.c_str(),&status);
491    status = 0;
492    fits_write_history(fptr, (char *)par.getImageFile().c_str(), &status);
493
494    if(par.getFlagSubsection()){
495      status = 0;
496      fits_write_comment(fptr,(char *)header_smoothSubsection_comment.c_str(),
497                         &status);
498      status = 0;
499      fits_write_key(fptr, TSTRING, (char *)keyword_subsection.c_str(),
500                     (char *)par.getSubsection().c_str(),
501                     (char *)comment_subsection.c_str(), &status);
502    }
503   
504    if(par.getSmoothType()=="spatial"){
505      // if kernMin is negative (not defined), make it equal to kernMaj
506      if(par.getKernMin() < 0) par.setKernMin(par.getKernMaj());
507
508      fits_write_key(fptr, TSTRING, (char *)keyword_smoothtype.c_str(),
509                     (char *)header_smoothSpatial.c_str(),
510                     (char *)comment_smoothtype.c_str(), &status);
511      float valf = par.getKernMaj();
512      fits_write_key(fptr, TFLOAT, (char *)keyword_kernmaj.c_str(), &valf,
513                     (char *)comment_kernmaj.c_str(), &status);
514      valf = par.getKernMin();
515      fits_write_key(fptr, TFLOAT, (char *)keyword_kernmin.c_str(), &valf,
516                     (char *)comment_kernmin.c_str(), &status);
517      valf = par.getKernPA();
518      fits_write_key(fptr, TFLOAT, (char *)keyword_kernpa.c_str(), &valf,
519                     (char *)comment_kernpa.c_str(), &status);
520    }
521    else if(par.getSmoothType()=="spectral"){
522      fits_write_key(fptr, TSTRING, (char *)keyword_smoothtype.c_str(),
523                     (char *)header_smoothSpectral.c_str(),
524                     (char *)comment_smoothtype.c_str(), &status);
525      int vali = par.getHanningWidth();
526      fits_write_key(fptr, TINT, (char *)keyword_hanningwidth.c_str(), &vali,
527                     (char *)comment_hanningwidth.c_str(), &status);
528    }
529  }
530
531  //---------------------------------------------------------------------------
532
533  void writeMaskHeaderInfo(fitsfile *fptr, Param &par)
534  {
535    /// @details
536    ///   A simple function that writes all the necessary keywords and comments
537    ///    to the FITS header pointed to by fptr.
538    ///   The keyword names and comments are taken from duchamp.hh
539
540    int status = 0;
541
542    fits_write_history(fptr, (char *)header_maskHistory.c_str(), &status);
543    status = 0;
544    fits_write_history(fptr, (char *)header_maskHistory_input.c_str(),&status);
545    status = 0;
546    fits_write_history(fptr, (char *)par.getImageFile().c_str(), &status);
547
548    if(par.getFlagSubsection()){
549      status = 0;
550      fits_write_comment(fptr,(char *)header_maskSubsection_comment.c_str(),
551                         &status);
552      status = 0;
553      fits_write_key(fptr, TSTRING, (char *)keyword_subsection.c_str(),
554                     (char *)par.getSubsection().c_str(),
555                     (char *)comment_subsection.c_str(), &status);
556    }
557  }
558
559}
Note: See TracBrowser for help on using the repository browser.