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 2 : 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 2 : if (offset+bin_size > num_in-1) return 0;
59 : // Can NOT offset (>0) when changing size (keepsize=false)
60 2 : AlwaysAssert((!keepsize || offset==0), AipsError);
61 2 : size_t num_bin = keepsize ? num_in : (num_in-offset)/bin_size;
62 2 : AlwaysAssert(num_bin <= num_out, AipsError);
63 2 : size_t num_loc_out = keepsize ? bin_size : 1;
64 2 : size_t out_idx = 0;
65 15 : for (size_t i=offset; (keepsize ? i : i+bin_size-1) < num_in && out_idx < num_bin;) {
66 13 : DataType data = static_cast<DataType>(0);
67 13 : bool mask = true;
68 13 : size_t iend=i+bin_size;
69 13 : size_t count=0;
70 51 : 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 38 : data += in_data[i];
74 38 : mask = mask && in_mask[i];
75 38 : ++count;
76 : }
77 13 : data = mask ? data/static_cast<DataType>(count) : static_cast<DataType>(0);
78 39 : for (size_t j = 0; j < num_loc_out && out_idx < num_bin; ++j){
79 26 : out_mask[out_idx] = mask;
80 26 : out_data[out_idx] = data;
81 26 : ++out_idx;
82 : }
83 : }
84 2 : AlwaysAssert(num_bin==out_idx, AipsError);
85 2 : return num_bin;
86 : }
87 :
88 1 : void LineFinderUtils::calculateMAD(size_t const num_data,
89 : float const in_data[],
90 : bool const in_mask[],
91 : float mad[])
92 : {
93 1 : float median_value = LineFinderUtils::maskedMedian(num_data, in_data, in_mask, 1.0);
94 : //cout << "median value for MAD = " << median_value << endl;
95 21 : for (size_t i = 0; i < num_data; ++i) {
96 20 : mad[i] = fabs(in_data[i]-median_value);
97 : }
98 1 : }
99 :
100 1 : 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 1 : LIBSAKURA_SYMBOL(Status) status = LIBSAKURA_SYMBOL(SetTrueIfGreaterThanOrEqualsFloat)(num_data, in_data, threshold, out_mask);
107 1 : AlwaysAssert(status == LIBSAKURA_SYMBOL(Status_kOK), AipsError);
108 21 : for (size_t i = 0; i < num_data; ++i) {
109 20 : out_mask[i] = out_mask[i] && in_mask[i];
110 : }
111 1 : }
112 :
113 1 : void LineFinderUtils::deBinRanges(size_t const bin_size, size_t const offset,
114 : list<pair<size_t,size_t>>& range_list)
115 : {
116 1 : for(list<pair<size_t,size_t>>::iterator iter=range_list.begin();
117 4 : iter!=range_list.end(); ++iter) {
118 3 : (*iter).first = (*iter).first*bin_size + offset;
119 3 : (*iter).second = (*iter).second*bin_size+(bin_size-1) + offset;
120 : }
121 1 : }
122 :
123 1 : void LineFinderUtils::extendRangeBySign(size_t num_sign, int8_t const* sign, bool const* mask,
124 : list<pair<size_t,size_t>>& range_list){
125 1 : for (list<pair<size_t,size_t>>::iterator iter = range_list.begin();
126 4 : iter!=range_list.end(); ++iter) {
127 3 : AlwaysAssert((*iter).first <= (*iter).second, AipsError);
128 : // extend left wing
129 3 : int8_t ref_sign = sign[(*iter).first];
130 5 : for (size_t i=(*iter).first-1; ; --i) {
131 5 : if (mask[i] && sign[i] == ref_sign) (*iter).first=i;
132 : else break;
133 2 : if (i==0) break;
134 : }
135 : // extend right wing
136 3 : ref_sign = sign[(*iter).second];
137 9 : for (size_t i=(*iter).second+1; i<num_sign; ++i) {
138 9 : if (mask[i] && sign[i] == ref_sign) (*iter).second=i;
139 : else break;
140 : }
141 : }
142 1 : }
143 :
144 :
145 2 : 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 2 : out_range.clear();
151 2 : if (num_mask==0) return;
152 2 : size_t istart=num_mask;
153 2 : if (mask[0]) istart=0;
154 40 : for (size_t i=1; i<num_mask; ++i) {
155 38 : int8_t diff = static_cast<int8_t>(mask[i]) - static_cast<int8_t>(mask[i-1]);
156 38 : if (diff==1) { // start of new line
157 10 : AlwaysAssert(istart==num_mask, AipsError);
158 10 : istart=i;
159 : }
160 28 : else if (diff==-1) { // end of line was i-1 chan
161 9 : AlwaysAssert(istart<num_mask, AipsError);
162 9 : out_range.push_back(pair<size_t,size_t>(istart, i-1));
163 9 : istart=num_mask;
164 : }
165 : }
166 2 : 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 1 : float LineFinderUtils::maskedMedian(size_t num_data, float const data[],
208 : bool const mask[], float fraction)
209 : {
210 :
211 1 : 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 1 : float *local_data_p = local_data.data();
216 1 : size_t num_valid(num_data+1);
217 1 : LIBSAKURA_SYMBOL(Status) status = LIBSAKURA_SYMBOL(SortValidValuesDenselyFloat)(num_data, mask, local_data_p, &num_valid);
218 1 : AlwaysAssert(status == LIBSAKURA_SYMBOL(Status_kOK), AipsError);
219 1 : AlwaysAssert(num_valid <= num_data, AipsError);
220 1 : if (fraction<1.0)
221 0 : num_valid = static_cast<size_t>(num_valid*fraction);
222 1 : float median_value = LineFinderUtils::getMedianOfSorted<float>(num_valid,
223 : local_data_p);
224 1 : return median_value;
225 :
226 1 : }
227 :
228 1 : void LineFinderUtils::mergeOverlappingRanges(list<pair<size_t,size_t>>& range_list)
229 : {
230 1 : if (range_list.size() < 2) return; // nothing to do
231 1 : list<pair<size_t,size_t>> temp_list;
232 1 : list<pair<size_t,size_t>>::iterator iter=range_list.begin();
233 1 : temp_list.push_back(*iter);
234 1 : ++iter;
235 6 : for ( ; iter!=range_list.end(); ++iter) {
236 : // always start searching from the begining of temp_list
237 : //size_t dummy_cursor =
238 5 : mergeARangeToList(temp_list, (*iter));
239 : }
240 1 : range_list.clear();
241 1 : range_list.splice(range_list.end(), temp_list);
242 1 : }
243 :
244 :
245 2 : void LineFinderUtils::mergeOverlapInTwoLists(list<pair<size_t,size_t>>& to,
246 : list<pair<size_t,size_t>>& from)
247 : {
248 2 : if (from.size()==0) return; // nothing to do
249 2 : if (to.size()==0) { // replace from with to
250 0 : to.splice(to.end(), from);
251 0 : return;
252 : }
253 2 : size_t cursor = 0;
254 2 : for (list<pair<size_t,size_t>>::iterator from_iter=from.begin();
255 12 : from_iter!=from.end(); ++from_iter) {
256 10 : cursor = mergeARangeToList(to, *from_iter, cursor);
257 : }
258 : }
259 :
260 15 : 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 15 : AlwaysAssert(new_range.first<=new_range.second, AipsError);
265 : // cursor should be zero if range_list is empty
266 15 : AlwaysAssert(range_list.size()>0 || cursor==0, AipsError);
267 15 : 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 15 : AlwaysAssert(cursor < range_list.size(), AipsError);
272 : // add range at the end
273 15 : if (new_range.first > range_list.back().second) {
274 3 : range_list.push_back(new_range);
275 3 : return range_list.size()-1;
276 : }
277 : // move iterator to the start point of the search
278 12 : list<pair<size_t,size_t>>::iterator start_cursor = range_list.begin();
279 12 : if (cursor > 0) {
280 19 : for (size_t i=1; i<cursor; ++i) {
281 13 : ++start_cursor;
282 : }
283 6 : AlwaysAssert(new_range.first > (*start_cursor).second, AipsError);
284 6 : ++start_cursor;
285 : }
286 : // add range to the begining
287 12 : if (new_range.second < (*start_cursor).first) {
288 1 : range_list.insert(start_cursor, new_range);
289 1 : return cursor+1;
290 : }
291 11 : 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 11 : list<pair<size_t,size_t>>::iterator insert_begin=range_list.end();
295 11 : bool any_overlap = false;
296 11 : for (list<pair<size_t,size_t>>::iterator iter=start_cursor;
297 22 : iter!=range_list.end(); ++iter, ++out_cursor) {
298 22 : 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 21 : } else if ((*iter).second < new_range.first) { //|| (*iter).first > new_range.second) {
303 : // still too small ranges
304 11 : continue;
305 : } else {//overlap started
306 10 : insert_begin = iter;
307 10 : any_overlap = true;
308 10 : break;
309 : }
310 : }
311 : // An overlap should be found.
312 10 : AlwaysAssert(any_overlap, AipsError);
313 : // find segment in range_list (final) list where overlap ends
314 10 : list<pair<size_t,size_t>>::iterator insert_end=range_list.end();
315 10 : for (list<pair<size_t,size_t>>::iterator iter=insert_begin;
316 23 : iter!=range_list.end(); ++iter) {
317 18 : if ((*iter).first <= new_range.second) {
318 13 : insert_end = iter;
319 : }
320 5 : else break; // the range exceed new_range.
321 : }
322 : // there was an overlap
323 10 : list<pair<size_t,size_t>> temp_list;
324 : // preceeding elements (including the first overlapped range)
325 10 : temp_list.splice(temp_list.end(), range_list, range_list.begin(), insert_begin);
326 : // overlaping segment
327 10 : pair<size_t,size_t> unified_range(std::min((*insert_begin).first, new_range.first),
328 20 : std::max((*insert_end).second, new_range.second));
329 10 : temp_list.push_back(unified_range);
330 : // trailing elements
331 10 : if (insert_end != range_list.end()){
332 10 : temp_list.splice(temp_list.end(),range_list,++insert_end, range_list.end());
333 : }
334 10 : range_list.clear();
335 10 : range_list.splice(range_list.end(), temp_list);
336 10 : return out_cursor;
337 10 : }
338 :
339 0 : void LineFinderUtils::mergeSmallGapByFraction(double const fraction,
340 : size_t const maxwidth,
341 : list<pair<size_t,size_t>>& range_list)
342 : {
343 0 : if (range_list.size() < 2) return; // nothing to do
344 0 : list<pair<size_t,size_t>> temp_list;
345 0 : list<pair<size_t,size_t>>::iterator iter=range_list.begin();
346 0 : temp_list.push_back(*iter);
347 0 : ++iter;
348 0 : for (; iter != range_list.end(); ++iter) {
349 0 : size_t curr_width = (*iter).second - (*iter).first;
350 0 : size_t prev_width = temp_list.back().second - temp_list.back().first;
351 0 : size_t gap = (*iter).first-temp_list.back().second;
352 0 : if (gap < static_cast<size_t>(std::min(prev_width, curr_width)*fraction) && \
353 0 : (*iter).second-temp_list.back().first < maxwidth) {
354 0 : temp_list.back().second = (*iter).second;
355 : }
356 : else {
357 0 : temp_list.push_back(*iter);
358 : }
359 : }
360 0 : range_list.clear();
361 0 : range_list.splice(range_list.end(), temp_list);
362 0 : }
363 :
364 1 : void LineFinderUtils::rejectWideRange(size_t const maxwidth,
365 : list<pair<size_t,size_t>>& range_list)
366 : {
367 1 : list<pair<size_t, size_t> > temp_list;
368 1 : for(list<pair<size_t,size_t>>::iterator iter=range_list.begin();
369 5 : iter!=range_list.end(); ++iter) {
370 4 : AlwaysAssert((*iter).second >= (*iter).first, AipsError);
371 4 : size_t width = (*iter).second - (*iter).first + 1;
372 4 : if (width <= maxwidth) {
373 3 : temp_list.push_back(*iter);
374 : }
375 : }
376 1 : range_list.clear();
377 1 : range_list.merge(temp_list);
378 1 : }
379 :
380 1 : void LineFinderUtils::rejectNarrowRange(size_t const minwidth,
381 : list<pair<size_t,size_t>>& range_list)
382 : {
383 1 : AlwaysAssert(minwidth>0, AipsError);
384 1 : list<pair<size_t, size_t> > temp_list;
385 1 : for(list<pair<size_t,size_t>>::iterator iter=range_list.begin();
386 4 : iter!=range_list.end(); ++iter) {
387 3 : AlwaysAssert((*iter).second >= (*iter).first, AipsError);
388 3 : size_t width = (*iter).second - (*iter).first + 1;
389 3 : if (width >= minwidth) {
390 3 : temp_list.push_back(*iter);
391 : }
392 : }
393 1 : range_list.clear();
394 1 : range_list.merge(temp_list);
395 1 : }
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
|
