Elaboradar  0.1
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 
28 namespace radarelab {
29 namespace algo {
30 
31 using namespace std;
32 using namespace radarelab;
33  //using std::filesystem::current_path;
34 
35 //--------------------------------------------------------------------------------
36 // These methods use only VRAD and WRAD values to clean the beam
37 //
38 std::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 
105 std::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 //
175 std::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 //
294 std::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
739 std::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 
801 void 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 
807 void 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 
849 void 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
897 void 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 
931 void 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 }
935 void 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 
1018 void 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 }
1022 void 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 // -------------------------------------------------------------------
1059 Matrix2D <unsigned char>img_tmp, z_clean;
1061 std::string ext;
1062 char pippo[200];
1063 if (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 
1117 void 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 
1123 double 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 
1133 vector<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:274
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
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
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 Z_missing
Valore dato mancante DBZH.
Definition: cleaner.h:26
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