Elaboradar 0.1
Caricamento in corso...
Ricerca in corso...
Nessun risultato
cleaner.cpp
1/*
2 * =================================================================================
3 *
4 * Filename: volume_cleaner.cpp
5 *
6 * Description:
7 *
8 * Version: 1.0
9 * Created: 18/02/2014 12:19:44
10 * Revision: none
11 * Compiler: gcc
12 *
13 * Author: YOUR NAME (),
14 * Organization:
15 *
16 * =================================================================================
17 */
18#include "cleaner.h"
19#include "elabora_volume.h"
20#include "radarelab/image.h"
21#include "radarelab/matrix.h"
22
23#include <string>
24#include <filesystem>
25#include <unistd.h>
26#include <cstring>
27
28namespace radarelab {
29namespace algo {
30
31using namespace std;
32using namespace radarelab;
33 //using std::filesystem::current_path;
34
35//--------------------------------------------------------------------------------
36// These methods use only VRAD and WRAD values to clean the beam
37//
38std::vector<bool> Cleaner::clean_beam(const Eigen::VectorXd& beam_z, const Eigen::VectorXd& beam_w, const Eigen::VectorXd& beam_v, int iray) const
39{
40 const unsigned beam_size = beam_z.rows();
41 vector<bool> res(beam_size, false);
42 bool in_a_segment = false;
43 unsigned start, end;
44 unsigned segment_length;
45 bool before, after;
46 unsigned counter = 0;
47
48 for (unsigned ibin = 0; ibin < beam_size; ++ibin)
49 {
50//printf(" %4d %4d %6.2f %6.2f %10.6f ",iray,ibin , beam_z(ibin),beam_v(ibin),beam_w(ibin));
51 if (!in_a_segment)
52 {
53 /* cerco la prima cella segmento da pulire*/
54 if (beam_w(ibin) == W_threshold && beam_v(ibin) == bin_wind_magic_number)
55 {
56//printf(" 1 ----- START SEGMENT ------");
57 in_a_segment = true;
58 start = ibin;
59 after = false;
60 before = false;
61 }
62// else printf(" 0 ");
63 } else {
64 /* cerco la fine segmento da pulire*/
65 if (beam_w(ibin) != W_threshold || beam_v(ibin) != bin_wind_magic_number || ibin == (beam_size - 1))
66 {
67 in_a_segment = false;
68 end = ibin - 1;
69 if (ibin == (beam_size - 1)) end = ibin; // caso particolare per fine raggio
70 /* Fine trovata ora procedo alla pulizia eventuale */
71 segment_length = end - start;
72 counter = counter + (unsigned)(segment_length);
73
74 unsigned c_b=0;
75 unsigned c_a=0;
76 /* Cerco dati validi in Z prima del segmento */
77 for (int ib = ibin - 12; ib < (signed)ibin; ++ib)
78 if (ib >= 0 && beam_z(ib) > Z_missing)
79 c_b++;
80 if (c_b > 0.25*12) before = true;
81
82 /* Cerco dati validi in Z dopo il segmento */
83 for (unsigned ia = ibin + 1; ia <= ibin + 12; ++ia)
84 if (ia < beam_size && beam_z(ia) >= Z_missing)
85 c_a++;
86 if (c_a > 0.25*12) after = true;
87
88//printf(" 0 ----- STOP SEGMENT ------ %4d -- %4d before %d %d after %d %d ",segment_length,counter, c_b,before, c_a, after);
89 if ((segment_length >= min_segment_length && !before && !after) ||
90 segment_length >= max_segment_length || counter > 100)
91 {
92 /* qui pulisco */
93 // printf (" pulisco %d %d %d \n",segment_length, min_segment_length, max_segment_length);
94 for (unsigned ib = start; ib <= end; ++ib)
95 if( beam_z(ib) > Z_missing)res[ib] = true;
96 }
97 }
98// else printf(" 1 ");
99 }
100//printf("\n");
101 }
102 return res;
103}
104
105std::vector<unsigned char> Cleaner::eval_clean_beam(const Eigen::VectorXd& beam_z, const Eigen::VectorXd& beam_w, const Eigen::VectorXd& beam_v, int iray) const
106{
107 const unsigned beam_size = beam_z.rows();
108 vector<unsigned char> res(beam_size, 0);
109 bool in_a_segment = false;
110 unsigned start = 0, end = 0;
111 unsigned segment_length;
112 bool before = false, after = false;
113 unsigned counter = 0;
114
115 for (unsigned ibin = 0; ibin < beam_size; ++ibin)
116 {
117//printf(" %4d %4d %6.2f %6.2f %10.6f ",iray,ibin , beam_z(ibin),beam_v(ibin),beam_w(ibin));
118 if (!in_a_segment)
119 {
120 /* search the first radar bin's segment to be cleaned*/
121 if (beam_w(ibin) == W_threshold && beam_v(ibin) == bin_wind_magic_number)
122 {
123//printf(" 1 ----- START SEGMENT ------");
124 in_a_segment = true;
125 start = ibin;
126 after = false;
127 before = false;
128 }
129 } else {
130 /* search the last radar bin's segment to be cleaned*/
131 if (beam_w(ibin) != W_threshold || beam_v(ibin) != bin_wind_magic_number || ibin == (beam_size - 1))
132 {
133 in_a_segment = false;
134 end = ibin - 1;
135 if (ibin == (beam_size - 1)) end = ibin; // beam ended
136 /* Fine trovata ora procedo alla pulizia eventuale */
137 segment_length = end - start;
138 counter = counter + (unsigned)(segment_length);
139
140 unsigned c_b=0;
141 unsigned c_a=0;
142 /* Cerco dati validi in Z prima del segmento */
143 for (int ib = ibin - 12; ib < (signed)ibin; ++ib)
144 if (ib >= 0 && beam_z(ib) > Z_missing)
145 c_b++;
146 if (c_b > 0.25*12) before = true;
147
148 /* Cerco dati validi in Z dopo il segmento */
149 for (unsigned ia = ibin + 1; ia <= ibin + 12; ++ia)
150 if (ia < beam_size && beam_z(ia) >= Z_missing)
151 c_a++;
152 if (c_a > 0.25*12) after = true;
153
154//printf(" 0 ----- STOP SEGMENT ------ %4d -- %4d before %d %d after %d %d ",segment_length,counter, c_b,before, c_a, after);
155 if ((segment_length >= min_segment_length && !before && !after) ||
156 segment_length >= max_segment_length || counter > 100)
157 {
158 /* qui pulisco */
159 // printf (" pulisco %d %d %d \n",segment_length, min_segment_length, max_segment_length);
160 for (unsigned ib = start; ib <= end; ++ib)
161 if( beam_z(ib) > Z_missing)res[ib] = 1;
162 }
163 }
164// else printf(" 1 ");
165 }
166//printf("\n");
167 }
168 return res;
169}
170//---------------------------------------------------------------------------------
171
172//---------------------------------------------------------------------------------
173// This method use VRAD, WRAD, sdDBZH, sdZDR values to clean the beam
174//
175std::vector<bool> Cleaner::clean_beam(const Eigen::VectorXd& beam_z, const Eigen::VectorXd& beam_w, const Eigen::VectorXd& beam_v, const Eigen::VectorXd& beam_sd, const Eigen::VectorXd& beam_sdzdr, PolarScan<double>& scan_z, PolarScan<double>& scan_w, PolarScan<double>& scan_v, PolarScan<double>& SD, int iray) const
176{
177 const unsigned beam_size = beam_z.rows();
178 vector<bool> res(beam_size, false);
179 bool in_a_segment = false;
180 unsigned start = 0, end;
181 unsigned segment_length;
182 bool before, after;
183 unsigned counter = 0;
184 unsigned counter_trash = 0;
185 unsigned counter_clutter =0;
186 for (unsigned ibin = 0; ibin < beam_size; ++ibin)
187 {
188 bool is_clutter = false;
189 bool is_trash = false;
190 unsigned flag = 0 ;
191// In our systems (ARPA ER) interferences and other non meteo echo are characterised by the following steps
192//
193// 1) Wind is not defined ad spectrumWidth is 0. with Z defined.
194 if ( beam_w(ibin) == W_threshold && beam_v(ibin) == bin_wind_magic_number && beam_z (ibin) != Z_missing ) {
195// 2) Std<Dev of ZDR coulb be close to 0
196 if( beam_sdzdr(ibin) <= 0.01 ){
197 is_trash = true;
198 flag=2;
199 } else {
200 if (beam_z (ibin) >= 45. ){
201// 2) inside thunderstorms (Z > 45) StdDev of Zdr and StdDev of Z are quite high)
202 if ((ibin >100 && double(counter_trash)/double(ibin) >=0.5 && ( beam_sdzdr(ibin) >1 || beam_sd (ibin) > 5. )) ||
203 (beam_sdzdr(ibin) >4.0 && beam_sd (ibin) > 20.) ) {
204 is_trash = true;
205 flag=2;
206 } else {
207 is_trash = false;
208 flag=0;
209 }
210 } else if ( (ibin >100 && double(counter_trash)/double(ibin) >=0.5 && ( beam_sdzdr(ibin) >1 || beam_sd (ibin) > 5. )) ||
211 (beam_sd (ibin) >2. && (beam_sdzdr(ibin) >2.0 || beam_sd (ibin) > 10. )) ) {
212// 2) outside thunderstorms (Z > 45) StdDev of Zdr and StdDev of Z are lower
213 is_trash = true;
214 flag=2;
215 }
216 }
217 } else {
218// 3) Clutter is characterised by low value of VRAD and WRAD
219 if ((beam_w(ibin) * fabs(beam_v(ibin))) <= 0.25 && beam_z (ibin) != Z_missing ) {
220 is_clutter = true;
221 flag = 1;
222 }
223 }
224 if( is_clutter) counter_clutter ++;
225 if( is_trash ) counter_trash ++;
226 //if(ibin <40 && false){
227 //printf(" %4d %4d %6.2f %6.2f %10.6f %6.2f %6.2f ",iray,ibin , beam_z(ibin),beam_v(ibin),beam_w(ibin), beam_sd(ibin),beam_sdzdr(ibin));
228 //printf(" ----- %2x %2x %2x %2x ",(unsigned char)((beam_z(ibin)-scan_z.offset)/scan_z.gain/256),
229 //(unsigned char)((beam_v(ibin)-scan_v.offset)/scan_v.gain/256),
230 //(unsigned char)((beam_w(ibin)-scan_w.offset)/scan_w.gain/256),
231 //(unsigned char)((beam_sd(ibin)-SD.offset)/SD.gain/256));
232 //}
233 if (!in_a_segment)
234 {
235 /* cerco la prima cella segmento da pulire*/
236 if ( is_clutter || is_trash )
237 {
238// if(ibin <40)printf(" %1d ----- START SEGMENT ------",flag);
239 in_a_segment = true;
240 start = ibin;
241 after = false;
242 before = false;
243 }
244// else if(ibin <40)printf(" %1d ",flag);
245 } else {
246 /* cerco la fine segmento da pulire*/
247 if ( ! (is_clutter || is_trash ) || ibin == (beam_size - 1))
248 {
249 in_a_segment = false;
250 end = ibin - 1;
251 if (ibin == (beam_size - 1)) end = ibin; // caso particolare per fine raggio
252 /* Fine trovata ora procedo alla pulizia eventuale */
253 segment_length = end - start+1;
254 counter = counter + (unsigned)(segment_length);
255
256/* il segmento è corto allora cerco nei dintorni dei dati validi, se li trovo non pulisco */
257 if (segment_length <= 2*min_segment_length ){
258 /* Cerco dati validi in Z prima del segmento */
259 int count=0;
260 for (int ib = ibin - 2*min_segment_length; ib < (signed)ibin; ++ib)
261 if (ib >= 0 && (beam_z(ib) > Z_missing && beam_w(ib) != W_threshold && ( beam_w(ib) > 0.5 || fabs(beam_v(ib)) > 0.5) ) )
262 count++;
263 if (double(count)/double(min(int(ibin),int(2*min_segment_length))) >=0.25) before = true;
264
265 /* Cerco dati validi in Z dopo il segmento */
266 count = 0;
267 for (unsigned ia = ibin + 1; ia <= ibin + 2*min_segment_length; ++ia)
268 if (ia < beam_size && (beam_z(ia) > Z_missing && (beam_w(ia) != W_threshold && ( beam_w(ia) > 0.5 || fabs(beam_v(ia)) > 0.5)) ))
269 count ++;
270 if (double(count)/double(min(int(beam_size - ibin),int(2*min_segment_length))) >=0.25) after = true;
271 }
272// if(ibin <40)printf(" %1d ----- STOP SEGMENT ------ %4d -- %4d before %d after %d ",flag, segment_length,counter, before,after);
273 if ((segment_length >= min_segment_length && (!before || !after) ) || segment_length >= max_segment_length)
274 // if ((segment_length >= min_segment_length ) || segment_length >= max_segment_length)
275 {
276 /* qui pulisco */
277// if(ibin <40)printf (" pulisco %d %d %d \n",segment_length, min_segment_length, max_segment_length);
278 for (unsigned ib = start; ib <= end; ++ib)
279 res[ib] = true;
280 }
281 }
282// else if(ibin <40)printf(" %1d ",flag);
283
284 }
285// if(ibin <40)printf(" %4d %4d \n",counter_clutter,counter_trash);
286 }
287 return res;
288}
289//---------------------------------------------------------------------------------
290
291//---------------------------------------------------------------------------------
292// These methods use VRAD, WRAD, sdDBZH, values to clean the beam
293//
294std::vector<bool> Cleaner::clean_beam(const Eigen::VectorXd& beam_z, const Eigen::VectorXd& beam_w, const Eigen::VectorXd& beam_v, const Eigen::VectorXd& beam_sd, PolarScan<double>& scan_z, PolarScan<double>& scan_w, PolarScan<double>& scan_v, PolarScan<double>& SD, int iray) const
295{
296 const unsigned beam_size = beam_z.rows();
297 vector<bool> res(beam_size, false);
298 bool in_a_segment = false;
299 unsigned start, end;
300 unsigned segment_length;
301 bool before, after;
302 unsigned counter = 0;
303
304 for (unsigned ibin = 0; ibin < beam_size; ++ibin)
305 {
306// printf(" %4d %4d %6.2f %6.2f %10.6f %6.2f ",iray,ibin , beam_z(ibin),beam_v(ibin),beam_w(ibin), beam_sd(ibin));
307//printf(" ----- %2x %2x %2x %2x ",(unsigned char)((beam_z(ibin)-scan_z.offset)/scan_z.gain/256),
308//(unsigned char)((beam_v(ibin)-scan_v.offset)/scan_v.gain/256),
309//(unsigned char)((beam_w(ibin)-scan_w.offset)/scan_w.gain/256),
310//(unsigned char)((beam_sd(ibin)-SD.offset)/SD.gain/256));
311 if (!in_a_segment)
312 {
313 /* cerco la prima cella segmento da pulire*/
314 if ( ((beam_w(ibin) == W_threshold && beam_v(ibin) == bin_wind_magic_number) ||(beam_w(ibin) * fabs(beam_v(ibin)) <= 0.25) ) && beam_z (ibin) != Z_missing && beam_sd(ibin) > sd_threshold )
315 {
316// printf(" 1 ----- START SEGMENT ------");
317 in_a_segment = true;
318 start = ibin;
319 after = false;
320 before = false;
321 }
322// else printf(" 0 ");
323 } else {
324 /* cerco la fine segmento da pulire*/
325 if ( ( ( beam_w(ibin) != W_threshold || beam_v(ibin) != bin_wind_magic_number) && (beam_w(ibin) * fabs(beam_v(ibin)) > 0.25) ) || ibin == (beam_size - 1) || beam_z(ibin) == Z_missing || beam_sd(ibin) <= sd_threshold )
326 {
327 in_a_segment = false;
328 end = ibin - 1;
329 if (ibin == (beam_size - 1)) end = ibin; // caso particolare per fine raggio
330 /* Fine trovata ora procedo alla pulizia eventuale */
331 segment_length = end - start+1;
332 counter = counter + (unsigned)(segment_length);
333
334/* il segmento è corto allora cerco nei dintorni dei dati validi, se li trovo non pulisco */
335 if (segment_length <= 2*min_segment_length ){
336 /* Cerco dati validi in Z prima del segmento */
337 int count=0;
338 for (int ib = ibin - 2*min_segment_length; ib < (signed)ibin; ++ib)
339 if (ib >= 0 && (beam_z(ib) > Z_missing && beam_w(ib) != W_threshold && ( beam_w(ib) > 0.5 || fabs(beam_v(ib)) > 0.5) ) )
340 count++;
341 if (double(count)/double(min(int(ibin),int(2*min_segment_length))) >=0.25) before = true;
342
343 /* Cerco dati validi in Z dopo il segmento */
344 count = 0;
345 for (unsigned ia = ibin + 1; ia <= ibin + 2*min_segment_length; ++ia)
346 if (ia < beam_size && (beam_z(ia) > Z_missing && (beam_w(ia) != W_threshold && ( beam_w(ia) > 0.5 || fabs(beam_v(ia)) > 0.5)) ))
347 count ++;
348 if (double(count)/double(min(int(beam_size - ibin),int(2*min_segment_length))) >=0.25) after = true;
349 }
350// printf(" 0 ----- STOP SEGMENT ------ %4d -- %4d before %d after %d ",segment_length,counter, before,after);
351 if ((segment_length >= min_segment_length && (!before || !after) ) || segment_length >= max_segment_length)
352 // if ((segment_length >= min_segment_length ) || segment_length >= max_segment_length)
353 {
354 /* qui pulisco */
355 // printf (" pulisco %d %d %d \n",segment_length, min_segment_length, max_segment_length);
356 for (unsigned ib = start; ib <= end; ++ib)
357 res[ib] = true;
358 }
359 }
360// else printf(" 1 ");
361
362 }
363// printf("\n");
364 }
365 return res;
366}
367
368// CC: fuzzy logic
369 tuple<std::vector<unsigned char>,std::vector<double>> Cleaner::eval_classID_beam(const Eigen::VectorXd& beam_z, const Eigen::VectorXd& beam_w, const Eigen::VectorXd& beam_v, const Eigen::VectorXd& beam_sd, const Eigen::VectorXd& beam_zdr, const Eigen::VectorXd& beam_rohv, const Eigen::VectorXd& beam_sqi, const Eigen::VectorXd& beam_snr, const Eigen::VectorXd& beam_zvd, const Eigen::VectorXd& beam_sdray, const Eigen::VectorXd& beam_sdaz, const Eigen::VectorXd& beam_zdr_sd, int iray, const string radar, double v_ny, const char* fuzzy_path, bool stamp, bool force_meteo) const
370{
371
372 const unsigned beam_size = beam_z.rows();
373 vector<unsigned char> res(beam_size, 0);
374 vector<double> diffs(beam_size, 0);
375 unsigned prevID = 4;
376 int Num_entries=0;
377 int Num_echoes = 5;
378 int Ntraps = 5; // 5 argomenti da passare a Trap : x1,x2,x3,x4,x5
379 //char f_dir[256];
380 //char *cwd = get_current_dir_name();
381 //strcpy(f_dir,cwd);
382
383 //string f_dir = cwd;
384 string f_dir = fuzzy_path;
385 //f_dir = f_dir +"/dati";
386
387 //leggo matrice dei pesi----------------------------------------------------
388 string fin_w = f_dir+"/matrix-"+radar+".txt";//strcat(f_dir,wname.c_str());
389 //cout<<"leggo pesi da "<<fin_w<<endl;
390 vector<string> w_vector;
391 w_vector = read_matrix_from_txt(fin_w);
392 Num_entries = w_vector.size()/Num_echoes;
393
394 Matrix2D<double> Wij(Num_echoes,Num_entries);
395 for(int i=0;i<Num_echoes;i++){ //itero colonna
396 for(int j=0;j<Num_entries;j++){ //itero rriga
397 Wij(i,j) = stod( w_vector[i*Num_entries+j]);
398 //cout<<" W["<<i<<","<<j<<"]="<<Wij(i,j);
399 }
400 //cout<<" "<<endl;
401 }
402
403 //--------------------------------------------------------------------------
404 //leggo matrice dei traps
405 string fin_t = f_dir+"/Trap-"+radar+".txt";
406 vector<string> t_vector;
407 t_vector = read_matrix_from_txt(fin_t);
408 double Traps[Num_entries][Num_echoes][Ntraps];
409 for(int i=0;i<Num_entries;i++){
410 for(int j=0;j<Num_echoes;j++){
411 for(int k=0;k<Ntraps;k++){
412 Traps[i][j][k] = stod( t_vector[i*(Num_echoes*Ntraps)+j*Num_echoes+k]);
413 //cout<<" Traps["<<i<<","<<j<<","<<k<<"]="<<Traps[i][j][k];
414 }
415 //cout<<endl;
416 }
417 //cout<<endl;
418 }
419 //cout<<"T shape"<<Traps[0].size()<<endl;
420 //------------------------------------------------------------------------
421
422 vector<unsigned> counter (Num_entries,0) ; // non sono sicura di cosa delle dimensioni di questo counter
423 for (unsigned ibin = 0; ibin < beam_size; ++ibin)
424 //for (unsigned ibin = 87*4; ibin < 87*4+12; ++ibin)
425 {
426 //printf(" %4d %4d %6.2f %6.2f %10.6f %10.6f %10.6f %10.6f %10.6f",iray,ibin , beam_z(ibin),beam_v(ibin),beam_w(ibin),beam_zdr(ibin),beam_sd(ibin),beam_sdray(ibin),beam_sdaz(ibin));
427 Matrix2D<double> Pij(Num_echoes,Num_entries);
428 Pij = Pij * 0.;
429 //cout<<Pij<<endl;
430 ArrayXd Class_WP(Num_echoes); // la dimensione di Class_WP deve essere Num_echoes, perchè alla fine ricavo un vettore con 6 valori, uno per ogni echo, passando per prodotto pesi prob
431 Class_WP.setZero();
432 if (beam_z(ibin) == Z_missing) {
433 unsigned ID=0;
434 res[ibin]=ID;
435 counter[ID]++;
436
437 continue;
438 }
439
440 //cout<<"V Nyquist: "<<v_ny<<endl;
441 //cout<<sizeof(Traps)<<" "<<Wij.size()<<" "<<Pij.size()<<endl;
442//eseguo un test unico per VRAD e WRAD e assegno tutte le prob assieme
443 if (beam_v(ibin) != bin_wind_magic_number ) { // VRAD
444 //Pij(0,1)=trap(v_ny*Traps[0][0][0], v_ny*Traps[0][0][1], v_ny*Traps[0][0][2], v_ny*Traps[0][0][3], beam_v(ibin),v_ny*Traps[0][0][4]); // METEO
445 if(radar=="GAT"){
446 double pv0 = trap(v_ny*Traps[0][0][0], v_ny*Traps[0][0][1], v_ny*Traps[0][0][2], v_ny*Traps[0][0][3], beam_v(ibin));
447 double pv3 = trap(v_ny*Traps[0][0][4],v_ny*Traps[0][0][4],v_ny*Traps[0][0][4],v_ny*Traps[0][0][4],beam_v(ibin));
448 Pij(0,1)=std::max(pv0,pv3);
449 }else Pij(0,1)=trap(v_ny*Traps[0][0][0], v_ny*Traps[0][0][1], v_ny*Traps[0][0][2], v_ny*Traps[0][0][3], beam_v(ibin),v_ny*Traps[0][0][4]);
450
451 double prob_v = trap (v_ny, v_ny, v_ny, v_ny, beam_v(ibin));
452 //cout<<"prob_v computed: "<<prob_v<<endl;
453 //for(int e=1;e<Num_echoes;e++){ Pij(e,1) = prob_v };
454 if(radar=="GAT"){
455 double pv0 = trap(v_ny*Traps[0][1][0], v_ny*Traps[0][1][1], v_ny*Traps[0][1][2], v_ny*Traps[0][1][3], beam_v(ibin));
456 double pv1 = trap(-8.,-8.,-7.,-7.,beam_v(ibin));
457 double pv2 = trap(5.5,5.5,6.5,6.5,beam_v(ibin));
458 double pv3 = trap(v_ny*Traps[0][1][4],v_ny*Traps[0][1][4],v_ny*Traps[0][1][4],v_ny*Traps[0][1][4],beam_v(ibin));
459 Pij(1,1)=std::max({pv0,pv1,pv2,pv3});
460 }else Pij(1,1)=trap(v_ny*Traps[0][1][0], v_ny*Traps[0][1][1], v_ny*Traps[0][1][2], v_ny*Traps[0][1][3], beam_v(ibin),v_ny*Traps[0][1][4]);
461 // CLUTTER
462 Pij(2,1)=prob_v; // INTERF. Strong
463 Pij(3,1)=prob_v; // INTERF. Med.
464 Pij(4,1)=prob_v; // INTERF. Weak
465 } else {
466 Pij(0,1)=0.3; // METEO
467 Pij(1,1)=1.; // CLUTTER
468 Pij(2,1)=1.; // INTERF. Strong
469 Pij(3,1)=1.; // INTERF. Med.
470 Pij(4,1)=1.; // INTERF. Weak
471 }
472
473 // WRAD
474 Pij(0,2)= trap(Traps[1][0][0],Traps[1][0][1],Traps[1][0][2],Traps[1][0][3],beam_w(ibin)); // METEO
475 Pij(1,2)= trap(Traps[1][1][0],Traps[1][1][1],Traps[1][1][2],Traps[1][1][3],beam_w(ibin)); // CLUTTER
476 double prob_w = trap(Traps[1][2][0],Traps[1][2][1],Traps[1][2][2],Traps[1][2][3],beam_w(ibin));
477 Pij(2,2) = prob_w; // INTERF MULTIPLE
478 Pij(3,2) = prob_w; // INTERF. Med.
479 Pij(4,2) = prob_w; // INTERF. Weak
480
481// METEO
482 Pij(0,0) = trap(Traps[2][0][0],Traps[2][0][1],Traps[2][0][2],Traps[2][0][3],beam_z(ibin), Traps[2][0][4]); // Z
483 Pij(0,3) = trap (Traps[3][0][0],Traps[3][0][1],Traps[3][0][2],Traps[3][0][3], beam_sd(ibin)); // SD_2D
484 Pij(0,4) = trap (Traps[4][0][0],Traps[4][0][1],Traps[4][0][2],Traps[4][0][3], beam_sdray(ibin)); // SD_RAY
485 Pij(0,5) = trap (Traps[5][0][0],Traps[5][0][1],Traps[5][0][2],Traps[5][0][3], beam_sdaz(ibin)); // SD_AZ
486 //if(Num_entries>6)
487 Pij(0,6) = trap(Traps[6][0][0],Traps[6][0][1],Traps[6][0][2],Traps[6][0][3],beam_zdr(ibin)); // ZDR
488 Pij(0,7) = trap(Traps[7][0][0],Traps[7][0][1],Traps[7][0][2],Traps[7][0][3],beam_rohv(ibin),Traps[7][0][4]); // ROHV
489 Pij(0,8) = trap(Traps[8][0][0],Traps[8][0][1],Traps[8][0][2],Traps[8][0][3],beam_zdr_sd(ibin)); // ZDR_SD2D
490 Pij(0,9) = trap(Traps[9][0][0],Traps[9][0][1],Traps[9][0][2],Traps[9][0][3], beam_sqi(ibin)); // SQI
491 Pij(0,10) = trap(Traps[10][0][0],Traps[10][0][1],Traps[10][0][2],Traps[10][0][3], beam_snr(ibin)); // SNR
492 if((radar=="GAT")&&(ibin<=120))
493 Pij(0,11) = trap(Traps[11][0][0],0.,Traps[11][0][2],Traps[11][0][3],beam_zvd(ibin));
494 else
495 Pij(0,11) = trap(Traps[11][0][0],Traps[11][0][1],Traps[11][0][2],Traps[11][0][3],beam_zvd(ibin)); // DBZH_VD
496
497// CLUTTER
498 Pij(1,0) = trap (Traps[2][1][0],Traps[2][1][1],Traps[2][1][2],Traps[2][1][3], beam_z(ibin), Traps[2][1][4]); // Z
499 Pij(1,3) = trap (Traps[3][1][0],Traps[3][1][1],Traps[3][1][2],Traps[3][1][3], beam_sd(ibin)); // SD_2D
500 Pij(1,4) = trap (Traps[4][1][0],Traps[4][1][1],Traps[4][1][2],Traps[4][1][3], beam_sdray(ibin)); // SD_RAY
501 Pij(1,5) = trap (Traps[5][1][0],Traps[5][1][1],Traps[5][1][2],Traps[5][1][3], beam_sdaz(ibin)); // SD_AZ
502 //if(Num_entries>6)
503 Pij(1,6) = trap(Traps[6][1][0],Traps[6][1][1],Traps[6][1][2],Traps[6][1][3],beam_zdr(ibin),Traps[6][1][4] ); // ZDR
504 Pij(1,7) = trap(Traps[7][1][0],Traps[7][1][1],Traps[7][1][2],Traps[7][1][3],beam_rohv(ibin),Traps[7][1][4]); // ROHV
505 Pij(1,8) = trap(Traps[8][1][0],Traps[8][1][1],Traps[8][1][2],Traps[8][1][3],beam_zdr_sd(ibin)); // ZDR_SD2D
506 Pij(1,9) = trap(Traps[9][1][0],Traps[9][1][1],Traps[9][1][2],Traps[9][1][3],beam_sqi(ibin)); // SQI
507 Pij(1,10) = trap(Traps[10][1][0],Traps[10][1][1],Traps[10][1][2],Traps[10][1][3], beam_snr(ibin)); // SNR
508 if(radar=="GAT"){
509 Pij(1,11) = std::max(trap(Traps[11][1][0],Traps[11][1][1],Traps[11][1][2],Traps[11][1][3],beam_zvd(ibin)),trap(35.,40.,60.,60.,beam_zvd(ibin)));
510 }else
511 Pij(1,11) = trap(Traps[11][1][0],Traps[11][1][1],Traps[11][1][2],Traps[11][1][3],beam_zvd(ibin));//, trap(70.,70.,100.,100.,beam_zvd(ibin)) ); // DBZH_VD
512 //if(radar=="GAT"&&beam_zvd(ibin)>=70.) Pij(1,11) = 1.0 ;
513// INTERF. Strong
514 Pij(2,0) = trap(Traps[2][2][0],Traps[2][2][1],Traps[2][2][2],Traps[2][2][3],beam_z(ibin), Traps[2][2][4]); // Z
515 Pij(2,3) = trap (Traps[3][2][0],Traps[3][2][1],Traps[3][2][2],Traps[3][2][3], beam_sd(ibin)); // SD_2D
516 Pij(2,4) = trap (Traps[4][2][0],Traps[4][2][1],Traps[4][2][2],Traps[4][2][3], beam_sdray(ibin)); // SD_RAY
517 Pij(2,5) = trap (Traps[5][2][0],Traps[5][2][1],Traps[5][2][2],Traps[5][2][3], beam_sdaz(ibin)); // SD_AZ
518 //if(Num_entries>6)
519 Pij(2,6) = trap(Traps[6][2][0],Traps[6][2][1],Traps[6][2][2],Traps[6][2][3],beam_zdr(ibin),Traps[6][2][4]); // ZDR
520 Pij(2,7) = trap(Traps[7][2][0],Traps[7][2][1],Traps[7][2][2],Traps[7][2][3],beam_rohv(ibin)); // ROHV
521 Pij(2,8) = trap(Traps[8][2][0],Traps[8][2][1],Traps[8][2][2],Traps[8][2][3],beam_zdr_sd(ibin)); // ZDR_SD2D
522 Pij(2,9) = trap(Traps[9][2][0],Traps[9][2][1],Traps[9][2][2],Traps[9][2][3],beam_sqi(ibin),Traps[9][2][4]); // SQI
523 Pij(2,10) = trap(Traps[10][2][0],Traps[10][2][1],Traps[10][2][2],Traps[10][2][3], beam_snr(ibin)); // SNR
524 if(radar=="GAT")
525 Pij(2,11) = std::max(trap(Traps[11][2][0],Traps[11][2][1],Traps[11][2][2],Traps[11][2][3],beam_zvd(ibin)),trap(10.,20.,40.,55.,beam_zvd(ibin)));
526 else
527 Pij(2,11) = trap(Traps[11][1][0],Traps[11][1][1],Traps[11][1][2],Traps[11][1][3],beam_zvd(ibin)); // DBZH_VD
528// INTERF. Med.
529 Pij(3,0) = trap(Traps[2][3][0],Traps[2][3][1],Traps[2][3][2],Traps[2][3][3],beam_z(ibin), Traps[2][3][4]); // Z
530 Pij(3,3) = trap (Traps[3][3][0],Traps[3][3][1],Traps[3][3][2],Traps[3][3][3], beam_sd(ibin)); // SD_2D
531 Pij(3,4) = trap (Traps[4][3][0],Traps[4][3][1],Traps[4][3][2],Traps[4][3][3], beam_sdray(ibin)); // SD_RAY
532 Pij(3,5) = trap (Traps[5][3][0],Traps[5][3][1],Traps[5][3][2],Traps[5][3][3], beam_sdaz(ibin)); // SD_AZ
533 //if(Num_entries>6)
534 Pij(3,6) = trap(Traps[6][3][0],Traps[6][3][1],Traps[6][3][2],Traps[6][3][3],beam_zdr(ibin),Traps[6][3][4]); // ZDR
535 Pij(3,7) = trap(Traps[7][3][0],Traps[7][3][1],Traps[7][3][2],Traps[7][3][3],beam_rohv(ibin)); // ROHV
536 Pij(3,8) = trap(Traps[8][3][0],Traps[8][3][1],Traps[8][3][2],Traps[8][3][3],beam_zdr_sd(ibin)); // ZDR_SD2D
537 Pij(3,9) = trap(Traps[9][3][0],Traps[9][3][1],Traps[9][3][2],Traps[9][3][3],beam_sqi(ibin),Traps[9][3][4]); // SQI
538 Pij(3,10) = trap(Traps[10][3][0],Traps[10][3][1],Traps[10][3][2],Traps[10][3][3], beam_snr(ibin)); // SNR
539 Pij(3,11) = trap(Traps[11][3][0],Traps[11][3][1],Traps[11][3][2],Traps[11][3][3],beam_zvd(ibin)); // DBZH_VD
540// INTERF. Weak
541 Pij(4,0) = trap(Traps[2][4][0],Traps[2][4][1],Traps[2][4][2],Traps[2][4][3],beam_z(ibin), Traps[2][4][4]); // Z
542 Pij(4,3) = trap (Traps[3][4][0],Traps[3][4][1],Traps[3][4][2],Traps[3][4][3], beam_sd(ibin)); // SD_2D
543 Pij(4,4) = trap (Traps[4][4][0],Traps[4][4][1],Traps[4][4][2],Traps[4][4][3], beam_sdray(ibin)); // SD_RAY
544 Pij(4,5) = trap (Traps[5][4][0],Traps[5][4][1],Traps[5][4][2],Traps[5][4][3], beam_sdaz(ibin)); // SD_AZ
545 //if(Num_entries>6)
546 Pij(4,6) = trap(Traps[6][4][0],Traps[6][4][1],Traps[6][4][2],Traps[6][4][3],beam_zdr(ibin),Traps[6][4][4]); // ZDR
547 Pij(4,7) = trap(Traps[7][4][0],Traps[7][4][1],Traps[7][4][2],Traps[7][4][3],beam_rohv(ibin),Traps[7][4][4]); // ROHV
548 Pij(4,8) = trap(Traps[8][4][0],Traps[8][4][1],Traps[8][4][2],Traps[8][4][3],beam_zdr_sd(ibin)); // ZDR_SD2D
549 Pij(4,9) = trap(Traps[9][4][0],Traps[9][4][1],Traps[9][4][2],Traps[9][4][3],beam_sqi(ibin),Traps[9][4][4]); // SQI
550 Pij(4,10) = trap(Traps[10][4][0],Traps[10][4][1],Traps[10][4][2],Traps[10][4][3], beam_snr(ibin)); // SNR
551 Pij(4,11) = trap(Traps[11][4][0],Traps[11][4][1],Traps[11][4][2],Traps[11][4][3],beam_zvd(ibin)); // DBZH_VD
552
553//---- fine calcolo probabilità
554// Calcolo classe appartenenza
555 Class_WP = ((Wij.array()*Pij.array()).matrix()*VectorXd::Ones(Num_entries)).array()/(Wij*VectorXd::Ones(Num_entries)).array();
556 unsigned i,ID;
557 Class_WP.maxCoeff(&i);
558 ID=i;
559 if (Class_WP(i) < 0.1 ) ID=5;
560 res[ibin]=ID;
561 double diff = Class_WP[i]-Class_WP[0];
562 diff = diff*100;
563 //printf("ID %d \n",ID);
564 //counter[ID]++;
565 if(stamp){
566 printf("bin %4d ",ibin);
567 printf("ID %d -- diff=%8.3f\n",ID, diff);
568 //printf("prevID %")
569 //lo aggiungi come array
570 for(unsigned c=0;c<5;c++) printf("%5.3f ",Class_WP(c));
571 printf(" ID %d \n",ID);
572 printf(" %8.3f %8.3f %8.3f %8.3f %8.3f %8.3f %8.3f %8.3f %8.3f %8.3f %8.3f %8.3f \n",beam_z(ibin), beam_v(ibin), beam_w(ibin), beam_sd(ibin), beam_sdray(ibin), beam_sdaz(ibin), beam_zdr(ibin),
573 beam_rohv(ibin), beam_zdr_sd(ibin), beam_sqi(ibin), beam_snr(ibin), beam_zvd(ibin));
574 for (unsigned c=0;c<5; c++){
575 for (unsigned d=0; d<12; d++) printf(" %8.3f", Pij(c,d));
576 printf("\n");
577 }
578 }
579 //if((force_meteo)&&(diff<=50.0)&revID==0)&&(ID==4)) res[ibin] = 0;
580 if(radar=="GAT"){
581 //if(ibin<=120){
582 //if((force_meteo)&&(diff<=8.0)&&(prevID==0)&&(ID!=1)&&(ID!=4)) res[ibin] = 0;
583 //}else{
584 if((force_meteo)&&(diff<=8.0)&&(prevID==0)&&(ID!=1)) res[ibin] = 0;
585 //}
586
587 //if((force_meteo)&&(diff<=10.0)&&(prevID==0)) res[ibin] = 0;
588 //if((force_meteo)&&(ID==0)&&(prevID==1)&&(ibin<=120)) res[ibin]= 1;
589 //if((force_meteo)&&(diff<=12.0)&&(prevID==0)&&(ID!=1)&&(ibin>120)) res[ibin] = 0;
590 }
591 else if(radar=="SPC"){
592 if((force_meteo)&&(diff<=10.)&&(prevID==0)) res[ibin] = 0;
593
594 if((force_meteo)&&(diff<=20.0)&&(prevID==0)&&(ID==4)) res[ibin] = 0;
595 }else cout<<"radar non specificato"<<endl;
596
597 diffs[ibin]=diff;
598 counter[ID]++;
599 prevID = res[ibin];
600 }
601
602 return {res, diffs};
603
604 }
605
606// CC: fuzzy logic senza zdr
607 std::vector<unsigned char> Cleaner::eval_classID_beam(const Eigen::VectorXd& beam_z, const Eigen::VectorXd& beam_w, const Eigen::VectorXd& beam_v, const Eigen::VectorXd& beam_sd, const Eigen::VectorXd& beam_sdray, const Eigen::VectorXd& beam_sdaz, int iray, const string radar, double v_ny, const char* fuzzy_path) const
608{
609
610 const unsigned beam_size = beam_z.rows();
611 vector<unsigned char> res(beam_size, 0);
612 int Num_entries=0;
613 int Num_echoes = 5;
614 int Ntraps = 5; // 5 argomenti da passare a Trap : x1,x2,x3,x4,x5
615 //char *cwd = get_current_dir_name();
616 //string f_dir = cwd;
617 string f_dir = fuzzy_path;
618
619 //leggo matrice dei pesi
620 string fin_w = f_dir+"/matrix-"+radar+"-nozdr.txt";
621 vector<string> w_vector;
622 w_vector = read_matrix_from_txt(fin_w);
623 Num_entries = w_vector.size()/Num_echoes;
624
625 Matrix2D<double> Wij(Num_echoes,Num_entries);
626 for(int i=0;i<Num_echoes;i++){ //itero colonna
627 for(int j=0;j<Num_entries;j++){ //itero rriga
628 Wij(i,j) = stod( w_vector[i*Num_entries+j]);
629 }
630 }
631
632 //leggo matrice dei traps
633 string fin_t = f_dir+"/Trap-"+radar+"-nozdr.txt";
634 vector<string> t_vector;
635 t_vector = read_matrix_from_txt(fin_t);
636 double Traps[Num_entries][Num_echoes][Ntraps];
637 for(int i=0;i<Num_entries;i++){
638 for(int j=0;j<Num_echoes;j++){
639 for(int k=0;k<Ntraps;k++){
640 Traps[i][j][k] = stod( t_vector[i*(Num_echoes*Ntraps)+j*Num_echoes+k]);
641 //cout<<" Traps["<<i<<","<<j<<","<<k<<"]="<<Traps[i][j][k];
642 }
643 //cout<<endl;
644 }
645 //cout<<endl;
646 }
647 //------------------------------------------------------------------------
648
649 vector<unsigned> counter (Num_entries,0) ; // non sono sicura di cosa delle dimensioni di questo counter
650 for (unsigned ibin = 0; ibin < beam_size; ++ibin)
651 //for (unsigned ibin = 87*4; ibin < 87*4+12; ++ibin)
652 {
653 //printf(" %4d %4d %6.2f %6.2f %10.6f %10.6f %10.6f %10.6f %10.6f",iray,ibin , beam_z(ibin),beam_v(ibin),beam_w(ibin),beam_zdr(ibin),beam_sd(ibin),beam_sdray(ibin),beam_sdaz(ibin));
654 Matrix2D<double> Pij(Num_echoes,Num_entries);
655 Pij = Pij * 0.;
656 //cout<<Pij<<endl;
657 ArrayXd Class_WP(Num_echoes); // la dimensione di Class_WP deve essere Num_echoes, perchè alla fine ricavo un vettore con 6 valori, uno per ogni echo, passando per prodotto pesi prob
658 Class_WP.setZero();
659 if (beam_z(ibin) == Z_missing) {
660 unsigned ID=0;
661 res[ibin]=ID;
662 counter[ID]++;
663
664 continue;
665 }
666
667//eseguo un test unico per VRAD e WRAD e assegno tutte le prob assieme
668 if (beam_v(ibin) != bin_wind_magic_number ) { // VRAD
669 Pij(0,1)=trap(v_ny*Traps[0][0][0], v_ny*Traps[0][0][1], v_ny*Traps[0][0][2], v_ny*Traps[0][0][3], beam_v(ibin)); // METEO
670 double prob_v = trap (v_ny, v_ny, v_ny, v_ny, beam_v(ibin),v_ny);
671 //cout<<"prob_v computed: "<<prob_v<<endl;
672 //for(int e=1;e<Num_echoes;e++){ Pij(e,1) = prob_v };
673 Pij(1,1)=trap(Traps[0][1][0],Traps[0][1][1],Traps[0][1][2],Traps[0][1][3], beam_v(ibin)); // CLUTTER
674 Pij(2,1)=prob_v; // INTERF. Strong
675 Pij(3,1)=prob_v; // INTERF. Med.
676 Pij(4,1)=prob_v; // INTERF. Weak
677 } else {
678 Pij(0,1)=0.3; // METEO
679 Pij(1,1)=1.; // CLUTTER
680 Pij(2,1)=1.; // INTERF. Strong
681 Pij(3,1)=1.; // INTERF. Med.
682 Pij(4,1)=1.; // INTERF. Weak
683 }
684
685 // WRAD
686 Pij(0,2)= trap(Traps[1][0][0],Traps[1][0][1],Traps[1][0][2],Traps[1][0][3],beam_w(ibin)); // METEO
687 Pij(1,2)= trap(Traps[1][1][0],Traps[1][1][1],Traps[1][1][2],Traps[1][1][3],beam_w(ibin)); // CLUTTER
688 double prob_w = trap(Traps[1][2][0],Traps[1][2][1],Traps[1][2][2],Traps[1][2][3],beam_w(ibin));
689 Pij(2,2) = prob_w; // INTERF MULTIPLE
690 Pij(3,2) = prob_w; // INTERF. Med.
691 Pij(4,2) = prob_w; // INTERF. Weak
692
693// METEO
694 Pij(0,0) = trap(Traps[2][0][0],Traps[2][0][1],Traps[2][0][2],Traps[2][0][3],beam_z(ibin), Traps[2][0][4]); // Z
695 Pij(0,3) = trap (Traps[3][0][0],Traps[3][0][1],Traps[3][0][2],Traps[3][0][3], beam_sd(ibin)); // SD_2D
696 Pij(0,4) = trap (Traps[4][0][0],Traps[4][0][1],Traps[4][0][2],Traps[4][0][3], beam_sdray(ibin)); // SD_RAY
697 Pij(0,5) = trap (Traps[5][0][0],Traps[5][0][1],Traps[5][0][2],Traps[5][0][3], beam_sdaz(ibin)); // SD_AZ
698
699// CLUTTER
700 Pij(1,0) = trap (Traps[2][1][0],Traps[2][1][1],Traps[2][1][2],Traps[2][1][3], beam_z(ibin), Traps[2][1][4]); // Z
701 Pij(1,3) = trap (Traps[3][1][0],Traps[3][1][1],Traps[3][1][2],Traps[3][1][3], beam_sd(ibin)); // SD_2D
702 Pij(1,4) = trap (Traps[4][1][0],Traps[4][1][1],Traps[4][1][2],Traps[4][1][3], beam_sdray(ibin)); // SD_RAY
703 Pij(1,5) = trap (Traps[5][1][0],Traps[5][1][1],Traps[5][1][2],Traps[5][1][3], beam_sdaz(ibin)); // SD_AZ
704// INTERF. Strong
705 Pij(2,0) = trap(Traps[2][2][0],Traps[2][2][1],Traps[2][2][2],Traps[2][2][3],beam_z(ibin), Traps[2][2][4]); // Z
706 Pij(2,3) = trap (Traps[3][2][0],Traps[3][2][1],Traps[3][2][2],Traps[3][2][3], beam_sd(ibin)); // SD_2D
707 Pij(2,4) = trap (Traps[4][2][0],Traps[4][2][1],Traps[4][2][2],Traps[4][2][3], beam_sdray(ibin)); // SD_RAY
708 Pij(2,5) = trap (Traps[5][2][0],Traps[5][2][1],Traps[5][2][2],Traps[5][2][3], beam_sdaz(ibin)); // SD_AZ
709// INTERF. Med.
710 Pij(3,0) = trap(Traps[2][3][0],Traps[2][3][1],Traps[2][3][2],Traps[2][3][3],beam_z(ibin), Traps[2][3][4]); // Z
711 Pij(3,3) = trap (Traps[3][3][0],Traps[3][3][1],Traps[3][3][2],Traps[3][3][3], beam_sd(ibin)); // SD_2D
712 Pij(3,4) = trap (Traps[4][3][0],Traps[4][3][1],Traps[4][3][2],Traps[4][3][3], beam_sdray(ibin)); // SD_RAY
713 Pij(3,5) = trap (Traps[5][3][0],Traps[5][3][1],Traps[5][3][2],Traps[5][3][3], beam_sdaz(ibin)); // SD_AZ
714// INTERF. Weak
715 Pij(4,0) = trap(Traps[2][4][0],Traps[2][4][1],Traps[2][4][2],Traps[2][4][3],beam_z(ibin), Traps[2][4][4]); // Z
716 Pij(4,3) = trap (Traps[3][4][0],Traps[3][4][1],Traps[3][4][2],Traps[3][4][3], beam_sd(ibin)); // SD_2D
717 Pij(4,4) = trap (Traps[4][4][0],Traps[4][4][1],Traps[4][4][2],Traps[4][4][3], beam_sdray(ibin)); // SD_RAY
718 Pij(4,5) = trap (Traps[5][4][0],Traps[5][4][1],Traps[5][4][2],Traps[5][4][3], beam_sdaz(ibin)); // SD_AZ
719
720//---- fine calcolo probabilità
721// Calcolo classe appartenenza
722
723 Class_WP = ((Wij.array()*Pij.array()).matrix()*VectorXd::Ones(Num_entries)).array()/(Wij*VectorXd::Ones(Num_entries)).array();
724 unsigned i,ID;
725 Class_WP.maxCoeff(&i);
726 ID=i;
727 if (Class_WP(i) < 0.1 ) ID=5;
728 res[ibin]=ID;
729 //printf("ID %d \n",ID);
730 counter[ID]++;
731
732 }
733
734 return res;
735
736 }
737
738// Senza ZDR - Basato su test
739std::vector<unsigned char> Cleaner::eval_clean_beam(const Eigen::VectorXd& beam_z, const Eigen::VectorXd& beam_w, const Eigen::VectorXd& beam_v, const Eigen::VectorXd& beam_sd, const Eigen::VectorXd& beam_sdray, const Eigen::VectorXd& beam_sdaz, int iray) const
740{
741
742 const unsigned beam_size = beam_z.rows();
743 vector<unsigned char> res(beam_size, 0);
744 unsigned counter = 0;
745 unsigned countIntWeak = 0;
746 unsigned countIntStrong = 0;
747 unsigned countIntMed = 0;
748 unsigned countClutter = 0;
749 unsigned countNoise = 0;
750
751 for (unsigned ibin = 0; ibin < beam_size; ++ibin)
752 {
753 //printf(" %4d %4d %6.2f %6.2f %10.6f %10.6f %10.6f %10.6f",iray,ibin , beam_z(ibin),beam_v(ibin),beam_w(ibin),beam_sd(ibin),beam_sdray(ibin),beam_sdaz(ibin));
754 if ((beam_w(ibin) > W_threshold && beam_v(ibin) != bin_wind_magic_number && beam_sd(ibin) >= 1. &&
755 beam_sd(ibin) <= 10. ) || beam_z(ibin) == Z_missing) {
756 // this should be a meteorological echo
757 ;
758 }else {
759 res [ibin]=1; // Not meteo but unclassified
760 counter++;
761 if (beam_w(ibin) == W_threshold && beam_v(ibin) == bin_wind_magic_number)
762 {
763 if (beam_sd(ibin) >= 5. && beam_sdray(ibin) >= 2 && beam_sdaz(ibin) > 4.){
764 // this should be clutter
765 res[ibin] = 2;
766 countClutter ++;
767 }
768 if (beam_sd(ibin) >= 1. && beam_sd(ibin) <= 5. && beam_sdray(ibin) < 2. && beam_sdaz(ibin) > 4.){
769 // this should be a strong Interference
770 res[ibin] = 3;
771 countIntStrong ++;
772 }
773 if (beam_sd(ibin) >= 3. && beam_sd(ibin) <= 7. && beam_sdray(ibin) < 2. && beam_sdaz(ibin) > 4.&& beam_sdaz(ibin) <7 ){
774 // this should be a medium Interference
775 res[ibin] = 4;
776 countIntMed ++;
777 }
778 //if (beam_sd(ibin) >= 5. && beam_sd(ibin) <= 3. && beam_sdray(ibin) < 3. && beam_sdaz(ibin) < 3.){
779 if (beam_sd(ibin) >= 5. && beam_sdray(ibin) < 3. && beam_sdaz(ibin) < 3.){
780 // this should be a weakInterference
781 res[ibin] = 5;
782 countIntWeak ++;
783 }
784 if (beam_sd(ibin) >= 10. && beam_sdray(ibin) < 2. && beam_sdaz(ibin) < 5. && beam_z(ibin) <10.){
785 // this should be a noise
786 res[ibin] = 6;
787 countNoise ++;
788 }
789 }
790 } // ELSE this should be not a meteo echo
791 //printf("%2d %4d %4d %4d %4d %4d\n",res[ibin],countClutter,countIntStrong, countIntMed, countIntWeak, countNoise);
792 }
793
794 return res;
795
796}
797//----------------------------------------------------------------------------------
798
799
800
801void Cleaner::evaluateCleanID(PolarScan<double>& scan_z, PolarScan<double>& scan_w, PolarScan<double>& scan_v,PolarScan<unsigned char>& scan_cleanID, unsigned iel)
802{
803 return evaluateCleanID(scan_z, scan_w, scan_v, scan_cleanID, scan_v.undetect,iel);
804 //return evaluateClassID(scan_z, scan_w, scan_v, scan_cleanID, scan_v.undetect, radar, iel);
805}
806
807void Cleaner::evaluateCleanID(PolarScan<double>& scan_z, PolarScan<double>& scan_w, PolarScan<double>& scan_v,PolarScan<unsigned char>& scan_cleanID, double bin_wind_magic_number, unsigned iel)
808{
809
810 if (scan_z.beam_count != scan_w.beam_count)
811 throw std::runtime_error("scan_z beam_count no equal to scan_w beam_count");
812 if (scan_z.beam_size != scan_w.beam_size)
813 throw std::runtime_error("scan_z beam_size no equal to scan_w beam_size");
814
815 if (scan_z.beam_count != scan_v.beam_count)
816 throw std::runtime_error("scan_z beam_count no equal to scan_v beam_count");
817 if (scan_z.beam_size != scan_v.beam_size)
818 throw std::runtime_error("scan_z beam_size no equal to scan_v beam_size");
819
820 Cleaner cleaner(scan_z.undetect, scan_w.undetect, scan_v.nodata, bin_wind_magic_number);
821
822 const unsigned beam_count = scan_z.beam_count;
823 const unsigned beam_size = scan_z.beam_size;
824
825 // fprintf(stderr, "NEWCLEANER zmis %f, wthr %f, vmis %f, mn %f\n",
826 // cleaner.Z_missing, cleaner.W_threshold, cleaner.V_missing, cleaner.bin_wind_magic_number);
827
828// radarelab::volume::Scans<double> Z_S, SD2D;
829// Z_S.push_back(scan_z);
830// radarelab::volume::textureSD( Z_S,SD2D, 1000. , 3,false);
831
832 radarelab::volume::Scans<double> Z_S, SD2D,SD_Ray,SD_Az;
833 Z_S.push_back(scan_z);
834 radarelab::volume::textureSD( Z_S,SD2D, 1000. , 3,false);
835
836 radarelab::volume::textureSD( Z_S,SD_Ray, scan_z.cell_size*9 , 360./scan_z.beam_count,true);
837 radarelab::volume::textureSD( Z_S,SD_Az, scan_z.cell_size , 5*360./scan_z.beam_count,true);
838
839 for (unsigned i = 0; i < beam_count; ++i)
840 {
841 //vector<unsigned char> corrected = cleaner.eval_clean_beam(scan_z.row(i), scan_w.row(i), scan_v.row(i),i);
842 vector<unsigned char> corrected = cleaner.eval_clean_beam(scan_z.row(i), scan_w.row(i), scan_v.row(i),SD2D[0].row(i), SD_Ray[0].row(i), SD_Az[0].row(i), i);
843 //vector<bool> corrected = cleaner.clean_beam(scan_z.row(i), scan_w.row(i), scan_v.row(i),SD2D[0].row(i), SDZDR2D[0].row(i), scan_z, scan_w, scan_v, SD2D[0],i);
844 for (unsigned ib = 0; ib < beam_size; ++ib)
845 scan_cleanID(i,ib)=corrected[ib];
846 }
847}
848
849void Cleaner::evaluateClassID(PolarScan<double>& scan_z, PolarScan<double>& scan_w, PolarScan<double>& scan_v, PolarScan<double>& scan_zdr, PolarScan<double>& scan_rohv, PolarScan<double>& scan_sqi, PolarScan<double>& scan_snr, PolarScan<double>& scan_zvd, PolarScan<unsigned char>& scan_cleanID, PolarScan<double>& scan_DiffProb, double bin_wind_magic_number,const string radar, const char* fuzzy_path, unsigned iel, bool force_meteo)
850{
851
852 if (scan_z.beam_count != scan_w.beam_count)
853 throw std::runtime_error("scan_z beam_count no equal to scan_w beam_count");
854 if (scan_z.beam_size != scan_w.beam_size)
855 throw std::runtime_error("scan_z beam_size no equal to scan_w beam_size");
856
857 if (scan_z.beam_count != scan_v.beam_count)
858 throw std::runtime_error("scan_z beam_count no equal to scan_v beam_count");
859 if (scan_z.beam_size != scan_v.beam_size)
860 throw std::runtime_error("scan_z beam_size no equal to scan_v beam_size");
861
862 Cleaner cleaner(scan_z.undetect, scan_w.undetect, scan_v.nodata, bin_wind_magic_number);
863
864 const unsigned beam_count = scan_z.beam_count;
865 const unsigned beam_size = scan_z.beam_size;
866
867 //cout<<"VRAD OFFSET="<<scan_v.offset<<" , GAIN="<<scan_v.gain<<endl;
868
869// compute texture volumes
870 radarelab::volume::Scans<double> Z_S, SD2D, SD_Ray, SD_Az, ZDR_S, ZDR_SD2D;
871 Z_S.push_back(scan_z);
872 ZDR_S.push_back(scan_zdr);
873 radarelab::volume::textureSD( Z_S,SD2D, 1000. , 3,false);
874 radarelab::volume::textureSD( Z_S,SD_Ray, scan_z.cell_size*21 , 360./scan_z.beam_count,true);
875 radarelab::volume::textureSD( Z_S,SD_Az, scan_z.cell_size , 5*360./scan_z.beam_count,true);
876
877 radarelab::volume::textureSD( ZDR_S,ZDR_SD2D, 1000. , 3,false);
878
879 //bool stamp=false;
880 for (unsigned i = 0; i <beam_count ; ++i)
881 {
882 bool stamp=false;
883 //if(i==100) stamp=true;//311
884 //vector<unsigned char> corrected = cleaner.eval_classID_beam(scan_z.row(i), scan_w.row(i), scan_v.row(i), SD2D[0].row(i), scan_zdr.row(i), scan_rohv.row(i), scan_sqi.row(i), scan_snr.row(i), scan_zvd.row(i), SD_Ray[0].row(i), SD_Az[0].row(i), ZDR_SD2D[0].row(i), i, radar, scan_v.offset, stamp);
885 auto [corrected, diff_prob] = cleaner.eval_classID_beam(scan_z.row(i), scan_w.row(i), scan_v.row(i), SD2D[0].row(i), scan_zdr.row(i), scan_rohv.row(i), scan_sqi.row(i), scan_snr.row(i), scan_zvd.row(i), SD_Ray[0].row(i), SD_Az[0].row(i), ZDR_SD2D[0].row(i), i, radar, scan_v.offset, fuzzy_path, stamp, force_meteo);
886 //vector<unsigned char> corrected = cleaner.eval_classID_beam(scan_z.row(i), scan_w.row(i), scan_v.row(i), SD2D[0].row(i), scan_zdr.row(i), scan_rohv.row(i), scan_sqi.row(i), scan_snr.row(i), scan_zvd.row(i), SD_Ray[0].row(i), SD_Az[0].row(i), ZDR_SD2D[0].row(i), i, radar, scan_v.offset, stamp, force_meteo);
887 for (unsigned ib = 0; ib < beam_size; ++ib){
888 scan_cleanID(i,ib)=corrected[ib];
889 scan_DiffProb(i,ib)=diff_prob[ib];
890 }
891 }
892 //cout<<"Scan_DiffProb(183,848)="<<scan_DiffProb(183,848)<<endl;
893 //cout << typeid(scan_DiffProb(0,100)).name() << endl;
894}
895
896// senza zdr
897void Cleaner::evaluateClassID(PolarScan<double>& scan_z, PolarScan<double>& scan_w, PolarScan<double>& scan_v, PolarScan<unsigned char>& scan_cleanID, double bin_wind_magic_number, const string radar, const char* fuzzy_path, unsigned iel)
898{
899
900 if (scan_z.beam_count != scan_w.beam_count)
901 throw std::runtime_error("scan_z beam_count no equal to scan_w beam_count");
902 if (scan_z.beam_size != scan_w.beam_size)
903 throw std::runtime_error("scan_z beam_size no equal to scan_w beam_size");
904
905 if (scan_z.beam_count != scan_v.beam_count)
906 throw std::runtime_error("scan_z beam_count no equal to scan_v beam_count");
907 if (scan_z.beam_size != scan_v.beam_size)
908 throw std::runtime_error("scan_z beam_size no equal to scan_v beam_size");
909
910 Cleaner cleaner(scan_z.undetect, scan_w.undetect, scan_v.nodata, bin_wind_magic_number);
911
912 const unsigned beam_count = scan_z.beam_count;
913 const unsigned beam_size = scan_z.beam_size;
914
915// compute texture volumes
916 radarelab::volume::Scans<double> Z_S, SD2D,SD_Ray,SD_Az;
917 Z_S.push_back(scan_z);
918 radarelab::volume::textureSD( Z_S,SD2D, 1000. , 3,false);
919 radarelab::volume::textureSD( Z_S,SD_Ray, scan_z.cell_size*21 , 360./scan_z.beam_count,true);
920 radarelab::volume::textureSD( Z_S,SD_Az, scan_z.cell_size , 5*360./scan_z.beam_count,true);
921
922
923 for (unsigned i = 0; i <beam_count ; ++i)
924 {
925 vector<unsigned char> corrected = cleaner.eval_classID_beam(scan_z.row(i), scan_w.row(i), scan_v.row(i), SD2D[0].row(i), SD_Ray[0].row(i), SD_Az[0].row(i), i, radar, scan_v.offset, fuzzy_path);
926 for (unsigned ib = 0; ib < beam_size; ++ib)
927 scan_cleanID(i,ib)=corrected[ib];
928 }
929}
930
931void Cleaner::clean(PolarScan<double>& scan_z, PolarScan<double>& scan_w, PolarScan<double>& scan_v, unsigned iel , bool set_undetect)
932{
933 return clean(scan_z, scan_w, scan_v, scan_v.undetect,iel);
934}
935void Cleaner::clean(PolarScan<double>& scan_z, PolarScan<double>& scan_w, PolarScan<double>& scan_v, double bin_wind_magic_number, unsigned iel, bool set_undetect)
936{
937
938 if (scan_z.beam_count != scan_w.beam_count)
939 throw std::runtime_error("scan_z beam_count no equal to scan_w beam_count");
940 if (scan_z.beam_size != scan_w.beam_size)
941 throw std::runtime_error("scan_z beam_size no equal to scan_w beam_size");
942
943 if (scan_z.beam_count != scan_v.beam_count)
944 throw std::runtime_error("scan_z beam_count no equal to scan_v beam_count");
945 if (scan_z.beam_size != scan_v.beam_size)
946 throw std::runtime_error("scan_z beam_size no equal to scan_v beam_size");
947
948 Cleaner cleaner(scan_z.undetect, scan_w.undetect, scan_v.nodata, bin_wind_magic_number);
949
950 const unsigned beam_count = scan_z.beam_count;
951 const unsigned beam_size = scan_z.beam_size;
952
953 // fprintf(stderr, "NEWCLEANER zmis %f, wthr %f, vmis %f, mn %f\n",
954 // cleaner.Z_missing, cleaner.W_threshold, cleaner.V_missing, cleaner.bin_wind_magic_number);
955
956 radarelab::volume::Scans<double> Z_S, SD2D,SD_Ray,SD_Az;
957 Z_S.push_back(scan_z);
958 radarelab::volume::textureSD( Z_S,SD2D, 1000. , 3,false);
959 radarelab::volume::textureSD( Z_S,SD_Ray, scan_z.cell_size*9 , 360./scan_z.beam_count,false);
960 radarelab::volume::textureSD( Z_S,SD_Az, scan_z.cell_size , 5*360./scan_z.beam_count,false);
961
962//radarelab::gdal_init_once();
963//printf("scrivo Z ");
964//Matrix2D <double>img;
965//img = (scan_z.array() - scan_z.offset )/ scan_z.gain /256 ;
966//Matrix2D <unsigned char>img_tmp, z_clean;
967//std::string ext;
968//char pippo[200];
969//sprintf(pippo, "_%02d.png",iel);
970//ext=pippo;
971
972//img_tmp=img.cast<unsigned char>();
973//z_clean=img_tmp;
974//radarelab::write_image(img_tmp,"/ponte/rad_svn/proc_operative/test_arch/rev_actual/radar/immagini/Cleaner/PPI_Z"+ext, "PNG");
975
976//printf("V ");
977//img = (scan_v.array()-scan_v.offset)/scan_v.gain/256 ;
978//img_tmp=img.cast<unsigned char>();
979//radarelab::write_image(img_tmp,"/ponte/rad_svn/proc_operative/test_arch/rev_actual/radar/immagini/Cleaner/PPI_V"+ext,"PNG");
980//printf("W ");
981//img = (scan_w.array()-scan_w.offset)/scan_w.gain/256 ;
982//img_tmp=img.cast<unsigned char>();
983//radarelab::write_image(img_tmp,"/ponte/rad_svn/proc_operative/test_arch/rev_actual/radar/immagini/Cleaner/PPI_W"+ext,"PNG");
984//printf("SD2d ");
985//img = (SD2D[0].array()-SD2D[0].offset)/SD2D[0].gain/256 ;
986//img_tmp=img.cast<unsigned char>();
987//radarelab::write_image(img_tmp,"/ponte/rad_svn/proc_operative/test_arch/rev_actual/radar/immagini/Cleaner/PPI_SD2d"+ext,"PNG");
988//printf("\n");
989
990 double new_val=cleaner.Z_missing;
991 if (set_undetect) new_val=scan_z.nodata;
992
993 for (unsigned i = 0; i < beam_count; ++i)
994 {
995 // Compute which elements need to be cleaned
996 // vector<bool> corrected = cleaner.clean_beam(scan_z.row(i), scan_w.row(i), scan_v.row(i),SD2D[0].row(i), scan_z, scan_w, scan_v, SD2D[0],i);
997 vector<bool> corrected = cleaner.clean_beam(scan_z.row(i), scan_w.row(i), scan_v.row(i),i);
998
999 for (unsigned ib = 0; ib < beam_size; ++ib)
1000 if (corrected[ib])
1001 {
1002 //scan_z(i, ib) = cleaner.Z_missing;
1003 scan_z(i, ib) = new_val;
1004 // scan_w(i, ib) = cleaner.W_threshold;
1005 // scan_v(i, ib) = cleaner.V_missing;
1006 }
1007// img_tmp(i,ib)=255;
1008// z_clean(i,ib)=0;
1009// } else img_tmp(i,ib)= 0 ;
1010
1011 }
1012//radarelab::write_image(img_tmp,"/ponte/rad_svn/proc_operative/test_arch/rev_actual/radar/immagini/Cleaner/PPI_clean"+ext,"PNG");
1013//radarelab::write_image(z_clean,"/ponte/rad_svn/proc_operative/test_arch/rev_actual/radar/immagini/Cleaner/PPI_Zclean"+ext,"PNG");
1014}
1015
1016
1017
1018void Cleaner::clean(PolarScan<double>& scan_z, PolarScan<double>& scan_w, PolarScan<double>& scan_v, PolarScan<double>& scan_zdr, unsigned iel, bool set_undetect )
1019{
1020 return clean(scan_z, scan_w, scan_v, scan_zdr, scan_v.undetect,iel,set_undetect);
1021}
1022void Cleaner::clean(PolarScan<double>& scan_z, PolarScan<double>& scan_w, PolarScan<double>& scan_v, PolarScan<double>& scan_zdr, double bin_wind_magic_number, unsigned iel, bool set_undetect )
1023{
1024 cout<<"Chiamato cleaner "<<set_undetect<<endl;
1025 if (scan_z.beam_count != scan_w.beam_count)
1026 throw std::runtime_error("scan_z beam_count no equal to scan_w beam_count");
1027 if (scan_z.beam_size != scan_w.beam_size)
1028 throw std::runtime_error("scan_z beam_size no equal to scan_w beam_size");
1029
1030 if (scan_z.beam_count != scan_v.beam_count)
1031 throw std::runtime_error("scan_z beam_count no equal to scan_v beam_count");
1032 if (scan_z.beam_size != scan_v.beam_size)
1033 throw std::runtime_error("scan_z beam_size no equal to scan_v beam_size");
1034 double z_val=scan_z.nodata;
1035 if(set_undetect) z_val=scan_z.undetect;
1036//Cleaner cleaner(scan_z.undetect, scan_w.undetect, scan_v.nodata, bin_wind_magic_number);
1037 Cleaner cleaner(z_val, scan_w.undetect, scan_v.nodata, bin_wind_magic_number);
1038
1039 const unsigned beam_count = scan_z.beam_count;
1040 const unsigned beam_size = scan_z.beam_size;
1041
1042 // fprintf(stderr, "NEWCLEANER zmis %f, wthr %f, vmis %f, mn %f\n",
1043 // cleaner.Z_missing, cleaner.W_threshold, cleaner.V_missing, cleaner.bin_wind_magic_number);
1044
1046 Z_S.push_back(scan_z);
1047 radarelab::volume::textureSD( Z_S,SD2D, 1000. , 3,false);
1048 radarelab::volume::Scans<double> ZDR_S, SDZDR2D;
1049 ZDR_S.push_back(scan_zdr);
1050 radarelab::volume::textureSD( ZDR_S,SDZDR2D, 1000. , 3,false);
1051
1052//----------------------------------------------------------------------------------
1053//----------------------------------------------------------------------------------
1054//----------------------------------------------------------------------------------
1055// Mettere a true per fare grafica per debug o false per non fare grafica
1056//
1057// RICORDARSI DI TOGLIERE/METTERE COMMENTI DOPO CLEAN_BEAM
1058// -------------------------------------------------------------------
1059Matrix2D <unsigned char>img_tmp, z_clean;
1060Matrix2D <double>img;
1061std::string ext;
1062char pippo[200];
1063if (false){
1065
1066 printf("scrivo Z ");
1067 img = (scan_z.array() - scan_z.offset )/ scan_z.gain /256 ;
1068 sprintf(pippo, "_%02d.png",iel);
1069 ext=pippo;
1070 img_tmp=img.cast<unsigned char>();
1071 z_clean=img_tmp;
1072 radarelab::write_image(img_tmp,"/ponte/rad_svn/proc_operative/test_arch/rev_actual/radar/immagini/Cleaner/PPI_Z"+ext, "PNG");
1073
1074// printf("V ");
1075// img = (scan_v.array()-scan_v.offset)/scan_v.gain/256 ;
1076// img_tmp=img.cast<unsigned char>();
1077// radarelab::write_image(img_tmp,"/ponte/rad_svn/proc_operative/test_arch/rev_actual/radar/immagini/Cleaner/PPI_V"+ext,"PNG");
1078// printf("W ");
1079// img = (scan_w.array()-scan_w.offset)/scan_w.gain/256 ;
1080// img_tmp=img.cast<unsigned char>();
1081// radarelab::write_image(img_tmp,"/ponte/rad_svn/proc_operative/test_arch/rev_actual/radar/immagini/Cleaner/PPI_W"+ext,"PNG");
1082 printf("SD2d ");
1083 img = (SD2D[0].array()-SD2D[0].offset)/SD2D[0].gain/256 ;
1084 img_tmp=img.cast<unsigned char>();
1085 radarelab::write_image(img_tmp,"/ponte/rad_svn/proc_operative/test_arch/rev_actual/radar/immagini/Cleaner/PPI_SD2d"+ext,"PNG");
1086 printf("SDZDR2d ");
1087 img = (SDZDR2D[0].array()-SDZDR2D[0].offset)/SDZDR2D[0].gain/256 ;
1088 img_tmp=img.cast<unsigned char>();
1089 radarelab::write_image(img_tmp,"/ponte/rad_svn/proc_operative/test_arch/rev_actual/radar/immagini/Cleaner/PPI_SDZDR2d"+ext,"PNG");
1090 printf("\n");
1091}
1092 double new_val=cleaner.Z_missing;
1093 if (set_undetect) new_val=scan_z.undetect;
1094//printf("valore new_val ---> %f\n",new_val);
1095 for (unsigned i = 0; i < beam_count; ++i)
1096 {
1097 // Compute which elements need to be cleaned
1098 vector<bool> corrected = cleaner.clean_beam(scan_z.row(i), scan_w.row(i), scan_v.row(i),SD2D[0].row(i), SDZDR2D[0].row(i), scan_z, scan_w, scan_v, SD2D[0],i);
1099
1100 for (unsigned ib = 0; ib < beam_size; ++ib)
1101 if (corrected[ib])
1102 {
1103 //scan_z(i, ib) = cleaner.Z_missing;
1104 scan_z(i, ib) = new_val;
1105 // scan_w(i, ib) = cleaner.W_threshold;
1106 // scan_v(i, ib) = cleaner.V_missing;
1107 }
1108// img_tmp(i,ib)=255;
1109// z_clean(i,ib)=0;
1110// } else img_tmp(i,ib)= 0 ;
1111
1112 }
1113//radarelab::write_image(img_tmp,"/ponte/rad_svn/proc_operative/test_arch/rev_actual/radar/immagini/Cleaner/PPI_clean"+ext,"PNG");
1114//radarelab::write_image(z_clean,"/ponte/rad_svn/proc_operative/test_arch/rev_actual/radar/immagini/Cleaner/PPI_Zclean"+ext,"PNG");
1115}
1116
1117void Cleaner::clean( radarelab::volume::Loader load_structure, double bin_wind_magic_number,unsigned iel, bool set_undetect)
1118{
1119 std::string Z_Quantity;
1120
1121}
1122
1123double Cleaner::trap(double x1, double x2, double x3, double x4, double val, double x5) const
1124{
1125 if((val<=x3&&val>=x2)) return 1.;
1126 else if(val<x2&&val>x1) return val/(x2-x1)-x1/(x2-x1);
1127 else if (val<x4&&val>x3) return val/(x3-x4)-x4/(x3-x4);
1128 else if(val<=x5) return 1.;
1129 else return 0.; // (val<=x1||val>=x4)
1130
1131}
1132
1133vector<string> Cleaner::read_matrix_from_txt(string fin) const
1134{
1135
1136 /*
1137 legge file txt e mette il contenuto in un vettore
1138 */
1139 vector<string> myVector;
1140 ifstream f(fin, ifstream::in);
1141 string line;
1142
1143 if(f.is_open()){
1144 while(getline(f,line)){
1145 stringstream stream (line);
1146 while( getline(stream, line, ' ')){
1147 //inserisco controllo commenti che li gestisce se attaccati tipo //commento
1148 if(line[0]=='\\'){ continue;}
1149 else{
1150 myVector.push_back(line);
1151 }
1152 }
1153 }
1154 }
1155 return myVector;
1156}
1157
1158
1159}
1160}
double nodata
Value used as 'no data' value.
Definition volume.h:116
double undetect
Minimum amount that can be measured.
Definition volume.h:118
double offset
Conversion factor.
Definition volume.h:122
double gain
Conversion factor.
Definition volume.h:120
PolarScan - structure to describe a polarScan containing a matrix of data and conversion factors.
Definition volume.h:113
T offset
Data Offset.
Definition volume.h:276
Sequence of PolarScans which can have a different beam count for each elevation.
Definition volume.h:264
void gdal_init_once()
Initialize the GDAL library when called for the first time; does nothing all other times.
Definition image.cpp:12
String functions.
Definition cart.cpp:4
Base for all matrices we use, since we rely on row-major data.
Definition matrix.h:37
unsigned beam_size
Number of samples in each beam.
Definition volume.h:33
unsigned beam_count
Count of beams in this scan.
Definition volume.h:30
double cell_size
Size of a beam cell in meters.
Definition volume.h:48
std::vector< unsigned char > eval_clean_beam(const Eigen::VectorXd &beam_z, const Eigen::VectorXd &beam_w, const Eigen::VectorXd &beam_v, int i) const
Funzione per ripulire raggio.Utilizza (sigmaV, V) Analoga a clean_beam(const Eigen::VectorXd& beam_z,...
Definition cleaner.cpp:105
double trap(double x1, double x2, double x3, double x4, double val, double x5=-9999.) const
Definition cleaner.cpp:1123
Cleaner(double Z_missing, double W_threshold, double V_missing, double bin_wind_magic_number)
Constructor.
Definition cleaner.h:33
std::vector< bool > clean_beam(const Eigen::VectorXd &beam_z, const Eigen::VectorXd &beam_w, const Eigen::VectorXd &beam_v, int i) const
Funzione per ripulire raggio.Utilizza (sigmaV, V)
Definition cleaner.cpp:38
const double sd_threshold
Soglia per devizione standard DBZH.
Definition cleaner.h:30
const double Z_missing
Valore dato mancante DBZH.
Definition cleaner.h:26
const double bin_wind_magic_number
valore magico per dati in formato SP20
Definition cleaner.h:29
const unsigned max_segment_length
lunghezza massima segmento in celle se più lungo pulisce in ogni caso
Definition cleaner.h:24
const double W_threshold
Soglia per WRAD.
Definition cleaner.h:27
const unsigned min_segment_length
lunghezza minima segmento in celle
Definition cleaner.h:23
Struttura per cleaner dati grezzi sulla base dei valori di V, W e la deviazione standard di Z.
Definition cleaner.h:22
Struttura che contiene mappa per caricamento dati.
Definition loader.h:25