source: tags/release-1.1.2b/src/Cubes/detectionIO.cc @ 1391

Last change on this file since 1391 was 392, checked in by MatthewWhiting, 17 years ago

Fixed up headers -- should have been done before tagging.
File size: 25.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    output << "\nFull stats:\n" << this->Stats;
429
430    output<<"--------------------\n";
431    output.close();
432  }
433
434  void Cube::outputDetectionList()
435  {
436    /**
437     *  A front-end to writing the full list of detected objects to a results
438     *   file and to cout.
439     *  Uses outputDetectionTextWCS for each objects.
440     *  Leaves the testing of whether the WCS parameters for each object
441     *   have been calculated to outputDetectionTextWCS.
442     */
443
444    std::string outfile;
445    if(this->par.getFlagSeparateHeader()) outfile = this->par.getHeaderFile();
446    else outfile = this->par.getOutFile();
447    std::ofstream output(outfile.c_str(),std::ios::app);
448    output<<"Total number of detections = "<<this->objectList->size()<<endl;
449    output<<"--------------------\n";
450    output.close();
451
452    if(this->par.getFlagSeparateHeader())
453      output.open(this->par.getOutFile().c_str());
454    else
455      output.open(this->par.getOutFile().c_str(),std::ios::app);
456
457    if(this->objectList->size()>0){
458      this->setupColumns();
459      this->objectList->at(0).outputDetectionTextHeaderFull(output,this->fullCols);
460      this->objectList->at(0).outputDetectionTextHeader(std::cout,this->fullCols);
461      for(int i=0;i<this->objectList->size();i++){
462        this->objectList->at(i).outputDetectionTextWCSFull(output,this->fullCols);
463        this->objectList->at(i).outputDetectionTextWCS(std::cout,this->fullCols);
464      }
465    }
466
467    output.close();
468  }
469
470  void Cube::logDetectionList()
471  {
472    /**
473     *  A front-end to writing a list of detected objects to the log file.
474     *  Does not assume WCS, so uses outputDetectionText.
475     *  Designed to be used by searching routines before returning their
476     *   final list.
477     */
478
479    if(this->objectList->size()>0){
480
481      long left = this->par.getBorderLeft();
482      long bottom = this->par.getBorderBottom();
483
484      std::ofstream fout(this->par.getLogFile().c_str(),std::ios::app);
485      this->setupColumns();
486      this->objectList->at(0).outputDetectionTextHeader(fout,this->logCols);
487
488      if(this->par.getFlagBaseline()){
489        for(int i=0;i<this->axisDim[0]*this->axisDim[1]*this->axisDim[2];i++)
490          this->array[i] += this->baseline[i];
491      }
492
493      for(int objCtr=0;objCtr<this->objectList->size();objCtr++){
494        Detection *obj = new Detection;
495        *obj = objectList->at(objCtr);
496        obj->setOffsets(par);
497        if(this->par.getFlagCubeTrimmed()){
498          obj->pixels().addOffsets(left,bottom,0);
499        }
500        obj->calcFluxes(this->array, this->axisDim);
501        obj->outputDetectionText(fout,this->logCols,objCtr+1);
502        delete obj;
503      }
504
505      if(this->par.getFlagBaseline()){
506        for(int i=0;i<this->axisDim[0]*this->axisDim[1]*this->axisDim[2];i++)
507          this->array[i] -= this->baseline[i];
508      }
509      fout.close();
510 
511    }
512
513  }
514
515  void Cube::logDetection(Detection obj, int counter)
516  {
517    /**
518     *  A front-end to writing a detected object to the log file.
519     *  Does not assume WCS, so uses outputDetectionText.
520     *  Corrects for changes to positions of pixels and removal of baselines.
521     *  Designed to be used by searching routines before returning their final
522     *   list.
523     *  \param obj Detection object to be written : passed by value, as we want
524     *    to potentially change positions etc, but not for the object in the
525     *    calling function.
526     *  \param counter The number to assign to the object : ideally its number
527     *    in a list of some kind.
528     */
529
530    std::ofstream fout(this->par.getLogFile().c_str(),std::ios::app);
531    // Need to deal with possibility of trimmed array
532    long left = this->par.getBorderLeft();
533    long right = this->par.getBorderRight();
534    long top = this->par.getBorderTop();
535    long bottom = this->par.getBorderBottom();
536    long *tempDim = new long[3];
537    tempDim[0] = (this->axisDim[0] + left + right);
538    tempDim[1] = (this->axisDim[1] + bottom + top);
539    tempDim[2] = this->axisDim[2];
540    long tempsize = tempDim[0] * tempDim[1] * tempDim[2];
541    float *temparray = new float[tempsize];
542    //  for(int i=0;i<this->numPixels;i++){ // loop over this->array
543    for(int z=0;z<tempDim[2];z++){
544      for(int y=0;y<tempDim[1];y++){
545        for(int x=0;x<tempDim[0];x++){
546
547          bool isDud = (x<left) || (x>=this->axisDim[0]+left) ||
548            (y<bottom) || (y>=this->axisDim[1]+bottom);
549       
550          int temppos = x + tempDim[0]*y + tempDim[1]*tempDim[0]*z;
551
552          int pos = (x-left) + (y-bottom)*this->axisDim[0] +
553            z*this->axisDim[0]*this->axisDim[1];
554
555          if(isDud) temparray[temppos] = this->par.getBlankPixVal();
556          else temparray[temppos] = this->array[pos];
557 
558          if(this->par.getFlagBaseline() && !isDud)
559            temparray[temppos] += this->baseline[pos];
560
561        }
562      }
563    }
564
565    if(this->par.getFlagCubeTrimmed()){
566      obj.pixels().addOffsets(left,bottom,0);
567    }
568    obj.calcFluxes(temparray, this->axisDim);
569    obj.outputDetectionText(fout,this->logCols,counter);
570    delete [] temparray;
571    delete [] tempDim;
572    fout.close();
573  }
574
575}
Note: See TracBrowser for help on using the repository browser.