Context Navigation

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

Visit:

Last change on this file since 1455 was 528, checked in by MatthewWhiting, 15 years ago
Changing the documentation comments to match the askapsoft style. Also have split ChanMap? and Object3D into separate files.
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
44	namespace 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(int o=0;o<this->objectList->size();o++){
152	std::vector<PixelInfo::Voxel> voxlist =
153	this->objectList->at(o).pixels().getPixelSet();
154	for(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.

Download in other formats: