38std::vector<bool>
Cleaner::clean_beam(
const Eigen::VectorXd& beam_z,
const Eigen::VectorXd& beam_w,
const Eigen::VectorXd& beam_v,
int iray)
const
40 const unsigned beam_size = beam_z.rows();
41 vector<bool> res(beam_size,
false);
42 bool in_a_segment =
false;
44 unsigned segment_length;
48 for (
unsigned ibin = 0; ibin < beam_size; ++ibin)
69 if (ibin == (beam_size - 1)) end = ibin;
71 segment_length = end - start;
72 counter = counter + (unsigned)(segment_length);
77 for (
int ib = ibin - 12; ib < (signed)ibin; ++ib)
80 if (c_b > 0.25*12) before =
true;
83 for (
unsigned ia = ibin + 1; ia <= ibin + 12; ++ia)
84 if (ia < beam_size && beam_z(ia) >=
Z_missing)
86 if (c_a > 0.25*12) after =
true;
94 for (
unsigned ib = start; ib <= end; ++ib)
95 if( beam_z(ib) >
Z_missing)res[ib] =
true;
105std::vector<unsigned char>
Cleaner::eval_clean_beam(
const Eigen::VectorXd& beam_z,
const Eigen::VectorXd& beam_w,
const Eigen::VectorXd& beam_v,
int iray)
const
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;
115 for (
unsigned ibin = 0; ibin < beam_size; ++ibin)
133 in_a_segment =
false;
135 if (ibin == (beam_size - 1)) end = ibin;
137 segment_length = end - start;
138 counter = counter + (unsigned)(segment_length);
143 for (
int ib = ibin - 12; ib < (signed)ibin; ++ib)
146 if (c_b > 0.25*12) before =
true;
149 for (
unsigned ia = ibin + 1; ia <= ibin + 12; ++ia)
150 if (ia < beam_size && beam_z(ia) >=
Z_missing)
152 if (c_a > 0.25*12) after =
true;
160 for (
unsigned ib = start; ib <= end; ++ib)
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;
183 unsigned counter = 0;
184 unsigned counter_trash = 0;
185 unsigned counter_clutter =0;
186 for (
unsigned ibin = 0; ibin < beam_size; ++ibin)
188 bool is_clutter =
false;
189 bool is_trash =
false;
196 if( beam_sdzdr(ibin) <= 0.01 ){
200 if (beam_z (ibin) >= 45. ){
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.) ) {
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. )) ) {
219 if ((beam_w(ibin) * fabs(beam_v(ibin))) <= 0.25 && beam_z (ibin) !=
Z_missing ) {
224 if( is_clutter) counter_clutter ++;
225 if( is_trash ) counter_trash ++;
236 if ( is_clutter || is_trash )
247 if ( ! (is_clutter || is_trash ) || ibin == (beam_size - 1))
249 in_a_segment =
false;
251 if (ibin == (beam_size - 1)) end = ibin;
253 segment_length = end - start+1;
254 counter = counter + (unsigned)(segment_length);
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) ) )
263 if (
double(count)/double(min(
int(ibin),
int(2*
min_segment_length))) >=0.25) before =
true;
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)) ))
270 if (
double(count)/double(min(
int(beam_size - ibin),
int(2*
min_segment_length))) >=0.25) after =
true;
278 for (
unsigned ib = start; ib <= end; ++ib)
296 const unsigned beam_size = beam_z.rows();
297 vector<bool> res(beam_size,
false);
298 bool in_a_segment =
false;
300 unsigned segment_length;
302 unsigned counter = 0;
304 for (
unsigned ibin = 0; ibin < beam_size; ++ibin)
327 in_a_segment =
false;
329 if (ibin == (beam_size - 1)) end = ibin;
331 segment_length = end - start+1;
332 counter = counter + (unsigned)(segment_length);
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) ) )
341 if (
double(count)/double(min(
int(ibin),
int(2*
min_segment_length))) >=0.25) before =
true;
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)) ))
348 if (
double(count)/double(min(
int(beam_size - ibin),
int(2*
min_segment_length))) >=0.25) after =
true;
356 for (
unsigned ib = start; ib <= end; ++ib)
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
372 const unsigned beam_size = beam_z.rows();
373 vector<unsigned char> res(beam_size, 0);
374 vector<double> diffs(beam_size, 0);
384 string f_dir = fuzzy_path;
388 string fin_w = f_dir+
"/matrix-"+radar+
".txt";
390 vector<string> w_vector;
391 w_vector = read_matrix_from_txt(fin_w);
392 Num_entries = w_vector.size()/Num_echoes;
395 for(
int i=0;i<Num_echoes;i++){
396 for(
int j=0;j<Num_entries;j++){
397 Wij(i,j) = stod( w_vector[i*Num_entries+j]);
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]);
422 vector<unsigned> counter (Num_entries,0) ;
423 for (
unsigned ibin = 0; ibin < beam_size; ++ibin)
430 ArrayXd Class_WP(Num_echoes);
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]);
451 double prob_v =
trap (v_ny, v_ny, v_ny, v_ny, beam_v(ibin));
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]);
474 Pij(0,2)=
trap(Traps[1][0][0],Traps[1][0][1],Traps[1][0][2],Traps[1][0][3],beam_w(ibin));
475 Pij(1,2)=
trap(Traps[1][1][0],Traps[1][1][1],Traps[1][1][2],Traps[1][1][3],beam_w(ibin));
476 double prob_w =
trap(Traps[1][2][0],Traps[1][2][1],Traps[1][2][2],Traps[1][2][3],beam_w(ibin));
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]);
483 Pij(0,3) =
trap (Traps[3][0][0],Traps[3][0][1],Traps[3][0][2],Traps[3][0][3], beam_sd(ibin));
484 Pij(0,4) =
trap (Traps[4][0][0],Traps[4][0][1],Traps[4][0][2],Traps[4][0][3], beam_sdray(ibin));
485 Pij(0,5) =
trap (Traps[5][0][0],Traps[5][0][1],Traps[5][0][2],Traps[5][0][3], beam_sdaz(ibin));
487 Pij(0,6) =
trap(Traps[6][0][0],Traps[6][0][1],Traps[6][0][2],Traps[6][0][3],beam_zdr(ibin));
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]);
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));
490 Pij(0,9) =
trap(Traps[9][0][0],Traps[9][0][1],Traps[9][0][2],Traps[9][0][3], beam_sqi(ibin));
491 Pij(0,10) =
trap(Traps[10][0][0],Traps[10][0][1],Traps[10][0][2],Traps[10][0][3], beam_snr(ibin));
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));
495 Pij(0,11) =
trap(Traps[11][0][0],Traps[11][0][1],Traps[11][0][2],Traps[11][0][3],beam_zvd(ibin));
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]);
499 Pij(1,3) =
trap (Traps[3][1][0],Traps[3][1][1],Traps[3][1][2],Traps[3][1][3], beam_sd(ibin));
500 Pij(1,4) =
trap (Traps[4][1][0],Traps[4][1][1],Traps[4][1][2],Traps[4][1][3], beam_sdray(ibin));
501 Pij(1,5) =
trap (Traps[5][1][0],Traps[5][1][1],Traps[5][1][2],Traps[5][1][3], beam_sdaz(ibin));
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] );
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]);
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));
506 Pij(1,9) =
trap(Traps[9][1][0],Traps[9][1][1],Traps[9][1][2],Traps[9][1][3],beam_sqi(ibin));
507 Pij(1,10) =
trap(Traps[10][1][0],Traps[10][1][1],Traps[10][1][2],Traps[10][1][3], beam_snr(ibin));
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)));
511 Pij(1,11) =
trap(Traps[11][1][0],Traps[11][1][1],Traps[11][1][2],Traps[11][1][3],beam_zvd(ibin));
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]);
515 Pij(2,3) =
trap (Traps[3][2][0],Traps[3][2][1],Traps[3][2][2],Traps[3][2][3], beam_sd(ibin));
516 Pij(2,4) =
trap (Traps[4][2][0],Traps[4][2][1],Traps[4][2][2],Traps[4][2][3], beam_sdray(ibin));
517 Pij(2,5) =
trap (Traps[5][2][0],Traps[5][2][1],Traps[5][2][2],Traps[5][2][3], beam_sdaz(ibin));
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]);
520 Pij(2,7) =
trap(Traps[7][2][0],Traps[7][2][1],Traps[7][2][2],Traps[7][2][3],beam_rohv(ibin));
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));
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]);
523 Pij(2,10) =
trap(Traps[10][2][0],Traps[10][2][1],Traps[10][2][2],Traps[10][2][3], beam_snr(ibin));
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)));
527 Pij(2,11) =
trap(Traps[11][1][0],Traps[11][1][1],Traps[11][1][2],Traps[11][1][3],beam_zvd(ibin));
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]);
530 Pij(3,3) =
trap (Traps[3][3][0],Traps[3][3][1],Traps[3][3][2],Traps[3][3][3], beam_sd(ibin));
531 Pij(3,4) =
trap (Traps[4][3][0],Traps[4][3][1],Traps[4][3][2],Traps[4][3][3], beam_sdray(ibin));
532 Pij(3,5) =
trap (Traps[5][3][0],Traps[5][3][1],Traps[5][3][2],Traps[5][3][3], beam_sdaz(ibin));
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]);
535 Pij(3,7) =
trap(Traps[7][3][0],Traps[7][3][1],Traps[7][3][2],Traps[7][3][3],beam_rohv(ibin));
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));
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]);
538 Pij(3,10) =
trap(Traps[10][3][0],Traps[10][3][1],Traps[10][3][2],Traps[10][3][3], beam_snr(ibin));
539 Pij(3,11) =
trap(Traps[11][3][0],Traps[11][3][1],Traps[11][3][2],Traps[11][3][3],beam_zvd(ibin));
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]);
542 Pij(4,3) =
trap (Traps[3][4][0],Traps[3][4][1],Traps[3][4][2],Traps[3][4][3], beam_sd(ibin));
543 Pij(4,4) =
trap (Traps[4][4][0],Traps[4][4][1],Traps[4][4][2],Traps[4][4][3], beam_sdray(ibin));
544 Pij(4,5) =
trap (Traps[5][4][0],Traps[5][4][1],Traps[5][4][2],Traps[5][4][3], beam_sdaz(ibin));
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]);
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]);
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));
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]);
550 Pij(4,10) =
trap(Traps[10][4][0],Traps[10][4][1],Traps[10][4][2],Traps[10][4][3], beam_snr(ibin));
551 Pij(4,11) =
trap(Traps[11][4][0],Traps[11][4][1],Traps[11][4][2],Traps[11][4][3],beam_zvd(ibin));
555 Class_WP = ((Wij.array()*Pij.array()).matrix()*VectorXd::Ones(Num_entries)).array()/(Wij*VectorXd::Ones(Num_entries)).array();
557 Class_WP.maxCoeff(&i);
559 if (Class_WP(i) < 0.1 ) ID=5;
561 double diff = Class_WP[i]-Class_WP[0];
566 printf(
"bin %4d ",ibin);
567 printf(
"ID %d -- diff=%8.3f\n",ID, diff);
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));
584 if((force_meteo)&&(diff<=8.0)&&(prevID==0)&&(ID!=1)) res[ibin] = 0;
591 else if(radar==
"SPC"){
592 if((force_meteo)&&(diff<=10.)&&(prevID==0)) res[ibin] = 0;
594 if((force_meteo)&&(diff<=20.0)&&(prevID==0)&&(ID==4)) res[ibin] = 0;
595 }
else cout<<
"radar non specificato"<<endl;
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
610 const unsigned beam_size = beam_z.rows();
611 vector<unsigned char> res(beam_size, 0);
617 string f_dir = fuzzy_path;
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;
626 for(
int i=0;i<Num_echoes;i++){
627 for(
int j=0;j<Num_entries;j++){
628 Wij(i,j) = stod( w_vector[i*Num_entries+j]);
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]);
649 vector<unsigned> counter (Num_entries,0) ;
650 for (
unsigned ibin = 0; ibin < beam_size; ++ibin)
657 ArrayXd Class_WP(Num_echoes);
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));
670 double prob_v =
trap (v_ny, v_ny, v_ny, v_ny, beam_v(ibin),v_ny);
673 Pij(1,1)=
trap(Traps[0][1][0],Traps[0][1][1],Traps[0][1][2],Traps[0][1][3], beam_v(ibin));
686 Pij(0,2)=
trap(Traps[1][0][0],Traps[1][0][1],Traps[1][0][2],Traps[1][0][3],beam_w(ibin));
687 Pij(1,2)=
trap(Traps[1][1][0],Traps[1][1][1],Traps[1][1][2],Traps[1][1][3],beam_w(ibin));
688 double prob_w =
trap(Traps[1][2][0],Traps[1][2][1],Traps[1][2][2],Traps[1][2][3],beam_w(ibin));
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]);
695 Pij(0,3) =
trap (Traps[3][0][0],Traps[3][0][1],Traps[3][0][2],Traps[3][0][3], beam_sd(ibin));
696 Pij(0,4) =
trap (Traps[4][0][0],Traps[4][0][1],Traps[4][0][2],Traps[4][0][3], beam_sdray(ibin));
697 Pij(0,5) =
trap (Traps[5][0][0],Traps[5][0][1],Traps[5][0][2],Traps[5][0][3], beam_sdaz(ibin));
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]);
701 Pij(1,3) =
trap (Traps[3][1][0],Traps[3][1][1],Traps[3][1][2],Traps[3][1][3], beam_sd(ibin));
702 Pij(1,4) =
trap (Traps[4][1][0],Traps[4][1][1],Traps[4][1][2],Traps[4][1][3], beam_sdray(ibin));
703 Pij(1,5) =
trap (Traps[5][1][0],Traps[5][1][1],Traps[5][1][2],Traps[5][1][3], beam_sdaz(ibin));
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]);
706 Pij(2,3) =
trap (Traps[3][2][0],Traps[3][2][1],Traps[3][2][2],Traps[3][2][3], beam_sd(ibin));
707 Pij(2,4) =
trap (Traps[4][2][0],Traps[4][2][1],Traps[4][2][2],Traps[4][2][3], beam_sdray(ibin));
708 Pij(2,5) =
trap (Traps[5][2][0],Traps[5][2][1],Traps[5][2][2],Traps[5][2][3], beam_sdaz(ibin));
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]);
711 Pij(3,3) =
trap (Traps[3][3][0],Traps[3][3][1],Traps[3][3][2],Traps[3][3][3], beam_sd(ibin));
712 Pij(3,4) =
trap (Traps[4][3][0],Traps[4][3][1],Traps[4][3][2],Traps[4][3][3], beam_sdray(ibin));
713 Pij(3,5) =
trap (Traps[5][3][0],Traps[5][3][1],Traps[5][3][2],Traps[5][3][3], beam_sdaz(ibin));
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]);
716 Pij(4,3) =
trap (Traps[3][4][0],Traps[3][4][1],Traps[3][4][2],Traps[3][4][3], beam_sd(ibin));
717 Pij(4,4) =
trap (Traps[4][4][0],Traps[4][4][1],Traps[4][4][2],Traps[4][4][3], beam_sdray(ibin));
718 Pij(4,5) =
trap (Traps[5][4][0],Traps[5][4][1],Traps[5][4][2],Traps[5][4][3], beam_sdaz(ibin));
723 Class_WP = ((Wij.array()*Pij.array()).matrix()*VectorXd::Ones(Num_entries)).array()/(Wij*VectorXd::Ones(Num_entries)).array();
725 Class_WP.maxCoeff(&i);
727 if (Class_WP(i) < 0.1 ) ID=5;
739std::vector<unsigned char>
Cleaner::eval_clean_beam(
const Eigen::VectorXd& beam_z,
const Eigen::VectorXd& beam_w,
const Eigen::VectorXd& beam_v,
const Eigen::VectorXd& beam_sd,
const Eigen::VectorXd& beam_sdray,
const Eigen::VectorXd& beam_sdaz,
int iray)
const
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;
751 for (
unsigned ibin = 0; ibin < beam_size; ++ibin)
755 beam_sd(ibin) <= 10. ) || beam_z(ibin) ==
Z_missing) {
763 if (beam_sd(ibin) >= 5. && beam_sdray(ibin) >= 2 && beam_sdaz(ibin) > 4.){
768 if (beam_sd(ibin) >= 1. && beam_sd(ibin) <= 5. && beam_sdray(ibin) < 2. && beam_sdaz(ibin) > 4.){
773 if (beam_sd(ibin) >= 3. && beam_sd(ibin) <= 7. && beam_sdray(ibin) < 2. && beam_sdaz(ibin) > 4.&& beam_sdaz(ibin) <7 ){
779 if (beam_sd(ibin) >= 5. && beam_sdray(ibin) < 3. && beam_sdaz(ibin) < 3.){
784 if (beam_sd(ibin) >= 10. && beam_sdray(ibin) < 2. && beam_sdaz(ibin) < 5. && beam_z(ibin) <10.){
803 return evaluateCleanID(scan_z, scan_w, scan_v, scan_cleanID, scan_v.
undetect,iel);
811 throw std::runtime_error(
"scan_z beam_count no equal to scan_w beam_count");
813 throw std::runtime_error(
"scan_z beam_size no equal to scan_w beam_size");
816 throw std::runtime_error(
"scan_z beam_count no equal to scan_v beam_count");
818 throw std::runtime_error(
"scan_z beam_size no equal to scan_v beam_size");
822 const unsigned beam_count = scan_z.
beam_count;
823 const unsigned beam_size = scan_z.
beam_size;
833 Z_S.push_back(scan_z);
834 radarelab::volume::textureSD( Z_S,SD2D, 1000. , 3,
false);
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);
839 for (
unsigned i = 0; i < beam_count; ++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);
844 for (
unsigned ib = 0; ib < beam_size; ++ib)
845 scan_cleanID(i,ib)=corrected[ib];
849void Cleaner::evaluateClassID(
PolarScan<double>& scan_z,
PolarScan<double>& scan_w,
PolarScan<double>& scan_v,
PolarScan<double>& scan_zdr,
PolarScan<double>& scan_rohv,
PolarScan<double>& scan_sqi,
PolarScan<double>& scan_snr,
PolarScan<double>& scan_zvd,
PolarScan<unsigned char>& scan_cleanID,
PolarScan<double>& scan_DiffProb,
double bin_wind_magic_number,
const string radar,
const char* fuzzy_path,
unsigned iel,
bool force_meteo)
853 throw std::runtime_error(
"scan_z beam_count no equal to scan_w beam_count");
855 throw std::runtime_error(
"scan_z beam_size no equal to scan_w beam_size");
858 throw std::runtime_error(
"scan_z beam_count no equal to scan_v beam_count");
860 throw std::runtime_error(
"scan_z beam_size no equal to scan_v beam_size");
864 const unsigned beam_count = scan_z.
beam_count;
865 const unsigned beam_size = scan_z.
beam_size;
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);
877 radarelab::volume::textureSD( ZDR_S,ZDR_SD2D, 1000. , 3,
false);
880 for (
unsigned i = 0; i <beam_count ; ++i)
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);
887 for (
unsigned ib = 0; ib < beam_size; ++ib){
888 scan_cleanID(i,ib)=corrected[ib];
889 scan_DiffProb(i,ib)=diff_prob[ib];
901 throw std::runtime_error(
"scan_z beam_count no equal to scan_w beam_count");
903 throw std::runtime_error(
"scan_z beam_size no equal to scan_w beam_size");
906 throw std::runtime_error(
"scan_z beam_count no equal to scan_v beam_count");
908 throw std::runtime_error(
"scan_z beam_size no equal to scan_v beam_size");
912 const unsigned beam_count = scan_z.
beam_count;
913 const unsigned beam_size = scan_z.
beam_size;
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);
923 for (
unsigned i = 0; i <beam_count ; ++i)
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];
933 return clean(scan_z, scan_w, scan_v, scan_v.
undetect,iel);
939 throw std::runtime_error(
"scan_z beam_count no equal to scan_w beam_count");
941 throw std::runtime_error(
"scan_z beam_size no equal to scan_w beam_size");
944 throw std::runtime_error(
"scan_z beam_count no equal to scan_v beam_count");
946 throw std::runtime_error(
"scan_z beam_size no equal to scan_v beam_size");
950 const unsigned beam_count = scan_z.
beam_count;
951 const unsigned beam_size = scan_z.
beam_size;
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);
991 if (set_undetect) new_val=scan_z.
nodata;
993 for (
unsigned i = 0; i < beam_count; ++i)
997 vector<bool> corrected = cleaner.
clean_beam(scan_z.row(i), scan_w.row(i), scan_v.row(i),i);
999 for (
unsigned ib = 0; ib < beam_size; ++ib)
1003 scan_z(i, ib) = new_val;
1020 return clean(scan_z, scan_w, scan_v, scan_zdr, scan_v.
undetect,iel,set_undetect);
1024 cout<<
"Chiamato cleaner "<<set_undetect<<endl;
1026 throw std::runtime_error(
"scan_z beam_count no equal to scan_w beam_count");
1028 throw std::runtime_error(
"scan_z beam_size no equal to scan_w beam_size");
1031 throw std::runtime_error(
"scan_z beam_count no equal to scan_v beam_count");
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;
1039 const unsigned beam_count = scan_z.
beam_count;
1040 const unsigned beam_size = scan_z.
beam_size;
1046 Z_S.push_back(scan_z);
1047 radarelab::volume::textureSD( Z_S,SD2D, 1000. , 3,
false);
1049 ZDR_S.push_back(scan_zdr);
1050 radarelab::volume::textureSD( ZDR_S,SDZDR2D, 1000. , 3,
false);
1059Matrix2D <unsigned char>img_tmp, z_clean;
1060Matrix2D <double>img;
1066 printf(
"scrivo Z ");
1067 img = (scan_z.array() - scan_z.
offset )/ scan_z.
gain /256 ;
1068 sprintf(pippo,
"_%02d.png",iel);
1070 img_tmp=img.cast<
unsigned char>();
1072 radarelab::write_image(img_tmp,
"/ponte/rad_svn/proc_operative/test_arch/rev_actual/radar/immagini/Cleaner/PPI_Z"+ext,
"PNG");
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");
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");
1092 double new_val=cleaner.Z_missing;
1093 if (set_undetect) new_val=scan_z.
undetect;
1095 for (
unsigned i = 0; i < beam_count; ++i)
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);
1100 for (
unsigned ib = 0; ib < beam_size; ++ib)
1104 scan_z(i, ib) = new_val;
1119 std::string Z_Quantity;
1123double Cleaner::trap(
double x1,
double x2,
double x3,
double x4,
double val,
double x5)
const
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.;
1133vector<string> Cleaner::read_matrix_from_txt(
string fin)
const
1139 vector<string> myVector;
1140 ifstream f(fin, ifstream::in);
1144 while(getline(f,line)){
1145 stringstream stream (line);
1146 while( getline(stream, line,
' ')){
1148 if(line[0]==
'\\'){
continue;}
1150 myVector.push_back(line);
double nodata
Value used as 'no data' value.
double undetect
Minimum amount that can be measured.
double offset
Conversion factor.
double gain
Conversion factor.
PolarScan - structure to describe a polarScan containing a matrix of data and conversion factors.
Sequence of PolarScans which can have a different beam count for each elevation.
void gdal_init_once()
Initialize the GDAL library when called for the first time; does nothing all other times.
Base for all matrices we use, since we rely on row-major data.
unsigned beam_size
Number of samples in each beam.
unsigned beam_count
Count of beams in this scan.
double cell_size
Size of a beam cell in meters.
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,...
double trap(double x1, double x2, double x3, double x4, double val, double x5=-9999.) const
Cleaner(double Z_missing, double W_threshold, double V_missing, double bin_wind_magic_number)
Constructor.
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)
const double sd_threshold
Soglia per devizione standard DBZH.
const double Z_missing
Valore dato mancante DBZH.
const double bin_wind_magic_number
valore magico per dati in formato SP20
const unsigned max_segment_length
lunghezza massima segmento in celle se più lungo pulisce in ogni caso
const double W_threshold
Soglia per WRAD.
const unsigned min_segment_length
lunghezza minima segmento in celle
Struttura per cleaner dati grezzi sulla base dei valori di V, W e la deviazione standard di Z.
Struttura che contiene mappa per caricamento dati.