source: tags/release-1.1.4/src/Cubes/detectionIO.cc

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

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