Elaboradar  0.1
 Tutto Classi Namespace File Funzioni Variabili Tipi enumerati (enum) Gruppi
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 namespace radarelab {
24 namespace algo {
25 
26 using namespace std;
27 using namespace radarelab;
28 
29 //--------------------------------------------------------------------------------
30 // These methods use only VRAD and WRAD values to clean the beam
31 //
32 std::vector<bool> Cleaner::clean_beam(const Eigen::VectorXd& beam_z, const Eigen::VectorXd& beam_w, const Eigen::VectorXd& beam_v, int iray) const
33 {
34  const unsigned beam_size = beam_z.rows();
35  vector<bool> res(beam_size, false);
36  bool in_a_segment = false;
37  unsigned start, end;
38  unsigned segment_length;
39  bool before, after;
40  unsigned counter = 0;
41 
42  for (unsigned ibin = 0; ibin < beam_size; ++ibin)
43  {
44 //printf(" %4d %4d %6.2f %6.2f %10.6f ",iray,ibin , beam_z(ibin),beam_v(ibin),beam_w(ibin));
45  if (!in_a_segment)
46  {
47  /* cerco la prima cella segmento da pulire*/
48  if (beam_w(ibin) == W_threshold && beam_v(ibin) == bin_wind_magic_number)
49  {
50 //printf(" 1 ----- START SEGMENT ------");
51  in_a_segment = true;
52  start = ibin;
53  after = false;
54  before = false;
55  }
56 // else printf(" 0 ");
57  } else {
58  /* cerco la fine segmento da pulire*/
59  if (beam_w(ibin) != W_threshold || beam_v(ibin) != bin_wind_magic_number || ibin == (beam_size - 1))
60  {
61  in_a_segment = false;
62  end = ibin - 1;
63  if (ibin == (beam_size - 1)) end = ibin; // caso particolare per fine raggio
64  /* Fine trovata ora procedo alla pulizia eventuale */
65  segment_length = end - start;
66  counter = counter + (unsigned)(segment_length);
67 
68  unsigned c_b=0;
69  unsigned c_a=0;
70  /* Cerco dati validi in Z prima del segmento */
71  for (int ib = ibin - 12; ib < (signed)ibin; ++ib)
72  if (ib >= 0 && beam_z(ib) > Z_missing)
73  c_b++;
74  if (c_b > 0.25*12) before = true;
75 
76  /* Cerco dati validi in Z dopo il segmento */
77  for (unsigned ia = ibin + 1; ia <= ibin + 12; ++ia)
78  if (ia < beam_size && beam_z(ia) >= Z_missing)
79  c_a++;
80  if (c_a > 0.25*12) after = true;
81 
82 //printf(" 0 ----- STOP SEGMENT ------ %4d -- %4d before %d %d after %d %d ",segment_length,counter, c_b,before, c_a, after);
83  if ((segment_length >= min_segment_length && !before && !after) ||
84  segment_length >= max_segment_length || counter > 100)
85  {
86  /* qui pulisco */
87  // printf (" pulisco %d %d %d \n",segment_length, min_segment_length, max_segment_length);
88  for (unsigned ib = start; ib <= end; ++ib)
89  if( beam_z(ib) > Z_missing)res[ib] = true;
90  }
91  }
92 // else printf(" 1 ");
93  }
94 //printf("\n");
95  }
96  return res;
97 }
98 
99 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
100 {
101  const unsigned beam_size = beam_z.rows();
102  vector<unsigned char> res(beam_size, 0);
103  bool in_a_segment = false;
104  unsigned start = 0, end = 0;
105  unsigned segment_length;
106  bool before = false, after = false;
107  unsigned counter = 0;
108 
109  for (unsigned ibin = 0; ibin < beam_size; ++ibin)
110  {
111 //printf(" %4d %4d %6.2f %6.2f %10.6f ",iray,ibin , beam_z(ibin),beam_v(ibin),beam_w(ibin));
112  if (!in_a_segment)
113  {
114  /* search the first radar bin's segment to be cleaned*/
115  if (beam_w(ibin) == W_threshold && beam_v(ibin) == bin_wind_magic_number)
116  {
117 //printf(" 1 ----- START SEGMENT ------");
118  in_a_segment = true;
119  start = ibin;
120  after = false;
121  before = false;
122  }
123  } else {
124  /* search the last radar bin's segment to be cleaned*/
125  if (beam_w(ibin) != W_threshold || beam_v(ibin) != bin_wind_magic_number || ibin == (beam_size - 1))
126  {
127  in_a_segment = false;
128  end = ibin - 1;
129  if (ibin == (beam_size - 1)) end = ibin; // beam ended
130  /* Fine trovata ora procedo alla pulizia eventuale */
131  segment_length = end - start;
132  counter = counter + (unsigned)(segment_length);
133 
134  unsigned c_b=0;
135  unsigned c_a=0;
136  /* Cerco dati validi in Z prima del segmento */
137  for (int ib = ibin - 12; ib < (signed)ibin; ++ib)
138  if (ib >= 0 && beam_z(ib) > Z_missing)
139  c_b++;
140  if (c_b > 0.25*12) before = true;
141 
142  /* Cerco dati validi in Z dopo il segmento */
143  for (unsigned ia = ibin + 1; ia <= ibin + 12; ++ia)
144  if (ia < beam_size && beam_z(ia) >= Z_missing)
145  c_a++;
146  if (c_a > 0.25*12) after = true;
147 
148 //printf(" 0 ----- STOP SEGMENT ------ %4d -- %4d before %d %d after %d %d ",segment_length,counter, c_b,before, c_a, after);
149  if ((segment_length >= min_segment_length && !before && !after) ||
150  segment_length >= max_segment_length || counter > 100)
151  {
152  /* qui pulisco */
153  // printf (" pulisco %d %d %d \n",segment_length, min_segment_length, max_segment_length);
154  for (unsigned ib = start; ib <= end; ++ib)
155  if( beam_z(ib) > Z_missing)res[ib] = 1;
156  }
157  }
158 // else printf(" 1 ");
159  }
160 //printf("\n");
161  }
162  return res;
163 }
164 //---------------------------------------------------------------------------------
165 
166 //---------------------------------------------------------------------------------
167 // This method use VRAD, WRAD, sdDBZH, sdZDR values to clean the beam
168 //
169 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
170 {
171  const unsigned beam_size = beam_z.rows();
172  vector<bool> res(beam_size, false);
173  bool in_a_segment = false;
174  unsigned start = 0, end;
175  unsigned segment_length;
176  bool before, after;
177  unsigned counter = 0;
178  unsigned counter_trash = 0;
179  unsigned counter_clutter =0;
180  for (unsigned ibin = 0; ibin < beam_size; ++ibin)
181  {
182  bool is_clutter = false;
183  bool is_trash = false;
184  unsigned flag = 0 ;
185 // In our systems (ARPA ER) interferences and other non meteo echo are characterised by the following steps
186 //
187 // 1) Wind is not defined ad spectrumWidth is 0. with Z defined.
188  if ( beam_w(ibin) == W_threshold && beam_v(ibin) == bin_wind_magic_number && beam_z (ibin) != Z_missing ) {
189 // 2) Std<Dev of ZDR coulb be close to 0
190  if( beam_sdzdr(ibin) <= 0.01 ){
191  is_trash = true;
192  flag=2;
193  } else {
194  if (beam_z (ibin) >= 45. ){
195 // 2) inside thunderstorms (Z > 45) StdDev of Zdr and StdDev of Z are quite high)
196  if ((ibin >100 && double(counter_trash)/double(ibin) >=0.5 && ( beam_sdzdr(ibin) >1 || beam_sd (ibin) > 5. )) ||
197  (beam_sdzdr(ibin) >4.0 && beam_sd (ibin) > 20.) ) {
198  is_trash = true;
199  flag=2;
200  } else {
201  is_trash = false;
202  flag=0;
203  }
204  } else if ( (ibin >100 && double(counter_trash)/double(ibin) >=0.5 && ( beam_sdzdr(ibin) >1 || beam_sd (ibin) > 5. )) ||
205  (beam_sd (ibin) >2. && (beam_sdzdr(ibin) >2.0 || beam_sd (ibin) > 10. )) ) {
206 // 2) outside thunderstorms (Z > 45) StdDev of Zdr and StdDev of Z are lower
207  is_trash = true;
208  flag=2;
209  }
210  }
211  } else {
212 // 3) Clutter is characterised by low value of VRAD and WRAD
213  if ((beam_w(ibin) * fabs(beam_v(ibin))) <= 0.25 && beam_z (ibin) != Z_missing ) {
214  is_clutter = true;
215  flag = 1;
216  }
217  }
218  if( is_clutter) counter_clutter ++;
219  if( is_trash ) counter_trash ++;
220 if(ibin <40 && false){
221  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));
222  printf(" ----- %2x %2x %2x %2x ",(unsigned char)((beam_z(ibin)-scan_z.offset)/scan_z.gain/256),
223  (unsigned char)((beam_v(ibin)-scan_v.offset)/scan_v.gain/256),
224  (unsigned char)((beam_w(ibin)-scan_w.offset)/scan_w.gain/256),
225  (unsigned char)((beam_sd(ibin)-SD.offset)/SD.gain/256));
226 }
227  if (!in_a_segment)
228  {
229  /* cerco la prima cella segmento da pulire*/
230  if ( is_clutter || is_trash )
231  {
232 // if(ibin <40)printf(" %1d ----- START SEGMENT ------",flag);
233  in_a_segment = true;
234  start = ibin;
235  after = false;
236  before = false;
237  }
238 // else if(ibin <40)printf(" %1d ",flag);
239  } else {
240  /* cerco la fine segmento da pulire*/
241  if ( ! (is_clutter || is_trash ) || ibin == (beam_size - 1))
242  {
243  in_a_segment = false;
244  end = ibin - 1;
245  if (ibin == (beam_size - 1)) end = ibin; // caso particolare per fine raggio
246  /* Fine trovata ora procedo alla pulizia eventuale */
247  segment_length = end - start+1;
248  counter = counter + (unsigned)(segment_length);
249 
250 /* il segmento è corto allora cerco nei dintorni dei dati validi, se li trovo non pulisco */
251  if (segment_length <= 2*min_segment_length ){
252  /* Cerco dati validi in Z prima del segmento */
253  int count=0;
254  for (int ib = ibin - 2*min_segment_length; ib < (signed)ibin; ++ib)
255  if (ib >= 0 && (beam_z(ib) > Z_missing && beam_w(ib) != W_threshold && ( beam_w(ib) > 0.5 || fabs(beam_v(ib)) > 0.5) ) )
256  count++;
257  if (double(count)/double(min(int(ibin),int(2*min_segment_length))) >=0.25) before = true;
258 
259  /* Cerco dati validi in Z dopo il segmento */
260  count = 0;
261  for (unsigned ia = ibin + 1; ia <= ibin + 2*min_segment_length; ++ia)
262  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)) ))
263  count ++;
264  if (double(count)/double(min(int(beam_size - ibin),int(2*min_segment_length))) >=0.25) after = true;
265  }
266 // if(ibin <40)printf(" %1d ----- STOP SEGMENT ------ %4d -- %4d before %d after %d ",flag, segment_length,counter, before,after);
267  if ((segment_length >= min_segment_length && (!before || !after) ) || segment_length >= max_segment_length)
268  // if ((segment_length >= min_segment_length ) || segment_length >= max_segment_length)
269  {
270  /* qui pulisco */
271 // if(ibin <40)printf (" pulisco %d %d %d \n",segment_length, min_segment_length, max_segment_length);
272  for (unsigned ib = start; ib <= end; ++ib)
273  res[ib] = true;
274  }
275  }
276 // else if(ibin <40)printf(" %1d ",flag);
277 
278  }
279 // if(ibin <40)printf(" %4d %4d \n",counter_clutter,counter_trash);
280  }
281  return res;
282 }
283 //---------------------------------------------------------------------------------
284 
285 //---------------------------------------------------------------------------------
286 // These methods use VRAD, WRAD, sdDBZH, values to clean the beam
287 //
288 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
289 {
290  const unsigned beam_size = beam_z.rows();
291  vector<bool> res(beam_size, false);
292  bool in_a_segment = false;
293  unsigned start, end;
294  unsigned segment_length;
295  bool before, after;
296  unsigned counter = 0;
297 
298  for (unsigned ibin = 0; ibin < beam_size; ++ibin)
299  {
300 // 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));
301 //printf(" ----- %2x %2x %2x %2x ",(unsigned char)((beam_z(ibin)-scan_z.offset)/scan_z.gain/256),
302 //(unsigned char)((beam_v(ibin)-scan_v.offset)/scan_v.gain/256),
303 //(unsigned char)((beam_w(ibin)-scan_w.offset)/scan_w.gain/256),
304 //(unsigned char)((beam_sd(ibin)-SD.offset)/SD.gain/256));
305  if (!in_a_segment)
306  {
307  /* cerco la prima cella segmento da pulire*/
308  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 )
309  {
310 // printf(" 1 ----- START SEGMENT ------");
311  in_a_segment = true;
312  start = ibin;
313  after = false;
314  before = false;
315  }
316 // else printf(" 0 ");
317  } else {
318  /* cerco la fine segmento da pulire*/
319  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 )
320  {
321  in_a_segment = false;
322  end = ibin - 1;
323  if (ibin == (beam_size - 1)) end = ibin; // caso particolare per fine raggio
324  /* Fine trovata ora procedo alla pulizia eventuale */
325  segment_length = end - start+1;
326  counter = counter + (unsigned)(segment_length);
327 
328 /* il segmento è corto allora cerco nei dintorni dei dati validi, se li trovo non pulisco */
329  if (segment_length <= 2*min_segment_length ){
330  /* Cerco dati validi in Z prima del segmento */
331  int count=0;
332  for (int ib = ibin - 2*min_segment_length; ib < (signed)ibin; ++ib)
333  if (ib >= 0 && (beam_z(ib) > Z_missing && beam_w(ib) != W_threshold && ( beam_w(ib) > 0.5 || fabs(beam_v(ib)) > 0.5) ) )
334  count++;
335  if (double(count)/double(min(int(ibin),int(2*min_segment_length))) >=0.25) before = true;
336 
337  /* Cerco dati validi in Z dopo il segmento */
338  count = 0;
339  for (unsigned ia = ibin + 1; ia <= ibin + 2*min_segment_length; ++ia)
340  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)) ))
341  count ++;
342  if (double(count)/double(min(int(beam_size - ibin),int(2*min_segment_length))) >=0.25) after = true;
343  }
344 // printf(" 0 ----- STOP SEGMENT ------ %4d -- %4d before %d after %d ",segment_length,counter, before,after);
345  if ((segment_length >= min_segment_length && (!before || !after) ) || segment_length >= max_segment_length)
346  // if ((segment_length >= min_segment_length ) || segment_length >= max_segment_length)
347  {
348  /* qui pulisco */
349  // printf (" pulisco %d %d %d \n",segment_length, min_segment_length, max_segment_length);
350  for (unsigned ib = start; ib <= end; ++ib)
351  res[ib] = true;
352  }
353  }
354 // else printf(" 1 ");
355 
356  }
357 // printf("\n");
358  }
359  return res;
360 }
361 
362 
363 // Senza ZDR - Fuzzy logic
364 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
365 {
366 
367  const unsigned beam_size = beam_z.rows();
368  vector<unsigned char> res(beam_size, 0);
369  vector<unsigned> counter (6,0) ;
370 
371  Matrix2D<double> Wij(6,6);
372 // Z V Wrad SD2d SDRay SDAz
373  Wij << 1.0, 1.0, 1.0, 0.5, 0.3, 0.3, // METEO 4.1
374  0.8, 0.5, 0.5, 0.5, 0.5, 0.5, // CLUTTER 3.3
375  0.0, 0.2, 0.2, 0.6, 0.8, 0.5, // INTERF. Strong 2.3
376  0.0, 0.2, 0.2, 0.5, 0.6, 0.4, // INTERF. Med. 2.2
377  0.0, 0.2, 0.2, 0.6, 0.5, 0.5, // INTERF. Weak 1.6
378  0.6, 0.2, 0.2, 0.6, 0.5, 0.5; // NOISE 2.6
379 // TOT = 2.4 2.3 2.3 3.4 3.2 2.5
380 
381 
382 
383  for (unsigned ibin = 0; ibin < beam_size; ++ibin)
384  {
385 //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));
386  Matrix2D<double> Pij(6,6);
387  Pij = Pij * 0.;
388  ArrayXd Class_WP(6);
389  Class_WP.setZero();
390  if (beam_z(ibin) == Z_missing) {
391  unsigned ID=0;
392  res[ibin]=ID;
393  counter[ID]++;
394 //printf("%2d %5.3f %5.3f %5.3f %5.3f %5.3f %5.3f %4d %4d %4d %4d %4d\n",res[ibin],Class_WP(0),Class_WP(1),Class_WP(2),Class_WP(3),Class_WP(4),Class_WP(5),counter[1],counter[2], counter[3], counter[4], counter[5]);
395 //printf("%2d %5.3f %4d %4d %4d %4d %4d\n",res[ibin],Class_WP(ID),counter[1],counter[2], counter[3], counter[4], counter[5]);
396  continue;
397  }
398 
399 //eseguo un test unico per VRAD e WRAD e assegno tutte le prob assieme
400  if (beam_v(ibin) != bin_wind_magic_number ) { // VRAD
401  Pij(0,1)=1.; // METEO
402  double prob_v = trap (-0.5, -0.3, 0.3, 0.5, beam_v(ibin)); // CLUTTER
403  Pij(1,1)=prob_v; // INTERF. Strong
404  Pij(2,1)=prob_v; // INTERF. Med.
405  Pij(3,1)=prob_v; // INTERF. Weak
406  Pij(4,1)=prob_v; // NOISE
407  Pij(5,1)=prob_v;
408  } else {
409  Pij(0,1)=0.3; // METEO
410  Pij(1,1)=1.; // CLUTTER
411  Pij(2,1)=1.; // INTERF. Strong
412  Pij(3,1)=1.; // INTERF. Med.
413  Pij(4,1)=1.; // INTERF. Weak
414  Pij(5,1)=1.; // NOISE
415  }
416 
417  if (beam_w(ibin) > W_threshold) { // WRAD
418  Pij(0,2)=1.; // METEO
419  double prob_w = trap (0.0001, 0., 0.2, 0.3, beam_w(ibin)); // CLUTTER
420  Pij(1,2)=prob_w; // INTERF. Strong
421  Pij(2,2)=prob_w; // INTERF. Med.
422  Pij(3,2)=prob_w; // INTERF. Weak
423  Pij(4,2)=prob_w; // NOISE
424  Pij(5,2)=prob_w;
425  } else {
426  Pij(0,2)=0.3; // METEO
427  Pij(1,2)=1.; // CLUTTER
428  Pij(2,2)=1.; // INTERF. Strong
429  Pij(3,2)=1.; // INTERF. Med.
430  Pij(4,2)=1.; // INTERF. Weak
431  Pij(5,2)=1.; // NOISE
432  }
433 
434 // METEO
435  Pij(0,0) = 1.; // Z
436  Pij(0,3) = trap (0.5, 1., 10.,11., beam_sd(ibin)); // SD_2D
437  Pij(0,4) = trap (0.5, 1., 10.,11., beam_sdray(ibin)); // SD_RAY
438  Pij(0,5) = trap (0.5, 1., 10.,11., beam_sdaz(ibin)); // SD_AZ
439 
440 // CLUTTER
441  Pij(1,0) = trap (5., 15., 99., 99.9, beam_z(ibin)); // Z
442  Pij(1,3) = trap (4.5, 5., 99., 99.9, beam_sd(ibin)); // SD_2D
443  Pij(1,4) = trap (4., 4.5, 99., 99.9, beam_sdray(ibin)); // SD_RAY
444  Pij(1,5) = trap (4., 4.5, 99., 99.9, beam_sdaz(ibin)); // SD_AZ
445 // INTERF. Strong
446  Pij(2,0) = 1.; // Z
447  Pij(2,3) = trap (0.5, 1.5, 5., 7., beam_sd(ibin)); // SD_2D
448  Pij(2,4) = trap (0.0, 0.3, 2., 3., beam_sdray(ibin)); // SD_RAY
449  Pij(2,5) = trap (2.0, 4., 90.,90.9, beam_sdaz(ibin)); // SD_AZ
450 
451 // INTERF. Med.
452  Pij(3,0) = 1.; // Z
453  Pij(3,3) = trap (4.0, 5., 90., 90.9, beam_sd(ibin)); // SD_2D
454  Pij(3,4) = trap (0., 0.5, 10., 15.0, beam_sdray(ibin)); // SD_RAY
455  Pij(3,5) = trap (3., 4., 90., 90.9, beam_sdaz(ibin)); // SD_AZ
456 
457 // INTERF. Weak
458  Pij(4,0) = 1.; // Z
459  Pij(4,3) = trap (4., 7., 90., 90.9, beam_sd(ibin)); // SD_2D
460  Pij(4,4) = trap (0.5, 1.5, 3., 4., beam_sdray(ibin)); // SD_RAY
461  Pij(4,5) = trap (3., 4., 90., 90.9, beam_sdaz(ibin)); // SD_AZ
462 
463 // NOISE
464  double coeff = 1.;
465  if (ibin >= 40 ) {
466  if ( ibin <= 160) coeff = (160 - ibin)/ 120. ;
467  else coeff = 0.;
468  }
469  Pij(5,0) = trap(Z_missing-0.0001, Z_missing, 15., 20., beam_z(ibin)); // Z
470  Pij(5,3) = trap (5., 7., 99., 99.9, beam_sd(ibin)); // SD_2D
471  Pij(5,4) = trap (0., 0.001, 0.8 + coeff * 9.2 , 1. + coeff * 14.,beam_sdray(ibin)); // SD_RAY
472 // Pij(5,4) = trap (0., 0.001, 0.8, 1.,beam_sdray(ibin)); // SD_RAY
473 // Pij(5,4) = trap (0., 0.001, 10., 15.,beam_sdray(ibin)); // SD_RAY
474  Pij(5,5) = trap (1.5, 3., 99., 99.9,beam_sdaz(ibin)); // SD_AZ
475 
476 //---- fine calcolo probabilità
477 // Calcolo classe appartenenza
478 
479  Class_WP = ((Wij.array()*Pij.array()).matrix()*VectorXd::Ones(6)).array()/(Wij*VectorXd::Ones(6)).array();
480  unsigned i,ID;
481  Class_WP.maxCoeff(&i);
482  ID=i;
483  if (Class_WP(i) < 0.1 ) ID=6;
484  res[ibin]=ID;
485  counter[ID]++;
486 //printf("%2d %5.3f %5.3f %5.3f %5.3f %5.3f %5.3f %4d %4d %4d %4d %4d\n",res[ibin],Class_WP(0),Class_WP(1),Class_WP(2),Class_WP(3),Class_WP(4),Class_WP(5),counter[1],counter[2], counter[3], counter[4], counter[5]);
487 // if (iray == 225 and ibin <= 50) {
488 // std::cout <<"Valori per fuzzy - Z_missing "<<Z_missing<< "differenza " <<beam_z(ibin)-Z_missing << std::endl;
489 // std::cout << Pij<<std::endl;
490 // }
491  }
492 
493  return res;
494 }
495 
496 // Senza ZDR - Basato su test
497 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
498 {
499 
500  const unsigned beam_size = beam_z.rows();
501  vector<unsigned char> res(beam_size, 0);
502  unsigned counter = 0;
503  unsigned countIntWeak = 0;
504  unsigned countIntStrong = 0;
505  unsigned countIntMed = 0;
506  unsigned countClutter = 0;
507  unsigned countNoise = 0;
508 
509  for (unsigned ibin = 0; ibin < beam_size; ++ibin)
510  {
511 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));
512  if ((beam_w(ibin) > W_threshold && beam_v(ibin) != bin_wind_magic_number && beam_sd(ibin) >= 1. &&
513  beam_sd(ibin) <= 10. ) || beam_z(ibin) == Z_missing) {
514  // this should be a meteorological echo
515  ;
516  }else {
517  res [ibin]=1; // Not meteo but unclassified
518  counter++;
519  if (beam_w(ibin) == W_threshold && beam_v(ibin) == bin_wind_magic_number)
520  {
521  if (beam_sd(ibin) >= 5. && beam_sdray(ibin) >= 2 && beam_sdaz(ibin) > 4.){
522  // this should be clutter
523  res[ibin] = 2;
524  countClutter ++;
525  }
526  if (beam_sd(ibin) >= 1. && beam_sd(ibin) <= 5. && beam_sdray(ibin) < 2. && beam_sdaz(ibin) > 4.){
527  // this should be a strong Interference
528  res[ibin] = 3;
529  countIntStrong ++;
530  }
531  if (beam_sd(ibin) >= 3. && beam_sd(ibin) <= 7. && beam_sdray(ibin) < 2. && beam_sdaz(ibin) > 4.&& beam_sdaz(ibin) <7 ){
532  // this should be a medium Interference
533  res[ibin] = 4;
534  countIntMed ++;
535  }
536  //if (beam_sd(ibin) >= 5. && beam_sd(ibin) <= 3. && beam_sdray(ibin) < 3. && beam_sdaz(ibin) < 3.){
537  if (beam_sd(ibin) >= 5. && beam_sdray(ibin) < 3. && beam_sdaz(ibin) < 3.){
538  // this should be a weakInterference
539  res[ibin] = 5;
540  countIntWeak ++;
541  }
542  if (beam_sd(ibin) >= 10. && beam_sdray(ibin) < 2. && beam_sdaz(ibin) < 5. && beam_z(ibin) <10.){
543  // this should be a noise
544  res[ibin] = 6;
545  countNoise ++;
546  }
547  }
548  } // ELSE this should be not a meteo echo
549 printf("%2d %4d %4d %4d %4d %4d\n",res[ibin],countClutter,countIntStrong, countIntMed, countIntWeak, countNoise);
550  }
551 
552  return res;
553 
554 }
555 //----------------------------------------------------------------------------------
556 
557 
558 
559 void Cleaner::evaluateCleanID(PolarScan<double>& scan_z, PolarScan<double>& scan_w, PolarScan<double>& scan_v,PolarScan<unsigned char>& scan_cleanID, unsigned iel)
560 {
561  //return evaluateCleanID(scan_z, scan_w, scan_v, scan_cleanID, scan_v.undetect,iel);
562  return evaluateClassID(scan_z, scan_w, scan_v, scan_cleanID, scan_v.undetect,iel);
563 }
564 
565 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)
566 {
567 
568  if (scan_z.beam_count != scan_w.beam_count)
569  throw std::runtime_error("scan_z beam_count no equal to scan_w beam_count");
570  if (scan_z.beam_size != scan_w.beam_size)
571  throw std::runtime_error("scan_z beam_size no equal to scan_w beam_size");
572 
573  if (scan_z.beam_count != scan_v.beam_count)
574  throw std::runtime_error("scan_z beam_count no equal to scan_v beam_count");
575  if (scan_z.beam_size != scan_v.beam_size)
576  throw std::runtime_error("scan_z beam_size no equal to scan_v beam_size");
577 
578  Cleaner cleaner(scan_z.undetect, scan_w.undetect, scan_v.nodata, bin_wind_magic_number);
579 
580  const unsigned beam_count = scan_z.beam_count;
581  const unsigned beam_size = scan_z.beam_size;
582 
583  // fprintf(stderr, "NEWCLEANER zmis %f, wthr %f, vmis %f, mn %f\n",
584  // cleaner.Z_missing, cleaner.W_threshold, cleaner.V_missing, cleaner.bin_wind_magic_number);
585 
586 // radarelab::volume::Scans<double> Z_S, SD2D;
587 // Z_S.push_back(scan_z);
588 // radarelab::volume::textureSD( Z_S,SD2D, 1000. , 3,false);
589 
590  radarelab::volume::Scans<double> Z_S, SD2D,SD_Ray,SD_Az;
591  Z_S.push_back(scan_z);
592  radarelab::volume::textureSD( Z_S,SD2D, 1000. , 3,false);
593 
594  radarelab::volume::textureSD( Z_S,SD_Ray, scan_z.cell_size*9 , 360./scan_z.beam_count,true);
595  radarelab::volume::textureSD( Z_S,SD_Az, scan_z.cell_size , 5*360./scan_z.beam_count,true);
596 
597  for (unsigned i = 0; i < beam_count; ++i)
598  {
599  //vector<unsigned char> corrected = cleaner.eval_clean_beam(scan_z.row(i), scan_w.row(i), scan_v.row(i),i);
600  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);
601  //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);
602  for (unsigned ib = 0; ib < beam_size; ++ib)
603  scan_cleanID(i,ib)=corrected[ib];
604  }
605 }
606 
607 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, unsigned iel)
608 {
609 
610  if (scan_z.beam_count != scan_w.beam_count)
611  throw std::runtime_error("scan_z beam_count no equal to scan_w beam_count");
612  if (scan_z.beam_size != scan_w.beam_size)
613  throw std::runtime_error("scan_z beam_size no equal to scan_w beam_size");
614 
615  if (scan_z.beam_count != scan_v.beam_count)
616  throw std::runtime_error("scan_z beam_count no equal to scan_v beam_count");
617  if (scan_z.beam_size != scan_v.beam_size)
618  throw std::runtime_error("scan_z beam_size no equal to scan_v beam_size");
619 
620  Cleaner cleaner(scan_z.undetect, scan_w.undetect, scan_v.nodata, bin_wind_magic_number);
621 
622  const unsigned beam_count = scan_z.beam_count;
623  const unsigned beam_size = scan_z.beam_size;
624 
625 // compute texture volumes
626  radarelab::volume::Scans<double> Z_S, SD2D,SD_Ray,SD_Az;
627  Z_S.push_back(scan_z);
628  radarelab::volume::textureSD( Z_S,SD2D, 1000. , 3,false);
629  radarelab::volume::textureSD( Z_S,SD_Ray, scan_z.cell_size*21 , 360./scan_z.beam_count,true);
630  radarelab::volume::textureSD( Z_S,SD_Az, scan_z.cell_size , 5*360./scan_z.beam_count,true);
631 
632 
633  for (unsigned i = 0; i < beam_count; ++i)
634  {
635  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);
636  for (unsigned ib = 0; ib < beam_size; ++ib)
637  scan_cleanID(i,ib)=corrected[ib];
638  }
639 }
640 
641 void Cleaner::clean(PolarScan<double>& scan_z, PolarScan<double>& scan_w, PolarScan<double>& scan_v, unsigned iel , bool set_undetect)
642 {
643  return clean(scan_z, scan_w, scan_v, scan_v.undetect,iel);
644 }
645 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)
646 {
647 
648  if (scan_z.beam_count != scan_w.beam_count)
649  throw std::runtime_error("scan_z beam_count no equal to scan_w beam_count");
650  if (scan_z.beam_size != scan_w.beam_size)
651  throw std::runtime_error("scan_z beam_size no equal to scan_w beam_size");
652 
653  if (scan_z.beam_count != scan_v.beam_count)
654  throw std::runtime_error("scan_z beam_count no equal to scan_v beam_count");
655  if (scan_z.beam_size != scan_v.beam_size)
656  throw std::runtime_error("scan_z beam_size no equal to scan_v beam_size");
657 
658  Cleaner cleaner(scan_z.undetect, scan_w.undetect, scan_v.nodata, bin_wind_magic_number);
659 
660  const unsigned beam_count = scan_z.beam_count;
661  const unsigned beam_size = scan_z.beam_size;
662 
663  // fprintf(stderr, "NEWCLEANER zmis %f, wthr %f, vmis %f, mn %f\n",
664  // cleaner.Z_missing, cleaner.W_threshold, cleaner.V_missing, cleaner.bin_wind_magic_number);
665 
666  radarelab::volume::Scans<double> Z_S, SD2D,SD_Ray,SD_Az;
667  Z_S.push_back(scan_z);
668  radarelab::volume::textureSD( Z_S,SD2D, 1000. , 3,false);
669  radarelab::volume::textureSD( Z_S,SD_Ray, scan_z.cell_size*9 , 360./scan_z.beam_count,false);
670  radarelab::volume::textureSD( Z_S,SD_Az, scan_z.cell_size , 5*360./scan_z.beam_count,false);
671 
672 //radarelab::gdal_init_once();
673 //printf("scrivo Z ");
674 //Matrix2D <double>img;
675 //img = (scan_z.array() - scan_z.offset )/ scan_z.gain /256 ;
676 //Matrix2D <unsigned char>img_tmp, z_clean;
677 //std::string ext;
678 //char pippo[200];
679 //sprintf(pippo, "_%02d.png",iel);
680 //ext=pippo;
681 
682 //img_tmp=img.cast<unsigned char>();
683 //z_clean=img_tmp;
684 //radarelab::write_image(img_tmp,"/ponte/rad_svn/proc_operative/test_arch/rev_actual/radar/immagini/Cleaner/PPI_Z"+ext, "PNG");
685 
686 //printf("V ");
687 //img = (scan_v.array()-scan_v.offset)/scan_v.gain/256 ;
688 //img_tmp=img.cast<unsigned char>();
689 //radarelab::write_image(img_tmp,"/ponte/rad_svn/proc_operative/test_arch/rev_actual/radar/immagini/Cleaner/PPI_V"+ext,"PNG");
690 //printf("W ");
691 //img = (scan_w.array()-scan_w.offset)/scan_w.gain/256 ;
692 //img_tmp=img.cast<unsigned char>();
693 //radarelab::write_image(img_tmp,"/ponte/rad_svn/proc_operative/test_arch/rev_actual/radar/immagini/Cleaner/PPI_W"+ext,"PNG");
694 //printf("SD2d ");
695 //img = (SD2D[0].array()-SD2D[0].offset)/SD2D[0].gain/256 ;
696 //img_tmp=img.cast<unsigned char>();
697 //radarelab::write_image(img_tmp,"/ponte/rad_svn/proc_operative/test_arch/rev_actual/radar/immagini/Cleaner/PPI_SD2d"+ext,"PNG");
698 //printf("\n");
699 
700  double new_val=cleaner.Z_missing;
701  if (set_undetect) new_val=scan_z.nodata;
702 
703  for (unsigned i = 0; i < beam_count; ++i)
704  {
705  // Compute which elements need to be cleaned
706  // 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);
707  vector<bool> corrected = cleaner.clean_beam(scan_z.row(i), scan_w.row(i), scan_v.row(i),i);
708 
709  for (unsigned ib = 0; ib < beam_size; ++ib)
710  if (corrected[ib])
711  {
712  //scan_z(i, ib) = cleaner.Z_missing;
713  scan_z(i, ib) = new_val;
714  // scan_w(i, ib) = cleaner.W_threshold;
715  // scan_v(i, ib) = cleaner.V_missing;
716  }
717 // img_tmp(i,ib)=255;
718 // z_clean(i,ib)=0;
719 // } else img_tmp(i,ib)= 0 ;
720 
721  }
722 //radarelab::write_image(img_tmp,"/ponte/rad_svn/proc_operative/test_arch/rev_actual/radar/immagini/Cleaner/PPI_clean"+ext,"PNG");
723 //radarelab::write_image(z_clean,"/ponte/rad_svn/proc_operative/test_arch/rev_actual/radar/immagini/Cleaner/PPI_Zclean"+ext,"PNG");
724 }
725 
726 
727 
728 void Cleaner::clean(PolarScan<double>& scan_z, PolarScan<double>& scan_w, PolarScan<double>& scan_v, PolarScan<double>& scan_zdr, unsigned iel, bool set_undetect )
729 {
730  return clean(scan_z, scan_w, scan_v, scan_zdr, scan_v.undetect,iel);
731 }
732 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 )
733 {
734  if (scan_z.beam_count != scan_w.beam_count)
735  throw std::runtime_error("scan_z beam_count no equal to scan_w beam_count");
736  if (scan_z.beam_size != scan_w.beam_size)
737  throw std::runtime_error("scan_z beam_size no equal to scan_w beam_size");
738 
739  if (scan_z.beam_count != scan_v.beam_count)
740  throw std::runtime_error("scan_z beam_count no equal to scan_v beam_count");
741  if (scan_z.beam_size != scan_v.beam_size)
742  throw std::runtime_error("scan_z beam_size no equal to scan_v beam_size");
743 
744  Cleaner cleaner(scan_z.undetect, scan_w.undetect, scan_v.nodata, bin_wind_magic_number);
745 
746  const unsigned beam_count = scan_z.beam_count;
747  const unsigned beam_size = scan_z.beam_size;
748 
749  // fprintf(stderr, "NEWCLEANER zmis %f, wthr %f, vmis %f, mn %f\n",
750  // cleaner.Z_missing, cleaner.W_threshold, cleaner.V_missing, cleaner.bin_wind_magic_number);
751 
753  Z_S.push_back(scan_z);
754  radarelab::volume::textureSD( Z_S,SD2D, 1000. , 3,false);
755  radarelab::volume::Scans<double> ZDR_S, SDZDR2D;
756  ZDR_S.push_back(scan_zdr);
757  radarelab::volume::textureSD( ZDR_S,SDZDR2D, 1000. , 3,false);
758 
759 //----------------------------------------------------------------------------------
760 //----------------------------------------------------------------------------------
761 //----------------------------------------------------------------------------------
762 // Mettere a true per fare grafica per debug o false per non fare grafica
763 //
764 // RICORDARSI DI TOGLIERE/METTERE COMMENTI DOPO CLEAN_BEAM
765 // -------------------------------------------------------------------
766 Matrix2D <unsigned char>img_tmp, z_clean;
768 std::string ext;
769 char pippo[200];
770 if (false){
772 
773  printf("scrivo Z ");
774  img = (scan_z.array() - scan_z.offset )/ scan_z.gain /256 ;
775  sprintf(pippo, "_%02d.png",iel);
776  ext=pippo;
777  img_tmp=img.cast<unsigned char>();
778  z_clean=img_tmp;
779  radarelab::write_image(img_tmp,"/ponte/rad_svn/proc_operative/test_arch/rev_actual/radar/immagini/Cleaner/PPI_Z"+ext, "PNG");
780 
781 // printf("V ");
782 // img = (scan_v.array()-scan_v.offset)/scan_v.gain/256 ;
783 // img_tmp=img.cast<unsigned char>();
784 // radarelab::write_image(img_tmp,"/ponte/rad_svn/proc_operative/test_arch/rev_actual/radar/immagini/Cleaner/PPI_V"+ext,"PNG");
785 // printf("W ");
786 // img = (scan_w.array()-scan_w.offset)/scan_w.gain/256 ;
787 // img_tmp=img.cast<unsigned char>();
788 // radarelab::write_image(img_tmp,"/ponte/rad_svn/proc_operative/test_arch/rev_actual/radar/immagini/Cleaner/PPI_W"+ext,"PNG");
789  printf("SD2d ");
790  img = (SD2D[0].array()-SD2D[0].offset)/SD2D[0].gain/256 ;
791  img_tmp=img.cast<unsigned char>();
792  radarelab::write_image(img_tmp,"/ponte/rad_svn/proc_operative/test_arch/rev_actual/radar/immagini/Cleaner/PPI_SD2d"+ext,"PNG");
793  printf("SDZDR2d ");
794  img = (SDZDR2D[0].array()-SDZDR2D[0].offset)/SDZDR2D[0].gain/256 ;
795  img_tmp=img.cast<unsigned char>();
796  radarelab::write_image(img_tmp,"/ponte/rad_svn/proc_operative/test_arch/rev_actual/radar/immagini/Cleaner/PPI_SDZDR2d"+ext,"PNG");
797  printf("\n");
798 }
799  double new_val=cleaner.Z_missing;
800  if (set_undetect) new_val=scan_z.undetect;
801 
802  for (unsigned i = 0; i < beam_count; ++i)
803  {
804  // Compute which elements need to be cleaned
805  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);
806 
807  for (unsigned ib = 0; ib < beam_size; ++ib)
808  if (corrected[ib])
809  {
810  //scan_z(i, ib) = cleaner.Z_missing;
811  scan_z(i, ib) = new_val;
812  // scan_w(i, ib) = cleaner.W_threshold;
813  // scan_v(i, ib) = cleaner.V_missing;
814  }
815 // img_tmp(i,ib)=255;
816 // z_clean(i,ib)=0;
817 // } else img_tmp(i,ib)= 0 ;
818 
819  }
820 //radarelab::write_image(img_tmp,"/ponte/rad_svn/proc_operative/test_arch/rev_actual/radar/immagini/Cleaner/PPI_clean"+ext,"PNG");
821 //radarelab::write_image(z_clean,"/ponte/rad_svn/proc_operative/test_arch/rev_actual/radar/immagini/Cleaner/PPI_Zclean"+ext,"PNG");
822 }
823 
824 void Cleaner::clean( radarelab::volume::Loader load_structure, double bin_wind_magic_number,unsigned iel, bool set_undetect)
825 {
826  std::string Z_Quantity;
827 
828 }
829 
830 
831 double Cleaner::trap(double x1, double x2, double x3, double x4, double val) const
832 {
833  if(val<=x3&&val>=x2) return 1.;
834  else if(val<x2&&val>x1) return val/(x2-x1)-x1/(x2-x1);
835  else if (val<x4&&val>x3) return val/(x3-x4)-x4/(x3-x4);
836  else return 0.; // (val<=x1||val>=x4)
837 }
838 
839 
840 }
841 }
double gain
Conversion factor.
Definition: volume.h:120
double undetect
Minimum amount that can be measured.
Definition: volume.h:118
Struttura per cleaner dati grezzi sulla base dei valori di V, W e la deviazione standard di Z...
Definition: cleaner.h:17
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:32
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&amp; beam_z...
Definition: cleaner.cpp:99
PolarScan - structure to describe a polarScan containing a matrix of data and conversion factors...
Definition: volume.h:112
double cell_size
Size of a beam cell in meters.
Definition: volume.h:48
unsigned beam_count
Count of beams in this scan.
Definition: volume.h:30
double trap(double x1, double x2, double x3, double x4, double val) const
Definition: cleaner.cpp:831
double offset
Conversion factor.
Definition: volume.h:122
static void clean(PolarScan< double > &scan_Z, PolarScan< double > &scan_W, PolarScan< double > &scan_V, unsigned iel=0, bool set_undetect=false)
Funzione che crea l&#39;oggetto cleaner, lo inizializza, pulisce i dati e modifica il PolarScan di DBZH...
Definition: cleaner.cpp:641
unsigned beam_size
Number of samples in each beam.
Definition: volume.h:33
Struttura che contiene mappa per caricamento dati.
Definition: loader.h:24
const double Z_missing
Valore dato mancante DBZH.
Definition: cleaner.h:22
void gdal_init_once()
Initialize the GDAL library when called for the first time; does nothing all other times...
Definition: image.cpp:12
double nodata
Value used as &#39;no data&#39; value.
Definition: volume.h:116