Elaboradar 0.1
Caricamento in corso...
Ricerca in corso...
Nessun risultato
steiner.cpp
1#include "steiner.h"
4
5//#ifdef __cplusplus
6//extern "C" {
7//#endif
8//#include <func_Z_R.h>
9//#ifdef __cplusplus
10//}
11//#endif
12
13// valore mancante
14#define MISSING 0
15
16// fattore conversione gradi-radianti
17#define DTOR M_PI/180.
18
19//Definizioni geometriche
20static const double AMPLITUDE = 0.9; /* esternalizzo?*/ // ampiezza fascio radar
21
22using namespace std;
23
24namespace radarelab {
25namespace algo {
26
27namespace steiner {
28
29void Point::add_sample(double sample)
30{
31 if (sample <= MINVAL_DB) return;
32 Z_bckgr += DBZ::BYTEtoZ(DBZ::DBtoBYTE(sample));
33 bckgr += sample;
34 ++npoints;
35}
36
37void Point::finalize()
38{
39 if (npoints > 0){
40 Z_bckgr = Z_bckgr / npoints;
41 //bckgr[i]=bckgr[i]/npoints; //no
42 if (Z_bckgr > 0) bckgr = 10 * (log10(Z_bckgr));
43 }
44 //il valore del raggio convettivo varia a seconda del background, da 1 a 5 km
45 if (bckgr < 25.)
46 convective_radius = 1.;
47 else if (bckgr >= 25. && bckgr < 30.)
48 convective_radius = 2.;
49 else if (bckgr >= 30. && bckgr < 35.)
50 convective_radius = 3.;
51 else if (bckgr >= 35. && bckgr < 40.)
52 convective_radius = 4.;
53 else if (bckgr > 40.)
54 convective_radius = 5.;
55}
56
57}
58
59CalcoloSteiner::CalcoloSteiner(
60 const Volume<double>& volume,
61 const volume::ElevFin<double>& elev_fin,
62 unsigned max_bin)
63 // TODO: evaluate cell_size scan by scan?
64 : volume(volume), elev_fin(elev_fin), max_bin(max_bin), size_cell(volume.scan(0).cell_size),
65 conv_STEINER(Matrix2D<unsigned char>::Constant(volume.beam_count, max_bin, MISSING))
66{
67 using namespace steiner;
68
69 logging_category = log4c_category_get("radar.vpr");
70
71 // Traccio una lista dei punti che hanno valore non nullo e sotto base
72 // bright band (lista_bckg) contenente iaz e irange e conto i punti
73 for (unsigned i=0; i < volume.beam_count; ++i)
74 for (unsigned j=0; j < max_bin; ++j) // propongo max_bin visto che risoluzione è la stessa
75 //if ( volume.scan(0)[i][j] > 1 && (float)(quota[i][j])/1000. < hbbb ) //verifico che il dato usato per la ZLR cioè la Z al lowest level sia > soglia e la sua quota sia sotto bright band o sopra bright band
76 if (j < volume.scan(0).beam_size && volume.scan(0).get(i, j) > MINVAL_DB)
77 lista_bckg.push_back(Point(i, j)); // IAZIMUT, IRANGE
78}
79
80void CalcoloSteiner::calcolo_background() // sui punti precipitanti calcolo bckgr . nb LA CLASSIFICAZIONE DI STEINER NON HA BISOGNO DI RICAMPIONAMENTO CILINDRICO PERCIÒ uso direttamente la matrice polare
81 // definisco una lista di punti utili per l'analisi, cioè con valore non nullo e sotto la bright band o sopra (NB. per l'analisi considero i punti usati per la ZLR che si suppongono non affetti da clutter e beam blocking<50%. questo ha tutta una serie di implicazioni..... tra cui la proiezione, il fatto che nel confronto entrano dei 'buchi'...cioè il confronto è fatto su una matrice pseudo-orizzontale bucata e quote variabili tuttavia, considerato che i gradienti verticali in caso convettivo e fuori dalla bright band non sono altissimi spero funzioni ( andro' a verificare che l'ipotesi sia soddisfatta)
82 // per trovare il background uso uno pseudo quadrato anzichè cerchio, mantenendo come semi-lato il raggio di steiner (11KM); devo perciò definire una semi-ampiezza delle finestre in azimut e range corrispondente al raggio di 11 km
83 // quella in azimut dipende dal range
84
85{
86 using namespace steiner;
87
88 if (lista_bckg.size() < 2)
89 return;
90
91 // Per il calcolo della finestra range su cui calcolare il background
92 // divido il raggio di Steiner (11km) per la dimensione della cella.
93 // Definisco ampiezza semi-finestra range corrispondente al raggio di
94 // steiner (11km), in unità matrice polare.
95 unsigned delta_nr = round(STEINER_RADIUS * 1000. / size_cell);
96 LOG_DEBUG("delta_n range per analisi Steiner = %u --- dimensione lista_bckg %d", delta_nr,lista_bckg.size());
97
98 for (vector<Point>::iterator i = lista_bckg.begin(); i != lista_bckg.end(); ++i)
99 {
100 // estremi della finestra in range
101 int kmin = i->range - delta_nr;
102 unsigned kmax = min(i->range + delta_nr, max_bin);
103
104 if (kmin>0)
105 {
106 //definisco ampiezza semi finestra nazimut corrispondente al raggio di steiner (11km) (11/distanzacentrocella)(ampiezzaangoloscansione)
107 unsigned delta_naz=ceil(STEINER_RADIUS/((i->range * size_cell/1000. + size_cell/2000.)/(AMPLITUDE*DTOR)));
108 if (delta_naz > volume.beam_count / 2)
109 delta_naz = volume.beam_count / 2;
110
111 int jmin = i->azimut - delta_naz;
112 int jmax = i->azimut + delta_naz;
113
114 for (int j = jmin; j < jmax; ++j)
115 for (unsigned k = kmin; k < kmax; ++k)
116 i->add_sample(elev_fin.db_at_elev_preci((j + volume.beam_count) % volume.beam_count, k));
117 } else {
118 // FIXME: questo fa mezzo scan tra 0 e kmax, e mezzo scan tra 0 e
119 // -kmin. Sempre gli stessi mezzi scan a prescindere dalla
120 // posizione di i. Ha senso?
121 for (unsigned j=0 ; j<volume.beam_count/2 ; j++)
122 for (unsigned k=0 ; k<kmax ; k++)
123 i->add_sample(elev_fin.db_at_elev_preci(j, k));
124
125 for (unsigned j= volume.beam_count/2 ; j<volume.beam_count ; j++)
126 for (int k=0 ; k<-kmin ; k++)
127 i->add_sample(elev_fin.db_at_elev_preci(j, k));
128 }
129 i->finalize();
130 }
131}
132
133void CalcoloSteiner::ingrasso_nuclei(float cr,int ja,int kr)
134{
135 int dr=(int)(cr*1000./size_cell);//definisco ampiezza semi-finestra range corrispondente al raggio di steiner (11km), unità matrice polare
136
137 int kmin=kr-dr;
138 unsigned kmax=kr+dr;
139
140 int daz=ceil(cr/((kr*size_cell/1000.+size_cell/2000.)/(AMPLITUDE*DTOR)));
141 int jmin=ja-daz;
142 unsigned jmax=ja+daz;
143
144 LOG_DEBUG("dr cr kmin kmax %d %f %d %d %d %d", dr,cr, kmin,kmax,jmin,jmax);
145
146 if (kmin>0) {
147 if (kmax > max_bin) kmax = max_bin;
148
149 if (jmin<0) {
150 jmin=volume.beam_count-jmin%volume.beam_count;
151 for (unsigned j=jmin; j< volume.beam_count ; j++) {
152 for (unsigned k=kmin ; k<kmax ; k++) {
153 conv_STEINER(j, k)=CONV_VAL;
154 }
155 }
156 LOG_DEBUG("jmin %d", jmin);
157 jmin=0;
158
159 }
160
161 if (jmax>=volume.beam_count) {
162 jmax=jmax%volume.beam_count;
163 for (unsigned j=0; j<jmax ; j++) {
164 for (unsigned k=kmin; k<kmax ; k++) {
165 conv_STEINER(j, k)=CONV_VAL;
166 }
167 }
168 LOG_DEBUG("jmax %d", jmax);
169 jmax=volume.beam_count;
170 }
171 for (unsigned j=jmin; j<jmax ; j++) {
172 for (unsigned k=kmin; k<kmax ; k++) {
173 conv_STEINER(j%volume.beam_count, k)=CONV_VAL;
174 }
175 }
176 }
177 else
178 {
179 for (unsigned j=0 ; j<volume.beam_count/2 ; j++)
180 for (unsigned k=0 ; k<kmax ; k++){
181 conv_STEINER(j, k)=CONV_VAL;
182 }
183 for (unsigned j= volume.beam_count/2 ; j<volume.beam_count ; j++)
184 for (unsigned k=0 ; k < (unsigned)-kmin ; k++){
185 conv_STEINER(j, k)=CONV_VAL;
186 }
187 }
188 LOG_DEBUG ("Finita ingrasso nuclei");
189 return;
190}
191
192void CalcoloSteiner::classifico_STEINER()
193{
194 using namespace steiner;
195
196 for (vector<Point>::const_iterator i = lista_bckg.begin(); i != lista_bckg.end(); ++i)
197 {
198 int j = i->azimut;
199 int k = i->range;
200 if (j < 0 || k < 0) continue;
201
202 double db = elev_fin.db_at_elev_preci(j, k);
203 // calcolo diff col background
204 float diff_bckgr = db - i->bckgr;
205 // test su differenza con bckground , se soddisfatto e simultaneamente il VIZ non ha dato class convettiva (?)
206 if ((db > 40.) ||
207 (i->bckgr < 0 && diff_bckgr > 10) ||
208 (i->bckgr < 42.43 && i->bckgr > 0 && diff_bckgr > 10. - i->bckgr * i->bckgr / 180.) ||
209 (i->bckgr > 42.43 && diff_bckgr > 0))
210 {
211 // assegno il punto nucleo di Steiner
212 conv_STEINER(j, k) = CONV_VAL;
213
214 // ingrasso il nucleo
215 float cr = i->convective_radius;
216 LOG_DEBUG(" %f cr", cr);
217 ingrasso_nuclei(cr, j, k);
218 }
219 }
220 return;
221}
222
223}
224}
static unsigned char DBtoBYTE(double DB, double gain=80./255., double offset=-20.)
funzione che converte dB in valore intero tra 0 e 255
Definition: dbz.h:94
static double BYTEtoZ(unsigned char byte)
funzione che converte byte in Z
Definition: dbz.cpp:125
String functions.
Definition: cart.cpp:4