Elaboradar  0.1
cylindrical.cpp
2 #include <cmath>
3 
4 #define DTOR M_PI/180. /* esternalizzo?*/ //fattore conversione gradi-radianti
5 
6 using namespace std;
7 
8 namespace radarelab {
9 
10 CylindricalVolume::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 
36 void 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