[297] | 1 | //#--------------------------------------------------------------------------- |
---|
| 2 | //# SDLineFinder.cc: A class for automated spectral line search |
---|
| 3 | //#-------------------------------------------------------------------------- |
---|
| 4 | //# Copyright (C) 2004 |
---|
| 5 | //# ATNF |
---|
| 6 | //# |
---|
| 7 | //# This program is free software; you can redistribute it and/or modify it |
---|
| 8 | //# under the terms of the GNU General Public License as published by the Free |
---|
| 9 | //# Software Foundation; either version 2 of the License, or (at your option) |
---|
| 10 | //# any later version. |
---|
| 11 | //# |
---|
| 12 | //# This program is distributed in the hope that it will be useful, but |
---|
| 13 | //# WITHOUT ANY WARRANTY; without even the implied warranty of |
---|
| 14 | //# MERCHANTABILITY or FITNESS FOR A PARTICULAR PURPOSE. See the GNU General |
---|
| 15 | //# Public License for more details. |
---|
| 16 | //# |
---|
| 17 | //# You should have received a copy of the GNU General Public License along |
---|
| 18 | //# with this program; if not, write to the Free Software Foundation, Inc., |
---|
| 19 | //# 675 Massachusetts Ave, Cambridge, MA 02139, USA. |
---|
| 20 | //# |
---|
| 21 | //# Correspondence concerning this software should be addressed as follows: |
---|
| 22 | //# Internet email: Malte.Marquarding@csiro.au |
---|
| 23 | //# Postal address: Malte Marquarding, |
---|
| 24 | //# Australia Telescope National Facility, |
---|
| 25 | //# P.O. Box 76, |
---|
| 26 | //# Epping, NSW, 2121, |
---|
| 27 | //# AUSTRALIA |
---|
| 28 | //# |
---|
| 29 | //# $Id: |
---|
| 30 | //#--------------------------------------------------------------------------- |
---|
| 31 | |
---|
| 32 | |
---|
| 33 | // ASAP |
---|
| 34 | #include "SDLineFinder.h" |
---|
[369] | 35 | #include "SDFitter.h" |
---|
[297] | 36 | |
---|
| 37 | // STL |
---|
[343] | 38 | #include <functional> |
---|
| 39 | #include <algorithm> |
---|
[297] | 40 | #include <iostream> |
---|
[351] | 41 | #include <fstream> |
---|
[297] | 42 | |
---|
| 43 | using namespace asap; |
---|
| 44 | using namespace casa; |
---|
| 45 | using namespace std; |
---|
| 46 | using namespace boost::python; |
---|
| 47 | |
---|
[344] | 48 | namespace asap { |
---|
| 49 | |
---|
[343] | 50 | /////////////////////////////////////////////////////////////////////////////// |
---|
| 51 | // |
---|
[351] | 52 | // RunningBox - a running box calculator. This class implements |
---|
| 53 | // interations over the specified spectrum and calculates |
---|
| 54 | // running box filter statistics. |
---|
[343] | 55 | // |
---|
| 56 | |
---|
[351] | 57 | class RunningBox { |
---|
[331] | 58 | // The input data to work with. Use reference symantics to avoid |
---|
| 59 | // an unnecessary copying |
---|
| 60 | const casa::Vector<casa::Float> &spectrum; // a buffer for the spectrum |
---|
| 61 | const casa::Vector<casa::Bool> &mask; // associated mask |
---|
| 62 | const std::pair<int,int> &edge; // start and stop+1 channels |
---|
| 63 | // to work with |
---|
| 64 | |
---|
[351] | 65 | // statistics for running box filtering |
---|
| 66 | casa::Float sumf; // sum of fluxes |
---|
| 67 | casa::Float sumf2; // sum of squares of fluxes |
---|
| 68 | casa::Float sumch; // sum of channel numbers (for linear fit) |
---|
| 69 | casa::Float sumch2; // sum of squares of channel numbers (for linear fit) |
---|
| 70 | casa::Float sumfch; // sum of flux*(channel number) (for linear fit) |
---|
| 71 | |
---|
[331] | 72 | int box_chan_cntr; // actual number of channels in the box |
---|
| 73 | int max_box_nchan; // maximum allowed number of channels in the box |
---|
| 74 | // (calculated from boxsize and actual spectrum size) |
---|
[351] | 75 | // cache for derivative statistics |
---|
| 76 | mutable casa::Bool need2recalculate; // if true, values of the statistics |
---|
| 77 | // below are invalid |
---|
| 78 | mutable casa::Float linmean; // a value of the linear fit to the |
---|
| 79 | // points in the running box |
---|
| 80 | mutable casa::Float linvariance; // the same for variance |
---|
| 81 | int cur_channel; // the number of the current channel |
---|
| 82 | int start_advance; // number of channel from which the box can |
---|
| 83 | // be moved (the middle of the box, if there is no |
---|
| 84 | // masking) |
---|
| 85 | public: |
---|
| 86 | // set up the object with the references to actual data |
---|
| 87 | // as well as the number of channels in the running box |
---|
| 88 | RunningBox(const casa::Vector<casa::Float> &in_spectrum, |
---|
| 89 | const casa::Vector<casa::Bool> &in_mask, |
---|
| 90 | const std::pair<int,int> &in_edge, |
---|
| 91 | int in_max_box_nchan) throw(AipsError); |
---|
| 92 | |
---|
| 93 | // access to the statistics |
---|
| 94 | const casa::Float& getLinMean() const throw(AipsError); |
---|
| 95 | const casa::Float& getLinVariance() const throw(AipsError); |
---|
| 96 | const casa::Float aboveMean() const throw(AipsError); |
---|
| 97 | int getChannel() const throw(); |
---|
| 98 | |
---|
| 99 | // actual number of channels in the box (max_box_nchan, if no channels |
---|
| 100 | // are masked) |
---|
| 101 | int getNumberOfBoxPoints() const throw(); |
---|
[297] | 102 | |
---|
[351] | 103 | // next channel |
---|
| 104 | void next() throw(AipsError); |
---|
| 105 | |
---|
| 106 | // checking whether there are still elements |
---|
| 107 | casa::Bool haveMore() const throw(); |
---|
| 108 | |
---|
| 109 | // go to start |
---|
| 110 | void rewind() throw(AipsError); |
---|
| 111 | |
---|
| 112 | protected: |
---|
| 113 | // supplementary function to control running mean calculations. |
---|
| 114 | // It adds a specified channel to the running mean box and |
---|
| 115 | // removes (ch-maxboxnchan+1)'th channel from there |
---|
| 116 | // Channels, for which the mask is false or index is beyond the |
---|
| 117 | // allowed range, are ignored |
---|
| 118 | void advanceRunningBox(int ch) throw(casa::AipsError); |
---|
| 119 | |
---|
| 120 | // calculate derivative statistics. This function is const, because |
---|
| 121 | // it updates the cache only |
---|
| 122 | void updateDerivativeStatistics() const throw(AipsError); |
---|
| 123 | }; |
---|
| 124 | |
---|
| 125 | // |
---|
| 126 | /////////////////////////////////////////////////////////////////////////////// |
---|
| 127 | |
---|
| 128 | /////////////////////////////////////////////////////////////////////////////// |
---|
| 129 | // |
---|
| 130 | // LFAboveThreshold An algorithm for line detection using running box |
---|
| 131 | // statistics. Line is detected if it is above the |
---|
| 132 | // specified threshold at the specified number of |
---|
| 133 | // consequtive channels. Prefix LF stands for Line Finder |
---|
| 134 | // |
---|
[352] | 135 | class LFAboveThreshold : protected LFLineListOperations { |
---|
[331] | 136 | // temporary line edge channels and flag, which is True if the line |
---|
| 137 | // was detected in the previous channels. |
---|
| 138 | std::pair<int,int> cur_line; |
---|
| 139 | casa::Bool is_detected_before; |
---|
| 140 | int min_nchan; // A minimum number of consequtive |
---|
| 141 | // channels, which should satisfy |
---|
| 142 | // the detection criterion, to be |
---|
| 143 | // a detection |
---|
| 144 | casa::Float threshold; // detection threshold - the |
---|
| 145 | // minimal signal to noise ratio |
---|
[351] | 146 | std::list<pair<int,int> > &lines; // list where detections are saved |
---|
| 147 | // (pair: start and stop+1 channel) |
---|
| 148 | RunningBox *running_box; // running box filter |
---|
[551] | 149 | casa::Vector<Int> signs; // An array to store the signs of |
---|
| 150 | // the value - current mean |
---|
| 151 | // (used to search wings) |
---|
[331] | 152 | public: |
---|
[351] | 153 | |
---|
| 154 | // set up the detection criterion |
---|
| 155 | LFAboveThreshold(std::list<pair<int,int> > &in_lines, |
---|
| 156 | int in_min_nchan = 3, |
---|
| 157 | casa::Float in_threshold = 5) throw(); |
---|
| 158 | virtual ~LFAboveThreshold() throw(); |
---|
| 159 | |
---|
[331] | 160 | // replace the detection criterion |
---|
| 161 | void setCriterion(int in_min_nchan, casa::Float in_threshold) throw(); |
---|
[297] | 162 | |
---|
[551] | 163 | // return the array with signs of the value-current mean |
---|
| 164 | // An element is +1 if value>mean, -1 if less, 0 if equal. |
---|
| 165 | // This array is updated each time the findLines method is called and |
---|
| 166 | // is used to search the line wings |
---|
| 167 | const casa::Vector<Int>& getSigns() const throw(); |
---|
| 168 | |
---|
[331] | 169 | // find spectral lines and add them into list |
---|
[344] | 170 | // if statholder is not NULL, the accumulate function of it will be |
---|
| 171 | // called for each channel to save statistics |
---|
[351] | 172 | // spectrum, mask and edge - reference to the data |
---|
| 173 | // max_box_nchan - number of channels in the running box |
---|
| 174 | void findLines(const casa::Vector<casa::Float> &spectrum, |
---|
| 175 | const casa::Vector<casa::Bool> &mask, |
---|
| 176 | const std::pair<int,int> &edge, |
---|
[352] | 177 | int max_box_nchan) throw(casa::AipsError); |
---|
[351] | 178 | |
---|
[331] | 179 | protected: |
---|
[297] | 180 | |
---|
[331] | 181 | // process a channel: update curline and is_detected before and |
---|
| 182 | // add a new line to the list, if necessary using processCurLine() |
---|
[351] | 183 | // detect=true indicates that the current channel satisfies the criterion |
---|
| 184 | void processChannel(Bool detect, const casa::Vector<casa::Bool> &mask) |
---|
| 185 | throw(casa::AipsError); |
---|
[297] | 186 | |
---|
[331] | 187 | // process the interval of channels stored in curline |
---|
| 188 | // if it satisfies the criterion, add this interval as a new line |
---|
[351] | 189 | void processCurLine(const casa::Vector<casa::Bool> &mask) |
---|
| 190 | throw(casa::AipsError); |
---|
[331] | 191 | }; |
---|
[344] | 192 | |
---|
| 193 | // |
---|
| 194 | /////////////////////////////////////////////////////////////////////////////// |
---|
[351] | 195 | |
---|
[331] | 196 | } // namespace asap |
---|
[297] | 197 | |
---|
[344] | 198 | /////////////////////////////////////////////////////////////////////////////// |
---|
[343] | 199 | // |
---|
[351] | 200 | // RunningBox - a running box calculator. This class implements |
---|
| 201 | // interations over the specified spectrum and calculates |
---|
| 202 | // running box filter statistics. |
---|
[331] | 203 | // |
---|
[297] | 204 | |
---|
[331] | 205 | // set up the object with the references to actual data |
---|
| 206 | // and the number of channels in the running box |
---|
[351] | 207 | RunningBox::RunningBox(const casa::Vector<casa::Float> &in_spectrum, |
---|
| 208 | const casa::Vector<casa::Bool> &in_mask, |
---|
| 209 | const std::pair<int,int> &in_edge, |
---|
| 210 | int in_max_box_nchan) throw(AipsError) : |
---|
[331] | 211 | spectrum(in_spectrum), mask(in_mask), edge(in_edge), |
---|
[351] | 212 | max_box_nchan(in_max_box_nchan) |
---|
| 213 | { |
---|
| 214 | rewind(); |
---|
| 215 | } |
---|
[331] | 216 | |
---|
[351] | 217 | void RunningBox::rewind() throw(AipsError) { |
---|
| 218 | // fill statistics for initial box |
---|
| 219 | box_chan_cntr=0; // no channels are currently in the box |
---|
| 220 | sumf=0.; // initialize statistics |
---|
| 221 | sumf2=0.; |
---|
| 222 | sumch=0.; |
---|
| 223 | sumch2=0.; |
---|
| 224 | sumfch=0.; |
---|
| 225 | int initial_box_ch=edge.first; |
---|
| 226 | for (;initial_box_ch<edge.second && box_chan_cntr<max_box_nchan; |
---|
| 227 | ++initial_box_ch) |
---|
| 228 | advanceRunningBox(initial_box_ch); |
---|
| 229 | |
---|
| 230 | if (initial_box_ch==edge.second) |
---|
| 231 | throw AipsError("RunningBox::rewind - too much channels are masked"); |
---|
| 232 | |
---|
| 233 | cur_channel=edge.first; |
---|
| 234 | start_advance=initial_box_ch-max_box_nchan/2; |
---|
| 235 | } |
---|
| 236 | |
---|
| 237 | // access to the statistics |
---|
| 238 | const casa::Float& RunningBox::getLinMean() const throw(AipsError) |
---|
[331] | 239 | { |
---|
[351] | 240 | DebugAssert(cur_channel<edge.second, AipsError); |
---|
| 241 | if (need2recalculate) updateDerivativeStatistics(); |
---|
| 242 | return linmean; |
---|
[297] | 243 | } |
---|
| 244 | |
---|
[351] | 245 | const casa::Float& RunningBox::getLinVariance() const throw(AipsError) |
---|
| 246 | { |
---|
| 247 | DebugAssert(cur_channel<edge.second, AipsError); |
---|
| 248 | if (need2recalculate) updateDerivativeStatistics(); |
---|
| 249 | return linvariance; |
---|
| 250 | } |
---|
[331] | 251 | |
---|
[351] | 252 | const casa::Float RunningBox::aboveMean() const throw(AipsError) |
---|
| 253 | { |
---|
| 254 | DebugAssert(cur_channel<edge.second, AipsError); |
---|
| 255 | if (need2recalculate) updateDerivativeStatistics(); |
---|
| 256 | return spectrum[cur_channel]-linmean; |
---|
| 257 | } |
---|
| 258 | |
---|
| 259 | int RunningBox::getChannel() const throw() |
---|
| 260 | { |
---|
| 261 | return cur_channel; |
---|
| 262 | } |
---|
| 263 | |
---|
| 264 | // actual number of channels in the box (max_box_nchan, if no channels |
---|
| 265 | // are masked) |
---|
| 266 | int RunningBox::getNumberOfBoxPoints() const throw() |
---|
| 267 | { |
---|
| 268 | return box_chan_cntr; |
---|
| 269 | } |
---|
| 270 | |
---|
[297] | 271 | // supplementary function to control running mean calculations. |
---|
| 272 | // It adds a specified channel to the running mean box and |
---|
| 273 | // removes (ch-max_box_nchan+1)'th channel from there |
---|
| 274 | // Channels, for which the mask is false or index is beyond the |
---|
| 275 | // allowed range, are ignored |
---|
[351] | 276 | void RunningBox::advanceRunningBox(int ch) throw(AipsError) |
---|
[297] | 277 | { |
---|
| 278 | if (ch>=edge.first && ch<edge.second) |
---|
| 279 | if (mask[ch]) { // ch is a valid channel |
---|
| 280 | ++box_chan_cntr; |
---|
[351] | 281 | sumf+=spectrum[ch]; |
---|
| 282 | sumf2+=square(spectrum[ch]); |
---|
| 283 | sumch+=Float(ch); |
---|
| 284 | sumch2+=square(Float(ch)); |
---|
| 285 | sumfch+=spectrum[ch]*Float(ch); |
---|
| 286 | need2recalculate=True; |
---|
[297] | 287 | } |
---|
| 288 | int ch2remove=ch-max_box_nchan; |
---|
| 289 | if (ch2remove>=edge.first && ch2remove<edge.second) |
---|
| 290 | if (mask[ch2remove]) { // ch2remove is a valid channel |
---|
| 291 | --box_chan_cntr; |
---|
[351] | 292 | sumf-=spectrum[ch2remove]; |
---|
| 293 | sumf2-=square(spectrum[ch2remove]); |
---|
| 294 | sumch-=Float(ch2remove); |
---|
| 295 | sumch2-=square(Float(ch2remove)); |
---|
| 296 | sumfch-=spectrum[ch2remove]*Float(ch2remove); |
---|
| 297 | need2recalculate=True; |
---|
[297] | 298 | } |
---|
| 299 | } |
---|
| 300 | |
---|
[351] | 301 | // next channel |
---|
| 302 | void RunningBox::next() throw(AipsError) |
---|
[297] | 303 | { |
---|
[351] | 304 | AlwaysAssert(cur_channel<edge.second,AipsError); |
---|
| 305 | ++cur_channel; |
---|
| 306 | if (cur_channel+max_box_nchan/2<edge.second && cur_channel>=start_advance) |
---|
| 307 | advanceRunningBox(cur_channel+max_box_nchan/2); // update statistics |
---|
[297] | 308 | } |
---|
| 309 | |
---|
[351] | 310 | // checking whether there are still elements |
---|
| 311 | casa::Bool RunningBox::haveMore() const throw() |
---|
| 312 | { |
---|
| 313 | return cur_channel<edge.second; |
---|
| 314 | } |
---|
| 315 | |
---|
| 316 | // calculate derivative statistics. This function is const, because |
---|
| 317 | // it updates the cache only |
---|
| 318 | void RunningBox::updateDerivativeStatistics() const throw(AipsError) |
---|
| 319 | { |
---|
| 320 | AlwaysAssert(box_chan_cntr, AipsError); |
---|
| 321 | |
---|
| 322 | Float mean=sumf/Float(box_chan_cntr); |
---|
| 323 | |
---|
| 324 | // linear LSF formulae |
---|
| 325 | Float meanch=sumch/Float(box_chan_cntr); |
---|
| 326 | Float meanch2=sumch2/Float(box_chan_cntr); |
---|
| 327 | if (meanch==meanch2 || box_chan_cntr<3) { |
---|
| 328 | // vertical line in the spectrum, can't calculate linmean and linvariance |
---|
| 329 | linmean=0.; |
---|
| 330 | linvariance=0.; |
---|
| 331 | } else { |
---|
| 332 | Float coeff=(sumfch/Float(box_chan_cntr)-meanch*mean)/ |
---|
| 333 | (meanch2-square(meanch)); |
---|
| 334 | linmean=coeff*(Float(cur_channel)-meanch)+mean; |
---|
| 335 | linvariance=sqrt(sumf2/Float(box_chan_cntr)-square(mean)- |
---|
| 336 | square(coeff)*(meanch2-square(meanch))); |
---|
| 337 | } |
---|
| 338 | need2recalculate=False; |
---|
| 339 | } |
---|
| 340 | |
---|
| 341 | |
---|
| 342 | // |
---|
| 343 | /////////////////////////////////////////////////////////////////////////////// |
---|
| 344 | |
---|
| 345 | /////////////////////////////////////////////////////////////////////////////// |
---|
| 346 | // |
---|
| 347 | // LFAboveThreshold - a running mean algorithm for line detection |
---|
| 348 | // |
---|
| 349 | // |
---|
| 350 | |
---|
| 351 | |
---|
| 352 | // set up the detection criterion |
---|
| 353 | LFAboveThreshold::LFAboveThreshold(std::list<pair<int,int> > &in_lines, |
---|
| 354 | int in_min_nchan, |
---|
| 355 | casa::Float in_threshold) throw() : |
---|
| 356 | min_nchan(in_min_nchan), threshold(in_threshold), |
---|
| 357 | lines(in_lines), running_box(NULL) {} |
---|
| 358 | |
---|
| 359 | LFAboveThreshold::~LFAboveThreshold() throw() |
---|
| 360 | { |
---|
| 361 | if (running_box!=NULL) delete running_box; |
---|
| 362 | } |
---|
| 363 | |
---|
| 364 | // replace the detection criterion |
---|
| 365 | void LFAboveThreshold::setCriterion(int in_min_nchan, casa::Float in_threshold) |
---|
| 366 | throw() |
---|
| 367 | { |
---|
| 368 | min_nchan=in_min_nchan; |
---|
| 369 | threshold=in_threshold; |
---|
| 370 | } |
---|
| 371 | |
---|
| 372 | |
---|
[297] | 373 | // process a channel: update cur_line and is_detected before and |
---|
| 374 | // add a new line to the list, if necessary |
---|
[351] | 375 | void LFAboveThreshold::processChannel(Bool detect, |
---|
| 376 | const casa::Vector<casa::Bool> &mask) throw(casa::AipsError) |
---|
[297] | 377 | { |
---|
| 378 | try { |
---|
[351] | 379 | if (detect) { |
---|
[297] | 380 | if (is_detected_before) |
---|
[351] | 381 | cur_line.second=running_box->getChannel()+1; |
---|
[297] | 382 | else { |
---|
| 383 | is_detected_before=True; |
---|
[351] | 384 | cur_line.first=running_box->getChannel(); |
---|
| 385 | cur_line.second=running_box->getChannel()+1; |
---|
[297] | 386 | } |
---|
[351] | 387 | } else processCurLine(mask); |
---|
[297] | 388 | } |
---|
| 389 | catch (const AipsError &ae) { |
---|
| 390 | throw; |
---|
| 391 | } |
---|
| 392 | catch (const exception &ex) { |
---|
[351] | 393 | throw AipsError(String("LFAboveThreshold::processChannel - STL error: ")+ex.what()); |
---|
[297] | 394 | } |
---|
| 395 | } |
---|
| 396 | |
---|
| 397 | // process the interval of channels stored in cur_line |
---|
| 398 | // if it satisfies the criterion, add this interval as a new line |
---|
[351] | 399 | void LFAboveThreshold::processCurLine(const casa::Vector<casa::Bool> &mask) |
---|
[331] | 400 | throw(casa::AipsError) |
---|
[297] | 401 | { |
---|
| 402 | try { |
---|
| 403 | if (is_detected_before) { |
---|
| 404 | if (cur_line.second-cur_line.first>min_nchan) { |
---|
| 405 | // it was a detection. We need to change the list |
---|
| 406 | Bool add_new_line=False; |
---|
| 407 | if (lines.size()) { |
---|
| 408 | for (int i=lines.back().second;i<cur_line.first;++i) |
---|
| 409 | if (mask[i]) { // one valid channel in between |
---|
| 410 | // means that we deal with a separate line |
---|
| 411 | add_new_line=True; |
---|
| 412 | break; |
---|
| 413 | } |
---|
| 414 | } else add_new_line=True; |
---|
| 415 | if (add_new_line) |
---|
| 416 | lines.push_back(cur_line); |
---|
| 417 | else lines.back().second=cur_line.second; |
---|
| 418 | } |
---|
| 419 | is_detected_before=False; |
---|
| 420 | } |
---|
| 421 | } |
---|
| 422 | catch (const AipsError &ae) { |
---|
| 423 | throw; |
---|
| 424 | } |
---|
| 425 | catch (const exception &ex) { |
---|
[351] | 426 | throw AipsError(String("LFAboveThreshold::processCurLine - STL error: ")+ex.what()); |
---|
[297] | 427 | } |
---|
| 428 | } |
---|
| 429 | |
---|
[551] | 430 | // return the array with signs of the value-current mean |
---|
| 431 | // An element is +1 if value>mean, -1 if less, 0 if equal. |
---|
| 432 | // This array is updated each time the findLines method is called and |
---|
| 433 | // is used to search the line wings |
---|
| 434 | const casa::Vector<Int>& LFAboveThreshold::getSigns() const throw() |
---|
| 435 | { |
---|
| 436 | return signs; |
---|
| 437 | } |
---|
| 438 | |
---|
[331] | 439 | // find spectral lines and add them into list |
---|
[351] | 440 | void LFAboveThreshold::findLines(const casa::Vector<casa::Float> &spectrum, |
---|
| 441 | const casa::Vector<casa::Bool> &mask, |
---|
| 442 | const std::pair<int,int> &edge, |
---|
[352] | 443 | int max_box_nchan) |
---|
[331] | 444 | throw(casa::AipsError) |
---|
| 445 | { |
---|
| 446 | const int minboxnchan=4; |
---|
[351] | 447 | try { |
---|
[331] | 448 | |
---|
[351] | 449 | if (running_box!=NULL) delete running_box; |
---|
| 450 | running_box=new RunningBox(spectrum,mask,edge,max_box_nchan); |
---|
[368] | 451 | |
---|
| 452 | |
---|
| 453 | // determine the off-line variance first |
---|
| 454 | // an assumption made: lines occupy a small part of the spectrum |
---|
[369] | 455 | |
---|
[368] | 456 | std::vector<float> variances(edge.second-edge.first); |
---|
| 457 | DebugAssert(variances.size(),AipsError); |
---|
[369] | 458 | |
---|
[368] | 459 | for (;running_box->haveMore();running_box->next()) |
---|
| 460 | variances[running_box->getChannel()-edge.first]= |
---|
| 461 | running_box->getLinVariance(); |
---|
[369] | 462 | |
---|
[368] | 463 | // in the future we probably should do a proper Chi^2 estimation |
---|
| 464 | // now a simple 80% of smaller values will be used. |
---|
| 465 | // it may degrade the performance of the algorithm for weak lines |
---|
| 466 | // due to a bias of the Chi^2 distribution. |
---|
| 467 | stable_sort(variances.begin(),variances.end()); |
---|
[369] | 468 | |
---|
[368] | 469 | Float offline_variance=0; |
---|
| 470 | uInt offline_cnt=uInt(0.8*variances.size()); |
---|
| 471 | if (!offline_cnt) offline_cnt=variances.size(); // no much else left, |
---|
| 472 | // although it is very inaccurate |
---|
| 473 | for (uInt n=0;n<offline_cnt;++n) |
---|
| 474 | offline_variance+=variances[n]; |
---|
| 475 | offline_variance/=Float(offline_cnt); |
---|
| 476 | |
---|
[351] | 477 | // actual search algorithm |
---|
| 478 | is_detected_before=False; |
---|
[368] | 479 | |
---|
[551] | 480 | // initiate the signs array |
---|
| 481 | signs.resize(spectrum.nelements()); |
---|
| 482 | signs=Vector<Int>(spectrum.nelements(),0); |
---|
| 483 | |
---|
[369] | 484 | //ofstream os("dbg.dat"); |
---|
[368] | 485 | for (running_box->rewind();running_box->haveMore(); |
---|
| 486 | running_box->next()) { |
---|
[351] | 487 | const int ch=running_box->getChannel(); |
---|
| 488 | if (running_box->getNumberOfBoxPoints()>=minboxnchan) |
---|
[352] | 489 | processChannel(mask[ch] && (fabs(running_box->aboveMean()) >= |
---|
[368] | 490 | threshold*offline_variance), mask); |
---|
[351] | 491 | else processCurLine(mask); // just finish what was accumulated before |
---|
[352] | 492 | const Float buf=running_box->aboveMean(); |
---|
| 493 | if (buf>0) signs[ch]=1; |
---|
| 494 | else if (buf<0) signs[ch]=-1; |
---|
[368] | 495 | else if (buf==0) signs[ch]=0; |
---|
[369] | 496 | // os<<ch<<" "<<spectrum[ch]<<" "<<running_box->getLinMean()<<" "<< |
---|
| 497 | // threshold*offline_variance<<endl; |
---|
[351] | 498 | } |
---|
[352] | 499 | if (lines.size()) |
---|
| 500 | searchForWings(lines,signs,mask,edge); |
---|
[344] | 501 | } |
---|
[351] | 502 | catch (const AipsError &ae) { |
---|
| 503 | throw; |
---|
| 504 | } |
---|
| 505 | catch (const exception &ex) { |
---|
| 506 | throw AipsError(String("LFAboveThreshold::findLines - STL error: ")+ex.what()); |
---|
| 507 | } |
---|
[331] | 508 | } |
---|
| 509 | |
---|
| 510 | // |
---|
| 511 | /////////////////////////////////////////////////////////////////////////////// |
---|
| 512 | |
---|
[343] | 513 | /////////////////////////////////////////////////////////////////////////////// |
---|
| 514 | // |
---|
[352] | 515 | // LFLineListOperations::IntersectsWith - An auxiliary object function |
---|
| 516 | // to test whether two lines have a non-void intersection |
---|
[343] | 517 | // |
---|
[331] | 518 | |
---|
[343] | 519 | |
---|
| 520 | // line1 - range of the first line: start channel and stop+1 |
---|
[352] | 521 | LFLineListOperations::IntersectsWith::IntersectsWith(const std::pair<int,int> &in_line1) : |
---|
[343] | 522 | line1(in_line1) {} |
---|
| 523 | |
---|
| 524 | |
---|
| 525 | // return true if line2 intersects with line1 with at least one |
---|
| 526 | // common channel, and false otherwise |
---|
| 527 | // line2 - range of the second line: start channel and stop+1 |
---|
[352] | 528 | bool LFLineListOperations::IntersectsWith::operator()(const std::pair<int,int> &line2) |
---|
[343] | 529 | const throw() |
---|
| 530 | { |
---|
| 531 | if (line2.second<line1.first) return false; // line2 is at lower channels |
---|
| 532 | if (line2.first>line1.second) return false; // line2 is at upper channels |
---|
| 533 | return true; // line2 has an intersection or is adjacent to line1 |
---|
| 534 | } |
---|
| 535 | |
---|
| 536 | // |
---|
| 537 | /////////////////////////////////////////////////////////////////////////////// |
---|
| 538 | |
---|
| 539 | /////////////////////////////////////////////////////////////////////////////// |
---|
| 540 | // |
---|
[352] | 541 | // LFLineListOperations::BuildUnion - An auxiliary object function to build a union |
---|
[343] | 542 | // of several lines to account for a possibility of merging the nearby lines |
---|
| 543 | // |
---|
| 544 | |
---|
| 545 | // set an initial line (can be a first line in the sequence) |
---|
[352] | 546 | LFLineListOperations::BuildUnion::BuildUnion(const std::pair<int,int> &line1) : |
---|
[343] | 547 | temp_line(line1) {} |
---|
| 548 | |
---|
| 549 | // update temp_line with a union of temp_line and new_line |
---|
| 550 | // provided there is no gap between the lines |
---|
[352] | 551 | void LFLineListOperations::BuildUnion::operator()(const std::pair<int,int> &new_line) |
---|
[343] | 552 | throw() |
---|
| 553 | { |
---|
| 554 | if (new_line.first<temp_line.first) temp_line.first=new_line.first; |
---|
| 555 | if (new_line.second>temp_line.second) temp_line.second=new_line.second; |
---|
| 556 | } |
---|
| 557 | |
---|
| 558 | // return the result (temp_line) |
---|
[352] | 559 | const std::pair<int,int>& LFLineListOperations::BuildUnion::result() const throw() |
---|
[343] | 560 | { |
---|
| 561 | return temp_line; |
---|
| 562 | } |
---|
| 563 | |
---|
| 564 | // |
---|
| 565 | /////////////////////////////////////////////////////////////////////////////// |
---|
| 566 | |
---|
| 567 | /////////////////////////////////////////////////////////////////////////////// |
---|
| 568 | // |
---|
[352] | 569 | // LFLineListOperations::LaterThan - An auxiliary object function to test whether a |
---|
[343] | 570 | // specified line is at lower spectral channels (to preserve the order in |
---|
| 571 | // the line list) |
---|
| 572 | // |
---|
| 573 | |
---|
| 574 | // setup the line to compare with |
---|
[352] | 575 | LFLineListOperations::LaterThan::LaterThan(const std::pair<int,int> &in_line1) : |
---|
[343] | 576 | line1(in_line1) {} |
---|
| 577 | |
---|
| 578 | // return true if line2 should be placed later than line1 |
---|
| 579 | // in the ordered list (so, it is at greater channel numbers) |
---|
[352] | 580 | bool LFLineListOperations::LaterThan::operator()(const std::pair<int,int> &line2) |
---|
[343] | 581 | const throw() |
---|
| 582 | { |
---|
| 583 | if (line2.second<line1.first) return false; // line2 is at lower channels |
---|
| 584 | if (line2.first>line1.second) return true; // line2 is at upper channels |
---|
| 585 | |
---|
| 586 | // line2 intersects with line1. We should have no such situation in |
---|
| 587 | // practice |
---|
| 588 | return line2.first>line1.first; |
---|
| 589 | } |
---|
| 590 | |
---|
| 591 | // |
---|
| 592 | /////////////////////////////////////////////////////////////////////////////// |
---|
| 593 | |
---|
| 594 | |
---|
| 595 | /////////////////////////////////////////////////////////////////////////////// |
---|
| 596 | // |
---|
[331] | 597 | // SDLineFinder - a class for automated spectral line search |
---|
[343] | 598 | // |
---|
| 599 | // |
---|
[331] | 600 | |
---|
| 601 | SDLineFinder::SDLineFinder() throw() : edge(0,0) |
---|
| 602 | { |
---|
[369] | 603 | setOptions(); |
---|
[331] | 604 | } |
---|
| 605 | |
---|
[369] | 606 | // set the parameters controlling algorithm |
---|
| 607 | // in_threshold a single channel threshold default is sqrt(3), which |
---|
| 608 | // means together with 3 minimum channels at least 3 sigma |
---|
| 609 | // detection criterion |
---|
| 610 | // For bad baseline shape, in_threshold may need to be |
---|
| 611 | // increased |
---|
| 612 | // in_min_nchan minimum number of channels above the threshold to report |
---|
| 613 | // a detection, default is 3 |
---|
| 614 | // in_avg_limit perform the averaging of no more than in_avg_limit |
---|
| 615 | // adjacent channels to search for broad lines |
---|
| 616 | // Default is 8, but for a bad baseline shape this |
---|
| 617 | // parameter should be decreased (may be even down to a |
---|
| 618 | // minimum of 1 to disable this option) to avoid |
---|
| 619 | // confusing of baseline undulations with a real line. |
---|
| 620 | // Setting a very large value doesn't usually provide |
---|
| 621 | // valid detections. |
---|
| 622 | // in_box_size the box size for running mean calculation. Default is |
---|
| 623 | // 1./5. of the whole spectrum size |
---|
| 624 | void SDLineFinder::setOptions(const casa::Float &in_threshold, |
---|
| 625 | const casa::Int &in_min_nchan, |
---|
| 626 | const casa::Int &in_avg_limit, |
---|
| 627 | const casa::Float &in_box_size) throw() |
---|
| 628 | { |
---|
| 629 | threshold=in_threshold; |
---|
| 630 | min_nchan=in_min_nchan; |
---|
| 631 | avg_limit=in_avg_limit; |
---|
| 632 | box_size=in_box_size; |
---|
| 633 | } |
---|
| 634 | |
---|
[331] | 635 | SDLineFinder::~SDLineFinder() throw(AipsError) {} |
---|
| 636 | |
---|
| 637 | // set scan to work with (in_scan parameter), associated mask (in_mask |
---|
| 638 | // parameter) and the edge channel rejection (in_edge parameter) |
---|
| 639 | // if in_edge has zero length, all channels chosen by mask will be used |
---|
| 640 | // if in_edge has one element only, it represents the number of |
---|
| 641 | // channels to drop from both sides of the spectrum |
---|
| 642 | // in_edge is introduced for convinience, although all functionality |
---|
| 643 | // can be achieved using a spectrum mask only |
---|
| 644 | void SDLineFinder::setScan(const SDMemTableWrapper &in_scan, |
---|
| 645 | const std::vector<bool> &in_mask, |
---|
| 646 | const boost::python::tuple &in_edge) throw(AipsError) |
---|
| 647 | { |
---|
| 648 | try { |
---|
| 649 | scan=in_scan.getCP(); |
---|
| 650 | AlwaysAssert(!scan.null(),AipsError); |
---|
[516] | 651 | |
---|
[331] | 652 | mask=in_mask; |
---|
| 653 | if (mask.nelements()!=scan->nChan()) |
---|
| 654 | throw AipsError("SDLineFinder::setScan - in_scan and in_mask have different" |
---|
| 655 | "number of spectral channels."); |
---|
| 656 | |
---|
| 657 | // number of elements in the in_edge tuple |
---|
| 658 | int n=extract<int>(in_edge.attr("__len__")()); |
---|
| 659 | if (n>2 || n<0) |
---|
| 660 | throw AipsError("SDLineFinder::setScan - the length of the in_edge parameter" |
---|
| 661 | "should not exceed 2"); |
---|
| 662 | if (!n) { |
---|
| 663 | // all spectrum, no rejection |
---|
| 664 | edge.first=0; |
---|
| 665 | edge.second=scan->nChan(); |
---|
| 666 | } else { |
---|
| 667 | edge.first=extract<int>(in_edge.attr("__getitem__")(0)); |
---|
| 668 | if (edge.first<0) |
---|
| 669 | throw AipsError("SDLineFinder::setScan - the in_edge parameter has a negative" |
---|
| 670 | "number of channels to drop"); |
---|
| 671 | if (edge.first>=scan->nChan()) |
---|
| 672 | throw AipsError("SDLineFinder::setScan - all channels are rejected by the in_edge parameter"); |
---|
| 673 | if (n==2) { |
---|
| 674 | edge.second=extract<int>(in_edge.attr("__getitem__")(1)); |
---|
| 675 | if (edge.second<0) |
---|
| 676 | throw AipsError("SDLineFinder::setScan - the in_edge parameter has a negative" |
---|
| 677 | "number of channels to drop"); |
---|
| 678 | edge.second=scan->nChan()-edge.second; |
---|
| 679 | } else edge.second=scan->nChan()-edge.first; |
---|
[369] | 680 | if (edge.second<0 || (edge.first>=edge.second)) |
---|
[331] | 681 | throw AipsError("SDLineFinder::setScan - all channels are rejected by the in_edge parameter"); |
---|
| 682 | } |
---|
| 683 | } |
---|
| 684 | catch(const AipsError &ae) { |
---|
| 685 | // setScan is unsuccessfull, reset scan/mask/edge |
---|
| 686 | scan=CountedConstPtr<SDMemTable>(); // null pointer |
---|
| 687 | mask.resize(0); |
---|
| 688 | edge=pair<int,int>(0,0); |
---|
| 689 | throw; |
---|
| 690 | } |
---|
| 691 | } |
---|
| 692 | |
---|
| 693 | // search for spectral lines. Number of lines found is returned |
---|
[370] | 694 | int SDLineFinder::findLines(const casa::uInt &whichRow) throw(casa::AipsError) |
---|
[368] | 695 | { |
---|
[331] | 696 | const int minboxnchan=4; |
---|
| 697 | if (scan.null()) |
---|
| 698 | throw AipsError("SDLineFinder::findLines - a scan should be set first," |
---|
| 699 | " use set_scan"); |
---|
| 700 | DebugAssert(mask.nelements()==scan->nChan(), AipsError); |
---|
| 701 | int max_box_nchan=int(scan->nChan()*box_size); // number of channels in running |
---|
| 702 | // box |
---|
| 703 | if (max_box_nchan<2) |
---|
| 704 | throw AipsError("SDLineFinder::findLines - box_size is too small"); |
---|
| 705 | |
---|
[370] | 706 | scan->getSpectrum(spectrum, whichRow); |
---|
[331] | 707 | |
---|
| 708 | lines.resize(0); // search from the scratch |
---|
[370] | 709 | last_row_used=whichRow; |
---|
[331] | 710 | Vector<Bool> temp_mask(mask); |
---|
[351] | 711 | |
---|
| 712 | Bool first_pass=True; |
---|
[368] | 713 | Int avg_factor=1; // this number of adjacent channels is averaged together |
---|
| 714 | // the total number of the channels is not altered |
---|
| 715 | // instead, min_nchan is also scaled |
---|
| 716 | // it helps to search for broad lines |
---|
[551] | 717 | Vector<Int> signs; // a buffer for signs of the value - mean quantity |
---|
| 718 | // see LFAboveThreshold for details |
---|
| 719 | // We need only signs resulted from last iteration |
---|
| 720 | // because all previous values may be corrupted by the |
---|
| 721 | // presence of spectral lines |
---|
[344] | 722 | while (true) { |
---|
[351] | 723 | // a buffer for new lines found at this iteration |
---|
[368] | 724 | std::list<pair<int,int> > new_lines; |
---|
[351] | 725 | |
---|
| 726 | try { |
---|
[369] | 727 | // line find algorithm |
---|
| 728 | LFAboveThreshold lfalg(new_lines,avg_factor*min_nchan, threshold); |
---|
[352] | 729 | lfalg.findLines(spectrum,temp_mask,edge,max_box_nchan); |
---|
[551] | 730 | signs.resize(lfalg.getSigns().nelements()); |
---|
| 731 | signs=lfalg.getSigns(); |
---|
[368] | 732 | first_pass=False; |
---|
| 733 | if (!new_lines.size()) |
---|
| 734 | throw AipsError("spurious"); // nothing new - use the same |
---|
| 735 | // code as for a real exception |
---|
[351] | 736 | } |
---|
| 737 | catch(const AipsError &ae) { |
---|
| 738 | if (first_pass) throw; |
---|
[368] | 739 | // nothing new - proceed to the next step of averaging, if any |
---|
| 740 | // (to search for broad lines) |
---|
[369] | 741 | if (avg_factor>avg_limit) break; // averaging up to avg_limit |
---|
| 742 | // adjacent channels, |
---|
| 743 | // stop after that |
---|
[368] | 744 | avg_factor*=2; // twice as more averaging |
---|
[369] | 745 | subtractBaseline(temp_mask,9); |
---|
[368] | 746 | averageAdjacentChannels(temp_mask,avg_factor); |
---|
| 747 | continue; |
---|
[351] | 748 | } |
---|
[368] | 749 | keepStrongestOnly(temp_mask,new_lines,max_box_nchan); |
---|
[343] | 750 | // update the list (lines) merging intervals, if necessary |
---|
[344] | 751 | addNewSearchResult(new_lines,lines); |
---|
| 752 | // get a new mask |
---|
| 753 | temp_mask=getMask(); |
---|
[343] | 754 | } |
---|
[551] | 755 | |
---|
| 756 | // an additional search for wings because in the presence of very strong |
---|
| 757 | // lines temporary mean used at each iteration will be higher than |
---|
| 758 | // the true mean |
---|
| 759 | |
---|
| 760 | if (lines.size()) |
---|
| 761 | LFLineListOperations::searchForWings(lines,signs,mask,edge); |
---|
| 762 | |
---|
[331] | 763 | return int(lines.size()); |
---|
| 764 | } |
---|
| 765 | |
---|
[369] | 766 | // auxiliary function to fit and subtract a polynomial from the current |
---|
| 767 | // spectrum. It uses the SDFitter class. This action is required before |
---|
| 768 | // reducing the spectral resolution if the baseline shape is bad |
---|
| 769 | void SDLineFinder::subtractBaseline(const casa::Vector<casa::Bool> &temp_mask, |
---|
| 770 | const casa::Int &order) throw(casa::AipsError) |
---|
| 771 | { |
---|
| 772 | AlwaysAssert(spectrum.nelements(),AipsError); |
---|
| 773 | // use the fact that temp_mask excludes channels rejected at the edge |
---|
| 774 | SDFitter sdf; |
---|
| 775 | std::vector<float> absc(spectrum.nelements()); |
---|
| 776 | for (Int i=0;i<absc.size();++i) |
---|
| 777 | absc[i]=float(i)/float(spectrum.nelements()); |
---|
| 778 | std::vector<float> spec; |
---|
| 779 | spectrum.tovector(spec); |
---|
| 780 | std::vector<bool> std_mask; |
---|
| 781 | temp_mask.tovector(std_mask); |
---|
| 782 | sdf.setData(absc,spec,std_mask); |
---|
| 783 | sdf.setExpression("poly",order); |
---|
| 784 | if (!sdf.fit()) return; // fit failed, use old spectrum |
---|
| 785 | spectrum=casa::Vector<casa::Float>(sdf.getResidual()); |
---|
| 786 | } |
---|
| 787 | |
---|
[368] | 788 | // auxiliary function to average adjacent channels and update the mask |
---|
| 789 | // if at least one channel involved in summation is masked, all |
---|
| 790 | // output channels will be masked. This function works with the |
---|
| 791 | // spectrum and edge fields of this class, but updates the mask |
---|
| 792 | // array specified, rather than the field of this class |
---|
| 793 | // boxsize - a number of adjacent channels to average |
---|
| 794 | void SDLineFinder::averageAdjacentChannels(casa::Vector<casa::Bool> &mask2update, |
---|
| 795 | const casa::Int &boxsize) |
---|
| 796 | throw(casa::AipsError) |
---|
| 797 | { |
---|
| 798 | DebugAssert(mask2update.nelements()==spectrum.nelements(), AipsError); |
---|
| 799 | DebugAssert(boxsize!=0,AipsError); |
---|
| 800 | |
---|
| 801 | for (int n=edge.first;n<edge.second;n+=boxsize) { |
---|
| 802 | DebugAssert(n<spectrum.nelements(),AipsError); |
---|
| 803 | int nboxch=0; // number of channels currently in the box |
---|
| 804 | Float mean=0; // buffer for mean calculations |
---|
| 805 | for (int k=n;k<n+boxsize && k<edge.second;++k) |
---|
| 806 | if (mask2update[k]) { // k is a valid channel |
---|
| 807 | mean+=spectrum[k]; |
---|
| 808 | ++nboxch; |
---|
| 809 | } |
---|
| 810 | if (nboxch<boxsize) // mask these channels |
---|
| 811 | for (int k=n;k<n+boxsize && k<edge.second;++k) |
---|
| 812 | mask2update[k]=False; |
---|
| 813 | else { |
---|
| 814 | mean/=Float(boxsize); |
---|
| 815 | for (int k=n;k<n+boxsize && k<edge.second;++k) |
---|
| 816 | spectrum[k]=mean; |
---|
| 817 | } |
---|
| 818 | } |
---|
| 819 | } |
---|
[331] | 820 | |
---|
[368] | 821 | |
---|
[297] | 822 | // get the mask to mask out all lines that have been found (default) |
---|
| 823 | // if invert=true, only channels belong to lines will be unmasked |
---|
| 824 | // Note: all channels originally masked by the input mask (in_mask |
---|
| 825 | // in setScan) or dropped out by the edge parameter (in_edge |
---|
| 826 | // in setScan) are still excluded regardless on the invert option |
---|
| 827 | std::vector<bool> SDLineFinder::getMask(bool invert) |
---|
| 828 | const throw(casa::AipsError) |
---|
| 829 | { |
---|
| 830 | try { |
---|
| 831 | if (scan.null()) |
---|
| 832 | throw AipsError("SDLineFinder::getMask - a scan should be set first," |
---|
| 833 | " use set_scan followed by find_lines"); |
---|
| 834 | DebugAssert(mask.nelements()==scan->nChan(), AipsError); |
---|
| 835 | /* |
---|
| 836 | if (!lines.size()) |
---|
| 837 | throw AipsError("SDLineFinder::getMask - one have to search for " |
---|
| 838 | "lines first, use find_lines"); |
---|
| 839 | */ |
---|
| 840 | std::vector<bool> res_mask(mask.nelements()); |
---|
| 841 | // iterator through lines |
---|
| 842 | std::list<std::pair<int,int> >::const_iterator cli=lines.begin(); |
---|
| 843 | for (int ch=0;ch<res_mask.size();++ch) |
---|
| 844 | if (ch<edge.first || ch>=edge.second) res_mask[ch]=false; |
---|
| 845 | else if (!mask[ch]) res_mask[ch]=false; |
---|
| 846 | else { |
---|
| 847 | res_mask[ch]=!invert; // no line by default |
---|
| 848 | if (cli==lines.end()) continue; |
---|
| 849 | if (ch>=cli->first && ch<cli->second) |
---|
| 850 | res_mask[ch]=invert; // this is a line |
---|
| 851 | if (ch>=cli->second) |
---|
| 852 | ++cli; // next line in the list |
---|
| 853 | } |
---|
| 854 | |
---|
| 855 | return res_mask; |
---|
| 856 | } |
---|
| 857 | catch (const AipsError &ae) { |
---|
| 858 | throw; |
---|
| 859 | } |
---|
| 860 | catch (const exception &ex) { |
---|
| 861 | throw AipsError(String("SDLineFinder::getMask - STL error: ")+ex.what()); |
---|
| 862 | } |
---|
| 863 | } |
---|
| 864 | |
---|
[370] | 865 | // get range for all lines found. The same units as used in the scan |
---|
| 866 | // will be returned (e.g. velocity instead of channels). |
---|
| 867 | std::vector<double> SDLineFinder::getLineRanges() |
---|
[297] | 868 | const throw(casa::AipsError) |
---|
| 869 | { |
---|
[370] | 870 | // convert to required abscissa units |
---|
| 871 | std::vector<double> vel=scan->getAbcissa(last_row_used); |
---|
| 872 | std::vector<int> ranges=getLineRangesInChannels(); |
---|
| 873 | std::vector<double> res(ranges.size()); |
---|
| 874 | |
---|
| 875 | std::vector<int>::const_iterator cri=ranges.begin(); |
---|
| 876 | std::vector<double>::iterator outi=res.begin(); |
---|
| 877 | for (;cri!=ranges.end() && outi!=res.end();++cri,++outi) |
---|
| 878 | if (uInt(*cri)>=vel.size()) |
---|
| 879 | throw AipsError("SDLineFinder::getLineRanges - getAbcissa provided less channels than reqired"); |
---|
| 880 | else *outi=vel[*cri]; |
---|
| 881 | return res; |
---|
| 882 | } |
---|
| 883 | |
---|
| 884 | // The same as getLineRanges, but channels are always used to specify |
---|
| 885 | // the range |
---|
| 886 | std::vector<int> SDLineFinder::getLineRangesInChannels() |
---|
| 887 | const throw(casa::AipsError) |
---|
| 888 | { |
---|
[297] | 889 | try { |
---|
| 890 | if (scan.null()) |
---|
[370] | 891 | throw AipsError("SDLineFinder::getLineRangesInChannels - a scan should be set first," |
---|
[297] | 892 | " use set_scan followed by find_lines"); |
---|
| 893 | DebugAssert(mask.nelements()==scan->nChan(), AipsError); |
---|
| 894 | |
---|
| 895 | if (!lines.size()) |
---|
[370] | 896 | throw AipsError("SDLineFinder::getLineRangesInChannels - one have to search for " |
---|
[297] | 897 | "lines first, use find_lines"); |
---|
| 898 | |
---|
| 899 | std::vector<int> res(2*lines.size()); |
---|
| 900 | // iterator through lines & result |
---|
| 901 | std::list<std::pair<int,int> >::const_iterator cli=lines.begin(); |
---|
| 902 | std::vector<int>::iterator ri=res.begin(); |
---|
| 903 | for (;cli!=lines.end() && ri!=res.end();++cli,++ri) { |
---|
| 904 | *ri=cli->first; |
---|
| 905 | if (++ri!=res.end()) |
---|
| 906 | *ri=cli->second-1; |
---|
[370] | 907 | } |
---|
[297] | 908 | return res; |
---|
| 909 | } |
---|
| 910 | catch (const AipsError &ae) { |
---|
| 911 | throw; |
---|
| 912 | } |
---|
| 913 | catch (const exception &ex) { |
---|
| 914 | throw AipsError(String("SDLineFinder::getLineRanges - STL error: ")+ex.what()); |
---|
| 915 | } |
---|
| 916 | } |
---|
[331] | 917 | |
---|
[370] | 918 | |
---|
| 919 | |
---|
[368] | 920 | // an auxiliary function to remove all lines from the list, except the |
---|
| 921 | // strongest one (by absolute value). If the lines removed are real, |
---|
| 922 | // they will be find again at the next iteration. This approach |
---|
| 923 | // increases the number of iterations required, but is able to remove |
---|
| 924 | // the sidelobes likely to occur near strong lines. |
---|
| 925 | // Later a better criterion may be implemented, e.g. |
---|
| 926 | // taking into consideration the brightness of different lines. Now |
---|
| 927 | // use the simplest solution |
---|
| 928 | // temp_mask - mask to work with (may be different from original mask as |
---|
| 929 | // the lines previously found may be masked) |
---|
| 930 | // lines2update - a list of lines to work with |
---|
| 931 | // nothing will be done if it is empty |
---|
| 932 | // max_box_nchan - channels in the running box for baseline filtering |
---|
| 933 | void SDLineFinder::keepStrongestOnly(const casa::Vector<casa::Bool> &temp_mask, |
---|
| 934 | std::list<std::pair<int, int> > &lines2update, |
---|
| 935 | int max_box_nchan) |
---|
| 936 | throw (casa::AipsError) |
---|
| 937 | { |
---|
| 938 | try { |
---|
| 939 | if (!lines2update.size()) return; // ignore an empty list |
---|
| 940 | |
---|
| 941 | // current line |
---|
| 942 | std::list<std::pair<int,int> >::iterator li=lines2update.begin(); |
---|
| 943 | // strongest line |
---|
| 944 | std::list<std::pair<int,int> >::iterator strongli=lines2update.begin(); |
---|
| 945 | // the flux (absolute value) of the strongest line |
---|
| 946 | Float peak_flux=-1; // negative value - a flag showing uninitialized |
---|
| 947 | // value |
---|
| 948 | // the algorithm below relies on the list being ordered |
---|
| 949 | Float tmp_flux=-1; // a temporary peak |
---|
| 950 | for (RunningBox running_box(spectrum,temp_mask,edge,max_box_nchan); |
---|
| 951 | running_box.haveMore(); running_box.next()) { |
---|
| 952 | |
---|
| 953 | if (li==lines2update.end()) break; // no more lines |
---|
| 954 | const int ch=running_box.getChannel(); |
---|
| 955 | if (ch>=li->first && ch<li->second) |
---|
| 956 | if (temp_mask[ch] && tmp_flux<fabs(running_box.aboveMean())) |
---|
| 957 | tmp_flux=fabs(running_box.aboveMean()); |
---|
| 958 | if (ch==li->second-1) { |
---|
| 959 | if (peak_flux<tmp_flux) { // if peak_flux=-1, this condition |
---|
| 960 | peak_flux=tmp_flux; // will be satisfied |
---|
| 961 | strongli=li; |
---|
| 962 | } |
---|
| 963 | ++li; |
---|
| 964 | tmp_flux=-1; |
---|
| 965 | } |
---|
| 966 | } |
---|
| 967 | std::list<std::pair<int,int> > res; |
---|
| 968 | res.splice(res.end(),lines2update,strongli); |
---|
| 969 | lines2update.clear(); |
---|
| 970 | lines2update.splice(lines2update.end(),res); |
---|
| 971 | } |
---|
| 972 | catch (const AipsError &ae) { |
---|
| 973 | throw; |
---|
| 974 | } |
---|
| 975 | catch (const exception &ex) { |
---|
| 976 | throw AipsError(String("SDLineFinder::keepStrongestOnly - STL error: ")+ex.what()); |
---|
| 977 | } |
---|
| 978 | |
---|
| 979 | } |
---|
| 980 | |
---|
[352] | 981 | // |
---|
| 982 | /////////////////////////////////////////////////////////////////////////////// |
---|
| 983 | |
---|
| 984 | |
---|
| 985 | /////////////////////////////////////////////////////////////////////////////// |
---|
| 986 | // |
---|
| 987 | // LFLineListOperations - a class incapsulating operations with line lists |
---|
| 988 | // The LF prefix stands for Line Finder |
---|
| 989 | // |
---|
| 990 | |
---|
[331] | 991 | // concatenate two lists preserving the order. If two lines appear to |
---|
| 992 | // be adjacent, they are joined into the new one |
---|
[352] | 993 | void LFLineListOperations::addNewSearchResult(const std::list<pair<int, int> > &newlines, |
---|
[344] | 994 | std::list<std::pair<int, int> > &lines_list) |
---|
[331] | 995 | throw(AipsError) |
---|
| 996 | { |
---|
| 997 | try { |
---|
| 998 | for (std::list<pair<int,int> >::const_iterator cli=newlines.begin(); |
---|
| 999 | cli!=newlines.end();++cli) { |
---|
[343] | 1000 | |
---|
| 1001 | // the first item, which has a non-void intersection or touches |
---|
| 1002 | // the new line |
---|
[344] | 1003 | std::list<pair<int,int> >::iterator pos_beg=find_if(lines_list.begin(), |
---|
| 1004 | lines_list.end(), IntersectsWith(*cli)); |
---|
[343] | 1005 | // the last such item |
---|
| 1006 | std::list<pair<int,int> >::iterator pos_end=find_if(pos_beg, |
---|
[344] | 1007 | lines_list.end(), not1(IntersectsWith(*cli))); |
---|
[343] | 1008 | |
---|
| 1009 | // extract all lines which intersect or touch a new one into |
---|
| 1010 | // a temporary buffer. This may invalidate the iterators |
---|
| 1011 | // line_buffer may be empty, if no lines intersects with a new |
---|
| 1012 | // one. |
---|
| 1013 | std::list<pair<int,int> > lines_buffer; |
---|
[344] | 1014 | lines_buffer.splice(lines_buffer.end(),lines_list, pos_beg, pos_end); |
---|
[343] | 1015 | |
---|
| 1016 | // build a union of all intersecting lines |
---|
| 1017 | pair<int,int> union_line=for_each(lines_buffer.begin(), |
---|
| 1018 | lines_buffer.end(),BuildUnion(*cli)).result(); |
---|
| 1019 | |
---|
| 1020 | // search for a right place for the new line (union_line) and add |
---|
[344] | 1021 | std::list<pair<int,int> >::iterator pos2insert=find_if(lines_list.begin(), |
---|
| 1022 | lines_list.end(), LaterThan(union_line)); |
---|
| 1023 | lines_list.insert(pos2insert,union_line); |
---|
[331] | 1024 | } |
---|
| 1025 | } |
---|
| 1026 | catch (const AipsError &ae) { |
---|
| 1027 | throw; |
---|
| 1028 | } |
---|
| 1029 | catch (const exception &ex) { |
---|
[352] | 1030 | throw AipsError(String("LFLineListOperations::addNewSearchResult - STL error: ")+ex.what()); |
---|
[331] | 1031 | } |
---|
| 1032 | } |
---|
[344] | 1033 | |
---|
| 1034 | // extend all line ranges to the point where a value stored in the |
---|
| 1035 | // specified vector changes (e.g. value-mean change its sign) |
---|
| 1036 | // This operation is necessary to include line wings, which are below |
---|
| 1037 | // the detection threshold. If lines becomes adjacent, they are |
---|
| 1038 | // merged together. Any masked channel stops the extension |
---|
[352] | 1039 | void LFLineListOperations::searchForWings(std::list<std::pair<int, int> > &newlines, |
---|
| 1040 | const casa::Vector<casa::Int> &signs, |
---|
| 1041 | const casa::Vector<casa::Bool> &mask, |
---|
| 1042 | const std::pair<int,int> &edge) throw(casa::AipsError) |
---|
[344] | 1043 | { |
---|
| 1044 | try { |
---|
| 1045 | for (std::list<pair<int,int> >::iterator li=newlines.begin(); |
---|
| 1046 | li!=newlines.end();++li) { |
---|
| 1047 | // update the left hand side |
---|
| 1048 | for (int n=li->first-1;n>=edge.first;--n) { |
---|
| 1049 | if (!mask[n]) break; |
---|
| 1050 | if (signs[n]==signs[li->first] && signs[li->first]) |
---|
| 1051 | li->first=n; |
---|
| 1052 | else break; |
---|
| 1053 | } |
---|
| 1054 | // update the right hand side |
---|
| 1055 | for (int n=li->second;n<edge.second;++n) { |
---|
| 1056 | if (!mask[n]) break; |
---|
| 1057 | if (signs[n]==signs[li->second-1] && signs[li->second-1]) |
---|
| 1058 | li->second=n; |
---|
| 1059 | else break; |
---|
| 1060 | } |
---|
| 1061 | } |
---|
| 1062 | // need to search for possible mergers. |
---|
| 1063 | std::list<std::pair<int, int> > result_buffer; |
---|
| 1064 | addNewSearchResult(newlines,result_buffer); |
---|
| 1065 | newlines.clear(); |
---|
| 1066 | newlines.splice(newlines.end(),result_buffer); |
---|
| 1067 | } |
---|
| 1068 | catch (const AipsError &ae) { |
---|
| 1069 | throw; |
---|
| 1070 | } |
---|
| 1071 | catch (const exception &ex) { |
---|
[352] | 1072 | throw AipsError(String("LFLineListOperations::extendLines - STL error: ")+ex.what()); |
---|
[344] | 1073 | } |
---|
| 1074 | } |
---|
[352] | 1075 | |
---|
| 1076 | // |
---|
| 1077 | /////////////////////////////////////////////////////////////////////////////// |
---|