Line data Source code
1 : //# tLineFinder.cc: this defines tests of SingleDishMS
2 : //#
3 : //# Copyright (C) 2014,2015
4 : //# National Astronomical Observatory of Japan
5 : //#
6 : //# This library is free software; you can redistribute it and/or modify it
7 : //# under the terms of the GNU Library General Public License as published by
8 : //# the Free Software Foundation; either version 2 of the License, or (at your
9 : //# option) any later version.
10 : //#
11 : //# This library is distributed in the hope that it will be useful, but WITHOUT
12 : //# ANY WARRANTY; without even the implied warranty of MERCHANTABILITY or
13 : //# FITNESS FOR A PARTICULAR PURPOSE. See the GNU Library General Public
14 : //# License for more details.
15 : //#
16 : //# You should have received a copy of the GNU Library General Public License
17 : //# along with this library; if not, write to the Free Software Foundation,
18 : //# Inc., 675 Massachusetts Ave, Cambridge, MA 02139, USA.
19 : //#
20 : //# Correspondence concerning AIPS++ should be addressed as follows:
21 : //# Internet email: aips2-request@nrao.edu.
22 : //# Postal address: AIPS++ Project Office
23 : //# National Radio Astronomy Observatory
24 : //# 520 Edgemont Road
25 : //# Charlottesville, VA 22903-2475 USA
26 : //#
27 : //# $Id$
28 : #include <iostream>
29 : #include <sstream>
30 : #include <list>
31 : #include <cassert>
32 :
33 :
34 : #include <casacore/casa/Arrays/Vector.h>
35 : #include <casacore/tables/Tables/Table.h>
36 : #include <casacore/tables/Tables/ArrayColumn.h>
37 : #include <casacore/tables/Tables/ScalarColumn.h>
38 :
39 : #include <casa_sakura/SakuraAlignedArray.h>
40 : #include <singledish/SingleDish/LineFindingUtils.h>
41 : #include <singledish/SingleDish/LineFinder.h>
42 :
43 : using namespace casacore;
44 : using namespace casa;
45 : using namespace std;
46 :
47 : template <typename DataType>
48 20 : string print_array(size_t const num_data, DataType const* data) {
49 20 : cout << "got " << num_data << " elements" << endl;
50 40 : ostringstream oss;
51 20 : oss << "[";
52 20 : if (num_data > 0)
53 20 : oss << data[0] ;
54 372 : for (size_t i=1; i<num_data; ++i)
55 352 : oss << ", " << data[i];
56 20 : oss << "]";
57 40 : return oss.str();
58 : }
59 :
60 20 : void print_line(list<pair<size_t,size_t>> &line_list) {
61 20 : cout << "Number of Lines: " << line_list.size() << endl;
62 100 : for (list<pair<size_t,size_t>>::iterator iter = line_list.begin();
63 180 : iter!=line_list.end(); ++iter) {
64 80 : cout << "- [ " << (*iter).first << ", " << (*iter).second << " ]" << endl;
65 : }
66 20 : }
67 :
68 1 : int main(int argc, char *argv[]) {
69 1 : bool dobin(true), domask(true), doline(true), domaskline(true), dolinefinder(false);
70 1 : if (dobin)
71 : {// Test binning
72 1 : cout << "\n**********\nTest binning\n**********" << endl;
73 1 : size_t const num_data=20;
74 : float data[num_data];
75 : bool mask[num_data];
76 21 : for (size_t i = 0; i<num_data; ++i) {
77 20 : data[i] = static_cast<float>(i);
78 20 : mask[i] = (i%4 != 0) ? true : false;
79 : }
80 1 : cout << "data = " << print_array<float>(num_data, data) << endl;
81 1 : cout << "mask = " << print_array<bool>(num_data, mask) << endl;
82 1 : size_t const bin_size = 3;
83 1 : size_t offset = 1;
84 1 : size_t num_out = (num_data-offset)/bin_size;
85 1 : bool keepsize=false;
86 1 : float out_data1[num_out];
87 1 : bool out_mask1[num_out];
88 1 : cout << "[Binning]" << endl;
89 1 : cout << "Bin = " << bin_size << ", offset = " << offset
90 1 : << ", keepsize = " << (keepsize ? "true" : "false") << endl;
91 1 : size_t nout1 = LineFinderUtils::binDataAndMask<float>(num_data, data, mask, bin_size, num_out, out_data1, out_mask1, offset, keepsize);
92 1 : cout << "returned" << endl;
93 1 : cout << "out_data = " << print_array<float>(nout1, out_data1) << endl;
94 1 : cout << "out_mask = " << print_array<bool>(nout1, out_mask1) << endl;
95 1 : keepsize = true;
96 1 : offset = 0;
97 : float out_data2[num_data];
98 : bool out_mask2[num_data];
99 1 : cout << "\n[keepsize]" << endl;
100 1 : cout << "Bin = " << bin_size << ", offset = " << offset
101 1 : << ", keepsize = " << (keepsize ? "true" : "false") << endl;
102 1 : size_t nout2 = LineFinderUtils::binDataAndMask<float>(num_data, data, mask, bin_size, num_data, out_data2, out_mask2, offset, keepsize);
103 1 : cout << "out_data = " << print_array<float>(nout2, out_data2) << endl;
104 2 : cout << "out_mask = " << print_array<bool>(nout2, out_mask2) << endl;
105 : }
106 :
107 1 : if (domask)
108 : {// Test mask creation
109 1 : cout << "\n**********\nTest mask creation\n**********" << endl;
110 1 : size_t const num_data = 20;
111 2 : SakuraAlignedArray<float> data(num_data);
112 2 : SakuraAlignedArray<float> mad(num_data);
113 2 : SakuraAlignedArray<bool> in_mask(num_data);
114 2 : SakuraAlignedArray<bool> out_mask(num_data);
115 1 : float* dptr = data.data();
116 1 : bool* mptr = in_mask.data();
117 21 : for (size_t i = 0 ; i < num_data; ++i) {
118 20 : dptr[i] = float(i);
119 20 : mptr[i] = (i%4 != 0) ? true : false;
120 : }
121 1 : LineFinderUtils::calculateMAD(num_data,data.data(),in_mask.data(),mad.data());
122 1 : cout << "[MAD array]" << endl;
123 1 : cout << "data = " << print_array<float>(num_data, data.data()) << endl;
124 1 : cout << "in_mask = " << print_array<bool>(num_data, in_mask.data()) << endl;
125 1 : cout << "MAD = " << print_array<float>(num_data, mad.data()) << endl;
126 :
127 1 : float threshold = 5.0;
128 1 : cout << "\n[Mask creation]" << endl;
129 1 : cout << "MAD = " << print_array<float>(num_data, mad.data()) << endl;
130 1 : cout << "in_mask = " << print_array<bool>(num_data, in_mask.data()) << endl;
131 1 : cout << "threshold = " << threshold << endl;
132 1 : LineFinderUtils::createMaskByAThreshold(num_data, mad.data(), in_mask.data(), threshold, out_mask.data());
133 1 : cout << "out_mask = " << print_array<bool>(num_data, out_mask.data()) << endl;
134 1 : cout << "\n[Sign creation]" << endl;
135 2 : Vector<int8_t> signval(num_data);
136 1 : LineFinderUtils::createSignByAThreshold(num_data, mad.data(),
137 : threshold,signval.data());
138 1 : cout << "data = " << print_array<float>(num_data, mad.data()) << endl;
139 1 : cout << "threshold = " << threshold << endl;
140 1 : cout << "sign = " << print_array<int8_t>(num_data, signval.data()) << endl;
141 :
142 1 : cout << "\n[Mask to line list]" << endl;
143 2 : list<pair<size_t,size_t>> lines;
144 1 : cout << "in_mask = " << print_array<bool>(num_data, in_mask.data()) << endl;
145 1 : LineFinderUtils::maskToRangesList(num_data, in_mask.data(), lines);
146 1 : print_line(lines);
147 :
148 21 : for (size_t i = 0; i<num_data; ++i) {
149 20 : in_mask.data()[i] = (i%4 != 1) ? true : false;
150 : }
151 1 : cout << "\nin_mask = " << print_array<bool>(num_data, in_mask.data()) << endl;
152 1 : LineFinderUtils::maskToRangesList(num_data, in_mask.data(), lines);
153 1 : print_line(lines);
154 : }
155 1 : if (doline)
156 : { // Test handling of line_list
157 1 : cout << "\n**********\nTest line list handling\n**********" << endl;
158 2 : list<pair<size_t,size_t>> mylines;
159 1 : mylines.push_back(pair<size_t,size_t>(0,2));
160 1 : mylines.push_back(pair<size_t,size_t>(5,10));
161 1 : mylines.push_back(pair<size_t,size_t>(20,70));
162 1 : mylines.push_back(pair<size_t,size_t>(72,80));
163 1 : cout << "[ORIGINAL]" << endl;
164 1 : print_line(mylines);
165 1 : LineFinderUtils::rejectWideRange(20,mylines);
166 1 : cout << "[Regect wide (>20)]" << endl;
167 1 : print_line(mylines);
168 1 : LineFinderUtils::rejectNarrowRange(3,mylines);
169 1 : cout << "[Regect narrow (<3)]" << endl;
170 1 : print_line(mylines);
171 1 : size_t binsize = 4;
172 1 : size_t offset = 1;
173 1 : LineFinderUtils::deBinRanges(binsize, offset, mylines);
174 1 : cout << "[debin (bin=" << binsize << ", offset=" << offset << ")]" << endl;
175 1 : print_line(mylines);
176 : // Merge overlap in a list
177 1 : cout << "\nMerge overlap in a list" << endl;
178 1 : mylines.clear();
179 1 : mylines.push_back(pair<size_t,size_t>(5,8));
180 1 : mylines.push_back(pair<size_t,size_t>(4,6));
181 1 : mylines.push_back(pair<size_t,size_t>(12,14));
182 1 : mylines.push_back(pair<size_t,size_t>(11,15));
183 1 : mylines.push_back(pair<size_t,size_t>(20,23));
184 1 : mylines.push_back(pair<size_t,size_t>(23,25));
185 1 : cout << "[ORIGINAL]" << endl;
186 1 : print_line(mylines);
187 1 : LineFinderUtils::mergeOverlappingRanges(mylines);
188 1 : cout << "[Merged]" << endl;
189 1 : print_line(mylines);
190 : // Merge overlap in two lists
191 1 : mylines.clear();
192 1 : mylines.push_back(pair<size_t,size_t>(5,10));
193 1 : mylines.push_back(pair<size_t,size_t>(15,20));
194 1 : mylines.push_back(pair<size_t,size_t>(25,30));
195 1 : mylines.push_back(pair<size_t,size_t>(35,40));
196 1 : mylines.push_back(pair<size_t,size_t>(45,50));
197 1 : mylines.push_back(pair<size_t,size_t>(55,60));
198 2 : list<pair<size_t,size_t>> new_lines;
199 1 : new_lines.push_back(pair<size_t,size_t>(0,2));
200 1 : new_lines.push_back(pair<size_t,size_t>(12,13));
201 1 : new_lines.push_back(pair<size_t,size_t>(22,27));
202 1 : new_lines.push_back(pair<size_t,size_t>(37,39));
203 1 : new_lines.push_back(pair<size_t,size_t>(42,52));
204 1 : new_lines.push_back(pair<size_t,size_t>(57,62));
205 1 : new_lines.push_back(pair<size_t,size_t>(67,69));
206 1 : cout << "[ORIGINAL]" << endl;
207 1 : cout << "to = ";
208 1 : print_line(mylines);
209 1 : cout << "from = ";
210 1 : print_line(new_lines);
211 1 : LineFinderUtils::mergeOverlapInTwoLists(mylines, new_lines);
212 1 : cout << "out = ";
213 1 : print_line(mylines);
214 :
215 1 : cout << "\nMultiple overlap" << endl;
216 1 : mylines.clear();
217 1 : mylines.push_back(pair<size_t,size_t>(5,10));
218 1 : mylines.push_back(pair<size_t,size_t>(15,20));
219 1 : mylines.push_back(pair<size_t,size_t>(25,30));
220 1 : mylines.push_back(pair<size_t,size_t>(35,40));
221 1 : mylines.push_back(pair<size_t,size_t>(45,50));
222 1 : mylines.push_back(pair<size_t,size_t>(55,60));
223 1 : new_lines.clear();
224 1 : new_lines.push_back(pair<size_t,size_t>(2,17));
225 1 : new_lines.push_back(pair<size_t,size_t>(22,42));
226 1 : new_lines.push_back(pair<size_t,size_t>(47,62));
227 1 : cout << "to = ";
228 1 : print_line(mylines);
229 1 : cout << "from = ";
230 1 : print_line(new_lines);
231 1 : LineFinderUtils::mergeOverlapInTwoLists(mylines, new_lines);
232 1 : cout << "out = ";
233 1 : print_line(mylines);
234 : }
235 :
236 1 : if(domaskline)
237 : {
238 1 : cout << "\n**********\nTest line list handling with mask\n**********" << endl;
239 1 : size_t num_data(20);
240 2 : Vector<int8_t> sign(num_data);
241 2 : Vector<bool> mask(num_data);
242 21 : for (size_t i = 0; i<num_data; ++i) {
243 20 : sign[i] = ((i/6) % 2)==0 ? 1:-1;
244 20 : mask[i] = (i % 5)==0 ? false : true;
245 : }
246 1 : cout << "[Extend by sign]" << endl;
247 1 : cout << "sign = " << print_array<int8_t>(num_data, sign.data()) << endl;
248 1 : cout << "in_mask = " << print_array<bool>(num_data, mask.data()) << endl;
249 2 : list<pair<size_t,size_t>> mylines;
250 1 : mylines.push_back(pair<size_t,size_t>(2,2));
251 1 : mylines.push_back(pair<size_t,size_t>(7,7));
252 1 : mylines.push_back(pair<size_t,size_t>(12,12));
253 1 : cout << "input line:" << endl;
254 1 : print_line(mylines);
255 1 : LineFinderUtils::extendRangeBySign(num_data, sign.data(), mask.data(), mylines);
256 1 : cout << "extended line:" << endl;
257 1 : print_line(mylines);
258 :
259 1 : cout << "[Merge Lines by false]" << endl;
260 1 : mylines.clear();
261 1 : mylines.push_back(pair<size_t,size_t>(5,10));
262 1 : mylines.push_back(pair<size_t,size_t>(16,18));
263 21 : for (size_t i = 0; i < num_data; ++i) {
264 20 : mask[i] = (i>mylines.front().second && i<mylines.back().first)? false : true;
265 : }
266 1 : size_t maxgap = 5;
267 1 : cout << "max gap = " << maxgap << endl;
268 1 : cout << "mask = " << print_array<bool>(num_data, mask.data()) << endl;
269 1 : cout << "input line = ";
270 1 : print_line(mylines);
271 1 : LineFinderUtils::mergeGapByFalse(num_data, mask.data(), maxgap, mylines);
272 1 : cout << "output line: ";
273 1 : print_line(mylines);
274 : // change max gap to merge. This time lines would not be merged.
275 1 : maxgap = 4;
276 1 : mylines.clear();
277 1 : mylines.push_back(pair<size_t,size_t>(5,10));
278 1 : mylines.push_back(pair<size_t,size_t>(16,18));
279 1 : cout << "max gap = " << maxgap << endl;
280 1 : cout << "mask = " << print_array<bool>(num_data, mask.data()) << endl;
281 1 : cout << "input line = ";
282 1 : print_line(mylines);
283 1 : LineFinderUtils::mergeGapByFalse(num_data, mask.data(), maxgap, mylines);
284 1 : cout << "output line: ";
285 1 : print_line(mylines);
286 : }
287 1 : if (dolinefinder) {
288 0 : cout << "\n**********\nTest line finding\n**********" << endl;
289 0 : string table_name = "/almadev/work/reg_test/IRC+10216_rawACSmod_cal";
290 : //string table_name = "/almadev/work/OrionS_rawACSmod_calTPave.asap";
291 0 : size_t row_idx = 0;
292 0 : cout << "Table: " << table_name << endl;
293 0 : cout << "idx: " << row_idx << endl;
294 0 : Table mytab(table_name, Table::Old);
295 0 : assert(row_idx<mytab.nrow());
296 0 : ScalarColumn<unsigned int> flagRCol(mytab, "FLAGROW");
297 0 : assert(flagRCol.get(row_idx)==0);
298 0 : ArrayColumn<float> specCol(mytab, "SPECTRA");
299 0 : ArrayColumn<uint8_t> flagCol(mytab, "FLAGTRA");
300 0 : Vector<float> specvec(specCol.get(row_idx));
301 0 : Vector<uint8_t> flagvec(flagCol.get(row_idx));
302 0 : size_t num_data(specvec.nelements());
303 0 : cout << "nchan: " << num_data << endl;
304 0 : Vector<float> data(num_data);
305 0 : Vector<bool> mask(num_data);
306 0 : for (size_t i=0; i<num_data; ++i) {
307 0 : data[i] = specvec[i];
308 0 : mask[i] = (flagvec[i]==static_cast<uint8_t>(0));
309 : }
310 0 : pair<size_t,size_t> edge(500,500);
311 : list<pair<size_t,size_t>> line_list = \
312 0 : linefinder::MADLineFinder(num_data, data.data(), mask.data(), 5.0, 10.0, 3, 1000, 4, edge);
313 0 : cout << "[Line finding result]" << endl;
314 0 : print_line(line_list);
315 : }
316 :
317 1 : return 0;
318 : }
|
