source: trunk/src/Cubes/cubes.hh @ 221

Last change on this file since 221 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: 22.7 KB
Line 
1#ifndef CUBES_H
2#define CUBES_H
3
4#include <iostream>
5#include <string>
6#include <vector>
7
8#ifndef PARAM_H
9#include <param.hh>
10#endif
11#ifndef DETECTION_H
12#include <Detection/detection.hh>
13#endif
14#ifndef COLUMNS_H
15#include <Detection/columns.hh>
16#endif
17#ifndef PLOTS_H
18#include <Cubes/plots.hh>
19#endif
20#ifndef STATS_H
21#include <Utils/Statistics.hh>
22#endif
23
24using std::string;
25using std::vector;
26using namespace Column;
27using namespace Statistics;
28
29/****************************************************************/
30/**
31 * Base class for the image container.
32 *
33 * Definition of an n-dimensional data array:
34 *     array of pixel values, size & dimensions
35 *     array of Detection objects
36 */
37
38class DataArray
39{
40public:
41  DataArray();  ///< Basic DataArray constructor
42  DataArray(short int nDim); ///< Basic nDim-dimensional DataArray constructor
43  DataArray(short int nDim, long size);///< Basic nDim-dimensional DataArray constructor, specifying size.
44  DataArray(short int nDim, long *dimensions); ///< Basic nDim-dimensional DataArray constructor, specifying size of dimensions.
45  virtual ~DataArray(); ///< Basic DataArray constructor
46
47  // Size and Dimension related
48  long               getDimX();   ///< Return X dimension, defaulting to 0.
49  long               getDimY();   ///< Return Y dimension, defaulting to 1.
50  long               getDimZ();   ///< Return Z dimension, defaulting to 1.
51  long               getSize();   ///< Return number of voxels.
52  short int          getNumDim(); ///< Return number of dimensions.
53  void               getDim(long &x, long &y, long &z); ///< Return first three dimensional axes.
54
55  // Related to the various arrays
56  void               getDimArray(long *output); ///< Return array of dimensional sizes.
57  void               getArray(float *output); ///< Return pixel value array.
58  virtual void       saveArray(float *input, long size); ///< Save pixel value array.
59  float              getPixValue(long pos); ///< Returns the flux value at voxel number pos.
60  void               setPixValue(long pos, float f); ///< Sets the value of pixel pos to f.
61
62  // Related to the object lists
63  Detection          getObject(long number); ///< Returns the Detection at position 'number' in the list.
64  void               addObject(Detection object); ///< Adds a single detection to the object list.
65  vector <Detection> getObjectList(); ///< Returns the full list of Detections.
66  void               addObjectList(vector <Detection> newlist);  ///< Adds all objects in a detection list to the object list.
67  void               addObjectOffsets(); ///< Add pixel offsets to object coordinates.
68  long               getNumObj(); ///< Returns the number of detected objects.
69  void               clearDetectionList(); ///< Delete all Detections from the object list.
70
71  // Parameter list related.
72  int                readParam(string paramfile); ///< Read parameters.
73  void               showParam(std::ostream &stream); ///< Output the Param set.
74  Param              getParam(); ///< Return the Param set.
75  void               saveParam(Param newpar); ///< Save a Param set to the Cube.
76  Param&             pars(); ///< Provides a reference to the Param set.
77               
78  bool               isBlank(int vox); ///< Is the voxel number given by vox a BLANK value?
79
80  // Statistics
81  StatsContainer<float> getStats(); ///< Returns the StatsContainer.
82  StatsContainer<float>& stats(); ///< Provides a reference to the StatsContainer.
83  void saveStats(StatsContainer<float> newStats); ///< Save a StatsContainer to the Cube.
84
85  bool isDetection(float value); ///< A detection test for value.
86  bool isDetection(long voxel); ///< A detection test for pixel.
87
88  friend std::ostream& operator<< (std::ostream& theStream, DataArray &array);
89  ///< Output operator for DataArray.
90
91protected:
92  short int               numDim;     ///< Number of dimensions.
93  long                   *axisDim;    ///< Array of dimensions of cube (ie. how large in each direction).
94  long                    numPixels;  ///< Total number of pixels in cube.
95  float                  *array;      ///< Array of data.
96  vector <Detection>      objectList; ///< The list of detected objects.
97  Param                   par;        ///< A parameter list.
98  StatsContainer<float>   Stats;      ///< The statistics for the DataArray.
99};
100
101////////////
102//// DataArray inline function definitions
103////////////
104
105inline long      DataArray::getDimX(){ if(numDim>0) return this->axisDim[0];else return 0;};
106inline long      DataArray::getDimY(){ if(numDim>1) return this->axisDim[1];else return 1;};
107inline long      DataArray::getDimZ(){ if(numDim>2) return this->axisDim[2];else return 1;};
108inline long      DataArray::getSize(){ return this->numPixels; };
109inline short int DataArray::getNumDim(){ return this->numDim; };
110
111inline float     DataArray::getPixValue(long pos){ return array[pos]; };
112inline void      DataArray::setPixValue(long pos, float f){array[pos] = f;};
113inline Detection DataArray::getObject(long number){ return objectList[number]; };
114inline vector <Detection> DataArray::getObjectList(){ return objectList; };
115inline long      DataArray::getNumObj(){ return objectList.size(); };
116inline void      DataArray::clearDetectionList(){ this->objectList.clear(); };
117inline int       DataArray::readParam(string paramfile){
118    /**
119     *  Uses Param::readParams() to read parameters from a file.
120     *  \param paramfile The file to be read.
121     */
122  return par.readParams(paramfile); };
123inline void      DataArray::showParam(std::ostream &stream){ stream << this->par; };
124inline Param     DataArray::getParam(){ return this->par; };
125inline void      DataArray::saveParam(Param newpar){this->par = newpar;};
126inline Param&    DataArray::pars(){ Param &rpar = this->par; return rpar; };
127inline bool      DataArray::isBlank(int vox){ return this->par.isBlank(this->array[vox]); };
128inline StatsContainer<float> DataArray::getStats(){ return this->Stats; };
129inline StatsContainer<float>& DataArray::stats(){
130  StatsContainer<float> &rstats = this->Stats;
131  return rstats;
132};
133inline void      DataArray::saveStats(StatsContainer<float> newStats){ this->Stats = newStats;};
134
135
136
137
138/****************************************************************/
139/**
140 * Definition of an data-cube object (3D):
141 *     a DataArray object limited to dim=3
142 */
143
144class Cube : public DataArray
145{
146public:
147  Cube();                 ///< Basic Cube constructor.
148  Cube(long nPix);        ///< Alternative Cube constructor.
149  Cube(long *dimensions); ///< Alternative Cube constructor.
150  virtual ~Cube();        ///< Basic Cube destructor.
151
152  // INLINE functions -- definitions included after class declaration.
153  bool        isBlank(int vox);                               ///< Is the voxel number given by vox a BLANK value?
154  bool        isBlank(long x, long y, long z);                ///< Is the voxel at (x,y,z) a BLANK value?
155  float       getPixValue(long pos);                          ///< Returns the flux value at voxel number pos.
156  float       getPixValue(long x, long y, long z);            ///< Returns the flux value at voxel (x,y,z).
157  short       getDetectMapValue(long pos);                    ///< Returns the value of the detection map at spatial pixel pos.
158  short       getDetectMapValue(long x, long y);              ///< Returns the value of the detection map at spatial pixel (x,y).
159  bool        isRecon();                                      ///< Does the Cube::recon array exist?
160  float       getReconValue(long pos);                        ///< Returns the reconstructed flux value at voxel number pos.
161  float       getReconValue(long x, long y, long z);          ///< Returns the reconstructed flux value at voxel (x,y,z).
162  float       getBaselineValue(long pos);                     ///< Returns the baseline flux value at voxel number pos.
163  float       getBaselineValue(long x, long y, long z);       ///< Returns the baseline flux value at voxel (x,y,z).
164  void        setPixValue(long pos, float f);                 ///< Sets the value of voxel pos.
165  void        setPixValue(long x, long y, long z, float f);   ///< Sets the value of voxel (x,y,z).
166  void        setDetectMapValue(long pos, short f);           ///< Sets the value of the detection map at spatial pixel pos.
167  void        setDetectMapValue(long x, long y, short f);     ///< Sets the value of the detection map at spatial pixel (x,y).
168  void        setReconValue(long pos, float f);               ///< Sets the reconstructed flux value of voxel pos.
169  void        setReconValue(long x, long y, long z, float f); ///< Sets the reconstructed flux value of voxel (x,y,z).
170  void        setReconFlag(bool f);                           ///< Sets the value of the reconExists flag.
171  vector<Col> getLogCols();                                   ///< Return the vector of log file columns
172  void        setLogCols(vector<Col> C);                      ///< Set the vector of log file columns
173  vector<Col> getFullCols();                                  ///< Return the vector of the full output column set.
174  void        setFullCols(vector<Col> C);                     ///< Set the vector of the full output column set.
175
176  // additional functions -- in Cubes/cubes.cc
177  void        initialiseCube(long *dimensions);   ///< Allocate memory correctly, with WCS defining the correct axes.
178  int         getCube();                          ///< Read in a FITS file, with subsection correction.
179  int         getopts(int argc, char ** argv);    ///< Read in command-line options.
180  void        readSavedArrays();                  ///< Read in reconstructed & smoothed arrays from FITS files on disk.
181  void        saveArray(float *input, long size); ///< Save an external array to the Cube's pixel array
182  void        saveRecon(float *input, long size); ///< Save an external array to the Cube's reconstructed array.
183  void        getRecon(float *output);            ///< Save reconstructed array to an external array.
184  void        removeMW();                         ///< Set Milky Way channels to zero.
185  // Statistics for cube
186  void        setCubeStatsOld();                  ///< Calculate the statistics for the Cube (older version).
187  void        setCubeStats();                     ///< Calculate the statistics for the Cube.
188  int         setupFDR();                         ///< Set up thresholds for the False Discovery Rate routine.
189  bool        isDetection(long x, long y, long z);///< A detection test for a given voxel.
190  // Dealing with the detections
191  void        calcObjectWCSparams();              ///< Calculate the WCS parameters for each Cube Detection.
192  void        sortDetections();                   ///< Sort the list of detections.
193  void        updateDetectMap();                  ///< Update the map of detected pixels.
194  void        updateDetectMap(Detection obj);     ///< Update the map of detected pixels for a given Detection.
195  float       enclosedFlux(Detection obj);        ///< Find the flux enclosed by a Detection.
196  void        setupColumns();                     ///< Set up the output column definitions for the Cube and its Detection list.
197  bool        objAtSpatialEdge(Detection obj);    ///< Is the object at the edge of the image?
198  bool        objAtSpectralEdge(Detection obj);   ///< Is the object at an end of the spectrum?
199  void        setObjectFlags();                   ///< Set warning flags for the detections.
200  // Graphical plotting of the cube and the detections.
201  void        plotBlankEdges();                   ///< Draw blank edges of cube.
202  // Dealing with the WCS
203  FitsHeader  getHead();                          ///< Return the FitsHeader.
204  void        setHead(FitsHeader F);              ///< Set the FitsHeader.
205  FitsHeader& header();                           ///< Provides a reference to the FitsHeader.
206  int         wcsToPix(const double *world, double *pix); ///< Convert a point from WCS to Pixel coords.
207  int         wcsToPix(const double *world, double *pix, const int npts); ///< Convert a set of points from WCS to Pixel coords.
208  int         pixToWCS(const double *pix, double *world); ///< Convert a point from Pixel to WCS coords.
209  int         pixToWCS(const double *pix, double *world, const int npts); ///< Convert a set of points from Pixel to WCS coords.
210
211  // FITS-I/O related functions -- not in cubes.cc
212  // in Cubes/getImage.cc
213  int         getCube(string fname);     ///< Function to read in FITS file.
214  // in FitsIO/dataIO.cc
215  int         getFITSdata(string fname); ///< Function to retrieve FITS data array
216  // in Cubes/saveImage.cc
217  void        saveSmoothedCube();        ///< Save Hanning-smoothed array to disk.
218  void        saveReconstructedCube();   ///< Save Reconstructed array to disk.
219  // in Cubes/readRecon.cc
220  int         readReconCube();           ///< Read in reconstructed array from FITS file.
221  // in Cubes/readSmooth.cc 
222  int         readSmoothCube();          ///< Read in Hanning-smoothed array from FITS file.
223
224  // Functions that act on the cube
225  // in Cubes/trimImage.cc
226  void        trimCube();        ///< Remove excess BLANK pixels from spatial edge of cube.
227  void        unTrimCube();      ///< Replace BLANK pixels to spatial edge of cube.
228  // in Cubes/baseline.cc
229  void        removeBaseline();  ///< Removes the baselines from the spectra, and stores in Cube::baseline
230  void        replaceBaseline(); ///< Replace the baselines stored in Cube::baseline
231  // in Cubes/invertCube.cc
232  void        invert();          ///< Multiply all pixel values by -1.
233  void        reInvert();        ///< Undo the inversion, and invert fluxes of all detected objects.
234
235  // Reconstruction, Searching and Merging functions
236// in ATrous/ReconSearch.cc
237  void        ReconSearch();     ///< Front-end to reconstruction & searching functions.
238  void        ReconCube();       ///< Switcher to reconstruction functions
239  void        ReconCube1D();     ///< Performs 1-dimensional a trous reconstruction on the Cube.
240  void        ReconCube2D();     ///< Performs 2-dimensional a trous reconstruction on the Cube.
241  void        ReconCube3D();     ///< Performs 3-dimensional a trous reconstruction on the Cube.
242  // in Cubes/CubicSearch.cc
243  void        CubicSearch();     ///< Front-end to the cubic searching routine.
244  // in Cubes/smoothCube.cc
245  void        SmoothSearch();    ///< Front-end to the smoothing and searching functions.
246  void        SmoothCube();      ///< A function to Hanning-smooth the cube.
247  // in Cubes/Merger.cc
248  void        ObjectMerger();    ///< Merge all objects in the detection list so that only distinct ones are left.
249 
250  // Text outputting of detected objects.
251  //   in Cubes/detectionIO.cc
252  void        outputDetectionsKarma(std::ostream &stream);   ///< Output detections to a Karma annotation file.
253  void        outputDetectionsVOTable(std::ostream &stream); ///< Output detections to a VOTable.
254  void        outputDetectionList();     ///< Output detections to the output file and standard output.
255  void        logDetectionList();        ///< Output detections to the log file.
256  void        logDetection(Detection obj, int counter);   ///< Output a single detection to the log file     
257
258  //  in Cubes/plotting.cc
259  void        plotDetectionMap(string pgDestination);     ///< Plot a spatial map of detections based on number of detected channels.
260  void        plotMomentMap(string pgDestination);        ///< Plot a spatial map of detections based on 0th moment map of each object.     
261  void        plotMomentMap(vector<string> pgDestination);///< Plot a spatial map of detections based on 0th moment map of each object to a number of PGPLOT devices.
262  void        plotWCSaxes();                              ///< Draw WCS axes over a PGPLOT map.
263
264  //  in Cubes/outputSpectra.cc
265  void        outputSpectra();    ///< Print spectra of each detected object.
266  void        plotSpectrum(Detection obj,Plot::SpectralPlot &plot);  ///< Print spectrum of a single object
267  //  in Cubes/drawMomentCutout.cc
268  void        drawMomentCutout(Detection &object); ///< Draw the 0th moment map for a single object.
269  void        drawScale(float xstart,float ystart,float channel); ///< Draw a scale bar indicating angular size of the cutout.
270  void        drawFieldEdge(); ///< Draw a yellow line around the edge of the spatial extent of pixels.
271
272private:
273  float      *recon;            ///< reconstructed array - used when doing a trous reconstruction.
274  bool        reconExists;      ///< flag saying whether there is a reconstruction
275  short      *detectMap;        ///< "moment map" - x,y locations of detected pixels
276  float      *baseline;         ///< array of spectral baseline values.
277                             
278  bool        reconAllocated;   ///< have we allocated memory for the recon array?
279  bool        baselineAllocated;///< have we allocated memory for the baseline array?
280                             
281  FitsHeader  head;             ///< the WCS and other header information.
282  vector<Col> fullCols;         ///< the list of all columns as printed in the results file
283  vector<Col> logCols;          ///< the list of columns as printed in the log file
284
285};
286
287////////////
288//// Cube inline function definitions
289////////////
290
291inline bool  Cube::isBlank(int vox){ return par.isBlank(array[vox]); };
292inline bool  Cube::isBlank(long x, long y, long z){ return par.isBlank(array[z*axisDim[0]*axisDim[1] + y*axisDim[0] + x]); };
293inline float Cube::getPixValue(long pos){ return array[pos]; };
294inline float Cube::getPixValue(long x, long y, long z){ return array[z*axisDim[0]*axisDim[1] + y*axisDim[0] + x]; };
295inline short Cube::getDetectMapValue(long pos){ return detectMap[pos]; };
296inline short Cube::getDetectMapValue(long x, long y){ return detectMap[y*axisDim[0]+x]; };
297inline bool  Cube::isRecon(){ return reconExists; };
298inline float Cube::getReconValue(long pos){ return recon[pos]; };
299inline float Cube::getReconValue(long x, long y, long z){ return recon[z*axisDim[0]*axisDim[1] + y*axisDim[0] + x];};
300inline float Cube::getBaselineValue(long pos){ return baseline[pos]; };
301inline float Cube::getBaselineValue(long x, long y, long z){ return baseline[z*axisDim[0]*axisDim[1] + y*axisDim[0] + x]; };
302  // these next ones should have checks against array overflow...
303inline void  Cube::setPixValue(long pos, float f){array[pos] = f;};
304inline void  Cube::setPixValue(long x, long y, long z, float f){ array[z*axisDim[0]*axisDim[1] + y*axisDim[0] + x] = f;};
305inline void  Cube::setDetectMapValue(long pos, short f){ detectMap[pos] = f;};
306inline void  Cube::setDetectMapValue(long x, long y, short f){ detectMap[y*axisDim[0] + x] = f;};
307inline void  Cube::setReconValue(long pos, float f){recon[pos] = f;};
308inline void  Cube::setReconValue(long x, long y, long z, float f){ recon[z*axisDim[0]*axisDim[1] + y*axisDim[0] + x] = f; };
309inline void  Cube::setReconFlag(bool f){reconExists = f;};
310inline vector<Col> Cube::getLogCols(){return logCols;};
311inline void        Cube::setLogCols(vector<Col> C){logCols=C;};
312inline vector<Col> Cube::getFullCols(){return fullCols;};
313inline void        Cube::setFullCols(vector<Col> C){fullCols=C;};
314inline FitsHeader  Cube::getHead(){ return head; };
315inline void        Cube::setHead(FitsHeader F){};
316inline FitsHeader& Cube::header(){ FitsHeader &h = head; return h; };
317inline int Cube::wcsToPix(const double *world, double *pix)
318{
319  /**
320   *  Use the WCS in the FitsHeader to convert from WCS to pixel coords for
321   *   a single point.
322   *  \param world The world coordinates.
323   *  \param pix The returned pixel coordinates.
324   */
325  return this->head.wcsToPix(world,pix);
326}
327inline int Cube::wcsToPix(const double *world, double *pix, const int npts)
328{
329  /**
330   *  Use the WCS in the FitsHeader to convert from WCS to pixel coords for
331   *   a set of points.
332   *  \param world The world coordinates.
333   *  \param pix The returned pixel coordinates.
334   *  \param npts The number of points being converted.
335   */
336  return this->head.wcsToPix(world,pix,npts);
337}
338inline int Cube::pixToWCS(const double *pix, double *world)
339{
340  /**
341   *  Use the WCS in the FitsHeader to convert from pixel to WCS coords for
342   *   a single point.
343   *  \param pix The pixel coordinates.
344   *  \param world The returned world coordinates.
345   */
346  return this->head.pixToWCS(pix,world);
347}
348inline int Cube::pixToWCS(const double *pix, double *world, const int npts)
349{
350  /**
351   *  Use the WCS in the FitsHeader to convert from pixel to WCS coords for
352   *   a set of points.
353   *  \param pix The pixel coordinates.
354   *  \param world The returned world coordinates.
355   *  \param npts The number of points being converted.
356   */
357  return this->head.pixToWCS(pix,world,npts);
358}
359
360
361//------------------------------------------------------------------------
362
363/****************************************************************/
364/**
365 *  A DataArray object limited to two dimensions, and with some additional
366 *   special functions.
367 *
368 *  It is used primarily for searching a 1- or 2-D array with
369 *   lutz_detect() and spectrumDetect().
370 */
371
372class Image : public DataArray
373{
374public:
375  Image(){numPixels=0; numDim=2;};
376  Image(long nPix);
377  Image(long *dimensions);
378  virtual ~Image(){};
379
380  // Defining the array
381  void      saveArray(float *input, long size);
382  void      extractSpectrum(float *Array, long *dim, long pixel);
383  void      extractImage(float *Array, long *dim, long channel);
384  void      extractSpectrum(Cube &cube, long pixel);
385  void      extractImage(Cube &cube, long channel);
386  // Accessing the data.
387  float     getPixValue(long x, long y){return array[y*axisDim[0] + x];};
388  float     getPixValue(long pos){return array[pos];};
389  // the next few should have checks against array overflow...
390  void      setPixValue(long x, long y, float f){array[y*axisDim[0] + x] = f;};
391  void      setPixValue(long pos, float f){array[pos] = f;};
392  bool      isBlank(int vox){return par.isBlank(array[vox]);};
393  bool      isBlank(long x,long y){return par.isBlank(array[y*axisDim[0]+x]);};
394
395  // Detection-related
396  void      lutz_detect();            // in Detection/lutz_detect.cc
397  void      spectrumDetect();         // in Detection/spectrumDetect.cc
398  int       getMinSize(){return minSize;};
399  void      setMinSize(int i){minSize=i;};
400  // the rest are in Detection/thresholding_functions.cc
401  bool      isDetection(long x, long y){
402    long voxel = y*axisDim[0] + x;
403    if(isBlank(x,y)) return false;
404    else return Stats.isDetection(array[voxel]);
405  }; 
406
407  void      removeMW();
408 
409private:
410  int        minSize;    ///< the minimum number of pixels for a detection to be accepted.
411};
412
413
414/****************************************************************/
415//////////////////////////////////////////////////////
416// Prototypes for functions that use above classes
417//////////////////////////////////////////////////////
418
419void findSources(Image &image);
420void findSources(Image &image, float mean, float sigma);
421
422vector <Detection> searchReconArray(long *dim, float *originalArray,
423                                    float *reconArray, Param &par,
424                                    StatsContainer<float> &stats);
425vector <Detection> search3DArray(long *dim, float *Array, Param &par,
426                                 StatsContainer<float> &stats);
427
428void growObject(Detection &object, Cube &cube);
429
430#endif
Note: See TracBrowser for help on using the repository browser.