27 using namespace radarelab;
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
34 const unsigned beam_size = beam_z.rows();
35 vector<bool> res(beam_size,
false);
36 bool in_a_segment =
false;
38 unsigned segment_length;
42 for (
unsigned ibin = 0; ibin < beam_size; ++ibin)
48 if (beam_w(ibin) == W_threshold && beam_v(ibin) == bin_wind_magic_number)
59 if (beam_w(ibin) != W_threshold || beam_v(ibin) != bin_wind_magic_number || ibin == (beam_size - 1))
63 if (ibin == (beam_size - 1)) end = ibin;
65 segment_length = end - start;
66 counter = counter + (unsigned)(segment_length);
71 for (
int ib = ibin - 12; ib < (signed)ibin; ++ib)
72 if (ib >= 0 && beam_z(ib) > Z_missing)
74 if (c_b > 0.25*12) before =
true;
77 for (
unsigned ia = ibin + 1; ia <= ibin + 12; ++ia)
78 if (ia < beam_size && beam_z(ia) >= Z_missing)
80 if (c_a > 0.25*12) after =
true;
83 if ((segment_length >= min_segment_length && !before && !after) ||
84 segment_length >= max_segment_length || counter > 100)
88 for (
unsigned ib = start; ib <= end; ++ib)
89 if( beam_z(ib) > Z_missing)res[ib] =
true;
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
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;
109 for (
unsigned ibin = 0; ibin < beam_size; ++ibin)
115 if (beam_w(ibin) == W_threshold && beam_v(ibin) == bin_wind_magic_number)
125 if (beam_w(ibin) != W_threshold || beam_v(ibin) != bin_wind_magic_number || ibin == (beam_size - 1))
127 in_a_segment =
false;
129 if (ibin == (beam_size - 1)) end = ibin;
131 segment_length = end - start;
132 counter = counter + (unsigned)(segment_length);
137 for (
int ib = ibin - 12; ib < (signed)ibin; ++ib)
138 if (ib >= 0 && beam_z(ib) > Z_missing)
140 if (c_b > 0.25*12) before =
true;
143 for (
unsigned ia = ibin + 1; ia <= ibin + 12; ++ia)
144 if (ia < beam_size && beam_z(ia) >= Z_missing)
146 if (c_a > 0.25*12) after =
true;
149 if ((segment_length >= min_segment_length && !before && !after) ||
150 segment_length >= max_segment_length || counter > 100)
154 for (
unsigned ib = start; ib <= end; ++ib)
155 if( beam_z(ib) > Z_missing)res[ib] = 1;
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;
177 unsigned counter = 0;
178 unsigned counter_trash = 0;
179 unsigned counter_clutter =0;
180 for (
unsigned ibin = 0; ibin < beam_size; ++ibin)
182 bool is_clutter =
false;
183 bool is_trash =
false;
188 if ( beam_w(ibin) == W_threshold && beam_v(ibin) == bin_wind_magic_number && beam_z (ibin) != Z_missing ) {
190 if( beam_sdzdr(ibin) <= 0.01 ){
194 if (beam_z (ibin) >= 45. ){
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.) ) {
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. )) ) {
213 if ((beam_w(ibin) * fabs(beam_v(ibin))) <= 0.25 && beam_z (ibin) != Z_missing ) {
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));
230 if ( is_clutter || is_trash )
241 if ( ! (is_clutter || is_trash ) || ibin == (beam_size - 1))
243 in_a_segment =
false;
245 if (ibin == (beam_size - 1)) end = ibin;
247 segment_length = end - start+1;
248 counter = counter + (unsigned)(segment_length);
251 if (segment_length <= 2*min_segment_length ){
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) ) )
257 if (
double(count)/double(min(
int(ibin),
int(2*min_segment_length))) >=0.25) before =
true;
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)) ))
264 if (
double(count)/double(min(
int(beam_size - ibin),
int(2*min_segment_length))) >=0.25) after =
true;
267 if ((segment_length >= min_segment_length && (!before || !after) ) || segment_length >= max_segment_length)
272 for (
unsigned ib = start; ib <= end; ++ib)
290 const unsigned beam_size = beam_z.rows();
291 vector<bool> res(beam_size,
false);
292 bool in_a_segment =
false;
294 unsigned segment_length;
296 unsigned counter = 0;
298 for (
unsigned ibin = 0; ibin < beam_size; ++ibin)
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 )
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 )
321 in_a_segment =
false;
323 if (ibin == (beam_size - 1)) end = ibin;
325 segment_length = end - start+1;
326 counter = counter + (unsigned)(segment_length);
329 if (segment_length <= 2*min_segment_length ){
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) ) )
335 if (
double(count)/double(min(
int(ibin),
int(2*min_segment_length))) >=0.25) before =
true;
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)) ))
342 if (
double(count)/double(min(
int(beam_size - ibin),
int(2*min_segment_length))) >=0.25) after =
true;
345 if ((segment_length >= min_segment_length && (!before || !after) ) || segment_length >= max_segment_length)
350 for (
unsigned ib = start; ib <= end; ++ib)
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
367 const unsigned beam_size = beam_z.rows();
368 vector<unsigned char> res(beam_size, 0);
369 vector<unsigned> counter (6,0) ;
373 Wij << 1.0, 1.0, 1.0, 0.5, 0.3, 0.3,
374 0.8, 0.5, 0.5, 0.5, 0.5, 0.5,
375 0.0, 0.2, 0.2, 0.6, 0.8, 0.5,
376 0.0, 0.2, 0.2, 0.5, 0.6, 0.4,
377 0.0, 0.2, 0.2, 0.6, 0.5, 0.5,
378 0.6, 0.2, 0.2, 0.6, 0.5, 0.5;
383 for (
unsigned ibin = 0; ibin < beam_size; ++ibin)
390 if (beam_z(ibin) == Z_missing) {
400 if (beam_v(ibin) != bin_wind_magic_number ) {
402 double prob_v = trap (-0.5, -0.3, 0.3, 0.5, beam_v(ibin));
417 if (beam_w(ibin) > W_threshold) {
419 double prob_w = trap (0.0001, 0., 0.2, 0.3, beam_w(ibin));
436 Pij(0,3) = trap (0.5, 1., 10.,11., beam_sd(ibin));
437 Pij(0,4) = trap (0.5, 1., 10.,11., beam_sdray(ibin));
438 Pij(0,5) = trap (0.5, 1., 10.,11., beam_sdaz(ibin));
441 Pij(1,0) = trap (5., 15., 99., 99.9, beam_z(ibin));
442 Pij(1,3) = trap (4.5, 5., 99., 99.9, beam_sd(ibin));
443 Pij(1,4) = trap (4., 4.5, 99., 99.9, beam_sdray(ibin));
444 Pij(1,5) = trap (4., 4.5, 99., 99.9, beam_sdaz(ibin));
447 Pij(2,3) = trap (0.5, 1.5, 5., 7., beam_sd(ibin));
448 Pij(2,4) = trap (0.0, 0.3, 2., 3., beam_sdray(ibin));
449 Pij(2,5) = trap (2.0, 4., 90.,90.9, beam_sdaz(ibin));
453 Pij(3,3) = trap (4.0, 5., 90., 90.9, beam_sd(ibin));
454 Pij(3,4) = trap (0., 0.5, 10., 15.0, beam_sdray(ibin));
455 Pij(3,5) = trap (3., 4., 90., 90.9, beam_sdaz(ibin));
459 Pij(4,3) = trap (4., 7., 90., 90.9, beam_sd(ibin));
460 Pij(4,4) = trap (0.5, 1.5, 3., 4., beam_sdray(ibin));
461 Pij(4,5) = trap (3., 4., 90., 90.9, beam_sdaz(ibin));
466 if ( ibin <= 160) coeff = (160 - ibin)/ 120. ;
469 Pij(5,0) = trap(Z_missing-0.0001, Z_missing, 15., 20., beam_z(ibin));
470 Pij(5,3) = trap (5., 7., 99., 99.9, beam_sd(ibin));
471 Pij(5,4) = trap (0., 0.001, 0.8 + coeff * 9.2 , 1. + coeff * 14.,beam_sdray(ibin));
474 Pij(5,5) = trap (1.5, 3., 99., 99.9,beam_sdaz(ibin));
479 Class_WP = ((Wij.array()*Pij.array()).matrix()*VectorXd::Ones(6)).array()/(Wij*VectorXd::Ones(6)).array();
481 Class_WP.maxCoeff(&i);
483 if (Class_WP(i) < 0.1 ) ID=6;
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
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;
509 for (
unsigned ibin = 0; ibin < beam_size; ++ibin)
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) {
519 if (beam_w(ibin) == W_threshold && beam_v(ibin) == bin_wind_magic_number)
521 if (beam_sd(ibin) >= 5. && beam_sdray(ibin) >= 2 && beam_sdaz(ibin) > 4.){
526 if (beam_sd(ibin) >= 1. && beam_sd(ibin) <= 5. && beam_sdray(ibin) < 2. && beam_sdaz(ibin) > 4.){
531 if (beam_sd(ibin) >= 3. && beam_sd(ibin) <= 7. && beam_sdray(ibin) < 2. && beam_sdaz(ibin) > 4.&& beam_sdaz(ibin) <7 ){
537 if (beam_sd(ibin) >= 5. && beam_sdray(ibin) < 3. && beam_sdaz(ibin) < 3.){
542 if (beam_sd(ibin) >= 10. && beam_sdray(ibin) < 2. && beam_sdaz(ibin) < 5. && beam_z(ibin) <10.){
549 printf(
"%2d %4d %4d %4d %4d %4d\n",res[ibin],countClutter,countIntStrong, countIntMed, countIntWeak, countNoise);
562 return evaluateClassID(scan_z, scan_w, scan_v, scan_cleanID, scan_v.
undetect,iel);
569 throw std::runtime_error(
"scan_z beam_count no equal to scan_w beam_count");
571 throw std::runtime_error(
"scan_z beam_size no equal to scan_w beam_size");
574 throw std::runtime_error(
"scan_z beam_count no equal to scan_v beam_count");
576 throw std::runtime_error(
"scan_z beam_size no equal to scan_v beam_size");
580 const unsigned beam_count = scan_z.
beam_count;
581 const unsigned beam_size = scan_z.
beam_size;
591 Z_S.push_back(scan_z);
592 radarelab::volume::textureSD( Z_S,SD2D, 1000. , 3,
false);
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);
597 for (
unsigned i = 0; i < beam_count; ++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);
602 for (
unsigned ib = 0; ib < beam_size; ++ib)
603 scan_cleanID(i,ib)=corrected[ib];
611 throw std::runtime_error(
"scan_z beam_count no equal to scan_w beam_count");
613 throw std::runtime_error(
"scan_z beam_size no equal to scan_w beam_size");
616 throw std::runtime_error(
"scan_z beam_count no equal to scan_v beam_count");
618 throw std::runtime_error(
"scan_z beam_size no equal to scan_v beam_size");
622 const unsigned beam_count = scan_z.
beam_count;
623 const unsigned beam_size = scan_z.
beam_size;
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);
633 for (
unsigned i = 0; i < beam_count; ++i)
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];
643 return clean(scan_z, scan_w, scan_v, scan_v.
undetect,iel);
649 throw std::runtime_error(
"scan_z beam_count no equal to scan_w beam_count");
651 throw std::runtime_error(
"scan_z beam_size no equal to scan_w beam_size");
654 throw std::runtime_error(
"scan_z beam_count no equal to scan_v beam_count");
656 throw std::runtime_error(
"scan_z beam_size no equal to scan_v beam_size");
660 const unsigned beam_count = scan_z.
beam_count;
661 const unsigned beam_size = scan_z.
beam_size;
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);
701 if (set_undetect) new_val=scan_z.
nodata;
703 for (
unsigned i = 0; i < beam_count; ++i)
707 vector<bool> corrected = cleaner.
clean_beam(scan_z.row(i), scan_w.row(i), scan_v.row(i),i);
709 for (
unsigned ib = 0; ib < beam_size; ++ib)
713 scan_z(i, ib) = new_val;
730 return clean(scan_z, scan_w, scan_v, scan_zdr, scan_v.
undetect,iel);
735 throw std::runtime_error(
"scan_z beam_count no equal to scan_w beam_count");
737 throw std::runtime_error(
"scan_z beam_size no equal to scan_w beam_size");
740 throw std::runtime_error(
"scan_z beam_count no equal to scan_v beam_count");
742 throw std::runtime_error(
"scan_z beam_size no equal to scan_v beam_size");
746 const unsigned beam_count = scan_z.
beam_count;
747 const unsigned beam_size = scan_z.
beam_size;
753 Z_S.push_back(scan_z);
754 radarelab::volume::textureSD( Z_S,SD2D, 1000. , 3,
false);
756 ZDR_S.push_back(scan_zdr);
757 radarelab::volume::textureSD( ZDR_S,SDZDR2D, 1000. , 3,
false);
774 img = (scan_z.array() - scan_z.
offset )/ scan_z.
gain /256 ;
775 sprintf(pippo,
"_%02d.png",iel);
777 img_tmp=img.cast<
unsigned char>();
779 radarelab::write_image(img_tmp,
"/ponte/rad_svn/proc_operative/test_arch/rev_actual/radar/immagini/Cleaner/PPI_Z"+ext,
"PNG");
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");
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");
800 if (set_undetect) new_val=scan_z.
undetect;
802 for (
unsigned i = 0; i < beam_count; ++i)
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);
807 for (
unsigned ib = 0; ib < beam_size; ++ib)
811 scan_z(i, ib) = new_val;
826 std::string Z_Quantity;
831 double Cleaner::trap(
double x1,
double x2,
double x3,
double x4,
double val)
const
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);
double gain
Conversion factor.
double undetect
Minimum amount that can be measured.
Struttura per cleaner dati grezzi sulla base dei valori di V, W e la deviazione standard di Z...
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)
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...
PolarScan - structure to describe a polarScan containing a matrix of data and conversion factors...
double cell_size
Size of a beam cell in meters.
unsigned beam_count
Count of beams in this scan.
double trap(double x1, double x2, double x3, double x4, double val) const
double offset
Conversion factor.
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'oggetto cleaner, lo inizializza, pulisce i dati e modifica il PolarScan di DBZH...
unsigned beam_size
Number of samples in each beam.
Struttura che contiene mappa per caricamento dati.
const double Z_missing
Valore dato mancante DBZH.
void gdal_init_once()
Initialize the GDAL library when called for the first time; does nothing all other times...
double nodata
Value used as 'no data' value.