1 | #include <iostream> |
---|
2 | #include <iomanip> |
---|
3 | #include <fstream> |
---|
4 | #include <sstream> |
---|
5 | #include <string> |
---|
6 | #include <math.h> |
---|
7 | #include <param.hh> |
---|
8 | |
---|
9 | using std::setw; |
---|
10 | using std::endl; |
---|
11 | |
---|
12 | /****************************************************************/ |
---|
13 | /////////////////////////////////////////////////// |
---|
14 | //// Accessor Functions for Parameter class: |
---|
15 | /////////////////////////////////////////////////// |
---|
16 | Param::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 | |
---|
97 | string 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 | |
---|
107 | void Param::readParams(string ¶mfile) |
---|
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 | |
---|
176 | std::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 | |
---|
250 | void 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 | |
---|
314 | bool 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 | } |
---|