source: tags/release-1.0.5/src/Cubes/saveImage.cc @ 1455

Last change on this file since 1455 was 146, checked in by Matthew Whiting, 18 years ago

A few new things:

  • Made a mycpgplot.hh file, to keep PGPlot-related constants (enum types)

in a standard namespace (so one can do cpgsci(RED) rather than cpgsci(2)...).
Incorporated this into all code that uses pgplot.

  • Improved the outputting of the number of detected objects in the search

functions. Now shows how many have been detected, before they are merged
into the list.

  • Fixed a bug in columns.cc that was incorrectly updating the precision for

negative velocities.

  • Additional text in the Guide, and general clean up of code.
File size: 7.4 KB
Line 
1#include <iostream>
2#include <sstream>
3#include <string>
4#include <wcs.h>
5#include <wcshdr.h>
6#define WCSLIB_GETWCSTAB
7// define this so that we don't try and redefine wtbarr
8// (this is a problem when using cfitsio v.3 and g++ v.4)
9
10#include <fitsio.h>
11#include <duchamp.hh>
12#include <Cubes/cubes.hh>
13
14void write_header_info(fitsfile *fptr, Param &par, string nature);
15
16void Cube::saveReconstructedCube()
17{
18  /**
19   *  Cube::saveReconstructedCube()
20   *
21   *   A function to save the reconstructed and/or residual arrays.
22   *   A number of header keywords are written as well, indicating the
23   *    nature of the reconstruction that has been done.
24   *   The file is always written -- if the filename (as calculated
25   *    based on the recon parameters) exists, then it is overwritten.
26   */
27 
28  int newbitpix = FLOAT_IMG;
29
30  float *resid = new float[this->numPixels];
31  for(int i=0;i<this->numPixels;i++)
32    resid[i] = this->array[i] - this->recon[i];
33  float blankval = this->par.getBlankPixVal();
34
35  long *fpixel = new long[this->numDim];
36  for(int i=0;i<this->numDim;i++) fpixel[i]=1;
37  int status = 0;  /* MUST initialize status */
38  fitsfile *fptrOld, *fptrNew;         
39  fits_open_file(&fptrOld,this->par.getImageFile().c_str(),READONLY,&status);
40  if (status) fits_report_error(stderr, status);
41 
42  if(this->par.getFlagOutputRecon()){
43    string fileout = "!" + this->par.outputReconFile();
44    // the ! is there so that it writes over an existing file.
45
46    status = 0;
47    fits_create_file(&fptrNew,fileout.c_str(),&status);
48    if (status)
49      fits_report_error(stderr, status);
50    else
51      {
52        status = 0;
53        fits_movabs_hdu(fptrOld, 1, NULL, &status);
54        if (status) fits_report_error(stderr, status);
55        status = 0;
56        fits_copy_header(fptrOld, fptrNew, &status);
57        if (status) fits_report_error(stderr, status);
58        char *comment = new char[80];
59        long dud;
60        status = 0;
61        fits_update_key(fptrNew, TINT, "BITPIX", &newbitpix,
62                        "number of bits per data pixel", &status);
63        if (status) fits_report_error(stderr, status);
64        status = 0;
65        float bscale=1., bzero=0.;
66        fits_update_key(fptrNew, TFLOAT, "BSCALE", &bscale,
67                        "PHYSICAL = PIXEL*BSCALE + BZERO", &status);
68        fits_update_key(fptrNew, TFLOAT, "BZERO", &bzero, "", &status);
69        fits_set_bscale(fptrNew, 1., 0., &status);
70        if (status) fits_report_error(stderr, status);
71        // Need to correct the dimensions, if we have subsectioned the image
72        if(this->par.getFlagSubsection()){
73          fits_read_key(fptrOld, TLONG, "NAXIS1", &dud, comment, &status);
74          fits_update_key(fptrNew, TLONG, "NAXIS1",
75                          &(this->axisDim[0]), comment, &status);
76          fits_read_key(fptrOld, TLONG, "NAXIS2", &dud, comment, &status);
77          fits_update_key(fptrNew, TLONG, "NAXIS2",
78                          &(this->axisDim[1]), comment, &status);
79          fits_read_key(fptrOld, TLONG, "NAXIS3", &dud, comment, &status);
80          fits_update_key(fptrNew, TLONG, "NAXIS3",
81                          &(this->axisDim[2]), comment, &status);
82        }
83
84        write_header_info(fptrNew, this->par, "recon");
85
86        if(this->par.getFlagBlankPix())
87          fits_write_pixnull(fptrNew, TFLOAT, fpixel, this->numPixels,
88                             this->recon, &blankval, &status);
89        else  fits_write_pix(fptrNew, TFLOAT, fpixel, this->numPixels,
90                             this->recon, &status);
91
92        status = 0;
93        fits_close_file(fptrNew, &status);
94        if (status) fits_report_error(stderr, status);
95      }
96  }
97
98
99  if(this->par.getFlagOutputResid()){
100    string fileout = "!" + this->par.outputResidFile();
101    // the ! is there so that it writes over an existing file.
102    status = 0;
103    fits_create_file(&fptrNew,fileout.c_str(),&status);
104    if (status)
105      fits_report_error(stderr, status);
106    else
107      {
108        status = 0;
109        fits_movabs_hdu(fptrOld, 1, NULL, &status);
110        if (status) fits_report_error(stderr, status);
111        status = 0;
112        fits_copy_header(fptrOld, fptrNew, &status);
113        if (status) fits_report_error(stderr, status);
114
115        status = 0;
116        fits_update_key(fptrNew, TINT, "BITPIX", &newbitpix,
117                        "number of bits per data pixel", &status);
118        if (status) fits_report_error(stderr, status);
119        status = 0;
120        float bscale=1., bzero=0.;
121        fits_update_key(fptrNew, TFLOAT, "BSCALE", &bscale,
122                        "PHYSICAL = PIXEL*BSCALE + BZERO", &status);
123        fits_update_key(fptrNew, TFLOAT, "BZERO", &bzero, "", &status);
124        fits_set_bscale(fptrNew, 1., 0., &status);
125        if (status) fits_report_error(stderr, status);
126
127        // Need to correct the dimensions, if we have subsectioned the image...
128        char *comment = new char[80];
129        long dud;
130        if(this->pars().getFlagSubsection()){
131          fits_read_key(fptrOld, TLONG, "NAXIS1", &dud, comment, &status);
132          fits_update_key(fptrNew, TLONG, "NAXIS1",
133                          &(this->axisDim[0]), comment, &status);
134          fits_read_key(fptrOld, TLONG, "NAXIS2", &dud, comment, &status);
135          fits_update_key(fptrNew, TLONG, "NAXIS2",
136                          &(this->axisDim[1]), comment, &status);
137          fits_read_key(fptrOld, TLONG, "NAXIS3", &dud, comment, &status);
138          fits_update_key(fptrNew, TLONG, "NAXIS3",
139                          &(this->axisDim[2]), comment, &status);
140        }
141
142        write_header_info(fptrNew, this->par, "resid");
143
144        if(this->par.getFlagBlankPix())
145          fits_write_pixnull(fptrNew, TFLOAT, fpixel, this->numPixels,
146                             resid, &blankval, &status);
147        else  fits_write_pix(fptrNew, TFLOAT, fpixel, this->numPixels,
148                             resid, &status);
149
150        fits_close_file(fptrNew, &status);
151      }
152  }
153
154  delete [] resid;
155  delete [] fpixel;
156
157}
158
159
160void write_header_info(fitsfile *fptr, Param &par, string nature)
161{
162  /**
163   *  write_header_info(fptr, par, nature)
164   *
165   *   A simple function that writes all the necessary keywords and comments
166   *    to the FITS header pointed to by fptr.
167   *   The keyword names and comments are taken from duchamp.hh
168   *   The parameter "nature" indicates what type of file is being written:
169   *    should be either "recon" or "resid".
170   */
171
172
173  int status = 0;
174  char *comment, *keyword;
175  string explanation = "",ReconResid="";
176
177  fits_write_history(fptr, (char *)header_history1.c_str(), &status);
178                                   
179  fits_write_history(fptr, (char *)header_history2.c_str(), &status);
180
181  fits_write_history(fptr, (char *)header_history_input.c_str(), &status);
182
183  fits_write_history(fptr, (char *)par.getImageFile().c_str(), &status);
184
185  if(par.getFlagSubsection()){
186    fits_write_comment(fptr,(char *)header_subsection_comment.c_str(),&status);
187    fits_write_key(fptr, TSTRING, (char *)keyword_subsection.c_str(),
188                   (char *)par.getSubsection().c_str(),
189                   (char *)comment_subsection.c_str(), &status);
190  }
191   
192  fits_write_comment(fptr, (char *)header_atrous_comment.c_str(), &status);
193
194  float valf = par.getAtrousCut();
195  fits_write_key(fptr, TFLOAT, (char *)keyword_snrRecon.c_str(), &valf,
196                 (char *)comment_snrRecon.c_str(), &status);
197
198  int vali = par.getReconDim();
199  fits_write_key(fptr, TINT, (char *)keyword_reconDim.c_str(), &vali,
200                 (char *)comment_reconDim.c_str(), &status);
201
202  vali = par.getMinScale();
203  fits_write_key(fptr, TINT, (char *)keyword_scaleMin.c_str(), &vali,
204                 (char *)comment_scaleMin.c_str(), &status);
205
206  vali = par.getFilterCode();
207  fits_write_key(fptr, TINT, (char *)keyword_filterCode.c_str(), &vali,
208                 (char *)comment_filterCode.c_str(), &status);
209
210  if(nature == "recon"){
211    explanation = "Duchamp: This is the RECONSTRUCTED cube";
212    ReconResid = "RECON";
213  }
214  else if(nature == "resid"){
215    explanation = "Duchamp: This is the RESIDUAL cube";
216    ReconResid = "RESID";
217  }
218  else duchampWarning("write_header_info","explanation not present!\n");
219  fits_write_comment(fptr, (char *)explanation.c_str(), &status);
220  fits_write_key(fptr, TSTRING, (char *)keyword_ReconResid.c_str(),
221                 (char *)ReconResid.c_str(),
222                 (char *)comment_ReconResid.c_str(), &status);
223
224}
Note: See TracBrowser for help on using the repository browser.