Context Navigation

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

Visit:

Last change on this file since 525 was 521, checked in by MatthewWhiting, 16 years ago
Changes for ticket #48, allowing the MASK image to record the object ID of detected pixels. Involves a new parameter flagMaskWithObjectNum.
File size: 20.1 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	/** Write FITS headers in correct format for reconstructed array output */
48	void writeReconHeaderInfo(fitsfile *fptr, Param &par, std::string nature);
49
50	/** Write FITS headers in correct format for smoothed array output */
51	void writeSmoothHeaderInfo(fitsfile *fptr, Param &par);
52
53	/** 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	/**
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
68	int newbitpix = SHORT_IMG;
69	long *fpixel = new long[this->numDim];
70	for(int i=0;i<this->header().WCS().naxis;i++) fpixel[i]=1;
71	int status = 0; /* MUST initialize status */
72	fitsfile fptrOld, fptrNew;
73	fits_open_file(&fptrOld,this->par.getImageFile().c_str(),READONLY,&status);
74	if (status){
75	duchampError("saveMask","Fits Error 1:");
76	fits_report_error(stderr, status);
77	}
78
79	std::string fileout = "!" + this->par.outputMaskFile();
80	// the ! is there so that it writes over an existing file.
81
82	status = 0;
83	fits_create_file(&fptrNew,fileout.c_str(),&status);
84	if (status){
85	duchampError("saveMask","Fits Error 2:");
86	fits_report_error(stderr, status);
87	}
88	else {
89	status = 0;
90	fits_movabs_hdu(fptrOld, 1, NULL, &status);
91	if (status){
92	duchampError("saveMask","Fits Error 3:");
93	fits_report_error(stderr, status);
94	}
95	status = 0;
96	fits_copy_header(fptrOld, fptrNew, &status);
97	if (status){
98	duchampError("saveMask","Fits Error 4:");
99	fits_report_error(stderr, status);
100	}
101	status = 0;
102	fits_update_key(fptrNew, TINT, "BITPIX", &newbitpix,
103	"number of bits per data pixel", &status);
104	if (status){
105	duchampError("saveMask","Fits Error 5:");
106	fits_report_error(stderr, status);
107	}
108	float bscale=1., bzero=0.;
109	fits_update_key(fptrNew, TFLOAT, "BSCALE", &bscale, "", &status);
110	fits_update_key(fptrNew, TFLOAT, "BZERO", &bzero, "", &status);
111	fits_set_bscale(fptrNew, 1, 0, &status);
112	if (status){
113	duchampError("saveMask","Fits Error 6:");
114	fits_report_error(stderr, status);
115	}
116	std::string newunits;
117	if(this->par.getFlagMaskWithObjectNum())
118	newunits = "Object ID";
119	else
120	newunits = "Detection flag";
121	fits_update_key(fptrNew, TSTRING, "BUNIT", (char *)newunits.c_str(), "", &status);
122	if (status){
123	duchampError("saveMask","Fits Error 7:");
124	fits_report_error(stderr, status);
125	}
126	char *comment = new char[80];
127	long dud;
128	// Need to correct the dimensions, if we have subsectioned the image
129	if(this->par.getFlagSubsection()){
130	std::stringstream naxis;
131	fits_read_key(fptrOld, TLONG, (char *)naxis.str().c_str(), &dud, comment, &status);
132	fits_update_key(fptrNew, TLONG, (char *)naxis.str().c_str(),
133	&(this->axisDim[0]), comment, &status);
134	naxis.str("");
135	naxis << "NAXIS" << this->head.WCS().lat;
136	fits_read_key(fptrOld, TLONG, (char *)naxis.str().c_str(), &dud, comment, &status);
137	fits_update_key(fptrNew, TLONG, (char *)naxis.str().c_str(),
138	&(this->axisDim[1]), comment, &status);
139	naxis.str("");
140	naxis << "NAXIS" << this->head.WCS().spec;
141	fits_read_key(fptrOld, TLONG, (char *)naxis.str().c_str(), &dud, comment, &status);
142	fits_update_key(fptrNew, TLONG, (char *)naxis.str().c_str(),
143	&(this->axisDim[2]), comment, &status);
144	}
145
146	delete [] comment;
147
148	writeMaskHeaderInfo(fptrNew, this->par);
149
150	short *mask = new short[this->numPixels];
151	for(int i=0;i<this->numPixels;i++) mask[i]=0;
152	for(int o=0;o<this->objectList->size();o++){
153	std::vector<PixelInfo::Voxel> voxlist =
154	this->objectList->at(o).pixels().getPixelSet();
155	for(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	/**
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
196	int newbitpix = FLOAT_IMG;
197
198	float blankval = this->par.getBlankPixVal();
199
200	long *fpixel = new long[this->numDim];
201	for(int i=0;i<this->numDim;i++) fpixel[i]=1;
202	int status = 0; /* MUST initialize status */
203	fitsfile fptrOld, fptrNew;
204	fits_open_file(&fptrOld,this->par.getImageFile().c_str(),READONLY,&status);
205	if (status) fits_report_error(stderr, status);
206
207	if(this->par.getFlagOutputSmooth()){
208	std::string fileout = "!" + this->par.outputSmoothFile();
209	// the ! is there so that it writes over an existing file.
210
211	status = 0;
212	fits_create_file(&fptrNew,fileout.c_str(),&status);
213	if (status)
214	fits_report_error(stderr, status);
215	else {
216	status = 0;
217	fits_movabs_hdu(fptrOld, 1, NULL, &status);
218	if (status) fits_report_error(stderr, status);
219	status = 0;
220	fits_copy_header(fptrOld, fptrNew, &status);
221	if (status) fits_report_error(stderr, status);
222	char *comment = new char[80];
223	long dud;
224	status = 0;
225	fits_update_key(fptrNew, TINT, "BITPIX", &newbitpix,
226	"number of bits per data pixel", &status);
227	if (status) fits_report_error(stderr, status);
228	status = 0;
229	float bscale=1., bzero=0.;
230	fits_update_key(fptrNew, TFLOAT, "BSCALE", &bscale,
231	"PHYSICAL = PIXEL*BSCALE + BZERO", &status);
232	fits_update_key(fptrNew, TFLOAT, "BZERO", &bzero, "", &status);
233	fits_set_bscale(fptrNew, 1., 0., &status);
234	if (status) fits_report_error(stderr, status);
235	// Need to correct the dimensions, if we have subsectioned the image
236	if(this->par.getFlagSubsection()){
237	fits_read_key(fptrOld, TLONG, "NAXIS1", &dud, comment, &status);
238	fits_update_key(fptrNew, TLONG, "NAXIS1",
239	&(this->axisDim[0]), comment, &status);
240	fits_read_key(fptrOld, TLONG, "NAXIS2", &dud, comment, &status);
241	fits_update_key(fptrNew, TLONG, "NAXIS2",
242	&(this->axisDim[1]), comment, &status);
243	fits_read_key(fptrOld, TLONG, "NAXIS3", &dud, comment, &status);
244	fits_update_key(fptrNew, TLONG, "NAXIS3",
245	&(this->axisDim[2]), comment, &status);
246	}
247
248	delete [] comment;
249
250	writeSmoothHeaderInfo(fptrNew, this->par);
251
252	if(this->par.getFlagBlankPix())
253	fits_write_pixnull(fptrNew, TFLOAT, fpixel, this->numPixels,
254	this->recon, &blankval, &status);
255	else fits_write_pix(fptrNew, TFLOAT, fpixel, this->numPixels,
256	this->recon, &status);
257
258	status = 0;
259	fits_close_file(fptrNew, &status);
260	if (status) fits_report_error(stderr, status);
261
262	}
263	}
264
265	delete [] fpixel;
266
267	}
268
269	//---------------------------------------------------------------------------
270
271	void Cube::saveReconstructedCube()
272	{
273	/**
274	* A function to save the reconstructed and/or residual arrays.
275	* A number of header keywords are written as well, indicating the
276	* nature of the reconstruction that has been done.
277	* The file is always written -- if the filename (as calculated
278	* based on the recon parameters) exists, then it is overwritten.
279	*/
280
281	int newbitpix = FLOAT_IMG;
282
283	float *resid = new float[this->numPixels];
284	for(int i=0;i<this->numPixels;i++)
285	resid[i] = this->array[i] - this->recon[i];
286	float blankval = this->par.getBlankPixVal();
287
288	long *fpixel = new long[this->numDim];
289	for(int i=0;i<this->numDim;i++) fpixel[i]=1;
290	int status = 0; /* MUST initialize status */
291	fitsfile fptrOld, fptrNew;
292	fits_open_file(&fptrOld,this->par.getImageFile().c_str(),READONLY,&status);
293	if (status) fits_report_error(stderr, status);
294
295	if(this->par.getFlagOutputRecon()){
296	std::string fileout = "!" + this->par.outputReconFile();
297	// the ! is there so that it writes over an existing file.
298
299	status = 0;
300	fits_create_file(&fptrNew,fileout.c_str(),&status);
301	if (status)
302	fits_report_error(stderr, status);
303	else
304	{
305	status = 0;
306	fits_movabs_hdu(fptrOld, 1, NULL, &status);
307	if (status) fits_report_error(stderr, status);
308	status = 0;
309	fits_copy_header(fptrOld, fptrNew, &status);
310	if (status) fits_report_error(stderr, status);
311	char *comment = new char[80];
312	long dud;
313	status = 0;
314	fits_update_key(fptrNew, TINT, "BITPIX", &newbitpix,
315	"number of bits per data pixel", &status);
316	if (status) fits_report_error(stderr, status);
317	status = 0;
318	float bscale=1., bzero=0.;
319	fits_update_key(fptrNew, TFLOAT, "BSCALE", &bscale,
320	"PHYSICAL = PIXEL*BSCALE + BZERO", &status);
321	fits_update_key(fptrNew, TFLOAT, "BZERO", &bzero, "", &status);
322	fits_set_bscale(fptrNew, 1., 0., &status);
323	if (status) fits_report_error(stderr, status);
324	// Need to correct the dimensions, if we have subsectioned the image
325	if(this->par.getFlagSubsection()){
326	fits_read_key(fptrOld, TLONG, "NAXIS1", &dud, comment, &status);
327	fits_update_key(fptrNew, TLONG, "NAXIS1",
328	&(this->axisDim[0]), comment, &status);
329	fits_read_key(fptrOld, TLONG, "NAXIS2", &dud, comment, &status);
330	fits_update_key(fptrNew, TLONG, "NAXIS2",
331	&(this->axisDim[1]), comment, &status);
332	fits_read_key(fptrOld, TLONG, "NAXIS3", &dud, comment, &status);
333	fits_update_key(fptrNew, TLONG, "NAXIS3",
334	&(this->axisDim[2]), comment, &status);
335	}
336
337	delete [] comment;
338
339	writeReconHeaderInfo(fptrNew, this->par, "recon");
340
341	if(this->par.getFlagBlankPix())
342	fits_write_pixnull(fptrNew, TFLOAT, fpixel, this->numPixels,
343	this->recon, &blankval, &status);
344	else fits_write_pix(fptrNew, TFLOAT, fpixel, this->numPixels,
345	this->recon, &status);
346
347	status = 0;
348	fits_close_file(fptrNew, &status);
349	if (status) fits_report_error(stderr, status);
350	}
351	}
352
353
354	if(this->par.getFlagOutputResid()){
355	std::string fileout = "!" + this->par.outputResidFile();
356	// the ! is there so that it writes over an existing file.
357	status = 0;
358	fits_create_file(&fptrNew,fileout.c_str(),&status);
359	if (status)
360	fits_report_error(stderr, status);
361	else
362	{
363	status = 0;
364	fits_movabs_hdu(fptrOld, 1, NULL, &status);
365	if (status) fits_report_error(stderr, status);
366	status = 0;
367	fits_copy_header(fptrOld, fptrNew, &status);
368	if (status) fits_report_error(stderr, status);
369
370	status = 0;
371	fits_update_key(fptrNew, TINT, "BITPIX", &newbitpix,
372	"number of bits per data pixel", &status);
373	if (status) fits_report_error(stderr, status);
374	status = 0;
375	float bscale=1., bzero=0.;
376	fits_update_key(fptrNew, TFLOAT, "BSCALE", &bscale,
377	"PHYSICAL = PIXEL*BSCALE + BZERO", &status);
378	fits_update_key(fptrNew, TFLOAT, "BZERO", &bzero, "", &status);
379	fits_set_bscale(fptrNew, 1., 0., &status);
380	if (status) fits_report_error(stderr, status);
381
382	// Need to correct the dimensions, if we have subsectioned the image...
383	char *comment = new char[80];
384	long dud;
385	if(this->pars().getFlagSubsection()){
386	fits_read_key(fptrOld, TLONG, "NAXIS1", &dud, comment, &status);
387	fits_update_key(fptrNew, TLONG, "NAXIS1",
388	&(this->axisDim[0]), comment, &status);
389	fits_read_key(fptrOld, TLONG, "NAXIS2", &dud, comment, &status);
390	fits_update_key(fptrNew, TLONG, "NAXIS2",
391	&(this->axisDim[1]), comment, &status);
392	fits_read_key(fptrOld, TLONG, "NAXIS3", &dud, comment, &status);
393	fits_update_key(fptrNew, TLONG, "NAXIS3",
394	&(this->axisDim[2]), comment, &status);
395	}
396
397	delete [] comment;
398
399	writeReconHeaderInfo(fptrNew, this->par, "resid");
400
401	if(this->par.getFlagBlankPix())
402	fits_write_pixnull(fptrNew, TFLOAT, fpixel, this->numPixels,
403	resid, &blankval, &status);
404	else fits_write_pix(fptrNew, TFLOAT, fpixel, this->numPixels,
405	resid, &status);
406
407	fits_close_file(fptrNew, &status);
408	}
409	}
410
411	delete [] resid;
412	delete [] fpixel;
413
414	}
415
416	//---------------------------------------------------------------------------
417
418	void writeReconHeaderInfo(fitsfile *fptr, Param &par, std::string nature)
419	{
420	/**
421	* A simple function that writes all the necessary keywords and comments
422	* to the FITS header pointed to by fptr.
423	* The keyword names and comments are taken from duchamp.hh
424	* The parameter "nature" indicates what type of file is being written:
425	* should be either "recon" or "resid".
426	*/
427
428
429	int status = 0;
430	std::string explanation = "",ReconResid="";
431
432	fits_write_history(fptr, (char *)header_reconHistory1.c_str(), &status);
433
434	fits_write_history(fptr, (char *)header_reconHistory2.c_str(), &status);
435
436	fits_write_history(fptr, (char *)header_reconHistory_input.c_str(), &status);
437
438	fits_write_history(fptr, (char *)par.getImageFile().c_str(), &status);
439
440	if(par.getFlagSubsection()){
441	fits_write_comment(fptr,(char *)header_reconSubsection_comment.c_str(),
442	&status);
443	fits_write_key(fptr, TSTRING, (char *)keyword_subsection.c_str(),
444	(char *)par.getSubsection().c_str(),
445	(char *)comment_subsection.c_str(), &status);
446	}
447
448	fits_write_comment(fptr, (char *)header_atrous_comment.c_str(), &status);
449
450	float valf = par.getAtrousCut();
451	fits_write_key(fptr, TFLOAT, (char *)keyword_snrRecon.c_str(), &valf,
452	(char *)comment_snrRecon.c_str(), &status);
453
454	int vali = par.getReconDim();
455	fits_write_key(fptr, TINT, (char *)keyword_reconDim.c_str(), &vali,
456	(char *)comment_reconDim.c_str(), &status);
457
458	vali = par.getMinScale();
459	fits_write_key(fptr, TINT, (char *)keyword_scaleMin.c_str(), &vali,
460	(char *)comment_scaleMin.c_str(), &status);
461
462	vali = par.getFilterCode();
463	fits_write_key(fptr, TINT, (char *)keyword_filterCode.c_str(), &vali,
464	(char *)comment_filterCode.c_str(), &status);
465
466	if(nature == "recon"){
467	explanation = "Duchamp: This is the RECONSTRUCTED cube";
468	ReconResid = "RECON";
469	}
470	else if(nature == "resid"){
471	explanation = "Duchamp: This is the RESIDUAL cube";
472	ReconResid = "RESID";
473	}
474	else duchampWarning("write_header_info","explanation not present!\n");
475	fits_write_comment(fptr, (char *)explanation.c_str(), &status);
476	fits_write_key(fptr, TSTRING, (char *)keyword_ReconResid.c_str(),
477	(char *)ReconResid.c_str(),
478	(char *)comment_ReconResid.c_str(), &status);
479
480	}
481
482	//---------------------------------------------------------------------------
483
484	void writeSmoothHeaderInfo(fitsfile *fptr, Param &par)
485	{
486	/**
487	* A simple function that writes all the necessary keywords and comments
488	* to the FITS header pointed to by fptr.
489	* The keyword names and comments are taken from duchamp.hh
490	*/
491
492
493	int status = 0;
494
495	fits_write_history(fptr, (char *)header_smoothHistory.c_str(), &status);
496	status = 0;
497	fits_write_history(fptr, (char *)header_smoothHistory_input.c_str(),&status);
498	status = 0;
499	fits_write_history(fptr, (char *)par.getImageFile().c_str(), &status);
500
501	if(par.getFlagSubsection()){
502	status = 0;
503	fits_write_comment(fptr,(char *)header_smoothSubsection_comment.c_str(),
504	&status);
505	status = 0;
506	fits_write_key(fptr, TSTRING, (char *)keyword_subsection.c_str(),
507	(char *)par.getSubsection().c_str(),
508	(char *)comment_subsection.c_str(), &status);
509	}
510
511	if(par.getSmoothType()=="spatial"){
512	// if kernMin is negative (not defined), make it equal to kernMaj
513	if(par.getKernMin() < 0) par.setKernMin(par.getKernMaj());
514
515	fits_write_key(fptr, TSTRING, (char *)keyword_smoothtype.c_str(),
516	(char *)header_smoothSpatial.c_str(),
517	(char *)comment_smoothtype.c_str(), &status);
518	float valf = par.getKernMaj();
519	fits_write_key(fptr, TFLOAT, (char *)keyword_kernmaj.c_str(), &valf,
520	(char *)comment_kernmaj.c_str(), &status);
521	valf = par.getKernMin();
522	fits_write_key(fptr, TFLOAT, (char *)keyword_kernmin.c_str(), &valf,
523	(char *)comment_kernmin.c_str(), &status);
524	valf = par.getKernPA();
525	fits_write_key(fptr, TFLOAT, (char *)keyword_kernpa.c_str(), &valf,
526	(char *)comment_kernpa.c_str(), &status);
527	}
528	else if(par.getSmoothType()=="spectral"){
529	fits_write_key(fptr, TSTRING, (char *)keyword_smoothtype.c_str(),
530	(char *)header_smoothSpectral.c_str(),
531	(char *)comment_smoothtype.c_str(), &status);
532	int vali = par.getHanningWidth();
533	fits_write_key(fptr, TINT, (char *)keyword_hanningwidth.c_str(), &vali,
534	(char *)comment_hanningwidth.c_str(), &status);
535	}
536	}
537
538	//---------------------------------------------------------------------------
539
540	void writeMaskHeaderInfo(fitsfile *fptr, Param &par)
541	{
542	/**
543	* A simple function that writes all the necessary keywords and comments
544	* to the FITS header pointed to by fptr.
545	* The keyword names and comments are taken from duchamp.hh
546	*/
547
548
549	int status = 0;
550
551	fits_write_history(fptr, (char *)header_maskHistory.c_str(), &status);
552	status = 0;
553	fits_write_history(fptr, (char *)header_maskHistory_input.c_str(),&status);
554	status = 0;
555	fits_write_history(fptr, (char *)par.getImageFile().c_str(), &status);
556
557	if(par.getFlagSubsection()){
558	status = 0;
559	fits_write_comment(fptr,(char *)header_maskSubsection_comment.c_str(),
560	&status);
561	status = 0;
562	fits_write_key(fptr, TSTRING, (char *)keyword_subsection.c_str(),
563	(char *)par.getSubsection().c_str(),
564	(char *)comment_subsection.c_str(), &status);
565	}
566	}
567
568	}

Note: See TracBrowser for help on using the repository browser.

Download in other formats: