source: tags/release-1.1/src/Cubes/detectionIO.cc @ 1323

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

Adding distribution text at the start of each file...
File size: 24.4 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  this->setupColumns();
436  this->objectList->at(0).outputDetectionTextHeaderFull(output,this->fullCols);
437  this->objectList->at(0).outputDetectionTextHeader(std::cout,this->fullCols);
438  for(int i=0;i<this->objectList->size();i++){
439    this->objectList->at(i).outputDetectionTextWCSFull(output,this->fullCols);
440    this->objectList->at(i).outputDetectionTextWCS(std::cout,this->fullCols);
441  }
442
443  output.close();
444}
445
446void Cube::logDetectionList()
447{
448  /**
449   *  A front-end to writing a list of detected objects to the log file.
450   *  Does not assume WCS, so uses outputDetectionText.
451   *  Designed to be used by searching routines before returning their
452   *   final list.
453   */
454
455  long left = this->par.getBorderLeft();
456  long bottom = this->par.getBorderBottom();
457
458  std::ofstream fout(this->par.getLogFile().c_str(),std::ios::app);
459  this->setupColumns();
460  this->objectList->at(0).outputDetectionTextHeader(fout,this->logCols);
461
462  if(this->par.getFlagBaseline()){
463    for(int i=0;i<this->axisDim[0]*this->axisDim[1]*this->axisDim[2];i++)
464      this->array[i] += this->baseline[i];
465  }
466
467  for(int objCtr=0;objCtr<this->objectList->size();objCtr++){
468    Detection *obj = new Detection;
469    *obj = objectList->at(objCtr);
470    obj->setOffsets(par);
471    if(this->par.getFlagCubeTrimmed()){
472      obj->pixels().addOffsets(left,bottom,0);
473    }
474    obj->calcFluxes(this->array, this->axisDim);
475    obj->outputDetectionText(fout,this->logCols,objCtr+1);
476    delete obj;
477  }
478
479  if(this->par.getFlagBaseline()){
480    for(int i=0;i<this->axisDim[0]*this->axisDim[1]*this->axisDim[2];i++)
481      this->array[i] -= this->baseline[i];
482  }
483  fout.close();
484}
485
486void Cube::logDetection(Detection obj, int counter)
487{
488  /**
489   *  A front-end to writing a detected object to the log file.
490   *  Does not assume WCS, so uses outputDetectionText.
491   *  Corrects for changes to positions of pixels and removal of baselines.
492   *  Designed to be used by searching routines before returning their final
493   *   list.
494   *  \param obj Detection object to be written : passed by value, as we want
495   *    to potentially change positions etc, but not for the object in the
496   *    calling function.
497   *  \param counter The number to assign to the object : ideally its number
498   *    in a list of some kind.
499   */
500
501  std::ofstream fout(this->par.getLogFile().c_str(),std::ios::app);
502  // Need to deal with possibility of trimmed array
503  long left = this->par.getBorderLeft();
504  long right = this->par.getBorderRight();
505  long top = this->par.getBorderTop();
506  long bottom = this->par.getBorderBottom();
507  long *tempDim = new long[3];
508  tempDim[0] = (this->axisDim[0] + left + right);
509  tempDim[1] = (this->axisDim[1] + bottom + top);
510  tempDim[2] = this->axisDim[2];
511  long tempsize = tempDim[0] * tempDim[1] * tempDim[2];
512  float *temparray = new float[tempsize];
513  //  for(int i=0;i<this->numPixels;i++){ // loop over this->array
514  for(int z=0;z<tempDim[2];z++){
515    for(int y=0;y<tempDim[1];y++){
516      for(int x=0;x<tempDim[0];x++){
517
518        bool isDud = (x<left) || (x>=this->axisDim[0]+left) ||
519          (y<bottom) || (y>=this->axisDim[1]+bottom);
520       
521        int temppos = x + tempDim[0]*y + tempDim[1]*tempDim[0]*z;
522
523        int pos = (x-left) + (y-bottom)*this->axisDim[0] +
524          z*this->axisDim[0]*this->axisDim[1];
525
526        if(isDud) temparray[temppos] = this->par.getBlankPixVal();
527        else temparray[temppos] = this->array[pos];
528 
529        if(this->par.getFlagBaseline() && !isDud)
530          temparray[temppos] += this->baseline[pos];
531
532      }
533    }
534  }
535
536  if(this->par.getFlagCubeTrimmed()){
537    obj.pixels().addOffsets(left,bottom,0);
538  }
539  obj.calcFluxes(temparray, this->axisDim);
540  obj.outputDetectionText(fout,this->logCols,counter);
541  delete [] temparray;
542  delete [] tempDim;
543  fout.close();
544}
Note: See TracBrowser for help on using the repository browser.