source: trunk/src/Cubes/detectionIO.cc @ 431

Last change on this file since 431 was 431, checked in by MatthewWhiting, 16 years ago

Changed the way the borders are calculated, and the way the Annotation file is written. Now shows outline of object in same manner as moment map plots.
File size: 26.9 KB
Line 
1// -----------------------------------------------------------------------
2// detectionIO.cc: Screen and File output of the detected objects.
3// -----------------------------------------------------------------------
4// Copyright (C) 2006, Matthew Whiting, ATNF
5//
6// This program is free software; you can redistribute it and/or modify it
7// under the terms of the GNU General Public License as published by the
8// Free Software Foundation; either version 2 of the License, or (at your
9// option) any later version.
10//
11// Duchamp is distributed in the hope that it will be useful, but WITHOUT
12// ANY WARRANTY; without even the implied warranty of MERCHANTABILITY or
13// FITNESS FOR A PARTICULAR PURPOSE.  See the GNU General Public License
14// for more details.
15//
16// You should have received a copy of the GNU General Public License
17// along with Duchamp; if not, write to the Free Software Foundation,
18// Inc., 59 Temple Place, Suite 330, Boston, MA 02111-1307, USA
19//
20// Correspondence concerning Duchamp may be directed to:
21//    Internet email: Matthew.Whiting [at] atnf.csiro.au
22//    Postal address: Dr. Matthew Whiting
23//                    Australia Telescope National Facility, CSIRO
24//                    PO Box 76
25//                    Epping NSW 1710
26//                    AUSTRALIA
27// -----------------------------------------------------------------------
28#include <iostream>
29#include <sstream>
30#include <fstream>
31#include <iomanip>
32#include <string>
33#include <vector>
34#include <time.h>
35#include <duchamp/param.hh>
36#include <duchamp/fitsHeader.hh>
37#include <duchamp/Cubes/cubes.hh>
38#include <duchamp/PixelMap/Object3D.hh>
39#include <duchamp/Detection/detection.hh>
40#include <duchamp/Detection/columns.hh>
41#include <duchamp/Utils/utils.hh>
42#include <duchamp/Utils/Statistics.hh>
43 
44using std::endl;
45using std::setw;
46using std::setprecision;
47using namespace PixelInfo;
48using namespace Statistics;
49
50namespace duchamp
51{
52
53  using namespace Column;
54
55  void Cube::outputDetectionsKarma(std::ostream &stream)
56  {
57    /**
58     *  Prints to a stream (provided) the list of detected objects in the cube
59     *   in the format of an annotation file for the Karma suite of programs.
60     *  Annotation file draws a box enclosing the detection, and writes the
61     *   ID number of the detection to the right of the box.
62     */
63
64    std::string fname = this->par.getImageFile();
65    if(this->par.getFlagSubsection()) fname+=this->par.getSubsection();
66    stream << "# Duchamp Source Finder results for FITS file: " << fname << endl;
67    if(this->par.getFlagFDR())
68      stream<<"# FDR Significance = " << this->par.getAlpha() << endl;
69    else
70      stream<<"# Threshold = " << this->par.getCut() << endl;
71    if(this->par.getFlagATrous()){
72      stream<<"# The a trous reconstruction method was used, with the following parameters." << endl;
73      stream<<"#  ATrous Dimension = " << this->par.getReconDim() << endl;
74      stream<<"#  ATrous Threshold = " << this->par.getAtrousCut() << endl;
75      stream<<"#  ATrous Minimum Scale =" << this->par.getMinScale() << endl;
76      stream<<"#  ATrous Filter = " << this->par.getFilterName() << endl;
77    }
78    else if(this->par.getFlagSmooth()){
79      stream<<"# The data was smoothed prior to searching, with the following parameters." << endl;
80      stream<<"#  Smoothing type = " << this->par.getSmoothType() << endl;
81      if(this->par.getSmoothType()=="spectral"){
82        stream << "#  Hanning width = " << this->par.getHanningWidth() << endl;
83      }
84      else{
85        stream << "#  Kernel Major axis = " << this->par.getKernMaj() << endl;
86        if(this->par.getKernMin()>0)
87          stream << "#  Kernel Minor axis = " << this->par.getKernMin() << endl;
88        else
89          stream << "#  Kernel Minor axis = " << this->par.getKernMaj() << endl;
90        stream << "#  Kernel Major axis = " << this->par.getKernPA() << endl;
91      }
92    }
93    else stream << "# No ATrous reconstruction done." << endl;
94    stream << "#\n";
95    stream << "COLOR RED" << endl;
96    stream << "COORD W" << endl;
97    stream << std::setprecision(6);
98    stream.setf(std::ios::fixed);
99    double *pix = new double[3];
100    double *wld = new double[3];
101    std::vector<Detection>::iterator obj;
102    for(obj=this->objectList->begin();obj<this->objectList->end();obj++){
103      std::vector<int> vertexSet = obj->getVertexSet();
104      for(int i=0;i<vertexSet.size()/4;i++){
105        pix[0] = vertexSet[i*4]-0.5;
106        pix[1] = vertexSet[i*4+1]-0.5;
107        pix[2] = obj->getZcentre();
108        this->head.pixToWCS(pix,wld);
109        stream << "LINE " << wld[0] << " " << wld[1];
110        pix[0] = vertexSet[i*4+2]-0.5;
111        pix[1] = vertexSet[i*4+3]-0.5;
112        this->head.pixToWCS(pix,wld);
113        stream << " " << wld[0] << " " << wld[1] << "\n";
114      }
115      stream << "TEXT "
116             << obj->getRA() << " "
117             << obj->getDec() << " "
118             << obj->getID() << "\n\n";
119    }
120  }
121
122  std::string fixUnitsVOT(std::string oldstring)
123  {
124    /** Fix a string containing units to make it acceptable for a VOTable.
125     *
126     * This function makes the provided units string acceptable
127     * according to the standard found at
128     * http://vizier.u-strasbg.fr/doc/catstd-3.2.htx
129     * This should then be able to convert the units used in the text
130     * table to units suitable for putting in a VOTable.
131     *
132     * Specifically, it removes any square brackets [] from the
133     * start/end of the string, and replaces blank spaces (representing
134     * multiplication) with a '.' (full stop).
135     */
136
137    std::string newstring;
138    for(int i=0;i<oldstring.size();i++){
139      if((oldstring[i]!='[')&&(oldstring[i]!=']')){
140        if(oldstring[i]==' ') newstring += '.';
141        else newstring += oldstring[i];
142      }
143    }
144    return newstring; 
145  }
146
147  void Cube::outputDetectionsVOTable(std::ostream &stream)
148  {
149    /**
150     *  Prints to a stream (provided) the list of detected objects in the cube
151     *   in a VOTable format.
152     *  Uses WCS information and assumes WCS parameters have been calculated for each
153     *   detected object.
154     *  If they have not (given by the isWCS() function), then those objects are not written...
155     */
156
157    // Set up Column definitions here
158    std::vector<Col> localCol;
159    std::string posUCD[4];
160    localCol.push_back(this->fullCols[NUM]);  // 0 = objID
161    localCol.push_back(this->fullCols[NAME]);  // 1 = name
162    localCol.push_back(this->fullCols[RA]);  // 2 = ra
163    if(makelower(localCol[2].getName())=="ra"){
164      posUCD[0] = "pos.eq.ra;meta.main";
165      posUCD[2] = "phys.angSize;pos.eq.ra";
166    }
167    else{
168      posUCD[0] = "pos.galactic.lat;meta.main";
169      posUCD[2] = "phys.angSize;pos.galactic.lat";
170    }
171    localCol.push_back(this->fullCols[DEC]);  // 3 = dec
172    if(makelower(localCol[2].getName())=="dec"){
173      posUCD[1] = "pos.eq.dec;meta.main";
174      posUCD[3] = "phys.angSize;pos.eq.dec";
175    }
176    else{
177      posUCD[1] = "pos.galactic.lon;meta.main";
178      posUCD[3] = "phys.angSize;pos.galactic.lon";
179    }
180    localCol.push_back(this->fullCols[WRA]);  // 4 = w_ra
181    localCol.push_back(this->fullCols[WDEC]);  // 5 = w_dec
182    localCol.push_back(this->fullCols[VEL]);  // 6 = vel
183    localCol.push_back(this->fullCols[WVEL]); // 7 = w_vel
184    if(this->head.isSpecOK())
185      localCol.push_back(this->fullCols[FINT]); // 8 = f_int
186    else
187      localCol.push_back(this->fullCols[FTOT]); // 8 = f_tot
188    localCol.push_back(this->fullCols[FPEAK]); // 9 = f_peak
189    localCol.push_back(this->fullCols[FLAG]); // 10 = flag
190    localCol.push_back(this->fullCols[XCENT]); // 11 = x_cent
191    localCol.push_back(this->fullCols[YCENT]); // 12 = y_cent
192    localCol.push_back(this->fullCols[ZCENT]); // 13 = z_cent
193    localCol.push_back(this->fullCols[XAV]); // 14 = x_av
194    localCol.push_back(this->fullCols[YAV]); // 15 = y_av
195    localCol.push_back(this->fullCols[ZAV]); // 16 = z_av
196    localCol.push_back(this->fullCols[XPEAK]); // 17 = x_peak
197    localCol.push_back(this->fullCols[YPEAK]); // 18 = y_peak
198    localCol.push_back(this->fullCols[ZPEAK]); // 19 = z_peak
199
200    stream<<"<?xml version=\"1.0\"?>\n";
201    stream<<"<VOTABLE version=\"1.1\" xmlns:xsi=\"http://www.w3.org/2001/XMLSchema-instance\"\n";
202    stream<<" xsi:noNamespaceSchemaLocation=\"http://www.ivoa.net/xml/VOTable/VOTable/v1.1\">\n";
203
204    stream<<"  <COOSYS ID=\"J2000\" equinox=\"J2000.\" epoch=\"J2000.\" system=\"eq_FK5\"/>\n";
205    stream<<"  <RESOURCE name=\"Duchamp Output\">\n";
206    stream<<"    <TABLE name=\"Detections\">\n";
207    stream<<"      <DESCRIPTION>Detected sources and parameters from running the Duchamp source finder.</DESCRIPTION>\n";
208
209    // PARAM section -- parts that are not entry-specific ie. apply to whole dataset
210    std::string fname = this->par.getImageFile();
211    if(this->par.getFlagSubsection()) fname+=this->par.getSubsection();
212    stream<<"      <PARAM name=\"FITS file\" datatype=\"char\" ucd=\"meta.file;meta.fits\" value=\""
213          << fname << "\" arraysize=\""<<fname.size()<<"\"/>\n";
214    if(this->par.getFlagFDR())
215      stream<<"      <PARAM name=\"FDR Significance\" datatype=\"float\" ucd=\"stat.param\" value=\""
216            << this->par.getAlpha() << "\"/>\n";
217    else
218      stream<<"      <PARAM name=\"Threshold\" datatype=\"float\" ucd=\"stat.snr\" value=\""
219            << this->par.getCut() << "\"/>\n";
220    if(this->par.getFlagATrous()){
221      std::string note = "The a trous reconstruction method was used, with the following parameters.";
222      stream<<"      <PARAM name=\"ATrous note\" datatype=\"char\" ucd=\"meta.note\" value=\""
223            <<note << "\" arraysize=\"" << note.size() << "\"/>\n";
224      stream<<"      <PARAM name=\"ATrous Dimension\" datatype=\"int\" ucd=\"meta.code;stat\" value=\""
225            << this->par.getReconDim() << "\"/>\n";
226      stream<<"      <PARAM name=\"ATrous Threshold\" datatype=\"float\" ucd=\"stat.snr\" value=\""
227            << this->par.getAtrousCut() << "\"/>\n";
228      stream<<"      <PARAM name=\"ATrous Minimum Scale\" datatype=\"int\" ucd=\"stat.param\" value=\""
229            << this->par.getMinScale() << "\"/>\n";
230      stream<<"      <PARAM name=\"ATrous Filter\" datatype=\"char\" ucd=\"meta.code;stat\" value=\""
231            << this->par.getFilterName() << "\" arraysize=\"" << this->par.getFilterName().size() << "\"/>\n";
232    }
233    if(this->par.getFlagSmooth()){
234      if(this->par.getSmoothType()=="spectral"){
235        std::string note = "The cube was smoothed spectrally with a Hanning filter, with the following parameters.";
236        stream<<"      <PARAM name=\"Smoothing note\" datatype=\"char\" ucd=\"meta.note\" value=\""
237              <<note << "\" arraysize=\"" << note.size() << "\"/>\n";
238        stream<<"      <PARAM name=\"Hanning filter width\" datatype=\"int\" ucd=\"meta.code;stat\" value=\""
239              << this->par.getHanningWidth() << "\"/>\n";
240      }
241      else if(this->par.getSmoothType()=="spatial"){
242        std::string note = "The cube was smoothed spatially with a Gaussian kernel, with the following parameters.";
243        stream<<"      <PARAM name=\"Smoothing note\" datatype=\"char\" ucd=\"meta.note\" value=\""
244              <<note << "\" arraysize=\"" << note.size() << "\"/>\n";
245        stream<<"      <PARAM name=\"Gaussian kernel major-axis FWHM\" datatype=\"int\" ucd=\"meta.code;stat\" value=\""
246              << this->par.getKernMaj() << "\"/>\n";
247        stream<<"      <PARAM name=\"Gaussian kernel minor-axis FWHM\" datatype=\"int\" ucd=\"meta.code;stat\" value=\""
248              << this->par.getKernMin() << "\"/>\n";
249        stream<<"      <PARAM name=\"Gaussian kernel position angle\" datatype=\"int\" ucd=\"meta.code;stat\" value=\""
250              << this->par.getKernPA() << "\"/>\n";
251      }   
252    }
253    // FIELD section -- names, titles and info for each column.
254    stream<<"      <FIELD name=\"ID\" ID=\"col01\" ucd=\"meta.id\" datatype=\"int\" width=\""
255          <<localCol[0].getWidth()<<"\" unit=\"--\"/>\n";
256    stream<<"      <FIELD name=\"Name\" ID=\"col02\" ucd=\"meta.id;meta.main\" datatype=\"char\" arraysize=\""
257          <<localCol[1].getWidth()<<"\" unit=\"--\"/>\n";
258    stream<<"      <FIELD name=\""<<localCol[2].getName()<<"\" ID=\"col03\" ucd=\""
259          <<posUCD[0]<<"\" ref=\"J2000\" datatype=\"float\" width=\""
260          <<localCol[2].getWidth()<<"\" precision=\""<<prec[prPOS]<<"\" unit=\"deg\"/>\n";
261    stream<<"      <FIELD name=\""<<localCol[3].getName()<<"\" ID=\"col04\" ucd=\""
262          <<posUCD[1]<<"\" ref=\"J2000\" datatype=\"float\" width=\""
263          <<localCol[3].getWidth()<<"\" precision=\""<<prec[prPOS]<<"\" unit=\"deg\"/>\n";
264    stream<<"      <FIELD name=\""<<localCol[4].getName()<<"\" ID=\"col05\" ucd=\""
265          <<posUCD[2]<<"\" ref=\"J2000\" datatype=\"float\" width=\""
266          <<localCol[4].getWidth()<<"\" precision=\""<<prec[prWPOS]<<"\" unit=\""
267          <<fixUnitsVOT(localCol[4].getUnits())<<"\"/>\n";
268    stream<<"      <FIELD name=\""<<localCol[5].getName()<<"\" ID=\"col06\" ucd=\""
269          <<posUCD[2]<<"\" ref=\"J2000\" datatype=\"float\" width=\""
270          <<localCol[5].getWidth()<<"\" precision=\""<<prec[prWPOS]<<"\" unit=\""
271          <<fixUnitsVOT(localCol[5].getUnits())<<"\"/>\n";
272    stream<<"      <FIELD name=\"Vel\" ID=\"col07\" ucd=\"phys.veloc;src.dopplerVeloc\" datatype=\"float\" width=\""
273          <<localCol[6].getWidth()<<"\" precision=\""<<prec[prVEL]<<"\" unit=\""
274          <<fixUnitsVOT(localCol[6].getUnits())<<"\"/>\n";
275    stream<<"      <FIELD name=\"w_Vel\" ID=\"col08\" ucd=\"phys.veloc;src.dopplerVeloc;spect.line.width\""
276          <<" datatype=\"float\" width=\""<<localCol[7].getWidth()<<"\" precision=\""
277          <<prec[prVEL]<<"\" unit=\""<<fixUnitsVOT(localCol[7].getUnits())<<"\"/>\n";
278    if(this->head.isSpecOK())
279      stream<<"      <FIELD name=\"Integrated_Flux\" ID=\"col09\" ucd=\"phot.flux;spect.line.intensity\""
280            <<" datatype=\"float\" width=\""<<localCol[8].getWidth()<<"\" precision=\""
281            <<prec[prFLUX]<<"\" unit=\""<<fixUnitsVOT(localCol[8].getUnits())<<"\"/>\n";
282    else
283      stream<<"      <FIELD name=\"Total_Flux\" ID=\"col09\" ucd=\"phot.flux;spect.line.intensity\""
284            <<" datatype=\"float\" width=\""<<localCol[8].getWidth()<<"\" precision=\""
285            <<prec[prFLUX]<<"\" unit=\""<<fixUnitsVOT(localCol[8].getUnits())<<"\"/>\n";
286    stream<<"      <FIELD name=\"Peak_Flux\" ID=\"col10\" ucd=\"phot.flux;spect.line.intensity\""
287          <<" datatype=\"float\" width=\""<<localCol[9].getWidth()<<"\" precision=\""
288          <<prec[prFLUX]<<"\" unit=\""<<fixUnitsVOT(localCol[9].getUnits())<<"\"/>\n";
289    stream<<"      <FIELD name=\"Flag\" ID=\"col11\" ucd=\"meta.code.qual\" datatype=\"char\" arraysize=\"3\" unit=\"--\"/>\n";
290    stream<<"      <FIELD name=\"X_Centroid\" ID=\"col12\" ucd=\"pos.cartesian.x\""
291          <<" datatype=\"float\" width=\""<<localCol[11].getWidth()<<"\" precision=\""<<localCol[11].getPrecision()<<"\" unit=\"\"/>\n";
292    stream<<"      <FIELD name=\"Y_Centroid\" ID=\"col13\" ucd=\"pos.cartesian.y\""
293          <<" datatype=\"float\" width=\""<<localCol[12].getWidth()<<"\" precision=\""<<localCol[12].getPrecision()<<"\" unit=\"\"/>\n";
294    stream<<"      <FIELD name=\"Z_Centroid\" ID=\"col14\" ucd=\"pos.cartesian.z\""
295          <<" datatype=\"float\" width=\""<<localCol[13].getWidth()<<"\" precision=\""<<localCol[13].getPrecision()<<"\" unit=\"\"/>\n";
296    stream<<"      <FIELD name=\"X_Av\" ID=\"col15\" ucd=\"pos.cartesian.x\""
297          <<" datatype=\"float\" width=\""<<localCol[14].getWidth()<<"\" precision=\""<<localCol[14].getPrecision()<<"\" unit=\"\"/>\n";
298    stream<<"      <FIELD name=\"Y_Av\" ID=\"col16\" ucd=\"pos.cartesian.y\""
299          <<" datatype=\"float\" width=\""<<localCol[15].getWidth()<<"\" precision=\""<<localCol[15].getPrecision()<<"\" unit=\"\"/>\n";
300    stream<<"      <FIELD name=\"Z_Av\" ID=\"col17\" ucd=\"pos.cartesian.z\""
301          <<" datatype=\"float\" width=\""<<localCol[16].getWidth()<<"\" precision=\""<<localCol[16].getPrecision()<<"\" unit=\"\"/>\n";
302    stream<<"      <FIELD name=\"X_Peak\" ID=\"col18\" ucd=\"pos.cartesian.x\""
303          <<" datatype=\"int\" width=\""<<localCol[17].getWidth()<<"\" precision=\""<<localCol[17].getPrecision()<<"\" unit=\"\"/>\n";
304    stream<<"      <FIELD name=\"Y_Peak\" ID=\"col19\" ucd=\"pos.cartesian.y\""
305          <<" datatype=\"int\" width=\""<<localCol[18].getWidth()<<"\" precision=\""<<localCol[18].getPrecision()<<"\" unit=\"\"/>\n";
306    stream<<"      <FIELD name=\"Z_Peak\" ID=\"col20\" ucd=\"pos.cartesian.z\""
307          <<" datatype=\"int\" width=\""<<localCol[19].getWidth()<<"\" precision=\""<<localCol[19].getPrecision()<<"\" unit=\"\"/>\n";
308
309
310    stream<<"      <DATA>\n"
311          <<"        <TABLEDATA>\n";
312
313    stream.setf(std::ios::fixed); 
314    for(int i=0;i<this->objectList->size();i++){
315      if(this->objectList->at(i).isWCS()){
316        stream<<"        <TR>\n";
317        stream<<"          <TD>"<<setw(localCol[0].getWidth())<<i+1<<"</TD>";
318        stream<<"<TD>" << setw(localCol[1].getWidth()) << this->objectList->at(i).getName()       <<"</TD>";
319        stream<<setprecision(prec[prPOS]);
320        stream<<"<TD>" << setw(localCol[2].getWidth()) << this->objectList->at(i).getRA()         <<"</TD>";
321        stream<<"<TD>" << setw(localCol[3].getWidth()) << this->objectList->at(i).getDec()        <<"</TD>";
322        stream<<setprecision(prec[prWPOS]);
323        stream<<"<TD>" << setw(localCol[4].getWidth()) << this->objectList->at(i).getRAWidth()    <<"</TD>";
324        stream<<"<TD>" << setw(localCol[5].getWidth()) << this->objectList->at(i).getDecWidth()   <<"</TD>";
325        stream<<setprecision(prec[prVEL]);
326        stream<<"<TD>" << setw(localCol[6].getWidth()) << this->objectList->at(i).getVel()        <<"</TD>";
327        stream<<"<TD>" << setw(localCol[7].getWidth()) << this->objectList->at(i).getVelWidth()   <<"</TD>";
328        stream<<setprecision(prec[prFLUX]);
329        if(this->head.isSpecOK())
330          stream<<"<TD>" << setw(localCol[8].getWidth()) << this->objectList->at(i).getIntegFlux() <<"</TD>";
331        else
332          stream<<"<TD>" << setw(localCol[8].getWidth()) << this->objectList->at(i).getTotalFlux() <<"</TD>";
333        stream<<"<TD>" << setw(localCol[9].getWidth()) << this->objectList->at(i).getPeakFlux()   <<"</TD>";
334        stream<<"<TD>" << setw(localCol[10].getWidth())<< this->objectList->at(i).getFlagText()   <<"</TD>";
335        stream<<"<TD>" << setw(localCol[11].getWidth())<< this->objectList->at(i).getXCentroid()  <<"</TD>";
336        stream<<"<TD>" << setw(localCol[12].getWidth())<< this->objectList->at(i).getYCentroid()  <<"</TD>";
337        stream<<"<TD>" << setw(localCol[13].getWidth())<< this->objectList->at(i).getZCentroid()  <<"</TD>";
338        stream<<"<TD>" << setw(localCol[14].getWidth())<< this->objectList->at(i).getXAverage()   <<"</TD>";
339        stream<<"<TD>" << setw(localCol[15].getWidth())<< this->objectList->at(i).getYAverage()   <<"</TD>";
340        stream<<"<TD>" << setw(localCol[16].getWidth())<< this->objectList->at(i).getZAverage()   <<"</TD>";
341        stream<<"<TD>" << setw(localCol[17].getWidth())<< this->objectList->at(i).getXPeak()      <<"</TD>";
342        stream<<"<TD>" << setw(localCol[18].getWidth())<< this->objectList->at(i).getYPeak()      <<"</TD>";
343        stream<<"<TD>" << setw(localCol[19].getWidth())<< this->objectList->at(i).getZPeak()      <<"</TD>";
344     
345        stream<<endl;
346        stream<<"        </TR>\n";
347      }
348    }
349    stream<<"        </TABLEDATA>\n";
350    stream<<"      </DATA>\n";
351    stream<<"    </TABLE>\n";
352    stream<<"  </RESOURCE>\n";
353    stream<<"</VOTABLE>\n";
354    resetiosflags(std::ios::fixed);
355 
356  }
357
358  void Cube::prepareOutputFile()
359  {
360    /**
361     *  A function to write the paramters, time of execution, and
362     *  statistical information to the output file.
363     */
364
365    std::string outfile;
366    if(this->par.getFlagSeparateHeader()) outfile = this->par.getHeaderFile();
367    else outfile = this->par.getOutFile();
368    std::ofstream output(outfile.c_str());
369    output<<"Results of the Duchamp source finder: ";
370    time_t now = time(NULL);
371    output << asctime( localtime(&now) );
372    this->showParam(output);
373    output<<"--------------------\n";
374    output.close();
375    this->outputStats();
376 
377  }
378
379  void Cube::outputStats()
380  {
381    /**
382     *  A function to write the statistical information to the output
383     *  file. This writes the threshold, its equivalent S/N ratio, and
384     *  the noise level and spread.
385     *
386     *  If smoothing has been done, the noise level & spread for the
387     *  original array are calculated and printed as well.
388     */
389
390    std::string outfile;
391    if(this->par.getFlagSeparateHeader()) outfile = this->par.getHeaderFile();
392    else outfile = this->par.getOutFile();
393    std::ofstream output(outfile.c_str(),std::ios::app);
394    output<<"Summary of statistics:\n";
395    output<<"Detection threshold = " << this->Stats.getThreshold()
396          <<" " << this->head.getFluxUnits();
397    if(this->par.getFlagFDR())
398      output<<" (or S/N=" << this->Stats.getThresholdSNR()<<")";
399    if(this->par.getFlagSmooth()){
400      output << " in smoothed cube.";
401      if(!this->par.getFlagUserThreshold())
402        output<<"\nNoise level = " << this->Stats.getMiddle()
403              <<", Noise spread = " << this->Stats.getSpread()
404              <<" in smoothed cube.";
405   
406      // calculate the stats for the original array, so that we can
407      // quote S/N values correctly.
408      this->par.setFlagSmooth(false);
409      bool verb=this->par.isVerbose();
410      bool fdrflag=this->par.getFlagFDR();
411      this->par.setVerbosity(false);
412      this->par.setFlagFDR(false);
413      this->setCubeStats();
414      this->par.setVerbosity(verb);
415      this->par.setFlagFDR(fdrflag);
416      this->par.setFlagSmooth(true);
417     
418      output << "\nNoise properties for the original cube are:";
419    }
420     
421    if(!this->par.getFlagUserThreshold())
422      output<<"\nNoise level = " << this->Stats.getMiddle()
423            <<", Noise spread = " << this->Stats.getSpread()
424            <<"\n";
425
426    if(this->par.getFlagGrowth()){
427      StatsContainer<float> growthStats = this->Stats;
428      growthStats.setThresholdSNR(this->par.getGrowthCut());
429      growthStats.setUseFDR(false);
430      output<<"  Detections grown down to threshold of "
431            << growthStats.getThreshold() << ".\n";
432    }
433
434    if(!this->par.getFlagUserThreshold())
435      output << "\nFull stats:\n" << this->Stats;
436    else
437      output << "\n\nNot calculating full stats since threshold given.\n";
438
439    output<<"--------------------\n";
440    output.close();
441  }
442
443  void Cube::outputDetectionList()
444  {
445    /**
446     *  A front-end to writing the full list of detected objects to a results
447     *   file and to cout.
448     *  Uses outputDetectionTextWCS for each objects.
449     *  Leaves the testing of whether the WCS parameters for each object
450     *   have been calculated to outputDetectionTextWCS.
451     */
452
453    std::string outfile;
454    if(this->par.getFlagSeparateHeader()) outfile = this->par.getHeaderFile();
455    else outfile = this->par.getOutFile();
456    std::ofstream output(outfile.c_str(),std::ios::app);
457    output<<"Total number of detections = "<<this->objectList->size()<<endl;
458    output<<"--------------------\n";
459    output.close();
460
461    if(this->par.getFlagSeparateHeader())
462      output.open(this->par.getOutFile().c_str());
463    else
464      output.open(this->par.getOutFile().c_str(),std::ios::app);
465
466    if(this->objectList->size()>0){
467      this->setupColumns();
468      this->objectList->at(0).outputDetectionTextHeaderFull(output,this->fullCols);
469      this->objectList->at(0).outputDetectionTextHeader(std::cout,this->fullCols);
470      for(int i=0;i<this->objectList->size();i++){
471        this->objectList->at(i).outputDetectionTextWCSFull(output,this->fullCols);
472        this->objectList->at(i).outputDetectionTextWCS(std::cout,this->fullCols);
473      }
474    }
475
476    output.close();
477  }
478
479  void Cube::prepareLogFile(int argc, char *argv[])
480  {
481    /**
482     *  Opens the log file so that it can be written to, and writes
483     *  the parameter summary as well as the time of execution to the
484     *  file.
485     *
486     *  It also writes the command-line statement, hence the need for
487     *  argv and argc.
488     */
489    // Open the logfile and write the time on the first line
490    std::ofstream logfile(this->par.getLogFile().c_str());
491    logfile << "New run of the Duchamp sourcefinder: ";
492    time_t now = time(NULL);
493    logfile << asctime( localtime(&now) );
494    // Write out the command-line statement
495    logfile << "Executing statement : ";
496    for(int i=0;i<argc;i++) logfile << argv[i] << " ";
497    logfile << std::endl;
498    logfile << this->par;
499    logfile.close();
500
501  }
502
503
504  void Cube::logDetectionList()
505  {
506    /**
507     *  A front-end to writing a list of detected objects to the log file.
508     *  Does not assume WCS, so uses outputDetectionText.
509     *  Designed to be used by searching routines before returning their
510     *   final list.
511     */
512
513    if(this->objectList->size()>0){
514
515      long left = this->par.getBorderLeft();
516      long bottom = this->par.getBorderBottom();
517
518      std::ofstream fout(this->par.getLogFile().c_str(),std::ios::app);
519      this->calcObjectFluxes();
520      this->setupColumns();
521      this->objectList->at(0).outputDetectionTextHeader(fout,this->logCols);
522
523      if(this->par.getFlagBaseline()){
524        for(int i=0;i<this->axisDim[0]*this->axisDim[1]*this->axisDim[2];i++)
525          this->array[i] += this->baseline[i];
526      }
527
528      for(int objCtr=0;objCtr<this->objectList->size();objCtr++){
529        Detection *obj = new Detection;
530        *obj = objectList->at(objCtr);
531        obj->setOffsets(par);
532        if(this->par.getFlagCubeTrimmed()){
533          obj->pixels().addOffsets(left,bottom,0);
534        }
535        obj->calcFluxes(this->array, this->axisDim);
536        obj->outputDetectionText(fout,this->logCols,objCtr+1);
537        delete obj;
538      }
539
540      if(this->par.getFlagBaseline()){
541        for(int i=0;i<this->axisDim[0]*this->axisDim[1]*this->axisDim[2];i++)
542          this->array[i] -= this->baseline[i];
543      }
544      fout.close();
545 
546    }
547
548  }
549
550  void Cube::logDetection(Detection obj, int counter)
551  {
552    /**
553     *  A front-end to writing a detected object to the log file.
554     *  Does not assume WCS, so uses outputDetectionText.
555     *  Corrects for changes to positions of pixels and removal of baselines.
556     *  Designed to be used by searching routines before returning their final
557     *   list.
558     *  \param obj Detection object to be written : passed by value, as we want
559     *    to potentially change positions etc, but not for the object in the
560     *    calling function.
561     *  \param counter The number to assign to the object : ideally its number
562     *    in a list of some kind.
563     */
564
565    std::ofstream fout(this->par.getLogFile().c_str(),std::ios::app);
566    // Need to deal with possibility of trimmed array
567    long left = this->par.getBorderLeft();
568    long right = this->par.getBorderRight();
569    long top = this->par.getBorderTop();
570    long bottom = this->par.getBorderBottom();
571    long *tempDim = new long[3];
572    tempDim[0] = (this->axisDim[0] + left + right);
573    tempDim[1] = (this->axisDim[1] + bottom + top);
574    tempDim[2] = this->axisDim[2];
575    long tempsize = tempDim[0] * tempDim[1] * tempDim[2];
576    float *temparray = new float[tempsize];
577    //  for(int i=0;i<this->numPixels;i++){ // loop over this->array
578    for(int z=0;z<tempDim[2];z++){
579      for(int y=0;y<tempDim[1];y++){
580        for(int x=0;x<tempDim[0];x++){
581
582          bool isDud = (x<left) || (x>=this->axisDim[0]+left) ||
583            (y<bottom) || (y>=this->axisDim[1]+bottom);
584       
585          int temppos = x + tempDim[0]*y + tempDim[1]*tempDim[0]*z;
586
587          int pos = (x-left) + (y-bottom)*this->axisDim[0] +
588            z*this->axisDim[0]*this->axisDim[1];
589
590          if(isDud) temparray[temppos] = this->par.getBlankPixVal();
591          else temparray[temppos] = this->array[pos];
592 
593          if(this->par.getFlagBaseline() && !isDud)
594            temparray[temppos] += this->baseline[pos];
595
596        }
597      }
598    }
599
600    if(this->par.getFlagCubeTrimmed()){
601      obj.pixels().addOffsets(left,bottom,0);
602    }
603    obj.calcFluxes(temparray, this->axisDim);
604    obj.outputDetectionText(fout,this->logCols,counter);
605    delete [] temparray;
606    delete [] tempDim;
607    fout.close();
608  }
609
610}
Note: See TracBrowser for help on using the repository browser.