source: tags/release-0.9.2/param.cc

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

Changed if statement in << for Params so that reconExists info is printed
out only when the previously reconstructed array is actually being used.
File size: 14.3 KB
Line 
1#include <iostream>
2#include <iomanip>
3#include <fstream>
4#include <sstream>
5#include <string>
6#include <math.h>
7#include <param.hh>
8
9using std::setw;
10using std::endl;
11
12/****************************************************************/
13///////////////////////////////////////////////////
14//// Accessor Functions for Parameter class:
15///////////////////////////////////////////////////
16Param::Param(){
17  /**
18   * Param()
19   *  Default intial values for the parameters.
20   * imageFile has no default value!
21   */
22  // Input files
23  this->imageFile       = "";
24  this->flagSubsection  = false;
25  this->subsection      = "[*,*,*]";
26  this->flagReconExists = false;
27  this->reconFile       = "";
28  // Output files
29  this->flagLog         = true;
30  this->logFile         = "duchamp-Logfile.txt";
31  this->outFile         = "duchamp-Results.txt";
32  this->spectraFile     = "duchamp-Spectra.ps";
33  this->flagOutputRecon = false;
34  this->flagOutputResid = false;
35  this->flagVOT         = false;
36  this->votFile         = "duchamp-Results.xml";
37  this->flagKarma       = false;
38  this->karmaFile       = "duchamp-Results.ann";
39  this->flagMaps        = true;
40  this->detectionMap    = "duchamp-DetectionMap.ps";
41  this->momentMap       = "duchamp-MomentMap.ps";
42  // Cube related parameters
43  this->flagBlankPix    = true;
44  this->blankPixValue   = -8.00061;
45  this->blankKeyword    = 1;
46  this->bscaleKeyword   = -8.00061;
47  this->bzeroKeyword    = 0.;
48  this->nanAsBlank      = false;
49  this->flagMW          = false;
50  this->maxMW           = 112;
51  this->minMW           = 75;
52  this->numPixBeam      = 0.;
53  // Trim-related         
54  this->flagTrimmed     = false;
55  this->borderLeft      = 0;
56  this->borderRight     = 0;
57  this->borderBottom    = 0;
58  this->borderTop       = 0;
59  // Subsection offsets
60  this->xSubOffset      = 0;
61  this->ySubOffset      = 0;
62  this->zSubOffset      = 0;
63  // Baseline related
64  this->flagBaseline    = false;
65  // Detection-related   
66  this->flagNegative    = false;
67  // Object growth       
68  this->flagGrowth      = false;
69  this->growthCut       = 2.;
70  // FDR analysis         
71  this->flagFDR         = false;
72  this->alphaFDR        = 0.01;
73  // Other detection     
74  this->snrCut          = 3.;
75  // A trous reconstruction parameters
76  this->flagATrous      = true;
77  this->scaleMin        = 1;
78  this->snrRecon        = 4.;
79  this->filterCode      = 2;
80  // Volume-merging parameters
81  this->flagAdjacent    = true;
82  this->threshSpatial   = 3.;
83  this->threshVelocity  = 7.;
84  this->minChannels     = 3;
85  this->minPix          = 2;
86  // Input-Output related
87  this->spectralMethod  = "peak";
88  this->borders         = true;
89  this->verbose         = true;
90};
91
92/****************************************************************/
93///////////////////////////////////////////////////
94//// Other Functions using the  Parameter class:
95///////////////////////////////////////////////////
96
97string makelower( string s )
98{
99  // "borrowed" from Matt Howlett's 'fred'
100  string out = "";
101  for( int i=0; i<s.size(); ++i ) {
102    out += tolower(s[i]);
103  }
104  return out;
105}
106
107void Param::readParams(string &paramfile)
108{
109
110  std::ifstream fin(paramfile.c_str());
111  string line;
112  string sval;
113  float fval;
114  int ival;
115  bool bval;
116  while( !std::getline(fin,line,'\n').eof()){
117
118    if(line[0]!='#'){
119      std::stringstream ss;
120      ss.str(line);
121      string arg;
122      ss >> arg;     
123      if(makelower(arg)=="imagefile"){       ss >> sval; this->imageFile = sval;}
124      if(makelower(arg)=="flagsubsection"){  ss >> bval; this->flagSubsection = bval; }
125      if(makelower(arg)=="subsection"){      ss >> sval; this->subsection = sval; }
126      if(makelower(arg)=="flagreconexists"){ ss >> bval; this->flagReconExists = bval; }
127      if(makelower(arg)=="reconfile"){       ss >> sval; this->reconFile = sval; }
128
129      if(makelower(arg)=="flaglog"){         ss >> bval; this->flagLog = bval; }
130      if(makelower(arg)=="logfile"){         ss >> sval; this->logFile = sval; }
131      if(makelower(arg)=="outfile"){         ss >> sval; this->outFile = sval; }
132      if(makelower(arg)=="spectrafile"){     ss >> sval; this->spectraFile = sval; }
133      if(makelower(arg)=="flagoutputrecon"){ ss >> bval; this->flagOutputRecon = bval; }
134      if(makelower(arg)=="flagoutputresid"){ ss >> bval; this->flagOutputResid = bval; }
135      if(makelower(arg)=="flagvot"){         ss >> bval; this->flagVOT = bval; }
136      if(makelower(arg)=="votfile"){         ss >> sval; this->votFile = sval; }
137      if(makelower(arg)=="flagkarma"){       ss >> bval; this->flagKarma = bval; }
138      if(makelower(arg)=="karmafile"){       ss >> sval; this->karmaFile = sval; }
139      if(makelower(arg)=="flagmaps"){        ss >> bval; this->flagMaps = bval; }
140      if(makelower(arg)=="detectionmap"){    ss >> sval; this->detectionMap = sval; }
141      if(makelower(arg)=="momentmap"){       ss >> sval; this->momentMap = sval; }
142
143      if(makelower(arg)=="flagnegative"){    ss >> bval; this->flagNegative = bval;}
144      if(makelower(arg)=="flagblankpix"){    ss >> bval; this->flagBlankPix = bval; }
145      if(makelower(arg)=="blankpixvalue"){   ss >> fval; this->blankPixValue = fval; }
146      if(makelower(arg)=="flagmw"){          ss >> bval; this->flagMW = bval; }
147      if(makelower(arg)=="maxmw"){           ss >> ival; this->maxMW = ival; }
148      if(makelower(arg)=="minmw"){           ss >> ival; this->minMW = ival; }
149
150      if(makelower(arg)=="flagbaseline"){    ss >> bval; this->flagBaseline = bval; }
151      if(makelower(arg)=="minpix"){          ss >> ival; this->minPix = ival; }
152      if(makelower(arg)=="flaggrowth"){      ss >> bval; this->flagGrowth = bval; }
153      if(makelower(arg)=="growthcut"){       ss >> fval; this->growthCut = fval; }
154
155      if(makelower(arg)=="flagfdr"){         ss >> bval; this->flagFDR = bval; }
156      if(makelower(arg)=="alphafdr"){        ss >> fval; this->alphaFDR = fval; }
157      if(makelower(arg)=="snrcut"){          ss >> fval; this->snrCut = fval; }
158      if(makelower(arg)=="flagatrous"){      ss >> bval; this->flagATrous = bval; }
159      if(makelower(arg)=="scalemin"){        ss >> ival; this->scaleMin = ival; }
160      if(makelower(arg)=="snrrecon"){        ss >> fval; this->snrRecon = fval; }
161      if(makelower(arg)=="filtercode"){      ss >> ival; this->filterCode = ival; }
162
163      if(makelower(arg)=="flagadjacent"){    ss >> bval; this->flagAdjacent = bval; }
164      if(makelower(arg)=="threshspatial"){   ss >> fval; this->threshSpatial = fval; }
165      if(makelower(arg)=="threshvelocity"){  ss >> fval; this->threshVelocity = fval; }
166      if(makelower(arg)=="minchannels"){     ss >> ival; this->minChannels = ival; }
167
168      if(makelower(arg)=="spectralMethod"){  ss >> sval; this->spectralMethod = makelower(sval); }
169      if(makelower(arg)=="drawborders"){     ss >> bval; this->borders = bval; }
170      if(makelower(arg)=="verbose"){         ss >> bval; this->verbose = bval; }
171    }
172  }
173}
174
175
176std::ostream& operator<< ( std::ostream& theStream, Param& par)
177{
178  theStream.setf(std::ios::boolalpha);
179  theStream.setf(std::ios::left);
180  theStream  <<"---- Parameters ----"<<endl;
181  if(par.getFlagSubsection())
182    theStream<<setw(40)<<"Image to be analysed"                 <<"= "<<(par.getImageFile()+
183                                                                         par.getSubsection())   <<endl;
184  else
185    theStream<<setw(40)<<"Image to be analysed"                 <<"= "<<par.getImageFile()      <<endl;
186  if(par.getFlagReconExists() && par.getFlagATrous()){
187    theStream<<setw(40)<<"Reconstructed array exists?"          <<"= "<<par.getFlagReconExists()<<endl;
188    theStream<<setw(40)<<"FITS file containing reconstruction"  <<"= "<<par.getReconFile()      <<endl;
189  }
190  theStream  <<setw(40)<<"Intermediate Logfile"                 <<"= "<<par.getLogFile()        <<endl;
191  theStream  <<setw(40)<<"Final Results file"                   <<"= "<<par.getOutFile()        <<endl;
192  theStream  <<setw(40)<<"Spectrum file"                        <<"= "<<par.getSpectraFile()    <<endl;
193  if(par.getFlagVOT()){
194    theStream<<setw(40)<<"VOTable file"                         <<"= "<<par.getVOTFile()        <<endl;
195  }
196  if(par.getFlagKarma()){
197    theStream<<setw(40)<<"Karma annotation file"                <<"= "<<par.getKarmaFile()      <<endl;
198  }
199  if(par.getFlagMaps()){
200    theStream<<setw(40)<<"0th Moment Map"                       <<"= "<<par.getMomentMap()      <<endl;
201    theStream<<setw(40)<<"Detection Map"                        <<"= "<<par.getDetectionMap()   <<endl;
202  }
203  if(par.getFlagATrous()){                             
204    theStream<<setw(40)<<"Saving reconstructed cube?"           <<"= "<<par.getFlagOutputRecon()<<endl;
205    theStream<<setw(40)<<"Saving residuals from reconstruction?"<<"= "<<par.getFlagOutputResid()<<endl;
206  }                                                   
207  theStream  <<"------"<<endl;
208  theStream  <<setw(40)<<"Searching for Negative features?"     <<"= "<<par.getFlagNegative()   <<endl;
209  theStream  <<setw(40)<<"Fixing Blank Pixels?"                 <<"= "<<par.getFlagBlankPix()   <<endl;
210  if(par.getFlagBlankPix()){
211    theStream<<setw(40)<<"Blank Pixel Value"                    <<"= "<<par.getBlankPixVal()    <<endl;
212  }
213  theStream  <<setw(40)<<"Removing Milky Way channels?"         <<"= "<<par.getFlagMW()         <<endl;
214  if(par.getFlagMW()){
215    theStream<<setw(40)<<"Milky Way Channels"                   <<"= "<<par.getMinMW()
216                       <<"-"                                          <<par.getMaxMW()          <<endl;
217  }
218  theStream  <<setw(40)<<"Beam Size (pixels)"                   <<"= "<<par.getBeamSize()       <<endl;
219  theStream  <<setw(40)<<"Removing baselines before search?"    <<"= "<<par.getFlagBaseline()   <<endl;
220  theStream  <<setw(40)<<"Minimum # Pixels in a detection"      <<"= "<<par.getMinPix()         <<endl;
221  theStream  <<setw(40)<<"Growing objects after detection?"     <<"= "<<par.getFlagGrowth()     <<endl;
222  if(par.getFlagGrowth()) {                           
223    theStream<<setw(40)<<"SNR Threshold for growth"             <<"= "<<par.getGrowthCut()      <<endl;
224  }
225  theStream  <<setw(40)<<"Using A Trous reconstruction?"        <<"= "<<par.getFlagATrous()     <<endl;
226  if(par.getFlagATrous()){                             
227    theStream<<setw(40)<<"Minimum scale in reconstruction"      <<"= "<<par.getMinScale()       <<endl;
228    theStream<<setw(40)<<"SNR Threshold within reconstruction"  <<"= "<<par.getAtrousCut()      <<endl;
229    theStream<<setw(40)<<"Filter being used for reconstruction" <<"= "<<par.getFilterName()     <<endl;
230  }                                                   
231  theStream  <<setw(40)<<"Using FDR analysis?"                  <<"= "<<par.getFlagFDR()        <<endl;
232  if(par.getFlagFDR()){                               
233    theStream<<setw(40)<<"Alpha value for FDR analysis"         <<"= "<<par.getAlpha()          <<endl;
234  }                                                   
235  else {
236    theStream<<setw(40)<<"SNR Threshold"                        <<"= "<<par.getCut()            <<endl;
237  }
238  theStream  <<setw(40)<<"Using Adjacent-pixel criterion?"      <<"= "<<par.getFlagAdjacent()   <<endl;
239  if(!par.getFlagAdjacent()){
240    theStream<<setw(40)<<"Max. spatial separation for merging"  <<"= "<<par.getThreshS()        <<endl;
241  }
242  theStream  <<setw(40)<<"Max. velocity separation for merging" <<"= "<<par.getThreshV()        <<endl;
243  theStream  <<setw(40)<<"Min. # channels for merging"          <<"= "<<par.getMinChannels()    <<endl;
244  theStream  <<setw(40)<<"Method of spectral plotting"          <<"= "<<par.getSpectralMethod() <<endl;
245  theStream  <<"--------------------"<<endl;
246  theStream.unsetf(std::ios::left);
247  theStream.unsetf(std::ios::boolalpha);
248}
249
250void Param::parseSubsection()
251{
252  std::stringstream ss;
253  ss.str(this->subsection);
254  bool removeStep = false;
255  string x,y,z;
256  getline(ss,x,'[');
257  getline(ss,x,',');
258  getline(ss,y,',');
259  getline(ss,z,']');
260  char *end;
261  if(x!="*") {
262    int x1,x2,dx;
263    int a = x.find(':');  // first occurence of ':' in x-subsection string
264    int b = x.find(':',a+1); // location of second ':' -- will be -1 if there is no second occurence.
265    x1 = atoi( x.substr(0,a).c_str() ); // minimum x in subsection
266    this->xSubOffset = x1 - 1;
267    if(b>0){ // there is a dx component
268      x2 = atoi( x.substr(a+1,b-a-1).c_str() ); // maximum x in subsection
269      dx = atoi( x.substr(b+1,x.size()).c_str() ); // delta x in subsection
270      x = x.substr(0,b); // rewrite x without the delta-x part.
271      removeStep = true;
272    }
273  }
274  if(y!="*") {
275    int y1,y2,dy;
276    int a = y.find(':');  // first occurence of ':' in y-subsection string
277    int b = y.find(':',a+1); // location of second ':' -- will be -1 if there is no second occurence.
278    y1 = atoi( y.substr(0,a).c_str() ); // minimum y in subsection
279    this->ySubOffset = y1 - 1;
280    if(b>0){ // there is a dy component
281      y2 = atoi( y.substr(a+1,b-a-1).c_str() ); // maximum y in subsection
282      dy = atoi( y.substr(b+1,y.size()).c_str() ); // delta y in subsection
283      y = y.substr(0,b); // rewrite y without the delta-y part.
284      removeStep = true;
285    }
286  }
287  if(z!="*") {
288    int z1,z2,dz;
289    int a = z.find(':');  // first occurence of ':' in z-subsection string
290    int b = z.find(':',a+1); // location of second ':' -- will be -1 if there is no second occurence.
291    z1 = atoi( z.substr(0,a).c_str() ); // minimum z in subsection
292    this->zSubOffset = z1 - 1;
293    if(b>0){ // there is a dz component
294      z2 = atoi( z.substr(a+1,b-a-1).c_str() ); // maximum z in subsection
295      dz = atoi( z.substr(b+1,z.size()).c_str() ); // delta z in subsection
296      z = z.substr(0,b); // rewrite z without the delta-z part.
297      removeStep = true;
298    }
299  }
300
301  if(removeStep){
302    std::cerr << "\a\tThe subsection given is " << this->subsection <<".\n";
303    std::cerr << "\tDuchamp is currently unable to deal with pixel steps in the subsection.\n";
304    std::cerr << "\tThese have been ignored, and so the subection used is ";   
305    this->subsection = "[" + x + "," + y + "," + z + "]";  // rewrite subsection without any step sizes.
306    std::cerr << this->subsection << std::endl;
307  }
308 
309//   if(x!="*") this->xSubOffset = strtol(x.c_str(),&end,10) - 1;   
310//   if(y!="*") this->ySubOffset = strtol(y.c_str(),&end,10) - 1;   
311//   if(z!="*") this->zSubOffset = strtol(z.c_str(),&end,10) - 1;   
312}
313
314bool Param::isBlank(float &val)
315{
316//   std::cerr << "val = " << val << " " << int((val-this->bzeroKeyword)/this->bscaleKeyword)
317//          <<" " <<this->blankKeyword << std::endl;
318  if(this->flagBlankPix){
319    if(this->nanAsBlank) return bool(isnan(val));
320    else//  return (val == this->blankPixValue);
321      return ( this->blankKeyword == int((val-this->bzeroKeyword)/this->bscaleKeyword) );
322  }
323  else return false;
324}
Note: See TracBrowser for help on using the repository browser.