source: trunk/src/param.hh @ 189

Last change on this file since 189 was 189, checked in by Matthew Whiting, 18 years ago

A suite of changes aimed at calculating a cube-wide threshold, rather than individual thresholds for each channel map & spectrum. Summary of changes:

  • New Cube::setCubeStats() function, that stores the cube's stats in Cube member variables (mean, median, stddev, madfm).
  • New Parameter threshold, which can be entered by the user.
  • Make use of the cube statistics in the search functions, so that the Image class receives the correct information. Some tidying up of the CubicSearch? and ReconSearch? functions as well.
  • Changing the thresholding functions to be able to use the new threshold parameter.

Warning though -- the FDR setup doesn't quite work at this stage!

File size: 16.8 KB
Line 
1#ifndef PARAM_H
2#define PARAM_H
3
4#include <iostream>
5#include <string>
6#include <vector>
7#include <math.h>
8#include <wcs.h>
9
10using std::string;
11using std::vector;
12
13/**
14 * Param class.
15 *   Used for storing parameters used by all functions.
16 */
17
18class FitsHeader; // foreshadow this so that Param knows it exists
19
20class Param
21{
22public:
23  Param();
24  virtual ~Param(){};
25  Param& operator= (const Param& p);
26  int    verifySubsection();              // in FitsIO/subsection.cc
27  void   setOffsets(wcsprm *wcs);         // in FitsIO/subsection.cc
28  int    readParams(string paramfile);    // in param.cc
29  void   copyHeaderInfo(FitsHeader &head);// in param.cc
30  bool   isBlank(float &value);           // in param.cc
31  bool   isInMW(int z){return ( flagMW && (z>=minMW) && (z<=maxMW) );};
32  string outputReconFile();             // in param.cc
33  string outputResidFile();             // in param.cc
34  friend std::ostream& operator<< ( std::ostream& theStream, Param& par);
35  friend class Image;
36 
37  //
38  string getImageFile(){return imageFile;};
39  void   setImageFile(string fname){imageFile = fname;};
40  string getFullImageFile(){
41    if(flagSubsection) return imageFile+subsection;
42    else return imageFile;
43  };
44  bool   getFlagSubsection(){return flagSubsection;};
45  void   setFlagSubsection(bool flag){flagSubsection=flag;};
46  string getSubsection(){return subsection;};
47  void   setSubsection(string range){subsection = range;};
48  bool   getFlagReconExists(){return flagReconExists;};
49  void   setFlagReconExists(bool flag){flagReconExists=flag;};
50  string getReconFile(){return reconFile;};
51  void   setReconFile(string file){reconFile = file;};
52  //
53  bool   getFlagLog(){return flagLog;};
54  void   setFlagLog(bool flag){flagLog=flag;};
55  string getLogFile(){return logFile;};
56  void   setLogFile(string fname){logFile = fname;};
57  string getOutFile(){return outFile;};
58  void   setOutFile(string fname){outFile = fname;};
59  string getSpectraFile(){return spectraFile;};
60  void   setSpectraFile(string fname){spectraFile = fname;};
61  bool   getFlagOutputRecon(){return flagOutputRecon;};
62  void   setFlagOutputRecon(bool flag){flagOutputRecon=flag;};
63  bool   getFlagOutputResid(){return flagOutputResid;};
64  void   setFlagOutputResid(bool flag){flagOutputResid=flag;};
65  bool   getFlagVOT(){return flagVOT;};
66  void   setFlagVOT(bool flag){flagVOT=flag;};
67  string getVOTFile(){return votFile;};
68  void   setVOTFile(string fname){votFile = fname;};
69  bool   getFlagKarma(){return flagKarma;};
70  void   setFlagKarma(bool flag){flagKarma=flag;};
71  string getKarmaFile(){return karmaFile;};
72  void   setKarmaFile(string fname){karmaFile = fname;};
73  bool   getFlagMaps(){return flagMaps;};
74  void   setFlagMaps(bool flag){flagMaps=flag;};
75  string getDetectionMap(){return detectionMap;};
76  void   setDetectionMap(string fname){detectionMap = fname;};
77  string getMomentMap(){return momentMap;};
78  void   setMomentMap(string fname){momentMap = fname;};
79  //
80  bool   getFlagNegative(){return flagNegative;};
81  void   setFlagNegative(bool flag){flagNegative=flag;};
82  bool   getFlagBlankPix(){return flagBlankPix;};
83  void   setFlagBlankPix(bool flag){flagBlankPix=flag;};
84  float  getBlankPixVal(){return blankPixValue;};
85  void   setBlankPixVal(float v){blankPixValue=v;};
86  int    getBlankKeyword(){return blankKeyword;};
87  void   setBlankKeyword(int v){blankKeyword=v;};
88  float  getBscaleKeyword(){return bscaleKeyword;};
89  void   setBscaleKeyword(float v){bscaleKeyword=v;};
90  float  getBzeroKeyword(){return bzeroKeyword;};
91  void   setBzeroKeyword(float v){bzeroKeyword=v;};
92  bool   getFlagUsingBlank(){return flagUsingBlank;};
93  void   setFlagUsingBlank(bool b){flagUsingBlank=b;};
94  bool   getFlagMW(){return flagMW;};
95  bool   setFlagMW(bool flag){flagMW=flag;};
96  int    getMaxMW(){return maxMW;};
97  void   setMaxMW(int m){maxMW=m;};
98  int    getMinMW(){return minMW;};
99  void   setMinMW(int m){minMW=m;};
100  void   setBeamSize(float s){numPixBeam = s;};
101  float  getBeamSize(){return numPixBeam;};
102  bool   getFlagUsingBeam(){return flagUsingBeam;};
103  void   setFlagUsingBeam(bool b){flagUsingBeam=b;};
104  //
105  bool   getFlagCubeTrimmed(){return flagTrimmed;};
106  void   setFlagCubeTrimmed(bool flag){flagTrimmed = flag;};
107  long   getBorderLeft(){return borderLeft;};
108  void   setBorderLeft(long b){borderLeft = b;};
109  long   getBorderRight(){return borderRight;};
110  void   setBorderRight(long b){borderRight = b;};
111  long   getBorderBottom(){return borderBottom;};
112  void   setBorderBottom(long b){borderBottom = b;};
113  long   getBorderTop(){return borderTop;};
114  void   setBorderTop(long b){borderTop = b;};
115  //
116  long   getXOffset(){return xSubOffset;};
117  void   setXOffset(long o){xSubOffset = o;};
118  long   getYOffset(){return ySubOffset;};
119  void   setYOffset(long o){ySubOffset = o;};
120  long   getZOffset(){return zSubOffset;};
121  void   setZOffset(long o){zSubOffset = o;};
122  //
123  int    getMinPix(){return minPix;};
124  void   setMinPix(int m){minPix=m;};
125  //     
126  bool   getFlagGrowth(){return flagGrowth;};
127  void   setFlagGrowth(bool flag){flagGrowth=flag;};
128  float  getGrowthCut(){return growthCut;};
129  void   setGrowthCut(float c){growthCut=c;};
130  //     
131  bool   getFlagFDR(){return flagFDR;};
132  void   setFlagFDR(bool flag){flagFDR=flag;};
133  float  getAlpha(){return alphaFDR;};
134  void   setAlpha(float a){alphaFDR=a;};
135  //
136  bool   getFlagBaseline(){return flagBaseline;};
137  void   setFlagBaseline(bool flag){flagBaseline = flag;};
138  //
139  float  getCut(){return snrCut;};
140  void   setCut(float c){snrCut=c;};
141  float  getThreshold(){return threshold;};
142  void   setThreshold(float f){threshold=f;};
143  //     
144  bool   getFlagATrous(){return flagATrous;};
145  void   setFlagATrous(bool flag){flagATrous=flag;};
146  int    getReconDim(){return reconDim;};
147  void   setReconDim(int i){reconDim=i;};
148  int    getMinScale(){return scaleMin;};
149  void   setMinScale(int s){scaleMin=s;};
150  float  getAtrousCut(){return snrRecon;};
151  void   setAtrousCut(float c){snrRecon=c;};
152  int    getFilterCode(){return filterCode;};
153  void   setFilterCode(int c){filterCode=c;};
154  string getFilterName(){return filterName;};
155  void   setFilterName(string s){filterName=s;};
156  //     
157  bool   getFlagAdjacent(){return flagAdjacent;};
158  void   setFlagAdjacent(bool flag){flagAdjacent=flag;};
159  float  getThreshS(){return threshSpatial;};
160  void   setThreshS(float t){threshSpatial=t;};
161  float  getThreshV(){return threshVelocity;};
162  void   setThreshV(float t){threshVelocity=t;};
163  int    getMinChannels(){return minChannels;};
164  void   setMinChannels(int n){minChannels=n;};
165  //
166  string getSpectralMethod(){return spectralMethod;};
167  void   setSpectralMethod(string s){spectralMethod=s;};
168  string getSpectralUnits(){return spectralUnits;};
169  void   setSpectralUnits(string s){spectralUnits=s;};
170  bool   drawBorders(){return borders;};
171  void   setDrawBorders(bool f){borders=f;};
172  bool   drawBlankEdge(){return blankEdge;};
173  void   setDrawBlankEdge(bool f){blankEdge=f;};
174  bool   isVerbose(){return verbose;};
175  void   setVerbosity(bool f){verbose=f;};
176 
177private:
178  // Input files
179  string imageFile;       // The image to be analysed.
180  bool   flagSubsection;  // Whether we just want a subsection of the image
181  string subsection;      // The subsection requested, taking the form
182                          //  [x1:x2,y1:y2,z1:z2]
183                          //  If you want the full range of one index, use *
184  bool   flagReconExists; // The reconstructed array is in a FITS file on disk.
185  string reconFile;       // The FITS file containing the reconstructed array.
186
187  // Output files
188  bool   flagLog;         // Should we do the intermediate logging?
189  string logFile;         // Where the intermediate logging goes.
190  string outFile;         // Where the final results get put.
191  string spectraFile;     // Where the spectra are displayed
192  bool   flagOutputRecon; // Should the reconstructed cube be written?
193  bool   flagOutputResid; // Should the reconstructed cube be written?
194  bool   flagVOT;         // Should we save results in VOTable format?
195  string votFile;         // Where the VOTable goes.
196  bool   flagKarma;       // Should we save results in Karma annotation format?
197  string karmaFile;       // Where the Karma annotation file goes.
198  bool   flagMaps;        // Should we produce detection and moment maps
199                          //  in postscript form?
200  string detectionMap;    // The name of the detection map (ps file).
201  string momentMap;       // The name of the 0th moment map (ps file).
202
203  // Cube related parameters
204  bool   flagNegative;    // Are we going to search for negative features?
205  bool   flagBlankPix;    // A flag that indicates whether there are pixels
206                          //   defined as BLANK and whether we need to remove
207                          //   & ignore them in processing.
208  float  blankPixValue;   // Pixel value that is considered BLANK.
209  int    blankKeyword;    // The FITS header keyword BLANK.
210  float  bscaleKeyword;   // The FITS header keyword BSCALE.
211  float  bzeroKeyword;    // The FITS header keyword BZERO.
212  bool   flagUsingBlank;  // If true, we are using the blankPixValue keyword,
213                          // otherwise we use the value in the FITS header.
214  bool   flagMW;          // A flag that indicates whether to ignore the
215                          //  Milky Way channels.
216  int    maxMW;           // Last  Galactic velocity plane for HIPASS cubes
217  int    minMW;           // First Galactic velocity plane for HIPASS cubes
218  float  numPixBeam;      // Size (area) of the beam in pixels.
219  bool   flagUsingBeam;   // If true, we are using the numPixBeam parameter,
220                          // otherwise we use the value in the FITS header.
221  // Trim-related         
222  bool   flagTrimmed;     // Has the cube been trimmed of excess BLANKs
223                          //  around the edge?
224  long   borderLeft;      // The number of BLANK pixels trimmed from the
225                          //   Left of the cube;
226  long   borderRight;     // The number trimmed from the Right of the cube;
227  long   borderBottom;    // The number trimmed from the Bottom of the cube;
228  long   borderTop;       // The number trimmed from the Top of the cube;
229  // Subsection offsets
230  long  *offsets;         // The array of offsets for each FITS axis.
231  long   xSubOffset;      // The offset in the x-direction from the subsection
232  long   ySubOffset;      // The offset in the y-direction from the subsection
233  long   zSubOffset;      // The offset in the z-direction from the subsection
234  // Baseline related;
235  bool   flagBaseline;    // Whether to do baseline subtraction before
236                          //  reconstruction and/or searching.
237  // Detection-related   
238  int    minPix;          // Minimum number of pixels for a detected object
239                          //   to be counted
240  // Object growth       
241  bool   flagGrowth;      // Are we growing objects once they are found?
242  float  growthCut;       // The SNR that we are growing objects down to.
243  // FDR analysis         
244  bool   flagFDR;         // Should the FDR method be used?
245  float  alphaFDR;        // Alpha value for FDR detection algorithm
246  // Basic detection     
247  float  snrCut;          // How many sigma above mean is a detection
248                          //   when sigma-clipping
249  float  threshold;       // What the threshold is (when sigma-clipping).
250  // A trous reconstruction parameters
251  bool   flagATrous;      // Are we using the a trous reconstruction?
252  int    reconDim;        // How many dimensions to use for the reconstruction?
253  int    scaleMin;        // Min scale used in a trous reconstruction
254  float  snrRecon;        // SNR cutoff used in a trous reconstruction
255                          //   (only wavelet coefficients that survive this
256                          //    threshold are kept)
257  int    filterCode;      // The code number for the filter to be used
258                          //  (saves having to parse names)
259  string filterName;      // The code number converted into a name,
260                          //  for outputting purposes.
261
262  // Volume-merging parameters
263  bool   flagAdjacent;    // Whether to use the adjacent criterion for
264                          //    judging if objects are to be merged.
265  float  threshSpatial;   // Maximum spatial separation between objects
266  float  threshVelocity;  // Maximum channels separation between objects
267  int    minChannels;     // Minimum no. of channels to make an object
268  // Input-Output related
269  string spectralMethod;  // A string indicating choice of spectral plotting
270                          //  method: choices are "peak" (default) or "sum"
271  string spectralUnits;   // A string indicating what units the spectral
272                          //  axis should be quoted in.
273  bool   borders;         // Whether to draw a border around the individual
274                          //  pixels of a detection in the spectral display
275  bool   blankEdge;       // Whether to draw a border around the BLANK pixel
276                          //  region in the moment maps and cutout images
277  bool   verbose;         // Whether to use maximum verbosity -- use progress
278                          //  indicators in the reconstruction & merging steps.
279
280};
281
282class FitsHeader
283{
284  /**
285   *  FitsHeader Class
286   *
287   *   Stores information from a FITS header, including WCS information
288   *    in the form of a wcsprm struct, as well as various keywords.
289   */
290
291public:
292  FitsHeader();
293  ~FitsHeader(){};
294  FitsHeader(const FitsHeader& h);
295  FitsHeader& operator= (const FitsHeader& h);
296
297  wcsprm *getWCS();             // in param.cc
298  void    setWCS(wcsprm *w);    // in param.cc
299  bool    isWCS(){return wcsIsGood;};
300  int     getNWCS(){return nwcs;};
301  void    setNWCS(int i){nwcs=i;};
302  int     readHeaderInfo(string fname, Param &par);
303  int     defineWCS(string fname, Param &par);
304  int     getBUNIT(string fname);
305  int     getBLANKinfo(string fname, Param &par);
306  int     getBeamInfo(string fname, Param &par);
307  string  getSpectralUnits(){return spectralUnits;};
308  void    setSpectralUnits(string s){spectralUnits=s;};
309  string  getSpectralDescription(){return spectralDescription;};
310  void    setSpectralDescription(string s){spectralDescription=s;};
311  string  getFluxUnits(){return fluxUnits;};
312  void    setFluxUnits(string s){fluxUnits=s;};
313  string  getIntFluxUnits(){return intFluxUnits;};
314  void    setIntFluxUnits(string s){intFluxUnits=s;};
315  float   getBeamSize(){return beamSize;};
316  void    setBeamSize(float f){beamSize=f;};
317  float   getBmajKeyword(){return bmajKeyword;};
318  void    setBmajKeyword(float f){bmajKeyword=f;};
319  float   getBminKeyword(){return bminKeyword;};
320  void    setBminKeyword(float f){bminKeyword=f;};
321  int     getBlankKeyword(){return blankKeyword;};
322  void    setBlankKeyword(int f){blankKeyword=f;};
323  float   getBzeroKeyword(){return bzeroKeyword;};
324  void    setBzeroKeyword(float f){bzeroKeyword=f;};
325  float   getBscaleKeyword(){return bscaleKeyword;};
326  void    setBscaleKeyword(float f){bscaleKeyword=f;};
327  float   getAvPixScale(){
328    return sqrt( fabs ( (wcs->pc[0]*wcs->cdelt[0])*
329                        (wcs->pc[wcs->naxis+1]*wcs->cdelt[1])));
330  };
331
332  // front ends to WCS functions
333  int     wcsToPix(const double *world, double *pix);
334  int     pixToWCS(const double *pix, double *world);
335  int     wcsToPix(const double *world, double *pix, const int npts);
336  int     pixToWCS(const double *pix, double *world, const int npts);
337  double  pixToVel(double &x, double &y, double &z);
338  double *pixToVel(double &x, double &y, double *zarray, int size);
339  double  specToVel(const double &z);
340  double  velToSpec(const float &vel);
341  string  getIAUName(double ra, double dec);
342
343  void    fixUnits(Param &par);
344 
345private:
346  wcsprm *wcs;                  // The WCS parameters for the cube in a
347                                //  struct from the wcslib library.
348  int     nwcs;                 // The number of WCS parameters
349  bool    wcsIsGood;            // A flag indicating whether there is a
350                                //  valid WCS present.
351  string  spectralUnits;        // The units of the spectral dimension
352  string  spectralDescription;  // The description of the spectral dimension
353                                //   (Frequency, Velocity, ...)
354  string  fluxUnits;            // The units of pixel flux (from header)
355  string  intFluxUnits;         // The units of pixel flux (from header)
356  float   beamSize;             // The calculated beam size in pixels.
357  float   bmajKeyword;          // The FITS header keyword BMAJ.
358  float   bminKeyword;          // The FITS header keyword BMIN.
359  int     blankKeyword;         // The FITS header keyword BLANK.
360  float   bzeroKeyword;         // The FITS header keyword BZERO.
361  float   bscaleKeyword;        // The FITS header keyword BSCALE.
362  double  scale;                // scale param for converting spectral coords
363  double  offset;               // offset param for converting spectral coords
364  double  power;                // power param for converting spectral coords
365};
366
367string makelower( string s );
368
369#endif
Note: See TracBrowser for help on using the repository browser.