[3] | 1 | #include <iostream> |
---|
| 2 | #include <sstream> |
---|
| 3 | #include <string> |
---|
| 4 | #include <wcs.h> |
---|
| 5 | #include <wcshdr.h> |
---|
[146] | 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 | |
---|
[66] | 10 | #include <fitsio.h> |
---|
[71] | 11 | #include <duchamp.hh> |
---|
[3] | 12 | #include <Cubes/cubes.hh> |
---|
| 13 | |
---|
[208] | 14 | void writeReconHeaderInfo(fitsfile *fptr, Param &par, string nature); |
---|
| 15 | void writeSmoothHeaderInfo(fitsfile *fptr, Param &par); |
---|
[86] | 16 | |
---|
[208] | 17 | void Cube::saveSmoothedCube() |
---|
| 18 | { |
---|
| 19 | /** |
---|
| 20 | * A function to save the Hanning-smoothed arrays to a FITS file. |
---|
| 21 | * Additional header keywords are written as well, indicating the |
---|
| 22 | * width of the Hanning filter. |
---|
| 23 | * The file is always written -- if the filename (as calculated |
---|
| 24 | * based on the parameters) exists, then it is overwritten. |
---|
| 25 | */ |
---|
| 26 | |
---|
| 27 | int newbitpix = FLOAT_IMG; |
---|
| 28 | |
---|
| 29 | float blankval = this->par.getBlankPixVal(); |
---|
| 30 | |
---|
| 31 | long *fpixel = new long[this->numDim]; |
---|
| 32 | for(int i=0;i<this->numDim;i++) fpixel[i]=1; |
---|
| 33 | int status = 0; /* MUST initialize status */ |
---|
| 34 | fitsfile *fptrOld, *fptrNew; |
---|
| 35 | fits_open_file(&fptrOld,this->par.getImageFile().c_str(),READONLY,&status); |
---|
| 36 | if (status) fits_report_error(stderr, status); |
---|
| 37 | |
---|
| 38 | if(this->par.getFlagOutputSmooth()){ |
---|
| 39 | string fileout = "!" + this->par.outputSmoothFile(); |
---|
| 40 | // the ! is there so that it writes over an existing file. |
---|
| 41 | |
---|
| 42 | status = 0; |
---|
| 43 | fits_create_file(&fptrNew,fileout.c_str(),&status); |
---|
| 44 | if (status) |
---|
| 45 | fits_report_error(stderr, status); |
---|
| 46 | else { |
---|
| 47 | status = 0; |
---|
| 48 | fits_movabs_hdu(fptrOld, 1, NULL, &status); |
---|
| 49 | if (status) fits_report_error(stderr, status); |
---|
| 50 | status = 0; |
---|
| 51 | fits_copy_header(fptrOld, fptrNew, &status); |
---|
| 52 | if (status) fits_report_error(stderr, status); |
---|
| 53 | char *comment = new char[80]; |
---|
| 54 | long dud; |
---|
| 55 | status = 0; |
---|
| 56 | fits_update_key(fptrNew, TINT, "BITPIX", &newbitpix, |
---|
| 57 | "number of bits per data pixel", &status); |
---|
| 58 | if (status) fits_report_error(stderr, status); |
---|
| 59 | status = 0; |
---|
| 60 | float bscale=1., bzero=0.; |
---|
| 61 | fits_update_key(fptrNew, TFLOAT, "BSCALE", &bscale, |
---|
| 62 | "PHYSICAL = PIXEL*BSCALE + BZERO", &status); |
---|
| 63 | fits_update_key(fptrNew, TFLOAT, "BZERO", &bzero, "", &status); |
---|
| 64 | fits_set_bscale(fptrNew, 1., 0., &status); |
---|
| 65 | if (status) fits_report_error(stderr, status); |
---|
| 66 | // Need to correct the dimensions, if we have subsectioned the image |
---|
| 67 | if(this->par.getFlagSubsection()){ |
---|
| 68 | fits_read_key(fptrOld, TLONG, "NAXIS1", &dud, comment, &status); |
---|
| 69 | fits_update_key(fptrNew, TLONG, "NAXIS1", |
---|
| 70 | &(this->axisDim[0]), comment, &status); |
---|
| 71 | fits_read_key(fptrOld, TLONG, "NAXIS2", &dud, comment, &status); |
---|
| 72 | fits_update_key(fptrNew, TLONG, "NAXIS2", |
---|
| 73 | &(this->axisDim[1]), comment, &status); |
---|
| 74 | fits_read_key(fptrOld, TLONG, "NAXIS3", &dud, comment, &status); |
---|
| 75 | fits_update_key(fptrNew, TLONG, "NAXIS3", |
---|
| 76 | &(this->axisDim[2]), comment, &status); |
---|
| 77 | } |
---|
| 78 | |
---|
| 79 | writeSmoothHeaderInfo(fptrNew, this->par); |
---|
| 80 | |
---|
| 81 | if(this->par.getFlagBlankPix()) |
---|
| 82 | fits_write_pixnull(fptrNew, TFLOAT, fpixel, this->numPixels, |
---|
| 83 | this->recon, &blankval, &status); |
---|
| 84 | else fits_write_pix(fptrNew, TFLOAT, fpixel, this->numPixels, |
---|
| 85 | this->recon, &status); |
---|
| 86 | |
---|
| 87 | status = 0; |
---|
| 88 | fits_close_file(fptrNew, &status); |
---|
| 89 | if (status) fits_report_error(stderr, status); |
---|
| 90 | } |
---|
| 91 | } |
---|
| 92 | |
---|
| 93 | } |
---|
| 94 | |
---|
| 95 | |
---|
[3] | 96 | void Cube::saveReconstructedCube() |
---|
| 97 | { |
---|
[86] | 98 | /** |
---|
[220] | 99 | * A function to save the reconstructed and/or residual arrays. |
---|
[146] | 100 | * A number of header keywords are written as well, indicating the |
---|
| 101 | * nature of the reconstruction that has been done. |
---|
| 102 | * The file is always written -- if the filename (as calculated |
---|
| 103 | * based on the recon parameters) exists, then it is overwritten. |
---|
[86] | 104 | */ |
---|
[3] | 105 | |
---|
[71] | 106 | int newbitpix = FLOAT_IMG; |
---|
| 107 | |
---|
[3] | 108 | float *resid = new float[this->numPixels]; |
---|
[146] | 109 | for(int i=0;i<this->numPixels;i++) |
---|
| 110 | resid[i] = this->array[i] - this->recon[i]; |
---|
[3] | 111 | float blankval = this->par.getBlankPixVal(); |
---|
| 112 | |
---|
[71] | 113 | long *fpixel = new long[this->numDim]; |
---|
| 114 | for(int i=0;i<this->numDim;i++) fpixel[i]=1; |
---|
| 115 | int status = 0; /* MUST initialize status */ |
---|
[3] | 116 | fitsfile *fptrOld, *fptrNew; |
---|
| 117 | fits_open_file(&fptrOld,this->par.getImageFile().c_str(),READONLY,&status); |
---|
| 118 | if (status) fits_report_error(stderr, status); |
---|
| 119 | |
---|
| 120 | if(this->par.getFlagOutputRecon()){ |
---|
[146] | 121 | string fileout = "!" + this->par.outputReconFile(); |
---|
| 122 | // the ! is there so that it writes over an existing file. |
---|
[3] | 123 | |
---|
| 124 | status = 0; |
---|
[71] | 125 | fits_create_file(&fptrNew,fileout.c_str(),&status); |
---|
[103] | 126 | if (status) |
---|
| 127 | fits_report_error(stderr, status); |
---|
| 128 | else |
---|
| 129 | { |
---|
| 130 | status = 0; |
---|
| 131 | fits_movabs_hdu(fptrOld, 1, NULL, &status); |
---|
| 132 | if (status) fits_report_error(stderr, status); |
---|
| 133 | status = 0; |
---|
| 134 | fits_copy_header(fptrOld, fptrNew, &status); |
---|
| 135 | if (status) fits_report_error(stderr, status); |
---|
| 136 | char *comment = new char[80]; |
---|
| 137 | long dud; |
---|
| 138 | status = 0; |
---|
[146] | 139 | fits_update_key(fptrNew, TINT, "BITPIX", &newbitpix, |
---|
| 140 | "number of bits per data pixel", &status); |
---|
[103] | 141 | if (status) fits_report_error(stderr, status); |
---|
| 142 | status = 0; |
---|
| 143 | float bscale=1., bzero=0.; |
---|
[146] | 144 | fits_update_key(fptrNew, TFLOAT, "BSCALE", &bscale, |
---|
| 145 | "PHYSICAL = PIXEL*BSCALE + BZERO", &status); |
---|
[103] | 146 | fits_update_key(fptrNew, TFLOAT, "BZERO", &bzero, "", &status); |
---|
| 147 | fits_set_bscale(fptrNew, 1., 0., &status); |
---|
| 148 | if (status) fits_report_error(stderr, status); |
---|
[146] | 149 | // Need to correct the dimensions, if we have subsectioned the image |
---|
[103] | 150 | if(this->par.getFlagSubsection()){ |
---|
| 151 | fits_read_key(fptrOld, TLONG, "NAXIS1", &dud, comment, &status); |
---|
[146] | 152 | fits_update_key(fptrNew, TLONG, "NAXIS1", |
---|
| 153 | &(this->axisDim[0]), comment, &status); |
---|
[103] | 154 | fits_read_key(fptrOld, TLONG, "NAXIS2", &dud, comment, &status); |
---|
[146] | 155 | fits_update_key(fptrNew, TLONG, "NAXIS2", |
---|
| 156 | &(this->axisDim[1]), comment, &status); |
---|
[103] | 157 | fits_read_key(fptrOld, TLONG, "NAXIS3", &dud, comment, &status); |
---|
[146] | 158 | fits_update_key(fptrNew, TLONG, "NAXIS3", |
---|
| 159 | &(this->axisDim[2]), comment, &status); |
---|
[103] | 160 | } |
---|
[46] | 161 | |
---|
[208] | 162 | writeReconHeaderInfo(fptrNew, this->par, "recon"); |
---|
[3] | 163 | |
---|
[103] | 164 | if(this->par.getFlagBlankPix()) |
---|
[146] | 165 | fits_write_pixnull(fptrNew, TFLOAT, fpixel, this->numPixels, |
---|
| 166 | this->recon, &blankval, &status); |
---|
| 167 | else fits_write_pix(fptrNew, TFLOAT, fpixel, this->numPixels, |
---|
| 168 | this->recon, &status); |
---|
[3] | 169 | |
---|
[103] | 170 | status = 0; |
---|
| 171 | fits_close_file(fptrNew, &status); |
---|
| 172 | if (status) fits_report_error(stderr, status); |
---|
| 173 | } |
---|
[3] | 174 | } |
---|
| 175 | |
---|
[103] | 176 | |
---|
[3] | 177 | if(this->par.getFlagOutputResid()){ |
---|
[146] | 178 | string fileout = "!" + this->par.outputResidFile(); |
---|
| 179 | // the ! is there so that it writes over an existing file. |
---|
[3] | 180 | status = 0; |
---|
[71] | 181 | fits_create_file(&fptrNew,fileout.c_str(),&status); |
---|
[103] | 182 | if (status) |
---|
| 183 | fits_report_error(stderr, status); |
---|
| 184 | else |
---|
| 185 | { |
---|
| 186 | status = 0; |
---|
| 187 | fits_movabs_hdu(fptrOld, 1, NULL, &status); |
---|
| 188 | if (status) fits_report_error(stderr, status); |
---|
| 189 | status = 0; |
---|
| 190 | fits_copy_header(fptrOld, fptrNew, &status); |
---|
| 191 | if (status) fits_report_error(stderr, status); |
---|
[46] | 192 | |
---|
[103] | 193 | status = 0; |
---|
[146] | 194 | fits_update_key(fptrNew, TINT, "BITPIX", &newbitpix, |
---|
| 195 | "number of bits per data pixel", &status); |
---|
[103] | 196 | if (status) fits_report_error(stderr, status); |
---|
| 197 | status = 0; |
---|
| 198 | float bscale=1., bzero=0.; |
---|
[146] | 199 | fits_update_key(fptrNew, TFLOAT, "BSCALE", &bscale, |
---|
| 200 | "PHYSICAL = PIXEL*BSCALE + BZERO", &status); |
---|
[103] | 201 | fits_update_key(fptrNew, TFLOAT, "BZERO", &bzero, "", &status); |
---|
| 202 | fits_set_bscale(fptrNew, 1., 0., &status); |
---|
| 203 | if (status) fits_report_error(stderr, status); |
---|
[71] | 204 | |
---|
[103] | 205 | // Need to correct the dimensions, if we have subsectioned the image... |
---|
| 206 | char *comment = new char[80]; |
---|
| 207 | long dud; |
---|
| 208 | if(this->pars().getFlagSubsection()){ |
---|
| 209 | fits_read_key(fptrOld, TLONG, "NAXIS1", &dud, comment, &status); |
---|
[146] | 210 | fits_update_key(fptrNew, TLONG, "NAXIS1", |
---|
| 211 | &(this->axisDim[0]), comment, &status); |
---|
[103] | 212 | fits_read_key(fptrOld, TLONG, "NAXIS2", &dud, comment, &status); |
---|
[146] | 213 | fits_update_key(fptrNew, TLONG, "NAXIS2", |
---|
| 214 | &(this->axisDim[1]), comment, &status); |
---|
[103] | 215 | fits_read_key(fptrOld, TLONG, "NAXIS3", &dud, comment, &status); |
---|
[146] | 216 | fits_update_key(fptrNew, TLONG, "NAXIS3", |
---|
| 217 | &(this->axisDim[2]), comment, &status); |
---|
[103] | 218 | } |
---|
[3] | 219 | |
---|
[208] | 220 | writeReconHeaderInfo(fptrNew, this->par, "resid"); |
---|
[71] | 221 | |
---|
[103] | 222 | if(this->par.getFlagBlankPix()) |
---|
[146] | 223 | fits_write_pixnull(fptrNew, TFLOAT, fpixel, this->numPixels, |
---|
| 224 | resid, &blankval, &status); |
---|
| 225 | else fits_write_pix(fptrNew, TFLOAT, fpixel, this->numPixels, |
---|
| 226 | resid, &status); |
---|
[3] | 227 | |
---|
[103] | 228 | fits_close_file(fptrNew, &status); |
---|
| 229 | } |
---|
[3] | 230 | } |
---|
| 231 | |
---|
| 232 | delete [] resid; |
---|
| 233 | delete [] fpixel; |
---|
| 234 | |
---|
| 235 | } |
---|
[71] | 236 | |
---|
| 237 | |
---|
[208] | 238 | void writeReconHeaderInfo(fitsfile *fptr, Param &par, string nature) |
---|
[71] | 239 | { |
---|
[86] | 240 | /** |
---|
[208] | 241 | * writeReconHeaderInfo(fptr, par, nature) |
---|
[86] | 242 | * |
---|
| 243 | * A simple function that writes all the necessary keywords and comments |
---|
| 244 | * to the FITS header pointed to by fptr. |
---|
| 245 | * The keyword names and comments are taken from duchamp.hh |
---|
| 246 | * The parameter "nature" indicates what type of file is being written: |
---|
| 247 | * should be either "recon" or "resid". |
---|
| 248 | */ |
---|
| 249 | |
---|
| 250 | |
---|
[71] | 251 | int status = 0; |
---|
| 252 | char *comment, *keyword; |
---|
| 253 | string explanation = "",ReconResid=""; |
---|
| 254 | |
---|
[208] | 255 | fits_write_history(fptr, (char *)header_reconHistory1.c_str(), &status); |
---|
[71] | 256 | |
---|
[208] | 257 | fits_write_history(fptr, (char *)header_reconHistory2.c_str(), &status); |
---|
[71] | 258 | |
---|
[208] | 259 | fits_write_history(fptr, (char *)header_reconHistory_input.c_str(), &status); |
---|
[71] | 260 | |
---|
[105] | 261 | fits_write_history(fptr, (char *)par.getImageFile().c_str(), &status); |
---|
| 262 | |
---|
| 263 | if(par.getFlagSubsection()){ |
---|
[208] | 264 | fits_write_comment(fptr,(char *)header_reconSubsection_comment.c_str(), |
---|
| 265 | &status); |
---|
[105] | 266 | fits_write_key(fptr, TSTRING, (char *)keyword_subsection.c_str(), |
---|
| 267 | (char *)par.getSubsection().c_str(), |
---|
| 268 | (char *)comment_subsection.c_str(), &status); |
---|
| 269 | } |
---|
| 270 | |
---|
| 271 | fits_write_comment(fptr, (char *)header_atrous_comment.c_str(), &status); |
---|
| 272 | |
---|
[71] | 273 | float valf = par.getAtrousCut(); |
---|
| 274 | fits_write_key(fptr, TFLOAT, (char *)keyword_snrRecon.c_str(), &valf, |
---|
| 275 | (char *)comment_snrRecon.c_str(), &status); |
---|
| 276 | |
---|
[103] | 277 | int vali = par.getReconDim(); |
---|
| 278 | fits_write_key(fptr, TINT, (char *)keyword_reconDim.c_str(), &vali, |
---|
| 279 | (char *)comment_reconDim.c_str(), &status); |
---|
| 280 | |
---|
| 281 | vali = par.getMinScale(); |
---|
[71] | 282 | fits_write_key(fptr, TINT, (char *)keyword_scaleMin.c_str(), &vali, |
---|
| 283 | (char *)comment_scaleMin.c_str(), &status); |
---|
| 284 | |
---|
| 285 | vali = par.getFilterCode(); |
---|
| 286 | fits_write_key(fptr, TINT, (char *)keyword_filterCode.c_str(), &vali, |
---|
| 287 | (char *)comment_filterCode.c_str(), &status); |
---|
| 288 | |
---|
| 289 | if(nature == "recon"){ |
---|
| 290 | explanation = "Duchamp: This is the RECONSTRUCTED cube"; |
---|
| 291 | ReconResid = "RECON"; |
---|
| 292 | } |
---|
| 293 | else if(nature == "resid"){ |
---|
| 294 | explanation = "Duchamp: This is the RESIDUAL cube"; |
---|
| 295 | ReconResid = "RESID"; |
---|
| 296 | } |
---|
[120] | 297 | else duchampWarning("write_header_info","explanation not present!\n"); |
---|
[71] | 298 | fits_write_comment(fptr, (char *)explanation.c_str(), &status); |
---|
| 299 | fits_write_key(fptr, TSTRING, (char *)keyword_ReconResid.c_str(), |
---|
[146] | 300 | (char *)ReconResid.c_str(), |
---|
| 301 | (char *)comment_ReconResid.c_str(), &status); |
---|
[71] | 302 | |
---|
| 303 | } |
---|
[208] | 304 | |
---|
| 305 | void writeSmoothHeaderInfo(fitsfile *fptr, Param &par) |
---|
| 306 | { |
---|
| 307 | /** |
---|
| 308 | * writeSmoothHeaderInfo(fptr, par) |
---|
| 309 | * |
---|
| 310 | * A simple function that writes all the necessary keywords and comments |
---|
| 311 | * to the FITS header pointed to by fptr. |
---|
| 312 | * The keyword names and comments are taken from duchamp.hh |
---|
| 313 | */ |
---|
| 314 | |
---|
| 315 | |
---|
| 316 | int status = 0; |
---|
| 317 | char *comment, *keyword; |
---|
| 318 | |
---|
| 319 | fits_write_history(fptr, (char *)header_smoothHistory.c_str(), &status); |
---|
| 320 | status = 0; |
---|
| 321 | fits_write_history(fptr, (char *)header_smoothHistory_input.c_str(),&status); |
---|
| 322 | status = 0; |
---|
| 323 | fits_write_history(fptr, (char *)par.getImageFile().c_str(), &status); |
---|
| 324 | |
---|
| 325 | if(par.getFlagSubsection()){ |
---|
| 326 | status = 0; |
---|
| 327 | fits_write_comment(fptr,(char *)header_smoothSubsection_comment.c_str(), |
---|
| 328 | &status); |
---|
| 329 | status = 0; |
---|
| 330 | fits_write_key(fptr, TSTRING, (char *)keyword_subsection.c_str(), |
---|
| 331 | (char *)par.getSubsection().c_str(), |
---|
| 332 | (char *)comment_subsection.c_str(), &status); |
---|
| 333 | } |
---|
| 334 | |
---|
| 335 | int val = par.getHanningWidth(); |
---|
| 336 | fits_write_key(fptr, TINT, (char *)keyword_hanningwidth.c_str(), &val, |
---|
| 337 | (char *)comment_hanningwidth.c_str(), &status); |
---|
| 338 | |
---|
| 339 | } |
---|