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

Last change on this file since 220 was 220, checked in by Matthew Whiting, 17 years ago
  • Two new files: plots.cc and feedback.cc. Introduced to separate the declarations and definitions for various classes.
  • Mostly just improving the documentation for use with Doxygen.
File size: 14.5 KB
Line 
1#include <iostream>
2#include <sstream>
3#include <fstream>
4#include <iomanip>
5#include <string>
6#include <time.h>
7#include <Cubes/cubes.hh>
8#include <Detection/detection.hh>
9#include <Detection/columns.hh>
10#include <Utils/utils.hh>
11 
12using std::endl;
13using std::setw;
14using std::setprecision;
15using namespace Column;
16
17void Cube::outputDetectionsKarma(std::ostream &stream)
18{
19  /**
20   *  Prints to a stream (provided) the list of detected objects in the cube
21   *   in the format of an annotation file for the Karma suite of programs.
22   *  Annotation file draws a box enclosing the detection, and writes the
23   *   ID number of the detection to the right of the box.
24   */
25
26  string fname = this->par.getImageFile();
27  if(this->par.getFlagSubsection()) fname+=this->par.getSubsection();
28  stream << "# Duchamp Source Finder results for FITS file: " << fname << endl;
29  if(this->par.getFlagFDR())
30    stream<<"# FDR Significance = " << this->par.getAlpha() << endl;
31  else
32    stream<<"# Threshold = " << this->par.getCut() << endl;
33  if(this->par.getFlagATrous()){
34    stream<<"# The a trous reconstruction method was used, with the following parameters." << endl;
35    stream<<"#  ATrous Dimension = " << this->par.getReconDim() << endl;
36    stream<<"#  ATrous Threshold = " << this->par.getAtrousCut() << endl;
37    stream<<"#  ATrous Minimum Scale =" << this->par.getMinScale() << endl;
38    stream<<"#  ATrous Filter = " << this->par.getFilterName() << endl;
39  }
40  else stream << "# No ATrous reconstruction done." << endl;
41  stream << "#\n";
42  stream << "COLOR RED" << endl;
43  stream << "COORD W" << endl;
44  for(int i=0;i<this->objectList.size();i++){
45    if(this->head.isWCS()){
46      float radius = this->objectList[i].getRAWidth()/120.;
47      if(this->objectList[i].getDecWidth()/120.>radius)
48        radius = this->objectList[i].getDecWidth()/120.;
49      stream << "CIRCLE "
50             << this->objectList[i].getRA() << " "
51             << this->objectList[i].getDec() << " "
52             << radius << endl;
53      stream << "TEXT "
54             << this->objectList[i].getRA() << " "
55             << this->objectList[i].getDec() << " "
56             << this->objectList[i].getID() << endl;
57    }
58    else{
59      float radius = this->objectList[i].getXmax() -
60        this->objectList[i].getXmin() + 1;
61      if(this->objectList[i].getYmax()-this->objectList[i].getYmin() + 1
62         > radius)
63        radius = this->objectList[i].getYmax() -
64          this->objectList[i].getYmin() + 1;
65      stream << "CIRCLE "
66             << this->objectList[i].getXcentre() << " "
67             << this->objectList[i].getYcentre() << " "
68             << radius << endl;
69      stream << "TEXT "
70             << this->objectList[i].getXcentre() << " "
71             << this->objectList[i].getYcentre() << " "
72             << this->objectList[i].getID() << endl;
73    }
74  }     
75
76}
77
78void Cube::outputDetectionsVOTable(std::ostream &stream)
79{
80  /**
81   *  Prints to a stream (provided) the list of detected objects in the cube
82   *   in a VOTable format.
83   *  Uses WCS information and assumes WCS parameters have been calculated for each
84   *   detected object.
85   *  If they have not (given by the isWCS() function), then those objects are not written...
86   */
87
88  // Set up Column definitions here
89  vector<Col> localCol;
90  string posUCD[4];
91  localCol.push_back(this->fullCols[NUM]);  // objID
92  localCol.push_back(this->fullCols[NAME]);  // name
93  localCol.push_back(this->fullCols[RA]);  // ra
94  if(makelower(localCol[2].getName())=="ra"){
95    posUCD[0] = "pos.eq.ra;meta.main";
96    posUCD[2] = "phys.angSize;pos.eq.ra";
97  }
98  else{
99    posUCD[0] = "pos.galactic.lat;meta.main";
100    posUCD[2] = "phys.angSize;pos.galactic.lat";
101  }
102  localCol.push_back(this->fullCols[DEC]);  // dec
103  if(makelower(localCol[2].getName())=="dec"){
104    posUCD[1] = "pos.eq.dec;meta.main";
105    posUCD[3] = "phys.angSize;pos.eq.dec";
106  }
107  else{
108    posUCD[1] = "pos.galactic.lon;meta.main";
109    posUCD[3] = "phys.angSize;pos.galactic.lon";
110  }
111  localCol.push_back(this->fullCols[WRA]);  // w_ra
112  localCol.push_back(this->fullCols[WDEC]);  // w_dec
113  localCol.push_back(this->fullCols[VEL]);  // vel
114  localCol.push_back(this->fullCols[WVEL]); // w_vel
115  localCol.push_back(this->fullCols[FINT]); // f_int
116  localCol.push_back(this->fullCols[FPEAK]); // f_peak
117  localCol.push_back(this->fullCols[FLAG]); // flag
118
119  stream<<"<?xml version=\"1.0\"?>"<<endl;
120  stream<<"<VOTABLE version=\"1.1\" xmlns:xsi=\"http://www.w3.org/2001/XMLSchema-instance\""<<endl;
121  stream<<" xsi:noNamespaceSchemaLocation=\"http://www.ivoa.net/xml/VOTable/VOTable/v1.1\">"<<endl;
122
123  stream<<"  <COOSYS ID=\"J2000\" equinox=\"J2000.\" epoch=\"J2000.\" system=\"eq_FK5\"/>"<<endl;
124  stream<<"  <RESOURCE name=\"Duchamp Output\">"<<endl;
125  stream<<"    <TABLE name=\"Detections\">"<<endl;
126  stream<<"      <DESCRIPTION>Detected sources and parameters from running the Duchamp source finder.</DESCRIPTION>"<<endl;
127
128  // PARAM section -- parts that are not entry-specific ie. apply to whole dataset
129  string fname = this->par.getImageFile();
130  if(this->par.getFlagSubsection()) fname+=this->par.getSubsection();
131  stream<<"      <PARAM name=\"FITS file\" datatype=\"char\" ucd=\"meta.file;meta.fits\" value=\"" << fname << "\"/>"<<endl;
132  if(this->par.getFlagFDR())
133    stream<<"      <PARAM name=\"FDR Significance\" datatype=\"float\" ucd=\"stat.param\" value=\"" << this->par.getAlpha() << "\">" << endl;
134  else
135    stream<<"      <PARAM name=\"Threshold\" datatype=\"float\" ucd=\"stat.snr\" value=\"" << this->par.getCut() << "\">" << endl;
136  if(this->par.getFlagATrous()){
137    stream<<"      <PARAM name=\"ATrous note\" datatype=\"char\" ucd=\"meta.note\" value=\"The a trous reconstruction method was used, with the following parameters.\">" << endl;
138    stream<<"      <PARAM name=\"ATrous Dimension\" datatype=\"int\" ucd=\"meta.code;stat\" value=\"" << this->par.getReconDim() << "\">" << endl;
139    stream<<"      <PARAM name=\"ATrous Threshold\" datatype=\"float\" ucd=\"stat.snr\" value=\"" << this->par.getAtrousCut() << "\">" << endl;
140    stream<<"      <PARAM name=\"ATrous Minimum Scale\" datatype=\"int\" ucd=\"stat.param\" value=\"" << this->par.getMinScale() << "\">" << endl;
141    stream<<"      <PARAM name=\"ATrous Filter\" datatype=\"char\" ucd=\"meta.code;stat\" value=\"" << this->par.getFilterName() << "\">" << endl;
142  }
143  // FIELD section -- names, titles and info for each column.
144  stream<<"      <FIELD name=\"ID\" ID=\"col01\" ucd=\"meta.id\" datatype=\"int\" width=\""<<localCol[0].getWidth()<<"\"/>"<<endl;
145  stream<<"      <FIELD name=\"Name\" ID=\"col02\" ucd=\"meta.id;meta.main\" datatype=\"char\" arraysize=\""<<localCol[1].getWidth()<<"\"/>"<<endl;
146  stream<<"      <FIELD name=\""<<localCol[2].getName()<<"\" ID=\"col03\" ucd=\""<<posUCD[0]<<"\" ref=\"J2000\" datatype=\"float\" width=\""<<localCol[2].getWidth()<<"\" precision=\""<<prec[prPOS]<<"\" unit=\"[deg]\"/>"<<endl;
147  stream<<"      <FIELD name=\""<<localCol[3].getName()<<"\" ID=\"col04\" ucd=\""<<posUCD[1]<<"\" ref=\"J2000\" datatype=\"float\" width=\""<<localCol[3].getWidth()<<"\" precision=\""<<prec[prPOS]<<"\" unit=\"[deg]\"/>"<<endl;
148  stream<<"      <FIELD name=\""<<localCol[4].getName()<<"\" ID=\"col05\" ucd=\""<<posUCD[2]<<"\" ref=\"J2000\" datatype=\"float\" width=\""<<localCol[4].getWidth()<<"\" precision=\""<<prec[prWPOS]<<"\" unit=\""<<localCol[4].getUnits()<<"\"/>"<<endl;
149  stream<<"      <FIELD name=\""<<localCol[5].getName()<<"\" ID=\"col06\" ucd=\""<<posUCD[2]<<"\" ref=\"J2000\" datatype=\"float\" width=\""<<localCol[5].getWidth()<<"\" precision=\""<<prec[prWPOS]<<"\" unit=\""<<localCol[5].getUnits()<<"\"/>"<<endl;
150  stream<<"      <FIELD name=\"Vel\" ID=\"col07\" ucd=\"phys.veloc;src.dopplerVeloc\" datatype=\"float\" width=\""<<localCol[6].getWidth()<<"\" precision=\""<<prec[prVEL]<<"\" unit=\""<<localCol[6].getUnits()<<"\"/>"<<endl;
151  stream<<"      <FIELD name=\"w_Vel\" ID=\"col08\" ucd=\"phys.veloc;src.dopplerVeloc;spect.line.width\" datatype=\"float\" width=\""<<localCol[7].getWidth()<<"\" precision=\""<<prec[prVEL]<<"\" unit=\""<<localCol[7].getUnits()<<"\"/>"<<endl;
152  stream<<"      <FIELD name=\"Integrated_Flux\" ID=\"col09\" ucd=\"phot.flux;spect.line.intensity\" datatype=\"float\" width=\""<<localCol[8].getWidth()<<"\" precision=\""<<prec[prFLUX]<<"\" unit=\""<<localCol[8].getUnits()<<"\"/>"<<endl;
153  stream<<"      <FIELD name=\"Peak_Flux\" ID=\"col10\" ucd=\"phot.flux;spect.line.intensity\" datatype=\"float\" width=\""<<localCol[9].getWidth()<<"\" precision=\""<<prec[prFLUX]<<"\" unit=\""<<localCol[9].getUnits()<<"\"/>"<<endl;
154  stream<<"      <FIELD name=\"Flag\" ID=\"col11\" ucd=\"meta.code.qual\" datatype=\"char\" arraysize=\"3\"/>"<<endl;
155
156
157  stream<<"      <DATA>"<<endl
158        <<"        <TABLEDATA>"<<endl;
159
160  stream.setf(std::ios::fixed); 
161  for(int i=0;i<this->objectList.size();i++){
162    if(this->objectList[i].isWCS()){
163      stream<<"        <TR>"<<endl;
164      stream<<"          <TD>"<<setw(localCol[0].getWidth())<<i+1<<"</TD>";
165      stream<<"<TD>" << setw(localCol[1].getWidth()) << this->objectList[i].getName()       <<"</TD>";
166      stream<<setprecision(prec[prPOS]);
167      stream<<"<TD>" << setw(localCol[2].getWidth()) << this->objectList[i].getRA()         <<"</TD>";
168      stream<<"<TD>" << setw(localCol[3].getWidth()) << this->objectList[i].getDec()        <<"</TD>";
169      stream<<setprecision(prec[prWPOS]);
170      stream<<"<TD>" << setw(localCol[4].getWidth()) << this->objectList[i].getRAWidth()    <<"</TD>";
171      stream<<"<TD>" << setw(localCol[5].getWidth()) << this->objectList[i].getDecWidth()   <<"</TD>";
172      stream<<setprecision(prec[prVEL]);
173      stream<<"<TD>" << setw(localCol[6].getWidth()) << this->objectList[i].getVel()        <<"</TD>";
174      stream<<"<TD>" << setw(localCol[7].getWidth()) << this->objectList[i].getVelWidth()   <<"</TD>";
175      stream<<setprecision(prec[prFLUX]);
176      stream<<"<TD>" << setw(localCol[8].getWidth()) << this->objectList[i].getIntegFlux()  <<"</TD>";
177      stream<<"<TD>" << setw(localCol[9].getWidth()) << this->objectList[i].getPeakFlux()   <<"</TD>";
178      stream<<"<TD>" << setw(localCol[10].getWidth())<< this->objectList[i].getFlagText()   <<"</TD>";
179     
180      stream<<endl;
181      stream<<"        </TR>"<<endl;
182    }
183  }
184  stream<<"        </TABLEDATA>"<<endl;
185  stream<<"      </DATA>"<<endl;
186  stream<<"    </TABLE>"<<endl;
187  stream<<"  </RESOURCE>"<<endl;
188  stream<<"</VOTABLE>"<<endl;
189  resetiosflags(std::ios::fixed);
190 
191}
192
193
194void Cube::outputDetectionList()
195{
196  /**
197   *  A front-end to writing the full list of detected objects to a results
198   *   file and to cout.
199   *  Uses outputDetectionTextWCS for each objects.
200   *  Leaves the testing of whether the WCS parameters for each object
201   *   have been calculated to outputDetectionTextWCS.
202   */
203
204  int count=0;
205  int flag;
206  std::ofstream output(this->par.getOutFile().c_str());
207  output<<"Results of the Duchamp source finder: ";
208  time_t now = time(NULL);
209  output << asctime( localtime(&now) );
210  this->showParam(output);
211  output<<"--------------------"<<endl;
212  output<<"Detection threshold = " << this->Stats.getThreshold()
213        <<" " << this->head.getFluxUnits();
214  if(this->par.getFlagFDR())
215    output<<" (or S/N=" << this->Stats.getThresholdSNR()<<")";
216  output<<endl;
217  output<<"  Noise level = " << this->Stats.getMiddle()
218        <<", Noise spread = " << this->Stats.getSpread() << endl;
219  output<<"Total number of detections = "<<this->objectList.size()<<endl;
220  output<<"--------------------"<<endl;
221  this->setupColumns();
222  this->objectList[0].outputDetectionTextHeader(output,this->fullCols);
223  this->objectList[0].outputDetectionTextHeader(std::cout,this->fullCols);
224  for(int i=0;i<this->objectList.size();i++){
225    this->objectList[i].outputDetectionTextWCS(output,this->fullCols);
226    this->objectList[i].outputDetectionTextWCS(std::cout,this->fullCols);
227  }
228}
229
230void Cube::logDetectionList()
231{
232  /**
233   * logDetectionList
234   *  A front-end to writing a list of detected objects to the log file.
235   *  Does not assume WCS, so uses outputDetectionText.
236   *  Designed to be used by searching routines before returning their final list.
237   */
238
239  std::ofstream fout(this->par.getLogFile().c_str(),std::ios::app);
240  this->setupColumns();
241  this->objectList[0].outputDetectionTextHeader(fout,this->logCols);
242  long pos;
243  bool baselineFlag = this->par.getFlagBaseline();
244  for(int objCtr=0;objCtr<this->objectList.size();objCtr++){
245    Detection *obj = new Detection;
246    *obj = objectList[objCtr];
247    if(this->par.getFlagCubeTrimmed()){
248      for(int pix=0;pix<obj->getSize();pix++){
249        // Need to correct the pixels first, as this hasn't been done yet.
250        // Corrections needed for positions (to account for trimmed region)
251        //  and for the baseline removal, if it has happened.
252        // Don't want to keep these changes, so just do it on a dummy variable.
253        pos = obj->getX(pix) + obj->getY(pix)*this->axisDim[0]
254          + obj->getZ(pix)*this->axisDim[0]*this->axisDim[1];
255        obj->setX(pix, obj->getX(pix) + this->par.getBorderLeft() );
256        obj->setY(pix, obj->getY(pix) + this->par.getBorderBottom() );
257        if(baselineFlag) obj->setF(pix, obj->getF(pix) + this->baseline[pos]);   // add in baseline
258      }
259      obj->calcParams();
260    }
261    obj->outputDetectionText(fout,this->logCols,objCtr+1);
262    delete obj;
263  }
264  fout.close();
265}
266
267void Cube::logDetection(Detection obj, int counter)
268{
269  /**
270   *  A front-end to writing a detected object to the log file.
271   *  Does not assume WCS, so uses outputDetectionText.
272   *  Corrects for changes to positions of pixels and removal of baselines.
273   *  Designed to be used by searching routines before returning their final
274   *   list.
275   *  \param obj Detection object to be written : passed by value, as we want
276   *    to potentially change positions etc, but not for the object in the
277   *    calling function.
278   *  \param counter The number to assign to the object : ideally its number
279   *    in a list of some kind.
280   */
281
282  std::ofstream fout(this->par.getLogFile().c_str(),std::ios::app);
283  bool baselineFlag = this->par.getFlagBaseline();
284  if(this->par.getFlagCubeTrimmed()){
285    for(int pix=0;pix<obj.getSize();pix++){
286      // Need to correct the pixels first, as this hasn't been done yet.
287      // Corrections needed for positions (to account for trimmed region)
288      //  and for the baseline removal, if it has happened.
289      // Don't want to keep these changes, so just do it on a dummy variable.
290      long pos = obj.getX(pix) + obj.getY(pix)*this->axisDim[0]
291        + obj.getZ(pix)*this->axisDim[0]*this->axisDim[1];
292      obj.setX(pix, obj.getX(pix) + this->par.getBorderLeft() );
293      obj.setY(pix, obj.getY(pix) + this->par.getBorderBottom() );
294      if(baselineFlag) obj.setF(pix, obj.getF(pix) + this->baseline[pos]);   // add in baseline
295    }
296    obj.calcParams();
297  }
298  obj.outputDetectionText(fout,this->logCols,counter);
299  fout.close();
300}
Note: See TracBrowser for help on using the repository browser.