10CylindricalVolume::CylindricalVolume(
const Volume<double>& volume,
double missing_value,
double x_res,
double z_res)
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.;
16 double xmin = floor(range_min*cos(volume.elevation_max()*DTOR));
17 double zmin = volume[0].sample_height(0) / 1000. + volume.radarSite.getTotalHeight();
18 double xmax = floor(range_maxLowestRay*cos(volume.elevation_min()*DTOR));
19 double zmax = volume.back().sample_height(volume.back().beam_size - 1) / 1000. + volume.radarSite.getTotalHeight();
22 x_size = (xmax - xmin) / x_res;
24 if (x_size > volume.max_beam_size()) x_size=volume.max_beam_size();
25 z_size = (zmax - zmin) / 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)));
36void CylindricalVolume::resample(
const Volume<double>& volume)
38 unsigned max_bin = x_size;
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.;
54 double xmin=floor(range_min*cos(volume.elevation_max()*DTOR));
55 double zmin=volume[0].sample_height(0) / 1000. + volume.radarSite.getTotalHeight();
59 const double w_size[2]={3.,0.3};
63 int w_x_size=ceil((w_size[0]/
resol[0])/2)*2+1;
64 int w_z_size=ceil((w_size[1]/
resol[1])/2)*2+1;
66 if (w_x_size < 3) w_x_size=3;
67 if (w_z_size < 3) w_z_size=3;
69 int w_x_size_2=w_x_size/2;
70 int w_z_size_2=w_z_size/2;
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.;
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();
84 double xx = range*cos(elev_rad);
85 int i_zz=floor((zz - zmin)/
resol[1]);
86 int i_xx=floor((xx - xmin)/
resol[0]);
91 if (i_xx - w_x_size_2 >= 0)
92 i_xx_min(i, k)= w_x_size_2;
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;
101 if (i_zz_min(i, k) - w_z_size_2 > 0)
102 i_zz_min(i, k) = w_z_size_2;
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;
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;
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.));
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.));
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];
155 Matrix2D<double> RHI_beam(volume.size(), max_bin);
156 for (
unsigned iaz=0; iaz<volume.beam_count; iaz++)
158 Matrix2D<double>& rhi_cart = *
slices[iaz];
159 Matrix2D<double> rhi_weight(Matrix2D<double>::Zero(x_size, z_size));
161 volume.read_vertical_slice(iaz, RHI_beam, MISSING_DB);
169 unsigned ray_size = volume.scan(0).beam_size;
170 if (ray_size > max_bin)
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];
177 for(
unsigned kx=0;kx<(unsigned)w_x_size;kx++){
178 for(
unsigned kz=0;kz<(unsigned)w_z_size;kz++){
183 beamXweight[kx][kz] = RHI_beam(iel, ibin) * w_tot(kx, kz);
186 unsigned imin = im(ibin, iel);
187 unsigned imax = ix(ibin, iel);
188 unsigned jmin = jm(ibin, iel);
189 unsigned jmax = jx(ibin, iel);
191 unsigned wimin=w_x_size_2-i_xx_min(ibin, iel);
193 unsigned wjmin=w_z_size_2-i_zz_min(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));
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);
208 rhi_cart(i, j)=MISSING_DB;
std::vector< Matrix2D< double > * > slices
Vertical rectangular x,z semi-slices of the cylinder, with one side resting on the cylinder axis.
double resol[2]
Resolution in x and z.