source: trunk/src/SDLineFinder.cc@ 364

Last change on this file since 364 was 352, checked in by vor010, 20 years ago

SDLineFinder - further interface changes to achieve a better structure of the program

  • Property svn:eol-style set to native
  • Property svn:keywords set to Author Date Id Revision
File size: 30.6 KB
Line 
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"
35
36// STL
37#include <functional>
38#include <algorithm>
39#include <iostream>
40#include <fstream>
41
42using namespace asap;
43using namespace casa;
44using namespace std;
45using namespace boost::python;
46
47namespace asap {
48
49///////////////////////////////////////////////////////////////////////////////
50//
51// RunningBox - a running box calculator. This class implements
52// interations over the specified spectrum and calculates
53// running box filter statistics.
54//
55
56class RunningBox {
57 // The input data to work with. Use reference symantics to avoid
58 // an unnecessary copying
59 const casa::Vector<casa::Float> &spectrum; // a buffer for the spectrum
60 const casa::Vector<casa::Bool> &mask; // associated mask
61 const std::pair<int,int> &edge; // start and stop+1 channels
62 // to work with
63
64 // statistics for running box filtering
65 casa::Float sumf; // sum of fluxes
66 casa::Float sumf2; // sum of squares of fluxes
67 casa::Float sumch; // sum of channel numbers (for linear fit)
68 casa::Float sumch2; // sum of squares of channel numbers (for linear fit)
69 casa::Float sumfch; // sum of flux*(channel number) (for linear fit)
70
71 int box_chan_cntr; // actual number of channels in the box
72 int max_box_nchan; // maximum allowed number of channels in the box
73 // (calculated from boxsize and actual spectrum size)
74 // cache for derivative statistics
75 mutable casa::Bool need2recalculate; // if true, values of the statistics
76 // below are invalid
77 mutable casa::Float linmean; // a value of the linear fit to the
78 // points in the running box
79 mutable casa::Float linvariance; // the same for variance
80 int cur_channel; // the number of the current channel
81 int start_advance; // number of channel from which the box can
82 // be moved (the middle of the box, if there is no
83 // masking)
84public:
85 // set up the object with the references to actual data
86 // as well as the number of channels in the running box
87 RunningBox(const casa::Vector<casa::Float> &in_spectrum,
88 const casa::Vector<casa::Bool> &in_mask,
89 const std::pair<int,int> &in_edge,
90 int in_max_box_nchan) throw(AipsError);
91
92 // access to the statistics
93 const casa::Float& getLinMean() const throw(AipsError);
94 const casa::Float& getLinVariance() const throw(AipsError);
95 const casa::Float aboveMean() const throw(AipsError);
96 int getChannel() const throw();
97
98 // actual number of channels in the box (max_box_nchan, if no channels
99 // are masked)
100 int getNumberOfBoxPoints() const throw();
101
102 // next channel
103 void next() throw(AipsError);
104
105 // checking whether there are still elements
106 casa::Bool haveMore() const throw();
107
108 // go to start
109 void rewind() throw(AipsError);
110
111protected:
112 // supplementary function to control running mean calculations.
113 // It adds a specified channel to the running mean box and
114 // removes (ch-maxboxnchan+1)'th channel from there
115 // Channels, for which the mask is false or index is beyond the
116 // allowed range, are ignored
117 void advanceRunningBox(int ch) throw(casa::AipsError);
118
119 // calculate derivative statistics. This function is const, because
120 // it updates the cache only
121 void updateDerivativeStatistics() const throw(AipsError);
122};
123
124//
125///////////////////////////////////////////////////////////////////////////////
126
127///////////////////////////////////////////////////////////////////////////////
128//
129// LFAboveThreshold An algorithm for line detection using running box
130// statistics. Line is detected if it is above the
131// specified threshold at the specified number of
132// consequtive channels. Prefix LF stands for Line Finder
133//
134class LFAboveThreshold : protected LFLineListOperations {
135 // temporary line edge channels and flag, which is True if the line
136 // was detected in the previous channels.
137 std::pair<int,int> cur_line;
138 casa::Bool is_detected_before;
139 int min_nchan; // A minimum number of consequtive
140 // channels, which should satisfy
141 // the detection criterion, to be
142 // a detection
143 casa::Float threshold; // detection threshold - the
144 // minimal signal to noise ratio
145 std::list<pair<int,int> > &lines; // list where detections are saved
146 // (pair: start and stop+1 channel)
147 RunningBox *running_box; // running box filter
148public:
149
150 // set up the detection criterion
151 LFAboveThreshold(std::list<pair<int,int> > &in_lines,
152 int in_min_nchan = 3,
153 casa::Float in_threshold = 5) throw();
154 virtual ~LFAboveThreshold() throw();
155
156 // replace the detection criterion
157 void setCriterion(int in_min_nchan, casa::Float in_threshold) throw();
158
159 // find spectral lines and add them into list
160 // if statholder is not NULL, the accumulate function of it will be
161 // called for each channel to save statistics
162 // spectrum, mask and edge - reference to the data
163 // max_box_nchan - number of channels in the running box
164 void findLines(const casa::Vector<casa::Float> &spectrum,
165 const casa::Vector<casa::Bool> &mask,
166 const std::pair<int,int> &edge,
167 int max_box_nchan) throw(casa::AipsError);
168
169protected:
170
171 // process a channel: update curline and is_detected before and
172 // add a new line to the list, if necessary using processCurLine()
173 // detect=true indicates that the current channel satisfies the criterion
174 void processChannel(Bool detect, const casa::Vector<casa::Bool> &mask)
175 throw(casa::AipsError);
176
177 // process the interval of channels stored in curline
178 // if it satisfies the criterion, add this interval as a new line
179 void processCurLine(const casa::Vector<casa::Bool> &mask)
180 throw(casa::AipsError);
181};
182
183//
184///////////////////////////////////////////////////////////////////////////////
185
186} // namespace asap
187
188///////////////////////////////////////////////////////////////////////////////
189//
190// RunningBox - a running box calculator. This class implements
191// interations over the specified spectrum and calculates
192// running box filter statistics.
193//
194
195// set up the object with the references to actual data
196// and the number of channels in the running box
197RunningBox::RunningBox(const casa::Vector<casa::Float> &in_spectrum,
198 const casa::Vector<casa::Bool> &in_mask,
199 const std::pair<int,int> &in_edge,
200 int in_max_box_nchan) throw(AipsError) :
201 spectrum(in_spectrum), mask(in_mask), edge(in_edge),
202 max_box_nchan(in_max_box_nchan)
203{
204 rewind();
205}
206
207void RunningBox::rewind() throw(AipsError) {
208 // fill statistics for initial box
209 box_chan_cntr=0; // no channels are currently in the box
210 sumf=0.; // initialize statistics
211 sumf2=0.;
212 sumch=0.;
213 sumch2=0.;
214 sumfch=0.;
215 int initial_box_ch=edge.first;
216 for (;initial_box_ch<edge.second && box_chan_cntr<max_box_nchan;
217 ++initial_box_ch)
218 advanceRunningBox(initial_box_ch);
219
220 if (initial_box_ch==edge.second)
221 throw AipsError("RunningBox::rewind - too much channels are masked");
222
223 cur_channel=edge.first;
224 start_advance=initial_box_ch-max_box_nchan/2;
225}
226
227// access to the statistics
228const casa::Float& RunningBox::getLinMean() const throw(AipsError)
229{
230 DebugAssert(cur_channel<edge.second, AipsError);
231 if (need2recalculate) updateDerivativeStatistics();
232 return linmean;
233}
234
235const casa::Float& RunningBox::getLinVariance() const throw(AipsError)
236{
237 DebugAssert(cur_channel<edge.second, AipsError);
238 if (need2recalculate) updateDerivativeStatistics();
239 return linvariance;
240}
241
242const casa::Float RunningBox::aboveMean() const throw(AipsError)
243{
244 DebugAssert(cur_channel<edge.second, AipsError);
245 if (need2recalculate) updateDerivativeStatistics();
246 return spectrum[cur_channel]-linmean;
247}
248
249int RunningBox::getChannel() const throw()
250{
251 return cur_channel;
252}
253
254// actual number of channels in the box (max_box_nchan, if no channels
255// are masked)
256int RunningBox::getNumberOfBoxPoints() const throw()
257{
258 return box_chan_cntr;
259}
260
261// supplementary function to control running mean calculations.
262// It adds a specified channel to the running mean box and
263// removes (ch-max_box_nchan+1)'th channel from there
264// Channels, for which the mask is false or index is beyond the
265// allowed range, are ignored
266void RunningBox::advanceRunningBox(int ch) throw(AipsError)
267{
268 if (ch>=edge.first && ch<edge.second)
269 if (mask[ch]) { // ch is a valid channel
270 ++box_chan_cntr;
271 sumf+=spectrum[ch];
272 sumf2+=square(spectrum[ch]);
273 sumch+=Float(ch);
274 sumch2+=square(Float(ch));
275 sumfch+=spectrum[ch]*Float(ch);
276 need2recalculate=True;
277 }
278 int ch2remove=ch-max_box_nchan;
279 if (ch2remove>=edge.first && ch2remove<edge.second)
280 if (mask[ch2remove]) { // ch2remove is a valid channel
281 --box_chan_cntr;
282 sumf-=spectrum[ch2remove];
283 sumf2-=square(spectrum[ch2remove]);
284 sumch-=Float(ch2remove);
285 sumch2-=square(Float(ch2remove));
286 sumfch-=spectrum[ch2remove]*Float(ch2remove);
287 need2recalculate=True;
288 }
289}
290
291// next channel
292void RunningBox::next() throw(AipsError)
293{
294 AlwaysAssert(cur_channel<edge.second,AipsError);
295 ++cur_channel;
296 if (cur_channel+max_box_nchan/2<edge.second && cur_channel>=start_advance)
297 advanceRunningBox(cur_channel+max_box_nchan/2); // update statistics
298}
299
300// checking whether there are still elements
301casa::Bool RunningBox::haveMore() const throw()
302{
303 return cur_channel<edge.second;
304}
305
306// calculate derivative statistics. This function is const, because
307// it updates the cache only
308void RunningBox::updateDerivativeStatistics() const throw(AipsError)
309{
310 AlwaysAssert(box_chan_cntr, AipsError);
311
312 Float mean=sumf/Float(box_chan_cntr);
313
314 // linear LSF formulae
315 Float meanch=sumch/Float(box_chan_cntr);
316 Float meanch2=sumch2/Float(box_chan_cntr);
317 if (meanch==meanch2 || box_chan_cntr<3) {
318 // vertical line in the spectrum, can't calculate linmean and linvariance
319 linmean=0.;
320 linvariance=0.;
321 } else {
322 Float coeff=(sumfch/Float(box_chan_cntr)-meanch*mean)/
323 (meanch2-square(meanch));
324 linmean=coeff*(Float(cur_channel)-meanch)+mean;
325 linvariance=sqrt(sumf2/Float(box_chan_cntr)-square(mean)-
326 square(coeff)*(meanch2-square(meanch)));
327 }
328 need2recalculate=False;
329}
330
331
332//
333///////////////////////////////////////////////////////////////////////////////
334
335///////////////////////////////////////////////////////////////////////////////
336//
337// LFAboveThreshold - a running mean algorithm for line detection
338//
339//
340
341
342// set up the detection criterion
343LFAboveThreshold::LFAboveThreshold(std::list<pair<int,int> > &in_lines,
344 int in_min_nchan,
345 casa::Float in_threshold) throw() :
346 min_nchan(in_min_nchan), threshold(in_threshold),
347 lines(in_lines), running_box(NULL) {}
348
349LFAboveThreshold::~LFAboveThreshold() throw()
350{
351 if (running_box!=NULL) delete running_box;
352}
353
354// replace the detection criterion
355void LFAboveThreshold::setCriterion(int in_min_nchan, casa::Float in_threshold)
356 throw()
357{
358 min_nchan=in_min_nchan;
359 threshold=in_threshold;
360}
361
362
363// process a channel: update cur_line and is_detected before and
364// add a new line to the list, if necessary
365void LFAboveThreshold::processChannel(Bool detect,
366 const casa::Vector<casa::Bool> &mask) throw(casa::AipsError)
367{
368 try {
369 if (detect) {
370 if (is_detected_before)
371 cur_line.second=running_box->getChannel()+1;
372 else {
373 is_detected_before=True;
374 cur_line.first=running_box->getChannel();
375 cur_line.second=running_box->getChannel()+1;
376 }
377 } else processCurLine(mask);
378 }
379 catch (const AipsError &ae) {
380 throw;
381 }
382 catch (const exception &ex) {
383 throw AipsError(String("LFAboveThreshold::processChannel - STL error: ")+ex.what());
384 }
385}
386
387// process the interval of channels stored in cur_line
388// if it satisfies the criterion, add this interval as a new line
389void LFAboveThreshold::processCurLine(const casa::Vector<casa::Bool> &mask)
390 throw(casa::AipsError)
391{
392 try {
393 if (is_detected_before) {
394 if (cur_line.second-cur_line.first>min_nchan) {
395 // it was a detection. We need to change the list
396 Bool add_new_line=False;
397 if (lines.size()) {
398 for (int i=lines.back().second;i<cur_line.first;++i)
399 if (mask[i]) { // one valid channel in between
400 // means that we deal with a separate line
401 add_new_line=True;
402 break;
403 }
404 } else add_new_line=True;
405 if (add_new_line)
406 lines.push_back(cur_line);
407 else lines.back().second=cur_line.second;
408 }
409 is_detected_before=False;
410 }
411 }
412 catch (const AipsError &ae) {
413 throw;
414 }
415 catch (const exception &ex) {
416 throw AipsError(String("LFAboveThreshold::processCurLine - STL error: ")+ex.what());
417 }
418}
419
420// find spectral lines and add them into list
421void LFAboveThreshold::findLines(const casa::Vector<casa::Float> &spectrum,
422 const casa::Vector<casa::Bool> &mask,
423 const std::pair<int,int> &edge,
424 int max_box_nchan)
425 throw(casa::AipsError)
426{
427 const int minboxnchan=4;
428 try {
429
430 if (running_box!=NULL) delete running_box;
431 running_box=new RunningBox(spectrum,mask,edge,max_box_nchan);
432
433 // actual search algorithm
434 is_detected_before=False;
435 Vector<Int> signs(spectrum.nelements(),0);
436
437 for (;running_box->haveMore();running_box->next()) {
438 const int ch=running_box->getChannel();
439 if (running_box->getNumberOfBoxPoints()>=minboxnchan)
440 processChannel(mask[ch] && (fabs(running_box->aboveMean()) >=
441 threshold*running_box->getLinVariance()), mask);
442 else processCurLine(mask); // just finish what was accumulated before
443 const Float buf=running_box->aboveMean();
444 if (buf>0) signs[ch]=1;
445 else if (buf<0) signs[ch]=-1;
446 else if (buf==0) signs[ch]=0;
447 }
448 if (lines.size())
449 searchForWings(lines,signs,mask,edge);
450 }
451 catch (const AipsError &ae) {
452 throw;
453 }
454 catch (const exception &ex) {
455 throw AipsError(String("LFAboveThreshold::findLines - STL error: ")+ex.what());
456 }
457}
458
459//
460///////////////////////////////////////////////////////////////////////////////
461
462///////////////////////////////////////////////////////////////////////////////
463//
464// LFLineListOperations::IntersectsWith - An auxiliary object function
465// to test whether two lines have a non-void intersection
466//
467
468
469// line1 - range of the first line: start channel and stop+1
470LFLineListOperations::IntersectsWith::IntersectsWith(const std::pair<int,int> &in_line1) :
471 line1(in_line1) {}
472
473
474// return true if line2 intersects with line1 with at least one
475// common channel, and false otherwise
476// line2 - range of the second line: start channel and stop+1
477bool LFLineListOperations::IntersectsWith::operator()(const std::pair<int,int> &line2)
478 const throw()
479{
480 if (line2.second<line1.first) return false; // line2 is at lower channels
481 if (line2.first>line1.second) return false; // line2 is at upper channels
482 return true; // line2 has an intersection or is adjacent to line1
483}
484
485//
486///////////////////////////////////////////////////////////////////////////////
487
488///////////////////////////////////////////////////////////////////////////////
489//
490// LFLineListOperations::BuildUnion - An auxiliary object function to build a union
491// of several lines to account for a possibility of merging the nearby lines
492//
493
494// set an initial line (can be a first line in the sequence)
495LFLineListOperations::BuildUnion::BuildUnion(const std::pair<int,int> &line1) :
496 temp_line(line1) {}
497
498// update temp_line with a union of temp_line and new_line
499// provided there is no gap between the lines
500void LFLineListOperations::BuildUnion::operator()(const std::pair<int,int> &new_line)
501 throw()
502{
503 if (new_line.first<temp_line.first) temp_line.first=new_line.first;
504 if (new_line.second>temp_line.second) temp_line.second=new_line.second;
505}
506
507// return the result (temp_line)
508const std::pair<int,int>& LFLineListOperations::BuildUnion::result() const throw()
509{
510 return temp_line;
511}
512
513//
514///////////////////////////////////////////////////////////////////////////////
515
516///////////////////////////////////////////////////////////////////////////////
517//
518// LFLineListOperations::LaterThan - An auxiliary object function to test whether a
519// specified line is at lower spectral channels (to preserve the order in
520// the line list)
521//
522
523// setup the line to compare with
524LFLineListOperations::LaterThan::LaterThan(const std::pair<int,int> &in_line1) :
525 line1(in_line1) {}
526
527// return true if line2 should be placed later than line1
528// in the ordered list (so, it is at greater channel numbers)
529bool LFLineListOperations::LaterThan::operator()(const std::pair<int,int> &line2)
530 const throw()
531{
532 if (line2.second<line1.first) return false; // line2 is at lower channels
533 if (line2.first>line1.second) return true; // line2 is at upper channels
534
535 // line2 intersects with line1. We should have no such situation in
536 // practice
537 return line2.first>line1.first;
538}
539
540//
541///////////////////////////////////////////////////////////////////////////////
542
543
544///////////////////////////////////////////////////////////////////////////////
545//
546// SDLineFinder - a class for automated spectral line search
547//
548//
549
550SDLineFinder::SDLineFinder() throw() : edge(0,0)
551{
552 // detection threshold - the minimal signal to noise ratio
553 threshold=3.; // 3 sigma is a default
554 box_size=1./5.; // default box size for running mean calculations is
555 // 1/5 of the whole spectrum
556 // A minimum number of consequtive channels, which should satisfy
557 // the detection criterion, to be a detection
558 min_nchan=3; // default is 3 channels
559}
560
561SDLineFinder::~SDLineFinder() throw(AipsError) {}
562
563// set scan to work with (in_scan parameter), associated mask (in_mask
564// parameter) and the edge channel rejection (in_edge parameter)
565// if in_edge has zero length, all channels chosen by mask will be used
566// if in_edge has one element only, it represents the number of
567// channels to drop from both sides of the spectrum
568// in_edge is introduced for convinience, although all functionality
569// can be achieved using a spectrum mask only
570void SDLineFinder::setScan(const SDMemTableWrapper &in_scan,
571 const std::vector<bool> &in_mask,
572 const boost::python::tuple &in_edge) throw(AipsError)
573{
574 try {
575 scan=in_scan.getCP();
576 AlwaysAssert(!scan.null(),AipsError);
577 if (scan->nRow()!=1)
578 throw AipsError("SDLineFinder::setScan - in_scan contains more than 1 row."
579 "Choose one first.");
580 mask=in_mask;
581 if (mask.nelements()!=scan->nChan())
582 throw AipsError("SDLineFinder::setScan - in_scan and in_mask have different"
583 "number of spectral channels.");
584
585 // number of elements in the in_edge tuple
586 int n=extract<int>(in_edge.attr("__len__")());
587 if (n>2 || n<0)
588 throw AipsError("SDLineFinder::setScan - the length of the in_edge parameter"
589 "should not exceed 2");
590 if (!n) {
591 // all spectrum, no rejection
592 edge.first=0;
593 edge.second=scan->nChan();
594 } else {
595 edge.first=extract<int>(in_edge.attr("__getitem__")(0));
596 if (edge.first<0)
597 throw AipsError("SDLineFinder::setScan - the in_edge parameter has a negative"
598 "number of channels to drop");
599 if (edge.first>=scan->nChan())
600 throw AipsError("SDLineFinder::setScan - all channels are rejected by the in_edge parameter");
601 if (n==2) {
602 edge.second=extract<int>(in_edge.attr("__getitem__")(1));
603 if (edge.second<0)
604 throw AipsError("SDLineFinder::setScan - the in_edge parameter has a negative"
605 "number of channels to drop");
606 edge.second=scan->nChan()-edge.second;
607 } else edge.second=scan->nChan()-edge.first;
608 if (edge.second<0 || (edge.second+edge.first)>scan->nChan())
609 throw AipsError("SDLineFinder::setScan - all channels are rejected by the in_edge parameter");
610 }
611 }
612 catch(const AipsError &ae) {
613 // setScan is unsuccessfull, reset scan/mask/edge
614 scan=CountedConstPtr<SDMemTable>(); // null pointer
615 mask.resize(0);
616 edge=pair<int,int>(0,0);
617 throw;
618 }
619}
620
621// search for spectral lines. Number of lines found is returned
622int SDLineFinder::findLines() throw(casa::AipsError)
623{
624 const int minboxnchan=4;
625 if (scan.null())
626 throw AipsError("SDLineFinder::findLines - a scan should be set first,"
627 " use set_scan");
628 DebugAssert(mask.nelements()==scan->nChan(), AipsError);
629 int max_box_nchan=int(scan->nChan()*box_size); // number of channels in running
630 // box
631 if (max_box_nchan<2)
632 throw AipsError("SDLineFinder::findLines - box_size is too small");
633
634 scan->getSpectrum(spectrum);
635
636 lines.resize(0); // search from the scratch
637 Vector<Bool> temp_mask(mask);
638
639 Bool first_pass=True;
640 while (true) {
641 // a buffer for new lines found at this iteration
642 std::list<pair<int,int> > new_lines;
643
644 // line find algorithm
645 LFAboveThreshold lfalg(new_lines,min_nchan, threshold);
646
647 try {
648 lfalg.findLines(spectrum,temp_mask,edge,max_box_nchan);
649 }
650 catch(const AipsError &ae) {
651 if (first_pass) throw;
652 break; // nothing new
653 }
654 first_pass=False;
655 if (!new_lines.size()) break; // nothing new
656
657 // update the list (lines) merging intervals, if necessary
658 addNewSearchResult(new_lines,lines);
659 // get a new mask
660 temp_mask=getMask();
661 }
662 return int(lines.size());
663}
664
665
666// get the mask to mask out all lines that have been found (default)
667// if invert=true, only channels belong to lines will be unmasked
668// Note: all channels originally masked by the input mask (in_mask
669// in setScan) or dropped out by the edge parameter (in_edge
670// in setScan) are still excluded regardless on the invert option
671std::vector<bool> SDLineFinder::getMask(bool invert)
672 const throw(casa::AipsError)
673{
674 try {
675 if (scan.null())
676 throw AipsError("SDLineFinder::getMask - a scan should be set first,"
677 " use set_scan followed by find_lines");
678 DebugAssert(mask.nelements()==scan->nChan(), AipsError);
679 /*
680 if (!lines.size())
681 throw AipsError("SDLineFinder::getMask - one have to search for "
682 "lines first, use find_lines");
683 */
684 std::vector<bool> res_mask(mask.nelements());
685 // iterator through lines
686 std::list<std::pair<int,int> >::const_iterator cli=lines.begin();
687 for (int ch=0;ch<res_mask.size();++ch)
688 if (ch<edge.first || ch>=edge.second) res_mask[ch]=false;
689 else if (!mask[ch]) res_mask[ch]=false;
690 else {
691 res_mask[ch]=!invert; // no line by default
692 if (cli==lines.end()) continue;
693 if (ch>=cli->first && ch<cli->second)
694 res_mask[ch]=invert; // this is a line
695 if (ch>=cli->second)
696 ++cli; // next line in the list
697 }
698
699 return res_mask;
700 }
701 catch (const AipsError &ae) {
702 throw;
703 }
704 catch (const exception &ex) {
705 throw AipsError(String("SDLineFinder::getMask - STL error: ")+ex.what());
706 }
707}
708
709// get range for all lines found. If defunits is true (default), the
710// same units as used in the scan will be returned (e.g. velocity
711// instead of channels). If defunits is false, channels will be returned
712std::vector<int> SDLineFinder::getLineRanges(bool defunits)
713 const throw(casa::AipsError)
714{
715 try {
716 if (scan.null())
717 throw AipsError("SDLineFinder::getLineRanges - a scan should be set first,"
718 " use set_scan followed by find_lines");
719 DebugAssert(mask.nelements()==scan->nChan(), AipsError);
720
721 if (!lines.size())
722 throw AipsError("SDLineFinder::getLineRanges - one have to search for "
723 "lines first, use find_lines");
724
725 // temporary
726 if (defunits)
727 throw AipsError("SDLineFinder::getLineRanges - sorry, defunits=true have not "
728 "yet been implemented");
729 //
730 std::vector<int> res(2*lines.size());
731 // iterator through lines & result
732 std::list<std::pair<int,int> >::const_iterator cli=lines.begin();
733 std::vector<int>::iterator ri=res.begin();
734 for (;cli!=lines.end() && ri!=res.end();++cli,++ri) {
735 *ri=cli->first;
736 if (++ri!=res.end())
737 *ri=cli->second-1;
738 }
739 return res;
740 }
741 catch (const AipsError &ae) {
742 throw;
743 }
744 catch (const exception &ex) {
745 throw AipsError(String("SDLineFinder::getLineRanges - STL error: ")+ex.what());
746 }
747}
748
749//
750///////////////////////////////////////////////////////////////////////////////
751
752
753///////////////////////////////////////////////////////////////////////////////
754//
755// LFLineListOperations - a class incapsulating operations with line lists
756// The LF prefix stands for Line Finder
757//
758
759// concatenate two lists preserving the order. If two lines appear to
760// be adjacent, they are joined into the new one
761void LFLineListOperations::addNewSearchResult(const std::list<pair<int, int> > &newlines,
762 std::list<std::pair<int, int> > &lines_list)
763 throw(AipsError)
764{
765 try {
766 for (std::list<pair<int,int> >::const_iterator cli=newlines.begin();
767 cli!=newlines.end();++cli) {
768
769 // the first item, which has a non-void intersection or touches
770 // the new line
771 std::list<pair<int,int> >::iterator pos_beg=find_if(lines_list.begin(),
772 lines_list.end(), IntersectsWith(*cli));
773 // the last such item
774 std::list<pair<int,int> >::iterator pos_end=find_if(pos_beg,
775 lines_list.end(), not1(IntersectsWith(*cli)));
776
777 // extract all lines which intersect or touch a new one into
778 // a temporary buffer. This may invalidate the iterators
779 // line_buffer may be empty, if no lines intersects with a new
780 // one.
781 std::list<pair<int,int> > lines_buffer;
782 lines_buffer.splice(lines_buffer.end(),lines_list, pos_beg, pos_end);
783
784 // build a union of all intersecting lines
785 pair<int,int> union_line=for_each(lines_buffer.begin(),
786 lines_buffer.end(),BuildUnion(*cli)).result();
787
788 // search for a right place for the new line (union_line) and add
789 std::list<pair<int,int> >::iterator pos2insert=find_if(lines_list.begin(),
790 lines_list.end(), LaterThan(union_line));
791 lines_list.insert(pos2insert,union_line);
792 }
793 }
794 catch (const AipsError &ae) {
795 throw;
796 }
797 catch (const exception &ex) {
798 throw AipsError(String("LFLineListOperations::addNewSearchResult - STL error: ")+ex.what());
799 }
800}
801
802// extend all line ranges to the point where a value stored in the
803// specified vector changes (e.g. value-mean change its sign)
804// This operation is necessary to include line wings, which are below
805// the detection threshold. If lines becomes adjacent, they are
806// merged together. Any masked channel stops the extension
807void LFLineListOperations::searchForWings(std::list<std::pair<int, int> > &newlines,
808 const casa::Vector<casa::Int> &signs,
809 const casa::Vector<casa::Bool> &mask,
810 const std::pair<int,int> &edge) throw(casa::AipsError)
811{
812 try {
813 for (std::list<pair<int,int> >::iterator li=newlines.begin();
814 li!=newlines.end();++li) {
815 // update the left hand side
816 for (int n=li->first-1;n>=edge.first;--n) {
817 if (!mask[n]) break;
818 if (signs[n]==signs[li->first] && signs[li->first])
819 li->first=n;
820 else break;
821 }
822 // update the right hand side
823 for (int n=li->second;n<edge.second;++n) {
824 if (!mask[n]) break;
825 if (signs[n]==signs[li->second-1] && signs[li->second-1])
826 li->second=n;
827 else break;
828 }
829 }
830 // need to search for possible mergers.
831 std::list<std::pair<int, int> > result_buffer;
832 addNewSearchResult(newlines,result_buffer);
833 newlines.clear();
834 newlines.splice(newlines.end(),result_buffer);
835 }
836 catch (const AipsError &ae) {
837 throw;
838 }
839 catch (const exception &ex) {
840 throw AipsError(String("LFLineListOperations::extendLines - STL error: ")+ex.what());
841 }
842}
843
844//
845///////////////////////////////////////////////////////////////////////////////
Note: See TracBrowser for help on using the repository browser.