Line data Source code
1 : //# --------------------------------------------------------------------
2 : //# LineFindingUtils.tcc: 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: casa-feedback@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 <casacore/casa/Utilities/Assert.h>
31 : #include <cmath>
32 :
33 : #include <libsakura/sakura.h>
34 : #include <singledish/SingleDish/LineFindingUtils.h>
35 :
36 : using namespace std;
37 :
38 : using namespace casacore;
39 : using namespace casacore;
40 : using namespace casacore;
41 : using namespace casacore;
42 : using namespace casacore;
43 : using namespace casacore;
44 : using namespace casacore;
45 : namespace casa { //# NAMESPACE CASA - BEGIN
46 :
47 : template <typename DataType>
48 7022 : size_t LineFinderUtils::binDataAndMask(size_t const num_in,
49 : DataType const* in_data,
50 : bool const* in_mask,
51 : size_t const bin_size,
52 : size_t const num_out,
53 : DataType* out_data,
54 : bool* out_mask,
55 : size_t const offset,
56 : bool const keepsize)
57 : {
58 7022 : if (offset+bin_size > num_in-1) return 0;
59 : // Can NOT offset (>0) when changing size (keepsize=false)
60 7022 : AlwaysAssert((!keepsize || offset==0), AipsError);
61 7022 : size_t num_bin = keepsize ? num_in : (num_in-offset)/bin_size;
62 7022 : AlwaysAssert(num_bin <= num_out, AipsError);
63 7022 : size_t num_loc_out = keepsize ? bin_size : 1;
64 7022 : size_t out_idx = 0;
65 15933583 : for (size_t i=offset; (keepsize ? i : i+bin_size-1) < num_in && out_idx < num_bin;) {
66 15926561 : DataType data = static_cast<DataType>(0);
67 15926561 : bool mask = true;
68 15926561 : size_t iend=i+bin_size;
69 15926561 : size_t count=0;
70 41488967 : for (;i<iend && i < num_in; ++i) {
71 : //Sum up regardless of in_mask because
72 : //it will be flagged anyway if any of mask=false.
73 25562406 : data += in_data[i];
74 25562406 : mask = mask && in_mask[i];
75 25562406 : ++count;
76 : }
77 15926561 : data = mask ? data/static_cast<DataType>(count) : static_cast<DataType>(0);
78 31853135 : for (size_t j = 0; j < num_loc_out && out_idx < num_bin; ++j){
79 15926574 : out_mask[out_idx] = mask;
80 15926574 : out_data[out_idx] = data;
81 15926574 : ++out_idx;
82 : }
83 : }
84 7022 : AlwaysAssert(num_bin==out_idx, AipsError);
85 7022 : return num_bin;
86 : }
87 :
88 28613 : void LineFinderUtils::calculateMAD(size_t const num_data,
89 : float const in_data[],
90 : bool const in_mask[],
91 : float mad[])
92 : {
93 28613 : float median_value = LineFinderUtils::maskedMedian(num_data, in_data, in_mask, 1.0);
94 : //cout << "median value for MAD = " << median_value << endl;
95 94550849 : for (size_t i = 0; i < num_data; ++i) {
96 94522236 : mad[i] = fabs(in_data[i]-median_value);
97 : }
98 28613 : }
99 :
100 25492 : void LineFinderUtils::createMaskByAThreshold(size_t const num_data,
101 : float const in_data[],
102 : bool const in_mask[],
103 : float const threshold,
104 : bool out_mask[])
105 : {
106 25492 : LIBSAKURA_SYMBOL(Status) status = LIBSAKURA_SYMBOL(SetTrueIfGreaterThanOrEqualsFloat)(num_data, in_data, threshold, out_mask);
107 25492 : AlwaysAssert(status == LIBSAKURA_SYMBOL(Status_kOK), AipsError);
108 91371195 : for (size_t i = 0; i < num_data; ++i) {
109 91345703 : out_mask[i] = out_mask[i] && in_mask[i];
110 : }
111 25492 : }
112 :
113 24690 : void LineFinderUtils::deBinRanges(size_t const bin_size, size_t const offset,
114 : list<pair<size_t,size_t>>& range_list)
115 : {
116 24690 : for(list<pair<size_t,size_t>>::iterator iter=range_list.begin();
117 98573 : iter!=range_list.end(); ++iter) {
118 73883 : (*iter).first = (*iter).first*bin_size + offset;
119 73883 : (*iter).second = (*iter).second*bin_size+(bin_size-1) + offset;
120 : }
121 24690 : }
122 :
123 24690 : void LineFinderUtils::extendRangeBySign(size_t num_sign, int8_t const* sign, bool const* mask,
124 : list<pair<size_t,size_t>>& range_list){
125 24690 : for (list<pair<size_t,size_t>>::iterator iter = range_list.begin();
126 184691 : iter!=range_list.end(); ++iter) {
127 160001 : AlwaysAssert((*iter).first <= (*iter).second, AipsError);
128 : // extend left wing
129 160001 : int8_t ref_sign = sign[(*iter).first];
130 48863750 : for (size_t i=(*iter).first-1; ; --i) {
131 48863750 : if (mask[i] && sign[i] == ref_sign) (*iter).first=i;
132 : else break;
133 48703749 : if (i==0) break;
134 : }
135 : // extend right wing
136 160001 : ref_sign = sign[(*iter).second];
137 40895214 : for (size_t i=(*iter).second+1; i<num_sign; ++i) {
138 40895214 : if (mask[i] && sign[i] == ref_sign) (*iter).second=i;
139 : else break;
140 : }
141 : }
142 24690 : }
143 :
144 :
145 25493 : void LineFinderUtils::maskToRangesList(size_t const num_mask, bool const* mask,
146 : list<pair<size_t,size_t>>& out_range)
147 : {
148 : static_assert(static_cast<int8_t>(true)==1, "cast of bool failed");
149 : static_assert(static_cast<int8_t>(false)==0, "cast of bool failed");
150 25493 : out_range.clear();
151 25493 : if (num_mask==0) return;
152 25493 : size_t istart=num_mask;
153 25493 : if (mask[0]) istart=0;
154 91345723 : for (size_t i=1; i<num_mask; ++i) {
155 91320230 : int8_t diff = static_cast<int8_t>(mask[i]) - static_cast<int8_t>(mask[i-1]);
156 91320230 : if (diff==1) { // start of new line
157 192014 : AlwaysAssert(istart==num_mask, AipsError);
158 192014 : istart=i;
159 : }
160 91128216 : else if (diff==-1) { // end of line was i-1 chan
161 192019 : AlwaysAssert(istart<num_mask, AipsError);
162 192019 : out_range.push_back(pair<size_t,size_t>(istart, i-1));
163 192019 : istart=num_mask;
164 : }
165 : }
166 25493 : if (mask[num_mask-1]) {
167 2 : AlwaysAssert(istart<num_mask, AipsError);
168 2 : out_range.push_back(pair<size_t,size_t>(istart, num_mask-1));
169 : }
170 : }
171 :
172 2 : void LineFinderUtils::mergeGapByFalse(size_t const num_mask, bool const* mask,
173 : size_t const maxgap,
174 : list<pair<size_t,size_t>>& range_list)
175 : {
176 2 : if (range_list.size() < 2) return; // nothing to do
177 2 : list<pair<size_t,size_t>>::iterator iter=range_list.begin();
178 2 : size_t from=(*iter).second+1;
179 2 : list<pair<size_t,size_t>> temp_list;
180 2 : temp_list.push_back(*iter);
181 2 : ++iter;
182 4 : for( ; iter!=range_list.end(); ++iter) {
183 2 : size_t to=(*iter).first;
184 2 : if (to-from > maxgap) { // do not merge if gap is too large
185 1 : temp_list.push_back(*iter);
186 1 : continue;
187 : }
188 1 : bool merge=true;
189 6 : for (size_t i=from; i<to && i < num_mask; ++i) {
190 5 : if (mask[i]) { // do not merge if any of channel in betwee the line
191 0 : merge = false;
192 0 : break;
193 : }
194 : }
195 1 : from=(*iter).second;
196 1 : if (merge) {
197 1 : temp_list.back().second = (*iter).second;
198 : }
199 : else {
200 0 : temp_list.push_back(*iter);
201 : }
202 : }
203 2 : range_list.clear();
204 2 : range_list.splice(range_list.end(), temp_list);
205 2 : }
206 :
207 57225 : float LineFinderUtils::maskedMedian(size_t num_data, float const data[],
208 : bool const mask[], float fraction)
209 : {
210 :
211 57225 : Vector<float> local_data(IPosition(1, num_data), const_cast<float *>(data), COPY);
212 : // for (size_t i = 0 ; i < num_data; ++i){
213 : // local_data.data[i] = data[i];
214 : // }
215 57225 : float *local_data_p = local_data.data();
216 57225 : size_t num_valid(num_data+1);
217 57225 : LIBSAKURA_SYMBOL(Status) status = LIBSAKURA_SYMBOL(SortValidValuesDenselyFloat)(num_data, mask, local_data_p, &num_valid);
218 57225 : AlwaysAssert(status == LIBSAKURA_SYMBOL(Status_kOK), AipsError);
219 57225 : AlwaysAssert(num_valid <= num_data, AipsError);
220 57225 : if (fraction<1.0)
221 28612 : num_valid = static_cast<size_t>(num_valid*fraction);
222 57225 : float median_value = LineFinderUtils::getMedianOfSorted<float>(num_valid,
223 : local_data_p);
224 57225 : return median_value;
225 :
226 57225 : }
227 :
228 24690 : void LineFinderUtils::mergeOverlappingRanges(list<pair<size_t,size_t>>& range_list)
229 : {
230 24690 : if (range_list.size() < 2) return; // nothing to do
231 24616 : list<pair<size_t,size_t>> temp_list;
232 24616 : list<pair<size_t,size_t>>::iterator iter=range_list.begin();
233 24616 : temp_list.push_back(*iter);
234 24616 : ++iter;
235 159930 : for ( ; iter!=range_list.end(); ++iter) {
236 : // always start searching from the begining of temp_list
237 : //size_t dummy_cursor =
238 135314 : mergeARangeToList(temp_list, (*iter));
239 : }
240 24616 : range_list.clear();
241 24616 : range_list.splice(range_list.end(), temp_list);
242 24616 : }
243 :
244 :
245 7022 : void LineFinderUtils::mergeOverlapInTwoLists(list<pair<size_t,size_t>>& to,
246 : list<pair<size_t,size_t>>& from)
247 : {
248 7022 : if (from.size()==0) return; // nothing to do
249 6220 : if (to.size()==0) { // replace from with to
250 3104 : to.splice(to.end(), from);
251 3104 : return;
252 : }
253 3116 : size_t cursor = 0;
254 3116 : for (list<pair<size_t,size_t>>::iterator from_iter=from.begin();
255 9322 : from_iter!=from.end(); ++from_iter) {
256 6206 : cursor = mergeARangeToList(to, *from_iter, cursor);
257 : }
258 : }
259 :
260 141520 : size_t LineFinderUtils::mergeARangeToList(list<pair<size_t,size_t>>& range_list,
261 : pair<size_t,size_t>& new_range,
262 : size_t const cursor)
263 : {
264 141520 : AlwaysAssert(new_range.first<=new_range.second, AipsError);
265 : // cursor should be zero if range_list is empty
266 141520 : AlwaysAssert(range_list.size()>0 || cursor==0, AipsError);
267 141520 : if (range_list.size()==0) {//just add new range to range_list;
268 0 : range_list.push_back(new_range);
269 0 : return cursor;
270 : }
271 141520 : AlwaysAssert(cursor < range_list.size(), AipsError);
272 : // add range at the end
273 141520 : if (new_range.first > range_list.back().second) {
274 58410 : range_list.push_back(new_range);
275 58410 : return range_list.size()-1;
276 : }
277 : // move iterator to the start point of the search
278 83110 : list<pair<size_t,size_t>>::iterator start_cursor = range_list.begin();
279 83110 : if (cursor > 0) {
280 20 : for (size_t i=1; i<cursor; ++i) {
281 13 : ++start_cursor;
282 : }
283 7 : AlwaysAssert(new_range.first > (*start_cursor).second, AipsError);
284 7 : ++start_cursor;
285 : }
286 : // add range to the begining
287 83110 : if (new_range.second < (*start_cursor).first) {
288 1 : range_list.insert(start_cursor, new_range);
289 1 : return cursor+1;
290 : }
291 83109 : size_t out_cursor = cursor;
292 : // new_range is in some where after cursor position of range_list.
293 : // find segment in range_list (final) list where overlap starts
294 83109 : list<pair<size_t,size_t>>::iterator insert_begin=range_list.end();
295 83109 : bool any_overlap = false;
296 83109 : for (list<pair<size_t,size_t>>::iterator iter=start_cursor;
297 163002 : iter!=range_list.end(); ++iter, ++out_cursor) {
298 163002 : if ((*iter).first > new_range.second) {
299 : // the range exceeded new_range. No overlap. Insert before this range.
300 1 : range_list.insert(iter, new_range);
301 1 : return out_cursor;
302 163001 : } else if ((*iter).second < new_range.first) { //|| (*iter).first > new_range.second) {
303 : // still too small ranges
304 79893 : continue;
305 : } else {//overlap started
306 83108 : insert_begin = iter;
307 83108 : any_overlap = true;
308 83108 : break;
309 : }
310 : }
311 : // An overlap should be found.
312 83108 : AlwaysAssert(any_overlap, AipsError);
313 : // find segment in range_list (final) list where overlap ends
314 83108 : list<pair<size_t,size_t>>::iterator insert_end=range_list.end();
315 83108 : for (list<pair<size_t,size_t>>::iterator iter=insert_begin;
316 166219 : iter!=range_list.end(); ++iter) {
317 89270 : if ((*iter).first <= new_range.second) {
318 83111 : insert_end = iter;
319 : }
320 6159 : else break; // the range exceed new_range.
321 : }
322 : // there was an overlap
323 83108 : list<pair<size_t,size_t>> temp_list;
324 : // preceeding elements (including the first overlapped range)
325 83108 : temp_list.splice(temp_list.end(), range_list, range_list.begin(), insert_begin);
326 : // overlaping segment
327 83108 : pair<size_t,size_t> unified_range(std::min((*insert_begin).first, new_range.first),
328 166216 : std::max((*insert_end).second, new_range.second));
329 83108 : temp_list.push_back(unified_range);
330 : // trailing elements
331 83108 : if (insert_end != range_list.end()){
332 83108 : temp_list.splice(temp_list.end(),range_list,++insert_end, range_list.end());
333 : }
334 83108 : range_list.clear();
335 83108 : range_list.splice(range_list.end(), temp_list);
336 83108 : return out_cursor;
337 83108 : }
338 :
339 24689 : void LineFinderUtils::mergeSmallGapByFraction(double const fraction,
340 : size_t const maxwidth,
341 : list<pair<size_t,size_t>>& range_list)
342 : {
343 24689 : if (range_list.size() < 2) return; // nothing to do
344 24611 : list<pair<size_t,size_t>> temp_list;
345 24611 : list<pair<size_t,size_t>>::iterator iter=range_list.begin();
346 24611 : temp_list.push_back(*iter);
347 24611 : ++iter;
348 83018 : for (; iter != range_list.end(); ++iter) {
349 58407 : size_t curr_width = (*iter).second - (*iter).first;
350 58407 : size_t prev_width = temp_list.back().second - temp_list.back().first;
351 58407 : size_t gap = (*iter).first-temp_list.back().second;
352 67623 : if (gap < static_cast<size_t>(std::min(prev_width, curr_width)*fraction) && \
353 9216 : (*iter).second-temp_list.back().first < maxwidth) {
354 9216 : temp_list.back().second = (*iter).second;
355 : }
356 : else {
357 49191 : temp_list.push_back(*iter);
358 : }
359 : }
360 24611 : range_list.clear();
361 24611 : range_list.splice(range_list.end(), temp_list);
362 24611 : }
363 :
364 49379 : void LineFinderUtils::rejectWideRange(size_t const maxwidth,
365 : list<pair<size_t,size_t>>& range_list)
366 : {
367 49379 : list<pair<size_t, size_t> > temp_list;
368 49379 : for(list<pair<size_t,size_t>>::iterator iter=range_list.begin();
369 206359 : iter!=range_list.end(); ++iter) {
370 156980 : AlwaysAssert((*iter).second >= (*iter).first, AipsError);
371 156980 : size_t width = (*iter).second - (*iter).first + 1;
372 156980 : if (width <= maxwidth) {
373 156979 : temp_list.push_back(*iter);
374 : }
375 : }
376 49379 : range_list.clear();
377 49379 : range_list.merge(temp_list);
378 49379 : }
379 :
380 50181 : void LineFinderUtils::rejectNarrowRange(size_t const minwidth,
381 : list<pair<size_t,size_t>>& range_list)
382 : {
383 50181 : AlwaysAssert(minwidth>0, AipsError);
384 50181 : list<pair<size_t, size_t> > temp_list;
385 50181 : for(list<pair<size_t,size_t>>::iterator iter=range_list.begin();
386 316074 : iter!=range_list.end(); ++iter) {
387 265893 : AlwaysAssert((*iter).second >= (*iter).first, AipsError);
388 265893 : size_t width = (*iter).second - (*iter).first + 1;
389 265893 : if (width >= minwidth) {
390 233881 : temp_list.push_back(*iter);
391 : }
392 : }
393 50181 : range_list.clear();
394 50181 : range_list.merge(temp_list);
395 50181 : }
396 :
397 : template size_t LineFinderUtils::binDataAndMask<float>(size_t const num_in,
398 : float const in_data[],
399 : bool const in_mask[],
400 : size_t const bin_size,
401 : size_t const num_out,
402 : float out_data[],
403 : bool out_mask[],
404 : size_t const offset,
405 : bool const keepsize);
406 :
407 : } //# NAMESPACE CASA - END
|
