Line data Source code
1 : //# --------------------------------------------------------------------
2 : //# LineFinder.cc: this defines utility functions of line finding
3 : //# --------------------------------------------------------------------
4 : //# Copyright (C) 2015
5 : //# National Astronomical Observatory of Japan
6 : //#
7 : //# This library is free software; you can redistribute it and/or modify it
8 : //# under the terms of the GNU Library General Public License as published by
9 : //# the Free Software Foundation; either version 2 of the License, or (at your
10 : //# option) any later version.
11 : //#
12 : //# This library is distributed in the hope that it will be useful, but WITHOUT
13 : //# ANY WARRANTY; without even the implied warranty of MERCHANTABILITY or
14 : //# FITNESS FOR A PARTICULAR PURPOSE. See the GNU Library General Public
15 : //# License for more details.
16 : //#
17 : //# You should have received a copy of the GNU Library General Public License
18 : //# along with this library; if not, write to the Free Software Foundation,
19 : //# Inc., 675 Massachusetts Ave, Cambridge, MA 02139, USA.
20 : //#
21 : //# Correspondence concerning AIPS++ should be addressed as follows:
22 : //# Internet email: aips2-request@nrao.edu.
23 : //# Postal address: AIPS++ Project Office
24 : //# National Radio Astronomy Observatory
25 : //# 520 Edgemont Road
26 : //# Charlottesville, VA 22903-2475 USA
27 : //#
28 : //# $Id$
29 : #include <algorithm>
30 : #include <cfloat>
31 : #include <iostream>
32 : #include <string>
33 :
34 : #include <libsakura/sakura.h>
35 :
36 : #include <casacore/casa/Logging/LogIO.h>
37 : #include <singledish/SingleDish/LineFinder.h>
38 : #include <singledish/SingleDish/LineFindingUtils.h>
39 : #include <casacore/casa/Utilities/Assert.h>
40 : #include <casa_sakura/SakuraAlignedArray.h>
41 :
42 : using namespace std;
43 :
44 : #define _ORIGIN LogOrigin("LineFinder", __func__, WHERE)
45 : //#define KS_DEBUG
46 :
47 : using namespace casacore;
48 : using namespace casacore;
49 : using namespace casacore;
50 : using namespace casacore;
51 : using namespace casacore;
52 : using namespace casacore;
53 : using namespace casacore;
54 : namespace casa { //# NAMESPACE CASA - BEGIN
55 :
56 : namespace linefinder {//# NAMESPACE LINEFINDER - BEGIN
57 :
58 3484 : void getMask(size_t const num_mask, bool *mask,
59 : list<pair<size_t,size_t>>& ranges,
60 : bool const invert, bool const initialize)
61 : {
62 3484 : bool const lineval = (!invert);
63 3484 : if (initialize) {
64 0 : for (size_t i = 0; i<num_mask; ++i) {
65 0 : mask[i] = invert;
66 : }
67 : }
68 12718 : for (list<pair<size_t,size_t>>::iterator iter=ranges.begin();
69 21952 : iter!=ranges.end(); ++iter) {
70 9234 : AlwaysAssert((*iter).second >= (*iter).first, AipsError);
71 5176266 : for (size_t i=(*iter).first; i<(*iter).second+1 && i<num_mask; ++i) {
72 5167032 : mask[i] = lineval;
73 : }
74 : }
75 3484 : }
76 :
77 3504 : list<pair<size_t,size_t>> MADLineFinder(size_t const num_data,
78 : float const* data,
79 : bool * mask,
80 : float const threshold,
81 : uint8_t max_iteration,
82 : size_t const minwidth,
83 : size_t const maxwidth,
84 : size_t const avg_limit,
85 : pair<size_t,size_t> edge)
86 : {
87 : #if defined(KS_DEBUG)
88 : LogIO os(_ORIGIN);
89 : os << LogIO::DEBUG1 << "Start MADLineFinder" << LogIO::POST;
90 : #endif
91 : // value check and edge handling
92 3504 : size_t min_valid_elem = 2;
93 3504 : AlwaysAssert(edge.first+edge.second<num_data-min_valid_elem, AipsError);
94 3504 : AlwaysAssert(maxwidth>minwidth, AipsError);
95 5044 : for (size_t i=0; i<edge.first; ++i){
96 1540 : mask[i] = false;
97 : }
98 5058 : for (size_t i=num_data-edge.second; i<num_data; ++i){
99 1554 : mask[i] = false;
100 : }
101 3504 : size_t num_valid_elem=LineFinderUtils::countTrue(num_data,mask);
102 3504 : AlwaysAssert(num_valid_elem >= min_valid_elem+edge.first+edge.second, AipsError);
103 : // factor for each iteration of binning
104 3504 : size_t const binfactor = 4;
105 3504 : size_t average_factor =1;
106 : //size_t maxgap = num_data;
107 3504 : list<pair<size_t,size_t>> line_list;
108 7008 : SakuraAlignedArray<float> binned_data(num_data);
109 7008 : SakuraAlignedArray<bool> binned_mask(num_data);
110 3504 : float *binned_data_p = binned_data.data();
111 3504 : bool *binned_mask_p = binned_mask.data();
112 : while (true) {
113 : // Bin spectrum and mask
114 7020 : size_t const maxwidth_bin = std::max(maxwidth/average_factor,static_cast<size_t>(1));
115 : //size_t const minwidth_bin = minwidth/average_factor;
116 7020 : size_t const offset = average_factor/2;
117 : size_t num_binned = \
118 7020 : LineFinderUtils::binDataAndMask<float>(num_data, data, mask,
119 : average_factor, num_data,
120 : binned_data_p,
121 : binned_mask_p,
122 : offset,false);
123 : #if defined(KS_DEBUG)
124 : os << LogIO::DEBUG1 << "Bin level: " << average_factor << " (offset: " << offset << ")" << LogIO::POST;
125 : os << LogIO::DEBUG1 << "---> Binned channel = " << num_binned << LogIO::POST;
126 : #endif
127 : // caluculate MAD array
128 7020 : SakuraAlignedArray<float> mad_data(num_binned);
129 7020 : SakuraAlignedArray<bool> line_mask(num_binned);
130 7020 : SakuraAlignedArray<bool> search_mask(num_binned);
131 7020 : float *mad_data_p = mad_data.data();
132 7020 : bool *line_mask_p = line_mask.data();
133 7020 : bool *search_mask_p = search_mask.data();
134 15933568 : for (size_t i=0; i<num_binned; ++i) {
135 15926548 : search_mask_p[i] = binned_mask_p[i];
136 : }
137 7020 : float prev_mad_threshold = FLT_MAX;
138 7020 : size_t prev_num_line_chan = 0;
139 7020 : list<pair<size_t,size_t>> prev_line_list;
140 7020 : prev_line_list.clear();
141 28612 : for (size_t iteration=0; iteration < max_iteration; ++iteration) {
142 : #if defined(KS_DEBUG)
143 : os << LogIO::DEBUG1 << "START iteration " << iteration << LogIO::POST;
144 : #endif
145 28612 : list<pair<size_t,size_t>> new_lines; // lines found in this iteration
146 : // working data array. sorted after median calculation.
147 28612 : LineFinderUtils::calculateMAD(num_binned, binned_data_p,
148 : search_mask_p, mad_data_p);
149 : // use lower 80% of mad_array and use median of the array as a criteria.
150 28612 : float mad80 = LineFinderUtils::maskedMedian(num_binned, mad_data_p,
151 : search_mask_p, 0.8);
152 28612 : float mad_threshold = threshold*sqrt(average_factor)* mad80;
153 : // cout << "mad_threshold = " << mad_threshold << endl;
154 : #if defined(KS_DEBUG)
155 : os << LogIO::DEBUG1 << "mad = " << mad80 << ", mad_threshold = " << mad_threshold << " (factor: " << threshold << ", bin=" << average_factor << ")" << LogIO::POST;
156 : #endif
157 28612 : if (mad_threshold >= prev_mad_threshold){
158 : #if defined(KS_DEBUG)
159 : os << LogIO::DEBUG1 << "threshold increased. stop iteration." << LogIO::POST;
160 : #endif
161 3121 : break; // stop iteration
162 : }
163 25491 : LineFinderUtils::createMaskByAThreshold(num_binned, mad_data_p,
164 : binned_mask_p, mad_threshold,
165 : line_mask_p);
166 : // channel mask -> mask range
167 25491 : LineFinderUtils::maskToRangesList(num_binned, line_mask_p, new_lines);
168 : #if defined(KS_DEBUG)
169 : os << LogIO::DEBUG1 << "raw detection:" << LineFinderUtils::FormatLineString(new_lines) << LogIO::POST;
170 : #endif
171 25491 : size_t const binned_min = max(minwidth/average_factor,1UL);
172 25491 : LineFinderUtils::rejectNarrowRange(binned_min, new_lines);
173 : #if defined(KS_DEBUG)
174 : os << LogIO::DEBUG1 << "rejection by min width (" << binned_min << " chan):" << LineFinderUtils::FormatLineString(new_lines) << LogIO::POST;
175 : #endif
176 : //rejectByRangeWidth(minwidth_bin,maxwidth_bin,new_lines);
177 25491 : if ( new_lines.size()>0 ) {
178 : // merge flagged gaps
179 : //LineFinderUtils::mergeGapByFalse(num_binned, binned_mask.data,
180 : // maxgap/average_factor, new_lines);
181 : // extend wing
182 24689 : SakuraAlignedArray<int8_t> sign(num_binned);
183 24689 : int8_t *sign_p = sign.data();
184 24689 : LineFinderUtils::createSignByAThreshold(num_binned, mad_data_p,
185 : //mad_threshold, sign_p);
186 : mad80, sign_p);
187 24689 : LineFinderUtils::extendRangeBySign(num_binned, sign_p,
188 : binned_mask_p,
189 : new_lines);
190 : #if defined(KS_DEBUG)
191 : os << LogIO::DEBUG1 << "extend for wing:" << LineFinderUtils::FormatLineString(new_lines) << LogIO::POST;
192 : #endif
193 : // merge overlap
194 24689 : LineFinderUtils::mergeOverlappingRanges(new_lines);
195 : // reject by max line width
196 24689 : LineFinderUtils::rejectWideRange(maxwidth_bin,new_lines);
197 : // merge lines with small gap
198 24689 : LineFinderUtils::mergeSmallGapByFraction(0.25,maxwidth_bin,new_lines);
199 : #if defined(KS_DEBUG)
200 : os << LogIO::DEBUG1 << "merge overlapping and nearby lines:" << LineFinderUtils::FormatLineString(new_lines) << LogIO::POST;
201 : #endif
202 : // de-bin line ranges
203 24689 : LineFinderUtils::deBinRanges(average_factor, offset, new_lines);
204 : #if defined(KS_DEBUG)
205 : os << LogIO::DEBUG1 << "debin lines: " << LineFinderUtils::FormatLineString(new_lines) << LogIO::POST;
206 : #endif
207 : // reject by min and max line width
208 24689 : LineFinderUtils::rejectWideRange(maxwidth,new_lines);
209 24689 : LineFinderUtils::rejectNarrowRange(minwidth, new_lines);
210 : #if defined(KS_DEBUG)
211 : os << LogIO::DEBUG1 << "check for line width: " << LineFinderUtils::FormatLineString(new_lines) << LogIO::POST;
212 : #endif
213 : // update search mask for the next iteration
214 98569 : for (list<pair<size_t,size_t>>::iterator iter=new_lines.begin();
215 172449 : iter!=new_lines.end(); ++iter) {
216 32668556 : for (size_t i=(*iter).first; i<=(*iter).second && i<num_binned; ++i) {
217 32594676 : search_mask_p[i] = false;
218 : }
219 : }
220 : // culculate the number of line channels
221 24689 : size_t curr_num_line_chan = 0;
222 98569 : for (list<pair<size_t,size_t>>::iterator iter=new_lines.begin();
223 98569 : iter!=new_lines.end(); ++iter) {
224 73880 : curr_num_line_chan += ((*iter).second - (*iter).first + 1);
225 : }
226 24689 : if (curr_num_line_chan <= prev_num_line_chan) {
227 : #if defined(KS_DEBUG)
228 : os << LogIO::DEBUG1 << "No new line found. stop at iteration " << iteration << LogIO::POST;
229 : #endif
230 3097 : break; // stop iteration
231 : }
232 : // save values of current iteration for next loop
233 21592 : prev_num_line_chan = curr_num_line_chan;
234 21592 : prev_mad_threshold = mad_threshold;
235 21592 : prev_line_list.clear();
236 21592 : prev_line_list.splice(prev_line_list.end(), new_lines);
237 : }
238 : else { // no lines found
239 : #if defined(KS_DEBUG)
240 : os << LogIO::DEBUG1 << "No line found" << LogIO::POST;
241 : #endif
242 802 : break;
243 : }
244 : } // end of iteration loop
245 : // merge line list from different bin level
246 : // cout << "Final line list for average_factor = " << average_factor << endl;
247 : // for (list<pair<size_t,size_t>>::iterator iter=prev_line_list.begin();
248 : // iter!=prev_line_list.end(); ++iter) {
249 : // cout << "[" << (*iter).first << ", " << (*iter).second << "], ";
250 : // }
251 : // cout << endl;
252 7020 : LineFinderUtils::mergeOverlapInTwoLists(line_list, prev_line_list);
253 7020 : average_factor *= binfactor;
254 7020 : if (average_factor <= avg_limit)
255 3516 : continue; // go to the next bin level.
256 : else
257 3504 : break;
258 3516 : } // end of bin level loop
259 7008 : return line_list;
260 : }
261 :
262 : } //# NAMESPACE LINEFINDER - END
263 :
264 : } //# NAMESPACE CASA - END
|
