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

Last change on this file since 1455 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.