Context Navigation

source: tags/release-1.1.2/src/Cubes/saveImage.cc @ 1441

Visit:

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

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

Download in other formats: