Elaboradar 0.1
Caricamento in corso...
Ricerca in corso...
Nessun risultato
cylindrical.cpp
2#include <cmath>
3
4#define DTOR M_PI/180. /* esternalizzo?*/ //fattore conversione gradi-radianti
5
6using namespace std;
7
8namespace radarelab {
9
10CylindricalVolume::CylindricalVolume(const Volume<double>& volume, double missing_value, double x_res, double z_res)
11{
12 const double size_cell = volume[0].cell_size;
13 double range_min=0.5 * size_cell/1000.;
14 double range_maxLowestRay = (volume[0].beam_size - 0.5) * size_cell / 1000.;
15
16 double xmin = floor(range_min*cos(volume.elevation_max()*DTOR)); // distanza orizzontale minima dal radar
17 double zmin = volume[0].sample_height(0) / 1000. + volume.radarSite.getTotalHeight(); // quota minima in prop standard
18 double xmax = floor(range_maxLowestRay*cos(volume.elevation_min()*DTOR)); // distanza orizzontale massima dal radar
19 double zmax = volume.back().sample_height(volume.back().beam_size - 1) / 1000. + volume.radarSite.getTotalHeight();//quota massima
20 //LOG_DEBUG(" Range min maxL maxU %7.3f %7.3f %7.3f -- xmin %7.3f xmax %7.3f zmin %7.3f zmax %7.3f", range_min, range_maxLowestRay, range_maxUpperRay, xmin,xmax,zmin,zmax);
21
22 x_size = (xmax - xmin) / x_res; //dimensione orizzontale
23 // FIXME: usiamo volume.max_beam_size invece di MyMAX_BIN?
24 if (x_size > volume.max_beam_size()) x_size=volume.max_beam_size();
25 z_size = (zmax - zmin) / z_res; //dimensione verticale
26
27 resol[0] = x_res;
28 resol[1] = z_res;
29 slices.reserve(volume.beam_count);
30 for (unsigned i = 0; i < volume.beam_count; ++i)
31 slices.push_back(new Matrix2D<double>(Matrix2D<double>::Constant(x_size, z_size, missing_value)));
32
33 resample(volume);
34}
35
36void CylindricalVolume::resample(const Volume<double>& volume)
37{
38 unsigned max_bin = x_size;
39 /* ---------------------------------- */
40 /* FASE 1 */
41 /* ---------------------------------- */
42 /* Costruzione matrice generale per indicizzare il caricamento dei dati */
43 /* da coordinate radar a coordinate (X,Z) */
44 /* Calcolo distanza per ogni punto sul raggio */
45 /* xx contiene la distanza sulla superfice terrestre in funzione del range radar e dell'elevazione */
46 /* Metodo 1 - -Calcolo le coordinate di ogni punto del RHI mediante utilizzo di cicli */
47
48 // estremi x e z (si procede per rhi)
49#warning to compute scan by scan?
50 double cell_size = volume.scan(0).cell_size;
51 double range_min=0.5 * cell_size / 1000.;
52 //double range_max=(max_bin-0.5) * size_cell/1000.;
53
54 double xmin=floor(range_min*cos(volume.elevation_max()*DTOR)); // distanza orizzontale minima dal radar
55 double zmin=volume[0].sample_height(0) / 1000. + volume.radarSite.getTotalHeight(); // quota minima in prop standard
56
57
58 //float w_size[2]={3.,1.5}; //dimensione della matrice pesi
59 const double w_size[2]={3.,0.3}; //dimensione della matrice pesi
60
61 //LOG_INFO("calcolati range_min e range_max , dimensione orizzontale e dimensione verticale range_min=%f range_max=%f x_size=%d z_size=%d",range_min,range_max,x_size,z_size);
62
63 int w_x_size=ceil((w_size[0]/resol[0])/2)*2+1; //dimensione x matrice pesi
64 int w_z_size=ceil((w_size[1]/resol[1])/2)*2+1; //dimensione z matrice pesi
65
66 if (w_x_size < 3) w_x_size=3;
67 if (w_z_size < 3) w_z_size=3;
68
69 int w_x_size_2=w_x_size/2;
70 int w_z_size_2=w_z_size/2;
71
72 Matrix2D<int> i_xx_min(Matrix2D<int>::Zero(max_bin, volume.size()));
73 Matrix2D<int> i_zz_min(Matrix2D<int>::Zero(max_bin, volume.size()));
74 Matrix2D<int> im(Matrix2D<int>::Zero(max_bin, volume.size()));
75 Matrix2D<int> ix(Matrix2D<int>::Zero(max_bin, volume.size()));
76 Matrix2D<int> jm(Matrix2D<int>::Zero(max_bin, volume.size()));
77 Matrix2D<int> jx(Matrix2D<int>::Zero(max_bin, volume.size()));
78 for (unsigned i = 0; i < max_bin; i++){
79 double range = (i + 0.5) * cell_size/1000.;
80
81 for (unsigned k=0; k < volume.size(); k++){
82 double elev_rad = volume.scan(k).elevation * DTOR;
83 double zz = volume.scan(k).sample_height(i) / 1000. + volume.radarSite.getTotalHeight();// quota
84 double xx = range*cos(elev_rad); // distanza
85 int i_zz=floor((zz - zmin)/resol[1]);// indice in z, nella proiezione cilindrica, del punto i,k PPA
86 int i_xx=floor((xx - xmin)/resol[0]);// indice in x, nella proiezione cilindrica, del punto i,k PPA
87 // Enrico RHI_ind[k][i]=i_xx+i_zz*x_size;
88 //shift orizzontale negativo del punto di indice i_xx per costruire la finestra in x
89 // se l'estremo minimo in x della finestra è negativo assegno come shift il massimo possibile e cioè la distanza del punto dall'origine
90 i_xx_min(i, k)=i_xx;
91 if (i_xx - w_x_size_2 >= 0)
92 i_xx_min(i, k)= w_x_size_2;
93
94 //shift orizzontale positivo attorno al punto di indice i_xx per costruire la finestra in x
95 int i_xx_max = x_size-i_xx-1;
96 if (i_xx+w_x_size_2 < (int) x_size)
97 i_xx_max = w_x_size_2;
98
99 //shift verticale negativo attorno al punto di indice i_zz per costruire la finestra in z
100 i_zz_min(i, k)=i_zz;
101 if (i_zz_min(i, k) - w_z_size_2 > 0)
102 i_zz_min(i, k) = w_z_size_2;
103
104 //shift verticale positivo attorno al punto di indice i_zz per costruire la finestra in z
105 int i_zz_max = z_size-i_zz-1;
106 if (i_zz+w_z_size_2 < z_size)
107 i_zz_max = w_z_size_2;
108
109 //indici minimo e massimo in x e z per definire la finestra sul punto
110 im(i, k)=i_xx-i_xx_min(i, k);
111 ix(i, k)=i_xx+i_xx_max;
112 jm(i, k)=i_zz-i_zz_min(i, k);
113 jx(i, k)=i_zz+i_zz_max;
114
115 }
116 }
117 /*
118 ;------------------------------------------------------------------------------
119 ; FASE 2
120 ;------------------------------------------------------------------------------
121 ; Costruzione matrice pesi
122 ; Questa matrice contiene i pesi (in funzione della distanza) per ogni punto.
123 ;-----------------------------------------------------------------------------*/
124
125 vector<double> w_x(w_x_size);
126 for (unsigned k=0;k<(unsigned) w_x_size;k++)
127 w_x[k]=exp(-pow(k-w_x_size_2,2.)/pow(w_x_size_2/2.,2.));
128
129 vector<double> w_z(w_z_size);
130 for (unsigned k=0;k<(unsigned)w_z_size;k++)
131 w_z[k]=exp(-pow(k-w_z_size_2,2.)/pow(w_z_size_2/2.,2.));
132
133 Matrix2D<double> w_tot(w_x_size, w_z_size);
134 for (unsigned i=0;i<(unsigned)w_x_size;i++){
135 for (unsigned j=0;j<(unsigned)w_z_size;j++){
136 w_tot(i, j)=w_x[i]*w_z[j];
137 }
138 }
139
140 /* ;----------------------------------------------------------- */
141 /* ; Matrici per puntare sul piano cartesiano velocemente */
142 /* ;---------------------------------- */
143 /* ; FASE 3 */
144 /* ;---------------------------------- */
145 /* ; Selezione dati per formare RHI */
146 /* ;---------------------------------- */
147
148/* for(k=0;k<MAX_BIN;k++){
149 beamXweight[k]=(float **) malloc(w_x_size*sizeof(float *));
150 for(i=0;i<w_x_size;i++){
151 beamXweight[k][i]=(float *) malloc(w_z_size*sizeof(float));
152 }
153 }
154*/
155 Matrix2D<double> RHI_beam(volume.size(), max_bin);
156 for (unsigned iaz=0; iaz<volume.beam_count; iaz++)
157 {
158 Matrix2D<double>& rhi_cart = *slices[iaz];
159 Matrix2D<double> rhi_weight(Matrix2D<double>::Zero(x_size, z_size));
160
161 volume.read_vertical_slice(iaz, RHI_beam, MISSING_DB);
162 /* ;---------------------------------- */
163 /* ; FASE 4 */
164 /* ;---------------------------------- */
165 /* ; Costruzione RHI */
166 /* ;---------------------------------- */
167
168 // Enrico: non sforare se il raggio è piú lungo di MAX_BIN
169 unsigned ray_size = volume.scan(0).beam_size;
170 if (ray_size > max_bin)
171 ray_size = max_bin;
172
173 for (unsigned iel=0;iel<volume.size();iel++){
174 for (unsigned ibin=0;ibin<ray_size;ibin++) {
175 double beamXweight[w_x_size][w_z_size];
176
177 for(unsigned kx=0;kx<(unsigned)w_x_size;kx++){
178 for(unsigned kz=0;kz<(unsigned)w_z_size;kz++){
179//std::cout<<"ibin , kx, kz "<<ibin<<" "<<kx<<" "<<kz<<" "<<w_x_size<< " "<<w_z_size<<" "<<MAX_BIN<<std::endl;
180//std::cout<<"beam "<< beamXweight[ibin][kx][kz]<<std::endl;
181//std::cout<<"RHI "<< RHI_beam[iel][ibin]<<std::endl;
182//std::cout<<"w_tot "<< w_tot[kx][kz]<<std::endl;
183 beamXweight[kx][kz] = RHI_beam(iel, ibin) * w_tot(kx, kz);
184 }
185 }
186 unsigned imin = im(ibin, iel);
187 unsigned imax = ix(ibin, iel);
188 unsigned jmin = jm(ibin, iel);
189 unsigned jmax = jx(ibin, iel);
190
191 unsigned wimin=w_x_size_2-i_xx_min(ibin, iel);
192 //wimax=w_x_size_2+i_xx_max[ibin][iel];
193 unsigned wjmin=w_z_size_2-i_zz_min(ibin, iel);
194 //wjmax=w_z_size_2+i_zz_max[ibin][iel];
195 for (unsigned i=imin;i<=imax;i++) {
196 for (unsigned j=jmin;j<=jmax;j++) {
197 rhi_cart(i, j) = rhi_cart(i, j) + beamXweight[wimin+(i-imin)][wjmin+(j-jmin)];
198 rhi_weight(i, j) = rhi_weight(i, j)+w_tot(wimin+(i-imin), wjmin+(j-jmin));
199 }
200 }
201 }
202 }
203 for (unsigned i=0;i<x_size;i++) {
204 for (unsigned j=0;j<z_size;j++) {
205 if (rhi_weight(i, j) > 0.0)
206 rhi_cart(i, j)=rhi_cart(i, j)/rhi_weight(i, j);
207 else {
208 rhi_cart(i, j)=MISSING_DB;
209
210 }
211 }
212 }
213 }
214}
215
216}
String functions.
Definition: cart.cpp:4
std::vector< Matrix2D< double > * > slices
Vertical rectangular x,z semi-slices of the cylinder, with one side resting on the cylinder axis.
Definition: cylindrical.h:26
double resol[2]
Resolution in x and z.
Definition: cylindrical.h:31