source: branches/pixel-map-branch/src/Cubes/detectionIO.cc

Last change on this file was 252, checked in by Matthew Whiting, 17 years ago
  • Have put all classes in the files in src/PixelMap/ into a PixelInfo? namespace.
  • Added "using namespace PixelInfo?" to all necessary files.
  • Removed "friend class Detection" from Voxel and Object3D classes -- not necessary and complicated now by them being in the namespace
  • Various minor changes, including fixing up commentary.
File size: 17.4 KB
Line 
1#include <iostream>
2#include <sstream>
3#include <fstream>
4#include <iomanip>
5#include <string>
6#include <vector>
7#include <time.h>
8#include <Cubes/cubes.hh>
9#include <PixelMap/Object3D.hh>
10#include <Detection/detection.hh>
11#include <Detection/columns.hh>
12#include <Utils/utils.hh>
13 
14using std::endl;
15using std::setw;
16using std::setprecision;
17using namespace Column;
18using namespace PixelInfo;
19
20void Cube::outputDetectionsKarma(std::ostream &stream)
21{
22  /**
23   *  Prints to a stream (provided) the list of detected objects in the cube
24   *   in the format of an annotation file for the Karma suite of programs.
25   *  Annotation file draws a box enclosing the detection, and writes the
26   *   ID number of the detection to the right of the box.
27   */
28
29  std::string fname = this->par.getImageFile();
30  if(this->par.getFlagSubsection()) fname+=this->par.getSubsection();
31  stream << "# Duchamp Source Finder results for FITS file: " << fname << endl;
32  if(this->par.getFlagFDR())
33    stream<<"# FDR Significance = " << this->par.getAlpha() << endl;
34  else
35    stream<<"# Threshold = " << this->par.getCut() << endl;
36  if(this->par.getFlagATrous()){
37    stream<<"# The a trous reconstruction method was used, with the following parameters." << endl;
38    stream<<"#  ATrous Dimension = " << this->par.getReconDim() << endl;
39    stream<<"#  ATrous Threshold = " << this->par.getAtrousCut() << endl;
40    stream<<"#  ATrous Minimum Scale =" << this->par.getMinScale() << endl;
41    stream<<"#  ATrous Filter = " << this->par.getFilterName() << endl;
42  }
43  else stream << "# No ATrous reconstruction done." << endl;
44  stream << "#\n";
45  stream << "COLOR RED" << endl;
46  stream << "COORD W" << endl;
47  for(int i=0;i<this->objectList.size();i++){
48    if(this->head.isWCS()){
49      float radius = this->objectList[i].getRAWidth()/120.;
50      if(this->objectList[i].getDecWidth()/120.>radius)
51        radius = this->objectList[i].getDecWidth()/120.;
52      stream << "CIRCLE "
53             << this->objectList[i].getRA() << " "
54             << this->objectList[i].getDec() << " "
55             << radius << endl;
56      stream << "TEXT "
57             << this->objectList[i].getRA() << " "
58             << this->objectList[i].getDec() << " "
59             << this->objectList[i].getID() << endl;
60    }
61    else{
62      float radius = this->objectList[i].getXmax() -
63        this->objectList[i].getXmin() + 1;
64      if(this->objectList[i].getYmax()-this->objectList[i].getYmin() + 1
65         > radius)
66        radius = this->objectList[i].getYmax() -
67          this->objectList[i].getYmin() + 1;
68      stream << "CIRCLE "
69             << this->objectList[i].getXcentre() << " "
70             << this->objectList[i].getYcentre() << " "
71             << radius << endl;
72      stream << "TEXT "
73             << this->objectList[i].getXcentre() << " "
74             << this->objectList[i].getYcentre() << " "
75             << this->objectList[i].getID() << endl;
76    }
77  }     
78
79}
80
81std::string fixUnitsVOT(std::string oldstring)
82{
83  /** Fix a string containing units to make it acceptable for a VOTable.
84   *
85   * This function makes the provided units string acceptable
86   * according to the standard found at
87   * http://vizier.u-strasbg.fr/doc/catstd-3.2.htx
88   * This should then be able to convert the units used in the text
89   * table to units suitable for putting in a VOTable.
90   *
91   * Specifically, it removes any square brackets [] from the
92   * start/end of the string, and replaces blank spaces (representing
93   * multiplication) with a '.' (full stop).
94   */
95
96  std::string newstring;
97  for(int i=0;i<oldstring.size();i++){
98    if((oldstring[i]!='[')&&(oldstring[i]!=']')){
99      if(oldstring[i]==' ') newstring += '.';
100      else newstring += oldstring[i];
101    }
102  }
103  return newstring; 
104}
105
106void Cube::outputDetectionsVOTable(std::ostream &stream)
107{
108  /**
109   *  Prints to a stream (provided) the list of detected objects in the cube
110   *   in a VOTable format.
111   *  Uses WCS information and assumes WCS parameters have been calculated for each
112   *   detected object.
113   *  If they have not (given by the isWCS() function), then those objects are not written...
114   */
115
116  // Set up Column definitions here
117  std::vector<Col> localCol;
118  std::string posUCD[4];
119  localCol.push_back(this->fullCols[NUM]);  // objID
120  localCol.push_back(this->fullCols[NAME]);  // name
121  localCol.push_back(this->fullCols[RA]);  // ra
122  if(makelower(localCol[2].getName())=="ra"){
123    posUCD[0] = "pos.eq.ra;meta.main";
124    posUCD[2] = "phys.angSize;pos.eq.ra";
125  }
126  else{
127    posUCD[0] = "pos.galactic.lat;meta.main";
128    posUCD[2] = "phys.angSize;pos.galactic.lat";
129  }
130  localCol.push_back(this->fullCols[DEC]);  // dec
131  if(makelower(localCol[2].getName())=="dec"){
132    posUCD[1] = "pos.eq.dec;meta.main";
133    posUCD[3] = "phys.angSize;pos.eq.dec";
134  }
135  else{
136    posUCD[1] = "pos.galactic.lon;meta.main";
137    posUCD[3] = "phys.angSize;pos.galactic.lon";
138  }
139  localCol.push_back(this->fullCols[WRA]);  // w_ra
140  localCol.push_back(this->fullCols[WDEC]);  // w_dec
141  localCol.push_back(this->fullCols[VEL]);  // vel
142  localCol.push_back(this->fullCols[WVEL]); // w_vel
143  localCol.push_back(this->fullCols[FINT]); // f_int
144  localCol.push_back(this->fullCols[FPEAK]); // f_peak
145  localCol.push_back(this->fullCols[FLAG]); // flag
146
147  stream<<"<?xml version=\"1.0\"?>\n";
148  stream<<"<VOTABLE version=\"1.1\" xmlns:xsi=\"http://www.w3.org/2001/XMLSchema-instance\"\n";
149  stream<<" xsi:noNamespaceSchemaLocation=\"http://www.ivoa.net/xml/VOTable/VOTable/v1.1\">\n";
150
151  stream<<"  <COOSYS ID=\"J2000\" equinox=\"J2000.\" epoch=\"J2000.\" system=\"eq_FK5\"/>\n";
152  stream<<"  <RESOURCE name=\"Duchamp Output\">\n";
153  stream<<"    <TABLE name=\"Detections\">\n";
154  stream<<"      <DESCRIPTION>Detected sources and parameters from running the Duchamp source finder.</DESCRIPTION>\n";
155
156  // PARAM section -- parts that are not entry-specific ie. apply to whole dataset
157  std::string fname = this->par.getImageFile();
158  if(this->par.getFlagSubsection()) fname+=this->par.getSubsection();
159  stream<<"      <PARAM name=\"FITS file\" datatype=\"char\" ucd=\"meta.file;meta.fits\" value=\""
160        << fname << "\" arraysize=\""<<fname.size()<<"\"/>\n";
161  if(this->par.getFlagFDR())
162    stream<<"      <PARAM name=\"FDR Significance\" datatype=\"float\" ucd=\"stat.param\" value=\""
163          << this->par.getAlpha() << "\"/>\n";
164  else
165    stream<<"      <PARAM name=\"Threshold\" datatype=\"float\" ucd=\"stat.snr\" value=\""
166          << this->par.getCut() << "\"/>\n";
167  if(this->par.getFlagATrous()){
168    std::string note = "The a trous reconstruction method was used, with the following parameters.";
169    stream<<"      <PARAM name=\"ATrous note\" datatype=\"char\" ucd=\"meta.note\" value=\""
170          <<note << "\" arraysize=\"" << note.size() << "\"/>\n";
171    stream<<"      <PARAM name=\"ATrous Dimension\" datatype=\"int\" ucd=\"meta.code;stat\" value=\""
172          << this->par.getReconDim() << "\"/>\n";
173    stream<<"      <PARAM name=\"ATrous Threshold\" datatype=\"float\" ucd=\"stat.snr\" value=\""
174          << this->par.getAtrousCut() << "\"/>\n";
175    stream<<"      <PARAM name=\"ATrous Minimum Scale\" datatype=\"int\" ucd=\"stat.param\" value=\""
176          << this->par.getMinScale() << "\"/>\n";
177    stream<<"      <PARAM name=\"ATrous Filter\" datatype=\"char\" ucd=\"meta.code;stat\" value=\""
178          << this->par.getFilterName() << "\" arraysize=\"" << this->par.getFilterName().size() << "\"/>\n";
179  }
180  // FIELD section -- names, titles and info for each column.
181  stream<<"      <FIELD name=\"ID\" ID=\"col01\" ucd=\"meta.id\" datatype=\"int\" width=\""
182        <<localCol[0].getWidth()<<"\" unit=\"--\"/>\n";
183  stream<<"      <FIELD name=\"Name\" ID=\"col02\" ucd=\"meta.id;meta.main\" datatype=\"char\" arraysize=\""
184        <<localCol[1].getWidth()<<"\" unit=\"--\"/>\n";
185  stream<<"      <FIELD name=\""<<localCol[2].getName()<<"\" ID=\"col03\" ucd=\""
186        <<posUCD[0]<<"\" ref=\"J2000\" datatype=\"float\" width=\""
187        <<localCol[2].getWidth()<<"\" precision=\""<<prec[prPOS]<<"\" unit=\"deg\"/>\n";
188  stream<<"      <FIELD name=\""<<localCol[3].getName()<<"\" ID=\"col04\" ucd=\""
189        <<posUCD[1]<<"\" ref=\"J2000\" datatype=\"float\" width=\""
190        <<localCol[3].getWidth()<<"\" precision=\""<<prec[prPOS]<<"\" unit=\"deg\"/>\n";
191  stream<<"      <FIELD name=\""<<localCol[4].getName()<<"\" ID=\"col05\" ucd=\""
192        <<posUCD[2]<<"\" ref=\"J2000\" datatype=\"float\" width=\""
193        <<localCol[4].getWidth()<<"\" precision=\""<<prec[prWPOS]<<"\" unit=\""
194        <<fixUnitsVOT(localCol[4].getUnits())<<"\"/>\n";
195  stream<<"      <FIELD name=\""<<localCol[5].getName()<<"\" ID=\"col06\" ucd=\""
196        <<posUCD[2]<<"\" ref=\"J2000\" datatype=\"float\" width=\""
197        <<localCol[5].getWidth()<<"\" precision=\""<<prec[prWPOS]<<"\" unit=\""
198        <<fixUnitsVOT(localCol[5].getUnits())<<"\"/>\n";
199  stream<<"      <FIELD name=\"Vel\" ID=\"col07\" ucd=\"phys.veloc;src.dopplerVeloc\" datatype=\"float\" width=\""
200        <<localCol[6].getWidth()<<"\" precision=\""<<prec[prVEL]<<"\" unit=\""
201        <<fixUnitsVOT(localCol[6].getUnits())<<"\"/>\n";
202  stream<<"      <FIELD name=\"w_Vel\" ID=\"col08\" ucd=\"phys.veloc;src.dopplerVeloc;spect.line.width\""
203        <<" datatype=\"float\" width=\""<<localCol[7].getWidth()<<"\" precision=\""
204        <<prec[prVEL]<<"\" unit=\""<<fixUnitsVOT(localCol[7].getUnits())<<"\"/>\n";
205  stream<<"      <FIELD name=\"Integrated_Flux\" ID=\"col09\" ucd=\"phot.flux;spect.line.intensity\""
206        <<" datatype=\"float\" width=\""<<localCol[8].getWidth()<<"\" precision=\""
207        <<prec[prFLUX]<<"\" unit=\""<<fixUnitsVOT(localCol[8].getUnits())<<"\"/>\n";
208  stream<<"      <FIELD name=\"Peak_Flux\" ID=\"col10\" ucd=\"phot.flux;spect.line.intensity\""
209        <<" datatype=\"float\" width=\""<<localCol[9].getWidth()<<"\" precision=\""
210        <<prec[prFLUX]<<"\" unit=\""<<fixUnitsVOT(localCol[9].getUnits())<<"\"/>\n";
211  stream<<"      <FIELD name=\"Flag\" ID=\"col11\" ucd=\"meta.code.qual\" datatype=\"char\" arraysize=\"3\" unit=\"--\"/>\n";
212
213
214  stream<<"      <DATA>\n"
215        <<"        <TABLEDATA>\n";
216
217  stream.setf(std::ios::fixed); 
218  for(int i=0;i<this->objectList.size();i++){
219    if(this->objectList[i].isWCS()){
220      stream<<"        <TR>\n";
221      stream<<"          <TD>"<<setw(localCol[0].getWidth())<<i+1<<"</TD>";
222      stream<<"<TD>" << setw(localCol[1].getWidth()) << this->objectList[i].getName()       <<"</TD>";
223      stream<<setprecision(prec[prPOS]);
224      stream<<"<TD>" << setw(localCol[2].getWidth()) << this->objectList[i].getRA()         <<"</TD>";
225      stream<<"<TD>" << setw(localCol[3].getWidth()) << this->objectList[i].getDec()        <<"</TD>";
226      stream<<setprecision(prec[prWPOS]);
227      stream<<"<TD>" << setw(localCol[4].getWidth()) << this->objectList[i].getRAWidth()    <<"</TD>";
228      stream<<"<TD>" << setw(localCol[5].getWidth()) << this->objectList[i].getDecWidth()   <<"</TD>";
229      stream<<setprecision(prec[prVEL]);
230      stream<<"<TD>" << setw(localCol[6].getWidth()) << this->objectList[i].getVel()        <<"</TD>";
231      stream<<"<TD>" << setw(localCol[7].getWidth()) << this->objectList[i].getVelWidth()   <<"</TD>";
232      stream<<setprecision(prec[prFLUX]);
233      stream<<"<TD>" << setw(localCol[8].getWidth()) << this->objectList[i].getIntegFlux()  <<"</TD>";
234      stream<<"<TD>" << setw(localCol[9].getWidth()) << this->objectList[i].getPeakFlux()   <<"</TD>";
235      stream<<"<TD>" << setw(localCol[10].getWidth())<< this->objectList[i].getFlagText()   <<"</TD>";
236     
237      stream<<endl;
238      stream<<"        </TR>\n";
239    }
240  }
241  stream<<"        </TABLEDATA>\n";
242  stream<<"      </DATA>\n";
243  stream<<"    </TABLE>\n";
244  stream<<"  </RESOURCE>\n";
245  stream<<"</VOTABLE>\n";
246  resetiosflags(std::ios::fixed);
247 
248}
249
250
251void Cube::outputDetectionList()
252{
253  /**
254   *  A front-end to writing the full list of detected objects to a results
255   *   file and to cout.
256   *  Uses outputDetectionTextWCS for each objects.
257   *  Leaves the testing of whether the WCS parameters for each object
258   *   have been calculated to outputDetectionTextWCS.
259   */
260
261  int count=0;
262  int flag;
263  std::ofstream output(this->par.getOutFile().c_str());
264  output<<"Results of the Duchamp source finder: ";
265  time_t now = time(NULL);
266  output << asctime( localtime(&now) );
267  this->showParam(output);
268  output<<"--------------------\n";
269  output<<"Detection threshold = " << this->Stats.getThreshold()
270        <<" " << this->head.getFluxUnits();
271  if(this->par.getFlagFDR())
272    output<<" (or S/N=" << this->Stats.getThresholdSNR()<<")";
273  output<<endl;
274  if(!this->par.getFlagUserThreshold())
275    output<<"  Noise level = " << this->Stats.getMiddle()
276          <<", Noise spread = " << this->Stats.getSpread() << endl;
277  output<<"Total number of detections = "<<this->objectList.size()<<endl;
278  output<<"--------------------\n";
279  this->setupColumns();
280  this->objectList[0].outputDetectionTextHeader(output,this->fullCols);
281  this->objectList[0].outputDetectionTextHeader(std::cout,this->fullCols);
282  for(int i=0;i<this->objectList.size();i++){
283    this->objectList[i].outputDetectionTextWCS(output,this->fullCols);
284    this->objectList[i].outputDetectionTextWCS(std::cout,this->fullCols);
285  }
286}
287
288void Cube::logDetectionList()
289{
290  /**
291   * logDetectionList
292   *  A front-end to writing a list of detected objects to the log file.
293   *  Does not assume WCS, so uses outputDetectionText.
294   *  Designed to be used by searching routines before returning their final list.
295   */
296
297  long left = this->par.getBorderLeft();
298  long right = this->par.getBorderRight();
299  long top = this->par.getBorderTop();
300  long bottom = this->par.getBorderBottom();
301
302  std::ofstream fout(this->par.getLogFile().c_str(),std::ios::app);
303  this->setupColumns();
304  this->objectList[0].outputDetectionTextHeader(fout,this->logCols);
305  long pos;
306
307//   std::ofstream ftemp("temp2.txt");
308//   for(int objCtr=0;objCtr<this->objectList.size();objCtr++){
309//     this->objectList[objCtr].pixels().order();
310//     ftemp << "Object #" << objCtr+1
311//        << ": size = " << this->objectList[objCtr].getSize()<<"\n";
312//     ftemp << this->objectList[objCtr];
313//   }
314//   ftemp.close();
315
316
317//   // Need to deal with possibility of trimmed array
318//   long *tempDim = new long[3];
319//   tempDim[0] = (this->axisDim[0] + left + right);
320//   tempDim[1] = (this->axisDim[1] + bottom + top);
321//   tempDim[2] = this->axisDim[2];
322//   long tempsize = tempDim[0] * tempDim[1] * tempDim[2];
323//   float *temparray = new float[tempsize];
324//   //  for(int i=0;i<this->numPixels;i++){ // loop over this->array
325//   for(int z=0;z<tempDim[2];z++){
326//     for(int y=0;y<tempDim[1];y++){
327//       for(int x=0;x<tempDim[0];x++){
328
329//      bool isDud = (x<left) || (x>=this->axisDim[0]+left) ||
330//        (y<bottom) || (y>=this->axisDim[1]+bottom);
331       
332//      int temppos = x + tempDim[0]*y + tempDim[1]*tempDim[0]*z;
333
334//      int pos = (x-left) + (y-bottom)*this->axisDim[0] +
335//        z*this->axisDim[0]*this->axisDim[1];
336
337//      if(isDud) temparray[temppos] = this->par.getBlankPixVal();
338//      else temparray[temppos] = this->array[pos];
339 
340//      if(this->par.getFlagBaseline() && !isDud)
341//        temparray[temppos] += this->baseline[pos];
342
343//       }
344//     }
345//   }
346
347  if(this->par.getFlagBaseline()){
348    for(int i=0;i<this->axisDim[0]*this->axisDim[1]*this->axisDim[2];i++)
349      this->array[i] += this->baseline[i];
350  }
351
352  for(int objCtr=0;objCtr<this->objectList.size();objCtr++){
353    Detection *obj = new Detection;
354    *obj = objectList[objCtr];
355    obj->setOffsets(par);
356    obj->calcFluxes(this->array, this->axisDim);
357    if(this->par.getFlagCubeTrimmed()){
358      obj->pixels().addOffsets(left,bottom,0);
359    }
360    obj->outputDetectionText(fout,this->logCols,objCtr+1);
361    delete obj;
362  }
363
364  if(this->par.getFlagBaseline()){
365    for(int i=0;i<this->axisDim[0]*this->axisDim[1]*this->axisDim[2];i++)
366      this->array[i] -= this->baseline[i];
367  }
368
369//   delete [] temparray;
370//   delete [] tempDim;
371  fout.close();
372}
373
374void Cube::logDetection(Detection obj, int counter)
375{
376  /**
377   *  A front-end to writing a detected object to the log file.
378   *  Does not assume WCS, so uses outputDetectionText.
379   *  Corrects for changes to positions of pixels and removal of baselines.
380   *  Designed to be used by searching routines before returning their final
381   *   list.
382   *  \param obj Detection object to be written : passed by value, as we want
383   *    to potentially change positions etc, but not for the object in the
384   *    calling function.
385   *  \param counter The number to assign to the object : ideally its number
386   *    in a list of some kind.
387   */
388
389  std::ofstream fout(this->par.getLogFile().c_str(),std::ios::app);
390  // Need to deal with possibility of trimmed array
391  long left = this->par.getBorderLeft();
392  long right = this->par.getBorderRight();
393  long top = this->par.getBorderTop();
394  long bottom = this->par.getBorderBottom();
395  long *tempDim = new long[3];
396  tempDim[0] = (this->axisDim[0] + left + right);
397  tempDim[1] = (this->axisDim[1] + bottom + top);
398  tempDim[2] = this->axisDim[2];
399  long tempsize = tempDim[0] * tempDim[1] * tempDim[2];
400  float *temparray = new float[tempsize];
401  //  for(int i=0;i<this->numPixels;i++){ // loop over this->array
402  for(int z=0;z<tempDim[2];z++){
403    for(int y=0;y<tempDim[1];y++){
404      for(int x=0;x<tempDim[0];x++){
405
406        bool isDud = (x<left) || (x>=this->axisDim[0]+left) ||
407          (y<bottom) || (y>=this->axisDim[1]+bottom);
408       
409        int temppos = x + tempDim[0]*y + tempDim[1]*tempDim[0]*z;
410
411        int pos = (x-left) + (y-bottom)*this->axisDim[0] +
412          z*this->axisDim[0]*this->axisDim[1];
413
414        if(isDud) temparray[temppos] = this->par.getBlankPixVal();
415        else temparray[temppos] = this->array[pos];
416 
417        if(this->par.getFlagBaseline() && !isDud)
418          temparray[temppos] += this->baseline[pos];
419
420      }
421    }
422  }
423
424  if(this->par.getFlagCubeTrimmed()){
425    obj.pixels().addOffsets(left,bottom,0);
426  }
427  obj.calcFluxes(temparray, this->axisDim);
428  obj.outputDetectionText(fout,this->logCols,counter);
429  delete [] temparray;
430  delete [] tempDim;
431  fout.close();
432}
Note: See TracBrowser for help on using the repository browser.