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

Last change on this file since 321 was 312, checked in by Matthew Whiting, 17 years ago

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