[297] | 1 | //#--------------------------------------------------------------------------- |
---|
[881] | 2 | //# STLineFinder.cc: A class for automated spectral line search |
---|
[297] | 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 | //# |
---|
[881] | 29 | //# $Id:$ |
---|
[297] | 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 | // |
---|
[881] | 52 | // RunningBox - a running box calculator. This class implements |
---|
[351] | 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 |
---|
[881] | 59 | // an unnecessary copying |
---|
[331] | 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 |
---|
[881] | 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) |
---|
[881] | 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); |
---|
[881] | 92 | |
---|
[351] | 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(); |
---|
[881] | 98 | |
---|
[351] | 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); |
---|
[881] | 111 | |
---|
[351] | 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 |
---|
[881] | 144 | casa::Float threshold; // detection threshold - the |
---|
[331] | 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(); |
---|
[881] | 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 | // |
---|
[881] | 200 | // RunningBox - a running box calculator. This class implements |
---|
[351] | 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); |
---|
[881] | 229 | |
---|
| 230 | if (initial_box_ch==edge.second) |
---|
[351] | 231 | throw AipsError("RunningBox::rewind - too much channels are masked"); |
---|
| 232 | |
---|
| 233 | cur_channel=edge.first; |
---|
[881] | 234 | start_advance=initial_box_ch-max_box_nchan/2; |
---|
[351] | 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]; |
---|
[881] | 293 | sumf2-=square(spectrum[ch2remove]); |
---|
[351] | 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); |
---|
[881] | 321 | |
---|
[351] | 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 | } |
---|
[881] | 387 | } else processCurLine(mask); |
---|
[297] | 388 | } |
---|
| 389 | catch (const AipsError &ae) { |
---|
| 390 | throw; |
---|
[881] | 391 | } |
---|
[297] | 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 { |
---|
[881] | 403 | if (is_detected_before) { |
---|
[297] | 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; |
---|
[881] | 407 | if (lines.size()) { |
---|
[297] | 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 | } |
---|
[881] | 414 | } else add_new_line=True; |
---|
| 415 | if (add_new_line) |
---|
[297] | 416 | lines.push_back(cur_line); |
---|
[881] | 417 | else lines.back().second=cur_line.second; |
---|
[297] | 418 | } |
---|
| 419 | is_detected_before=False; |
---|
[881] | 420 | } |
---|
[297] | 421 | } |
---|
| 422 | catch (const AipsError &ae) { |
---|
| 423 | throw; |
---|
[881] | 424 | } |
---|
[297] | 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 | |
---|
[881] | 452 | |
---|
[368] | 453 | // determine the off-line variance first |
---|
| 454 | // an assumption made: lines occupy a small part of the spectrum |
---|
[881] | 455 | |
---|
[368] | 456 | std::vector<float> variances(edge.second-edge.first); |
---|
| 457 | DebugAssert(variances.size(),AipsError); |
---|
[881] | 458 | |
---|
| 459 | for (;running_box->haveMore();running_box->next()) |
---|
[368] | 460 | variances[running_box->getChannel()-edge.first]= |
---|
| 461 | running_box->getLinVariance(); |
---|
[881] | 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()); |
---|
[881] | 468 | |
---|
[368] | 469 | Float offline_variance=0; |
---|
[881] | 470 | uInt offline_cnt=uInt(0.8*variances.size()); |
---|
[368] | 471 | if (!offline_cnt) offline_cnt=variances.size(); // no much else left, |
---|
| 472 | // although it is very inaccurate |
---|
[881] | 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; |
---|
[881] | 504 | } |
---|
[351] | 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 |
---|
[881] | 585 | |
---|
[343] | 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 | // |
---|
[881] | 597 | // STLineFinder - a class for automated spectral line search |
---|
[343] | 598 | // |
---|
| 599 | // |
---|
[331] | 600 | |
---|
[881] | 601 | STLineFinder::STLineFinder() throw() : edge(0,0) |
---|
[331] | 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 |
---|
[881] | 616 | // Default is 8, but for a bad baseline shape this |
---|
[369] | 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. |
---|
[881] | 620 | // Setting a very large value doesn't usually provide |
---|
| 621 | // valid detections. |
---|
[369] | 622 | // in_box_size the box size for running mean calculation. Default is |
---|
| 623 | // 1./5. of the whole spectrum size |
---|
[881] | 624 | void STLineFinder::setOptions(const casa::Float &in_threshold, |
---|
[369] | 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 | |
---|
[881] | 635 | STLineFinder::~STLineFinder() throw(AipsError) {} |
---|
[331] | 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 |
---|
[881] | 643 | // can be achieved using a spectrum mask only |
---|
| 644 | void STLineFinder::setScan(const ScantableWrapper &in_scan, |
---|
[331] | 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; |
---|
[881] | 653 | if (mask.nelements()!=scan->nchan()) |
---|
| 654 | throw AipsError("STLineFinder::setScan - in_scan and in_mask have different" |
---|
[331] | 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) |
---|
[881] | 660 | throw AipsError("STLineFinder::setScan - the length of the in_edge parameter" |
---|
[331] | 661 | "should not exceed 2"); |
---|
| 662 | if (!n) { |
---|
[881] | 663 | // all spectra, no rejection |
---|
[331] | 664 | edge.first=0; |
---|
[881] | 665 | edge.second=scan->nchan(); |
---|
[331] | 666 | } else { |
---|
| 667 | edge.first=extract<int>(in_edge.attr("__getitem__")(0)); |
---|
| 668 | if (edge.first<0) |
---|
[881] | 669 | throw AipsError("STLineFinder::setScan - the in_edge parameter has a negative" |
---|
[331] | 670 | "number of channels to drop"); |
---|
[881] | 671 | if (edge.first>=scan->nchan()) |
---|
| 672 | throw AipsError("STLineFinder::setScan - all channels are rejected by the in_edge parameter"); |
---|
[331] | 673 | if (n==2) { |
---|
| 674 | edge.second=extract<int>(in_edge.attr("__getitem__")(1)); |
---|
| 675 | if (edge.second<0) |
---|
[881] | 676 | throw AipsError("STLineFinder::setScan - the in_edge parameter has a negative" |
---|
[331] | 677 | "number of channels to drop"); |
---|
[881] | 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)) |
---|
[881] | 681 | throw AipsError("STLineFinder::setScan - all channels are rejected by the in_edge parameter"); |
---|
| 682 | } |
---|
[331] | 683 | } |
---|
| 684 | catch(const AipsError &ae) { |
---|
| 685 | // setScan is unsuccessfull, reset scan/mask/edge |
---|
[881] | 686 | scan=CountedConstPtr<Scantable>(); // null pointer |
---|
[331] | 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 |
---|
[881] | 694 | int STLineFinder::findLines(const casa::uInt &whichRow) throw(casa::AipsError) |
---|
| 695 | { |
---|
[331] | 696 | const int minboxnchan=4; |
---|
| 697 | if (scan.null()) |
---|
[881] | 698 | throw AipsError("STLineFinder::findLines - a scan should be set first," |
---|
[331] | 699 | " use set_scan"); |
---|
[881] | 700 | DebugAssert(mask.nelements()==scan->nchan(), AipsError); |
---|
| 701 | int max_box_nchan=int(scan->nchan()*box_size); // number of channels in running |
---|
[331] | 702 | // box |
---|
| 703 | if (max_box_nchan<2) |
---|
[881] | 704 | throw AipsError("STLineFinder::findLines - box_size is too small"); |
---|
[331] | 705 | |
---|
[881] | 706 | spectrum.resize(); |
---|
| 707 | spectrum = Vector<Float>(scan->getSpectrum(whichRow)); |
---|
[331] | 708 | |
---|
| 709 | lines.resize(0); // search from the scratch |
---|
[370] | 710 | last_row_used=whichRow; |
---|
[331] | 711 | Vector<Bool> temp_mask(mask); |
---|
[351] | 712 | |
---|
| 713 | Bool first_pass=True; |
---|
[368] | 714 | Int avg_factor=1; // this number of adjacent channels is averaged together |
---|
| 715 | // the total number of the channels is not altered |
---|
| 716 | // instead, min_nchan is also scaled |
---|
| 717 | // it helps to search for broad lines |
---|
[551] | 718 | Vector<Int> signs; // a buffer for signs of the value - mean quantity |
---|
| 719 | // see LFAboveThreshold for details |
---|
| 720 | // We need only signs resulted from last iteration |
---|
| 721 | // because all previous values may be corrupted by the |
---|
| 722 | // presence of spectral lines |
---|
[344] | 723 | while (true) { |
---|
[351] | 724 | // a buffer for new lines found at this iteration |
---|
[881] | 725 | std::list<pair<int,int> > new_lines; |
---|
[351] | 726 | |
---|
| 727 | try { |
---|
[369] | 728 | // line find algorithm |
---|
| 729 | LFAboveThreshold lfalg(new_lines,avg_factor*min_nchan, threshold); |
---|
[352] | 730 | lfalg.findLines(spectrum,temp_mask,edge,max_box_nchan); |
---|
[551] | 731 | signs.resize(lfalg.getSigns().nelements()); |
---|
| 732 | signs=lfalg.getSigns(); |
---|
[368] | 733 | first_pass=False; |
---|
| 734 | if (!new_lines.size()) |
---|
| 735 | throw AipsError("spurious"); // nothing new - use the same |
---|
| 736 | // code as for a real exception |
---|
[351] | 737 | } |
---|
| 738 | catch(const AipsError &ae) { |
---|
| 739 | if (first_pass) throw; |
---|
[368] | 740 | // nothing new - proceed to the next step of averaging, if any |
---|
| 741 | // (to search for broad lines) |
---|
[369] | 742 | if (avg_factor>avg_limit) break; // averaging up to avg_limit |
---|
| 743 | // adjacent channels, |
---|
| 744 | // stop after that |
---|
[368] | 745 | avg_factor*=2; // twice as more averaging |
---|
[369] | 746 | subtractBaseline(temp_mask,9); |
---|
[368] | 747 | averageAdjacentChannels(temp_mask,avg_factor); |
---|
[881] | 748 | continue; |
---|
[351] | 749 | } |
---|
[368] | 750 | keepStrongestOnly(temp_mask,new_lines,max_box_nchan); |
---|
[343] | 751 | // update the list (lines) merging intervals, if necessary |
---|
[344] | 752 | addNewSearchResult(new_lines,lines); |
---|
| 753 | // get a new mask |
---|
[881] | 754 | temp_mask=getMask(); |
---|
[343] | 755 | } |
---|
[881] | 756 | |
---|
[551] | 757 | // an additional search for wings because in the presence of very strong |
---|
| 758 | // lines temporary mean used at each iteration will be higher than |
---|
| 759 | // the true mean |
---|
[881] | 760 | |
---|
[551] | 761 | if (lines.size()) |
---|
| 762 | LFLineListOperations::searchForWings(lines,signs,mask,edge); |
---|
[881] | 763 | |
---|
[331] | 764 | return int(lines.size()); |
---|
| 765 | } |
---|
| 766 | |
---|
[369] | 767 | // auxiliary function to fit and subtract a polynomial from the current |
---|
| 768 | // spectrum. It uses the SDFitter class. This action is required before |
---|
| 769 | // reducing the spectral resolution if the baseline shape is bad |
---|
[881] | 770 | void STLineFinder::subtractBaseline(const casa::Vector<casa::Bool> &temp_mask, |
---|
[369] | 771 | const casa::Int &order) throw(casa::AipsError) |
---|
| 772 | { |
---|
| 773 | AlwaysAssert(spectrum.nelements(),AipsError); |
---|
| 774 | // use the fact that temp_mask excludes channels rejected at the edge |
---|
| 775 | SDFitter sdf; |
---|
| 776 | std::vector<float> absc(spectrum.nelements()); |
---|
| 777 | for (Int i=0;i<absc.size();++i) |
---|
| 778 | absc[i]=float(i)/float(spectrum.nelements()); |
---|
| 779 | std::vector<float> spec; |
---|
| 780 | spectrum.tovector(spec); |
---|
| 781 | std::vector<bool> std_mask; |
---|
| 782 | temp_mask.tovector(std_mask); |
---|
| 783 | sdf.setData(absc,spec,std_mask); |
---|
| 784 | sdf.setExpression("poly",order); |
---|
| 785 | if (!sdf.fit()) return; // fit failed, use old spectrum |
---|
[881] | 786 | spectrum=casa::Vector<casa::Float>(sdf.getResidual()); |
---|
[369] | 787 | } |
---|
| 788 | |
---|
[368] | 789 | // auxiliary function to average adjacent channels and update the mask |
---|
| 790 | // if at least one channel involved in summation is masked, all |
---|
| 791 | // output channels will be masked. This function works with the |
---|
| 792 | // spectrum and edge fields of this class, but updates the mask |
---|
| 793 | // array specified, rather than the field of this class |
---|
| 794 | // boxsize - a number of adjacent channels to average |
---|
[881] | 795 | void STLineFinder::averageAdjacentChannels(casa::Vector<casa::Bool> &mask2update, |
---|
[368] | 796 | const casa::Int &boxsize) |
---|
| 797 | throw(casa::AipsError) |
---|
| 798 | { |
---|
| 799 | DebugAssert(mask2update.nelements()==spectrum.nelements(), AipsError); |
---|
| 800 | DebugAssert(boxsize!=0,AipsError); |
---|
[881] | 801 | |
---|
[368] | 802 | for (int n=edge.first;n<edge.second;n+=boxsize) { |
---|
| 803 | DebugAssert(n<spectrum.nelements(),AipsError); |
---|
| 804 | int nboxch=0; // number of channels currently in the box |
---|
| 805 | Float mean=0; // buffer for mean calculations |
---|
| 806 | for (int k=n;k<n+boxsize && k<edge.second;++k) |
---|
| 807 | if (mask2update[k]) { // k is a valid channel |
---|
| 808 | mean+=spectrum[k]; |
---|
| 809 | ++nboxch; |
---|
[881] | 810 | } |
---|
[368] | 811 | if (nboxch<boxsize) // mask these channels |
---|
| 812 | for (int k=n;k<n+boxsize && k<edge.second;++k) |
---|
| 813 | mask2update[k]=False; |
---|
| 814 | else { |
---|
| 815 | mean/=Float(boxsize); |
---|
| 816 | for (int k=n;k<n+boxsize && k<edge.second;++k) |
---|
| 817 | spectrum[k]=mean; |
---|
| 818 | } |
---|
| 819 | } |
---|
| 820 | } |
---|
[331] | 821 | |
---|
[368] | 822 | |
---|
[297] | 823 | // get the mask to mask out all lines that have been found (default) |
---|
| 824 | // if invert=true, only channels belong to lines will be unmasked |
---|
| 825 | // Note: all channels originally masked by the input mask (in_mask |
---|
| 826 | // in setScan) or dropped out by the edge parameter (in_edge |
---|
| 827 | // in setScan) are still excluded regardless on the invert option |
---|
[881] | 828 | std::vector<bool> STLineFinder::getMask(bool invert) |
---|
[297] | 829 | const throw(casa::AipsError) |
---|
| 830 | { |
---|
| 831 | try { |
---|
| 832 | if (scan.null()) |
---|
[881] | 833 | throw AipsError("STLineFinder::getMask - a scan should be set first," |
---|
[297] | 834 | " use set_scan followed by find_lines"); |
---|
[881] | 835 | DebugAssert(mask.nelements()==scan->nchan(), AipsError); |
---|
[297] | 836 | /* |
---|
| 837 | if (!lines.size()) |
---|
[881] | 838 | throw AipsError("STLineFinder::getMask - one have to search for " |
---|
[297] | 839 | "lines first, use find_lines"); |
---|
[881] | 840 | */ |
---|
[297] | 841 | std::vector<bool> res_mask(mask.nelements()); |
---|
| 842 | // iterator through lines |
---|
| 843 | std::list<std::pair<int,int> >::const_iterator cli=lines.begin(); |
---|
[881] | 844 | for (int ch=0;ch<res_mask.size();++ch) |
---|
[297] | 845 | if (ch<edge.first || ch>=edge.second) res_mask[ch]=false; |
---|
| 846 | else if (!mask[ch]) res_mask[ch]=false; |
---|
[881] | 847 | else { |
---|
[297] | 848 | res_mask[ch]=!invert; // no line by default |
---|
| 849 | if (cli==lines.end()) continue; |
---|
| 850 | if (ch>=cli->first && ch<cli->second) |
---|
| 851 | res_mask[ch]=invert; // this is a line |
---|
| 852 | if (ch>=cli->second) |
---|
| 853 | ++cli; // next line in the list |
---|
| 854 | } |
---|
[881] | 855 | |
---|
[297] | 856 | return res_mask; |
---|
| 857 | } |
---|
| 858 | catch (const AipsError &ae) { |
---|
| 859 | throw; |
---|
[881] | 860 | } |
---|
[297] | 861 | catch (const exception &ex) { |
---|
[881] | 862 | throw AipsError(String("STLineFinder::getMask - STL error: ")+ex.what()); |
---|
[297] | 863 | } |
---|
| 864 | } |
---|
| 865 | |
---|
[370] | 866 | // get range for all lines found. The same units as used in the scan |
---|
| 867 | // will be returned (e.g. velocity instead of channels). |
---|
[881] | 868 | std::vector<double> STLineFinder::getLineRanges() |
---|
[297] | 869 | const throw(casa::AipsError) |
---|
| 870 | { |
---|
[370] | 871 | // convert to required abscissa units |
---|
| 872 | std::vector<double> vel=scan->getAbcissa(last_row_used); |
---|
| 873 | std::vector<int> ranges=getLineRangesInChannels(); |
---|
| 874 | std::vector<double> res(ranges.size()); |
---|
| 875 | |
---|
| 876 | std::vector<int>::const_iterator cri=ranges.begin(); |
---|
| 877 | std::vector<double>::iterator outi=res.begin(); |
---|
| 878 | for (;cri!=ranges.end() && outi!=res.end();++cri,++outi) |
---|
| 879 | if (uInt(*cri)>=vel.size()) |
---|
[881] | 880 | throw AipsError("STLineFinder::getLineRanges - getAbcissa provided less channels than reqired"); |
---|
[370] | 881 | else *outi=vel[*cri]; |
---|
| 882 | return res; |
---|
| 883 | } |
---|
| 884 | |
---|
| 885 | // The same as getLineRanges, but channels are always used to specify |
---|
| 886 | // the range |
---|
[881] | 887 | std::vector<int> STLineFinder::getLineRangesInChannels() |
---|
[370] | 888 | const throw(casa::AipsError) |
---|
| 889 | { |
---|
[297] | 890 | try { |
---|
| 891 | if (scan.null()) |
---|
[881] | 892 | throw AipsError("STLineFinder::getLineRangesInChannels - a scan should be set first," |
---|
[297] | 893 | " use set_scan followed by find_lines"); |
---|
[881] | 894 | DebugAssert(mask.nelements()==scan->nchan(), AipsError); |
---|
| 895 | |
---|
[297] | 896 | if (!lines.size()) |
---|
[881] | 897 | throw AipsError("STLineFinder::getLineRangesInChannels - one have to search for " |
---|
[297] | 898 | "lines first, use find_lines"); |
---|
[881] | 899 | |
---|
[297] | 900 | std::vector<int> res(2*lines.size()); |
---|
| 901 | // iterator through lines & result |
---|
| 902 | std::list<std::pair<int,int> >::const_iterator cli=lines.begin(); |
---|
| 903 | std::vector<int>::iterator ri=res.begin(); |
---|
[881] | 904 | for (;cli!=lines.end() && ri!=res.end();++cli,++ri) { |
---|
[297] | 905 | *ri=cli->first; |
---|
[881] | 906 | if (++ri!=res.end()) |
---|
| 907 | *ri=cli->second-1; |
---|
| 908 | } |
---|
[297] | 909 | return res; |
---|
| 910 | } |
---|
| 911 | catch (const AipsError &ae) { |
---|
| 912 | throw; |
---|
[881] | 913 | } |
---|
[297] | 914 | catch (const exception &ex) { |
---|
[881] | 915 | throw AipsError(String("STLineFinder::getLineRanges - STL error: ")+ex.what()); |
---|
[297] | 916 | } |
---|
| 917 | } |
---|
[331] | 918 | |
---|
[370] | 919 | |
---|
| 920 | |
---|
[368] | 921 | // an auxiliary function to remove all lines from the list, except the |
---|
| 922 | // strongest one (by absolute value). If the lines removed are real, |
---|
[881] | 923 | // they will be find again at the next iteration. This approach |
---|
| 924 | // increases the number of iterations required, but is able to remove |
---|
[368] | 925 | // the sidelobes likely to occur near strong lines. |
---|
| 926 | // Later a better criterion may be implemented, e.g. |
---|
| 927 | // taking into consideration the brightness of different lines. Now |
---|
[881] | 928 | // use the simplest solution |
---|
[368] | 929 | // temp_mask - mask to work with (may be different from original mask as |
---|
| 930 | // the lines previously found may be masked) |
---|
| 931 | // lines2update - a list of lines to work with |
---|
| 932 | // nothing will be done if it is empty |
---|
| 933 | // max_box_nchan - channels in the running box for baseline filtering |
---|
[881] | 934 | void STLineFinder::keepStrongestOnly(const casa::Vector<casa::Bool> &temp_mask, |
---|
[368] | 935 | std::list<std::pair<int, int> > &lines2update, |
---|
| 936 | int max_box_nchan) |
---|
| 937 | throw (casa::AipsError) |
---|
| 938 | { |
---|
| 939 | try { |
---|
| 940 | if (!lines2update.size()) return; // ignore an empty list |
---|
| 941 | |
---|
| 942 | // current line |
---|
| 943 | std::list<std::pair<int,int> >::iterator li=lines2update.begin(); |
---|
| 944 | // strongest line |
---|
| 945 | std::list<std::pair<int,int> >::iterator strongli=lines2update.begin(); |
---|
| 946 | // the flux (absolute value) of the strongest line |
---|
| 947 | Float peak_flux=-1; // negative value - a flag showing uninitialized |
---|
| 948 | // value |
---|
| 949 | // the algorithm below relies on the list being ordered |
---|
| 950 | Float tmp_flux=-1; // a temporary peak |
---|
| 951 | for (RunningBox running_box(spectrum,temp_mask,edge,max_box_nchan); |
---|
| 952 | running_box.haveMore(); running_box.next()) { |
---|
| 953 | |
---|
| 954 | if (li==lines2update.end()) break; // no more lines |
---|
[881] | 955 | const int ch=running_box.getChannel(); |
---|
[368] | 956 | if (ch>=li->first && ch<li->second) |
---|
| 957 | if (temp_mask[ch] && tmp_flux<fabs(running_box.aboveMean())) |
---|
| 958 | tmp_flux=fabs(running_box.aboveMean()); |
---|
| 959 | if (ch==li->second-1) { |
---|
| 960 | if (peak_flux<tmp_flux) { // if peak_flux=-1, this condition |
---|
| 961 | peak_flux=tmp_flux; // will be satisfied |
---|
| 962 | strongli=li; |
---|
| 963 | } |
---|
| 964 | ++li; |
---|
| 965 | tmp_flux=-1; |
---|
| 966 | } |
---|
[881] | 967 | } |
---|
[368] | 968 | std::list<std::pair<int,int> > res; |
---|
| 969 | res.splice(res.end(),lines2update,strongli); |
---|
| 970 | lines2update.clear(); |
---|
| 971 | lines2update.splice(lines2update.end(),res); |
---|
| 972 | } |
---|
| 973 | catch (const AipsError &ae) { |
---|
| 974 | throw; |
---|
[881] | 975 | } |
---|
[368] | 976 | catch (const exception &ex) { |
---|
[881] | 977 | throw AipsError(String("STLineFinder::keepStrongestOnly - STL error: ")+ex.what()); |
---|
[368] | 978 | } |
---|
| 979 | |
---|
| 980 | } |
---|
| 981 | |
---|
[352] | 982 | // |
---|
| 983 | /////////////////////////////////////////////////////////////////////////////// |
---|
| 984 | |
---|
| 985 | |
---|
| 986 | /////////////////////////////////////////////////////////////////////////////// |
---|
| 987 | // |
---|
| 988 | // LFLineListOperations - a class incapsulating operations with line lists |
---|
| 989 | // The LF prefix stands for Line Finder |
---|
| 990 | // |
---|
| 991 | |
---|
[331] | 992 | // concatenate two lists preserving the order. If two lines appear to |
---|
| 993 | // be adjacent, they are joined into the new one |
---|
[352] | 994 | void LFLineListOperations::addNewSearchResult(const std::list<pair<int, int> > &newlines, |
---|
[881] | 995 | std::list<std::pair<int, int> > &lines_list) |
---|
[331] | 996 | throw(AipsError) |
---|
| 997 | { |
---|
| 998 | try { |
---|
| 999 | for (std::list<pair<int,int> >::const_iterator cli=newlines.begin(); |
---|
| 1000 | cli!=newlines.end();++cli) { |
---|
[881] | 1001 | |
---|
[343] | 1002 | // the first item, which has a non-void intersection or touches |
---|
| 1003 | // the new line |
---|
[344] | 1004 | std::list<pair<int,int> >::iterator pos_beg=find_if(lines_list.begin(), |
---|
[881] | 1005 | lines_list.end(), IntersectsWith(*cli)); |
---|
| 1006 | // the last such item |
---|
[343] | 1007 | std::list<pair<int,int> >::iterator pos_end=find_if(pos_beg, |
---|
[344] | 1008 | lines_list.end(), not1(IntersectsWith(*cli))); |
---|
[343] | 1009 | |
---|
| 1010 | // extract all lines which intersect or touch a new one into |
---|
| 1011 | // a temporary buffer. This may invalidate the iterators |
---|
| 1012 | // line_buffer may be empty, if no lines intersects with a new |
---|
| 1013 | // one. |
---|
| 1014 | std::list<pair<int,int> > lines_buffer; |
---|
[344] | 1015 | lines_buffer.splice(lines_buffer.end(),lines_list, pos_beg, pos_end); |
---|
[343] | 1016 | |
---|
[881] | 1017 | // build a union of all intersecting lines |
---|
[343] | 1018 | pair<int,int> union_line=for_each(lines_buffer.begin(), |
---|
| 1019 | lines_buffer.end(),BuildUnion(*cli)).result(); |
---|
[881] | 1020 | |
---|
[343] | 1021 | // search for a right place for the new line (union_line) and add |
---|
[344] | 1022 | std::list<pair<int,int> >::iterator pos2insert=find_if(lines_list.begin(), |
---|
| 1023 | lines_list.end(), LaterThan(union_line)); |
---|
| 1024 | lines_list.insert(pos2insert,union_line); |
---|
[331] | 1025 | } |
---|
| 1026 | } |
---|
| 1027 | catch (const AipsError &ae) { |
---|
| 1028 | throw; |
---|
[881] | 1029 | } |
---|
[331] | 1030 | catch (const exception &ex) { |
---|
[352] | 1031 | throw AipsError(String("LFLineListOperations::addNewSearchResult - STL error: ")+ex.what()); |
---|
[331] | 1032 | } |
---|
| 1033 | } |
---|
[344] | 1034 | |
---|
| 1035 | // extend all line ranges to the point where a value stored in the |
---|
| 1036 | // specified vector changes (e.g. value-mean change its sign) |
---|
| 1037 | // This operation is necessary to include line wings, which are below |
---|
| 1038 | // the detection threshold. If lines becomes adjacent, they are |
---|
| 1039 | // merged together. Any masked channel stops the extension |
---|
[352] | 1040 | void LFLineListOperations::searchForWings(std::list<std::pair<int, int> > &newlines, |
---|
| 1041 | const casa::Vector<casa::Int> &signs, |
---|
| 1042 | const casa::Vector<casa::Bool> &mask, |
---|
| 1043 | const std::pair<int,int> &edge) throw(casa::AipsError) |
---|
[344] | 1044 | { |
---|
| 1045 | try { |
---|
| 1046 | for (std::list<pair<int,int> >::iterator li=newlines.begin(); |
---|
| 1047 | li!=newlines.end();++li) { |
---|
| 1048 | // update the left hand side |
---|
| 1049 | for (int n=li->first-1;n>=edge.first;--n) { |
---|
| 1050 | if (!mask[n]) break; |
---|
| 1051 | if (signs[n]==signs[li->first] && signs[li->first]) |
---|
| 1052 | li->first=n; |
---|
[881] | 1053 | else break; |
---|
[344] | 1054 | } |
---|
| 1055 | // update the right hand side |
---|
| 1056 | for (int n=li->second;n<edge.second;++n) { |
---|
| 1057 | if (!mask[n]) break; |
---|
| 1058 | if (signs[n]==signs[li->second-1] && signs[li->second-1]) |
---|
| 1059 | li->second=n; |
---|
[881] | 1060 | else break; |
---|
[344] | 1061 | } |
---|
| 1062 | } |
---|
| 1063 | // need to search for possible mergers. |
---|
| 1064 | std::list<std::pair<int, int> > result_buffer; |
---|
| 1065 | addNewSearchResult(newlines,result_buffer); |
---|
| 1066 | newlines.clear(); |
---|
| 1067 | newlines.splice(newlines.end(),result_buffer); |
---|
| 1068 | } |
---|
| 1069 | catch (const AipsError &ae) { |
---|
| 1070 | throw; |
---|
[881] | 1071 | } |
---|
[344] | 1072 | catch (const exception &ex) { |
---|
[352] | 1073 | throw AipsError(String("LFLineListOperations::extendLines - STL error: ")+ex.what()); |
---|
[344] | 1074 | } |
---|
| 1075 | } |
---|
[352] | 1076 | |
---|
| 1077 | // |
---|
| 1078 | /////////////////////////////////////////////////////////////////////////////// |
---|