\( \def\bold#1{\bf #1} \newcommand{\d}{\mathrm{d}} \) BTP: Manual and Source Code Documentation

Power Uphill

bike mass [kg]
body mass [kg]
altitude gain [m]
climb length [km]
gradient [%]
time [s]
speed [km/h]
power [W]
power/mass [W/kg]
climbrate [m/min]

average power on climb stage

BTP  3.0
Routing/ClimbAnalysis/PowerCalculation
HeightData.cpp
1 
2 #include "HeightData.h"
3 #define PI 3.14159265
4 #define RO 6371000
5 
6 HeightData::HeightData(){
7  nReserve = 2;
8  sx = sy = -1;
9  grid = NULL;
10 }
11 void HeightData::set_rect(double minlat,double maxlat,double minlon,double maxlon){
12  if(minlat > maxlat) maxlat = minlat = 0;
13  if(minlon > maxlon) maxlon = minlon = 0;
14  DEMRect.minlat = minlat;
15  DEMRect.maxlat = maxlat;
16  DEMRect.minlon = minlon;
17  DEMRect.maxlon = maxlon;
18 }
19 HeightData::~HeightData(void){
20  free_grid();
21 }
22 void HeightData::create_grid(){
23  free_grid();
24  int xminlon = floor(DEMRect.minlon*1200);
25  int xmaxlon = floor(DEMRect.maxlon*1200);
26  int yminlat = floor(DEMRect.minlat*1200);
27  int ymaxlat = floor(DEMRect.maxlat*1200);
28  GridRect.minlon = xminlon-nReserve;
29  GridRect.minlat = yminlat-nReserve;
30  GridRect.maxlon = xmaxlon+nReserve+1;
31  GridRect.maxlat = ymaxlat+nReserve+1;
32  sx = GridRect.maxlon-GridRect.minlon+1;
33  sy = GridRect.maxlat-GridRect.minlat+1;
34 
35  grid = new short*[sy];
36  for(int i = 0; i < sy; i++){
37  grid[i] = new short[sx];
38  for( int j = 0; j < sx; j++)
39  grid[i][j] = -10000;
40  }
41 
42  /*Gitter auslesen*/
43  for( int lon = floor(GridRect.minlon/1200.); lon <= floor(GridRect.maxlon/1200.);lon++)
44  for( int lat = floor(GridRect.minlat/1200.); lat <= floor(GridRect.maxlat/1200.);lat++){
45  char filename[12];
46  char path[1000];
47  if(lat < 0 && lon < 0) sprintf(&(filename[0]),"S%02iW%03i.hgt",-lat,-lon);
48  if(lat < 0 && lon >= 0) sprintf(&(filename[0]),"S%02iE%03i.hgt",-lat, lon);
49  if(lat >= 0 && lon < 0) sprintf(&(filename[0]),"N%02iW%03i.hgt", lat,-lon);
50  if(lat >= 0 && lon >= 0) sprintf(&(filename[0]),"N%02iE%03i.hgt", lat, lon);
51 
52  strcpy(path,SRTM3folder);
53  strcat(path,"\\");
54  strcat(path,filename);
55 
56  int ystart = (lat + 1) * 1200 - max(GridRect.minlat, lat * 1200);
57  int yend = (lat + 1) * 1200 - min(GridRect.maxlat,(lat+1) * 1200);
58  int yos = (lat + 1) * 1200 - GridRect.minlat;
59  int xstart = max(lon * 1200,GridRect.minlon) - lon * 1200;
60  int xend = min((lon+1) * 1200,GridRect.maxlon) - lon * 1200;
61  int xread = xend - xstart + 1;
62  int xos = max(lon * 1200,GridRect.minlon) - GridRect.minlon;
63 
64  FILE* f = fopen(path,"rb");
65  if(f != NULL){
66  /*Datei existiert*/
67  for(int y = ystart; y >= yend; y--){
68  fseek(f,2*(1201*y+xstart),SEEK_SET);
69  fread(&(grid[yos-y][xos]),2,xread,f);
70  }
71  fclose(f);
72  }
73  else
74  /*falls Datei nicht existiert alle Null setzen*/
75  for(int y = ystart; y >= yend; y--)
76  for(int x = xos; x < xos + xread; x++)
77  grid[yos-y][x]=0;
78  }
79 
80  /*Bytereihenfolge umkehren (Little-Big-Endian)*/
81  for(int y = 0; y < sy; y++)
82  for( int x = 0; x < sx; x++)
83  grid[y][x] = (((grid[y][x])>>8)&0xff)+(((grid[y][x]) << 8)&0xff00);
84 }
85 void HeightData::free_grid(){
86  if(grid != NULL){
87  for(int i = 0; i < sy; i++)
88  delete grid[i];
89  if(sy>0)
90  delete grid;
91  grid = NULL;
92  }
93 }
94 int HeightData::lty(double lat){
95  return floor(1200*lat) - GridRect.minlat;
96 }
97 int HeightData::ltx(double lon){
98  return floor(1200*lon) - GridRect.minlon;
99 }
100 double HeightData::ltdy(double lat){
101  return (1200*lat - floor(1200*lat));
102 }
103 double HeightData::ltdx(double lon){
104  return (1200*lon - floor(1200*lon));
105 }
106 void HeightData::interpolate_bicubic_conv3(int width,int height,double alpha,
107  double** h, double** dhdx, double** dhdy){
108  nReserve = 3;
109  create_grid();
110 
111  /*Rechenhilfen vorberechnen*/
112  double dlon = DEMRect.maxlon - DEMRect.minlon;
113  double dlat = DEMRect.maxlat - DEMRect.minlat;
114  CubicKernel3* kx = new CubicKernel3[width];
115  CubicKernel3* kxx = new CubicKernel3[width];
116  for(int x = 0; x < width; x++){
117  kxx[x].index[2] = kx[x].index[2] = ltx(1.*x/(width-1)*dlon + DEMRect.minlon);
118  kxx[x].index[0] = kx[x].index[0] = kx[x].index[2] - 2;
119  kxx[x].index[1] = kx[x].index[1] = kx[x].index[2] - 1;
120  kxx[x].index[3] = kx[x].index[3] = kx[x].index[2] + 1;
121  kxx[x].index[4] = kx[x].index[4] = kx[x].index[2] + 2;
122  kxx[x].index[5] = kx[x].index[5] = kx[x].index[2] + 3;
123  double dx = ltdx(1.*x/(width-1)*dlon + DEMRect.minlon);
124  kx[x].weights[0] = Cubic_conv_kernel3(dx+2);
125  kx[x].weights[1] = Cubic_conv_kernel3(dx+1);
126  kx[x].weights[2] = Cubic_conv_kernel3(dx+0);
127  kx[x].weights[3] = Cubic_conv_kernel3(dx-1);
128  kx[x].weights[4] = Cubic_conv_kernel3(dx-2);
129  kx[x].weights[5] = Cubic_conv_kernel3(dx-3);
130  kxx[x].weights[0] = Cubic_conv_kernel3_d(dx+2);
131  kxx[x].weights[1] = Cubic_conv_kernel3_d(dx+1);
132  kxx[x].weights[2] = Cubic_conv_kernel3_d(dx+0);
133  kxx[x].weights[3] = Cubic_conv_kernel3_d(dx-1);
134  kxx[x].weights[4] = Cubic_conv_kernel3_d(dx-2);
135  kxx[x].weights[5] = Cubic_conv_kernel3_d(dx-3);
136  }
137 
138  CubicKernel3* ky = new CubicKernel3[height];
139  CubicKernel3* kyy = new CubicKernel3[height];
140  for(int y = 0; y < height; y++){
141  kyy[y].index[2] = ky[y].index[2] = lty(1.*y/(height-1)*dlat + DEMRect.minlat);
142  kyy[y].index[0] = ky[y].index[0] = ky[y].index[2] - 2;
143  kyy[y].index[1] = ky[y].index[1] = ky[y].index[2] - 1;
144  kyy[y].index[3] = ky[y].index[3] = ky[y].index[2] + 1;
145  kyy[y].index[4] = ky[y].index[4] = ky[y].index[2] + 2;
146  kyy[y].index[5] = ky[y].index[5] = ky[y].index[2] + 3;
147  double dy = ltdy(1.*y/(height-1)*dlat + DEMRect.minlat);
148  ky[y].weights[0] = Cubic_conv_kernel3(dy+2);
149  ky[y].weights[1] = Cubic_conv_kernel3(dy+1);
150  ky[y].weights[2] = Cubic_conv_kernel3(dy+0);
151  ky[y].weights[3] = Cubic_conv_kernel3(dy-1);
152  ky[y].weights[4] = Cubic_conv_kernel3(dy-2);
153  ky[y].weights[5] = Cubic_conv_kernel3(dy-3);
154  kyy[y].weights[0] = Cubic_conv_kernel3_d(dy+2);
155  kyy[y].weights[1] = Cubic_conv_kernel3_d(dy+1);
156  kyy[y].weights[2] = Cubic_conv_kernel3_d(dy+0);
157  kyy[y].weights[3] = Cubic_conv_kernel3_d(dy-1);
158  kyy[y].weights[4] = Cubic_conv_kernel3_d(dy-2);
159  kyy[y].weights[5] = Cubic_conv_kernel3_d(dy-3);
160  }
161  if(h != NULL){
162  for(int x = 0; x < width; x++)
163  for(int y = 0; y < height; y++){
164  h [x][height-1-y] = Bicubic_conv3(kx [x],ky [y]);
165  /*fuer genaue aber halb so schnelle Berechnung folgende 2 Zeilen
166  aktivieren*/
167  //dhdx[x][height-1-y] = Bicubic_conv3(kxx[x],ky [y]);
168  //dhdy[x][height-1-y] = Bicubic_conv3(kx [x],kyy[y]);
169  }
170  }
171  free_grid();
172  delete(kx);
173  delete(kxx);
174  delete(ky);
175  delete(kyy);
176  /*hier die schnell-Berechnung des Gradienten*/
177  if(dhdx != NULL){
178  double pxperx = width/((DEMRect.maxlon-DEMRect.minlon)*1200.);
179  for(int x = 1; x < width-1; x++)
180  for(int y = 1; y < height-1; y++)
181  dhdx[x][y] = (h[x+1][y]-h[x-1][y])/2.*pxperx;
182  for(int i = 0 ; i < width; i++)
183  dhdx[i][0] = dhdx[i][height-1] = 0;
184  }
185  if(dhdy != NULL){
186  double pxpery = height/((DEMRect.maxlat-DEMRect.minlat)*1200.);
187  for(int x = 1; x < width-1; x++)
188  for(int y = 1; y < height-1; y++)
189  dhdy[x][y] = (h[x][y+1]-h[x][y-1])/2.*pxpery;
190  for(int i = 0 ; i < height; i++)
191  dhdy[0][i] = dhdy[width-1][i] = 0;
192  }
193 }
194 void HeightData::interpolate_bicubic_conv3_mercator(int width,int height,
195  double alpha,double** h,
196  double** dhdx, double** dhdy){
197  nReserve = 3;
198  create_grid();
199 
200  /*Rechenhilfen vorberechnen*/
201  double dlon = DEMRect.maxlon - DEMRect.minlon;
202  CubicKernel3* kx = new CubicKernel3[width];
203  CubicKernel3* kxx = new CubicKernel3[width];
204  for(int x = 0; x < width; x++){
205  kxx[x].index[2] = kx[x].index[2] = ltx((0.5+x)/width*dlon + DEMRect.minlon);
206  kxx[x].index[0] = kx[x].index[0] = kx[x].index[2] - 2;
207  kxx[x].index[1] = kx[x].index[1] = kx[x].index[2] - 1;
208  kxx[x].index[3] = kx[x].index[3] = kx[x].index[2] + 1;
209  kxx[x].index[4] = kx[x].index[4] = kx[x].index[2] + 2;
210  kxx[x].index[5] = kx[x].index[5] = kx[x].index[2] + 3;
211  double dx = ltdx((0.5+x)/width * dlon + DEMRect.minlon);
212  kx[x].weights[0] = Cubic_conv_kernel3(dx+2);
213  kx[x].weights[1] = Cubic_conv_kernel3(dx+1);
214  kx[x].weights[2] = Cubic_conv_kernel3(dx+0);
215  kx[x].weights[3] = Cubic_conv_kernel3(dx-1);
216  kx[x].weights[4] = Cubic_conv_kernel3(dx-2);
217  kx[x].weights[5] = Cubic_conv_kernel3(dx-3);
218  kxx[x].weights[0] = Cubic_conv_kernel3_d(dx+2);
219  kxx[x].weights[1] = Cubic_conv_kernel3_d(dx+1);
220  kxx[x].weights[2] = Cubic_conv_kernel3_d(dx+0);
221  kxx[x].weights[3] = Cubic_conv_kernel3_d(dx-1);
222  kxx[x].weights[4] = Cubic_conv_kernel3_d(dx-2);
223  kxx[x].weights[5] = Cubic_conv_kernel3_d(dx-3);
224  }
225 
226  double* ylat = new double[height];
227  CubicKernel3* ky = new CubicKernel3[height];
228  CubicKernel3* kyy = new CubicKernel3[height];
229  double mmax = lattomercator(DEMRect.maxlat);
230  double mmin = lattomercator(DEMRect.minlat);
231  double dlat = mmax - mmin;
232  for(int y = 0; y < height; y++){
233  ylat[y] = mercatortolat((0.5+y)/height*dlat + mmin);
234  kyy[y].index[2] = ky[y].index[2] = lty(ylat[y]);
235  kyy[y].index[0] = ky[y].index[0] = ky[y].index[2] - 2;
236  kyy[y].index[1] = ky[y].index[1] = ky[y].index[2] - 1;
237  kyy[y].index[3] = ky[y].index[3] = ky[y].index[2] + 1;
238  kyy[y].index[4] = ky[y].index[4] = ky[y].index[2] + 2;
239  kyy[y].index[5] = ky[y].index[5] = ky[y].index[2] + 3;
240  double dy = ltdy(ylat[y]);
241  ky[y].weights[0] = Cubic_conv_kernel3(dy+2);
242  ky[y].weights[1] = Cubic_conv_kernel3(dy+1);
243  ky[y].weights[2] = Cubic_conv_kernel3(dy+0);
244  ky[y].weights[3] = Cubic_conv_kernel3(dy-1);
245  ky[y].weights[4] = Cubic_conv_kernel3(dy-2);
246  ky[y].weights[5] = Cubic_conv_kernel3(dy-3);
247  kyy[y].weights[0] = Cubic_conv_kernel3_d(dy+2);
248  kyy[y].weights[1] = Cubic_conv_kernel3_d(dy+1);
249  kyy[y].weights[2] = Cubic_conv_kernel3_d(dy+0);
250  kyy[y].weights[3] = Cubic_conv_kernel3_d(dy-1);
251  kyy[y].weights[4] = Cubic_conv_kernel3_d(dy-2);
252  kyy[y].weights[5] = Cubic_conv_kernel3_d(dy-3);
253  }
254  if(h != NULL){
255  for(int x = 0; x < width; x++)
256  for(int y = 0; y < height; y++){
257  h [x][height-1-y] = Bicubic_conv3(kx [x],ky [y]);
258  /*fuer genaue aber halb so schnelle Berechnung folgende 2 Zeilen
259  aktivieren*/
260  //dhdx[x][height-1-y] = Bicubic_conv3(kxx[x],ky [y]);
261  //dhdy[x][height-1-y] = Bicubic_conv3(kx [x],kyy[y]);
262  }
263  }
264  free_grid();
265  delete(kx);
266  delete(kxx);
267  delete(ky);
268  delete(kyy);
269  /*hier die schnell-Berechnung des Gradienten*/
270  if(dhdx != NULL){
271  double pxperx = width/((DEMRect.maxlon-DEMRect.minlon)*1200.);
272  for(int x = 1; x < width-1; x++)
273  for(int y = 1; y < height-1; y++)
274  dhdx[x][y] = (h[x+1][y]-h[x-1][y])/2.*pxperx;
275  for(int i = 0 ; i < width; i++)
276  dhdx[i][0] = dhdx[i][height-1] = 0;
277  }
278  if(dhdy != NULL){
279  double pxpery;// = height/((DEMRect.maxlat-DEMRect.minlat)*1200.);
280 
281  for(int y = 1; y < height-1; y++){
282  pxpery = 2./((ylat[y+1]-ylat[y-1])*1200.);
283  for(int x = 1; x < width-1; x++)
284  dhdy[x][y] = (h[x][y+1]-h[x][y-1])/2.*pxpery;
285  }
286  for(int i = 0 ; i < height; i++)
287  dhdy[0][i] = dhdy[width-1][i] = 0;
288  }
289  delete ylat;
290 }
291 inline double HeightData::Cubic_conv_kernel3(double x){
292  x = qAbs(x);
293  if(x >= 3) return 0;
294  if(x > 2) return x*x*x/12.-2.*x*x/3.+21.*x/12.-3./2.;
295  if(x > 1) return -7.*x*x*x/12.+3*x*x-59.*x/12+15./6.;
296  if(x >=0) return 4.*x*x*x/3.-7.*x*x/3.+1.;
297  else return 0;
298 }
299 inline double HeightData::Cubic_conv_kernel3_d(double x){
300  int IsGreaterZero = 1;
301  if(x < 0){
302  IsGreaterZero = 0;
303  x = -x;
304  }
305  double result = 0;
306  if(x > 3) result = 0;
307  if(x < 3 && x >= 2) result = x*x/4.-4.*x/3.+7./4.;
308  if(x < 2 && x >= 1) result = -7.*x*x/4.+6*x-59./12.;
309  if(x < 1 && x >= 0) result = 4.*x*x-14.*x/3.;
310 
311  if(IsGreaterZero)
312  return result;
313  else
314  return -result;
315 }
316 inline double HeightData::Bicubic_conv3(CubicKernel3 kx, CubicKernel3 ky){
317  double result = 0;
318  for( int x = 0; x < 6; x++)
319  for(int y = 0; y < 6; y++)
320  result = result + grid[ky.index[y]][kx.index[x]]
321  * kx.weights[x] * ky.weights[y];
322  return result;
323 }
324 void HeightData::set_SRTM3_folder(char* folder){
325  strcpy(SRTM3folder,folder);
326 }
327 void HeightData::reset_t(double minlat,double maxlat,double minlon,double maxlon){
328  nReserve = 3;
329  set_rect(minlat,maxlat,minlon,maxlon);
330  create_grid();
331 }
332 int HeightData::max_t(){
333  int result = -10000;
334  if(grid != NULL){
335  result = int(grid[0][0]);
336  for(int y = 0; y < sy; y++)
337  for( int x = 0; x < sx; x++)
338  if(grid[y][x] > result)
339  result = grid[y][x];
340  }
341  return result;
342 }
343 double HeightData::height_t(double lon,double lat,void* hc,short mode){
344  CubicKernel3 kx;
345  kx.index[2] = ltx(lon);
346  kx.index[0] = kx.index[2] - 2;
347  kx.index[1] = kx.index[2] - 1;
348  kx.index[3] = kx.index[2] + 1;
349  kx.index[4] = kx.index[2] + 2;
350  kx.index[5] = kx.index[2] + 3;
351  CubicKernel3 ky;
352  ky.index[2] = lty(lat);
353  ky.index[0] = ky.index[2] - 2;
354  ky.index[1] = ky.index[2] - 1;
355  ky.index[3] = ky.index[2] + 1;
356  ky.index[4] = ky.index[2] + 2;
357  ky.index[5] = ky.index[2] + 3;
358 
359  double dx = ltdx(lon);
360  double dy = ltdy(lat);
361 
362  /*kontrolliere, dass indeces im grid enthalten sind*/
363  if(ky.index[0] < 0 || ky.index[5] >= sy ||
364  kx.index[0] < 0 || kx.index[5] >= sx)
365  return 0;
366 
367  double h = 0,dhdx=0,dhdy=0;
368  if(mode == 0 || hc != NULL){
369  /*einfache Hoehe berechnen*/
370  for(int i = 0; i < 6; i++){
371  kx.weights[i] = Cubic_conv_kernel3(dx+2-i);
372  ky.weights[i] = Cubic_conv_kernel3(dy+2-i);
373  }
374  h = Bicubic_conv3(kx,ky);
375  }
376  if(mode == 1 || mode == 3 || mode == 4){
377  /*Ableitung in x berechnen*/
378  for(int i = 0; i < 6; i++){
379  kx.weights[i] = Cubic_conv_kernel3_d(dx+2-i);
380  ky.weights[i] = Cubic_conv_kernel3(dy+2-i);
381  }
382  dhdx = Bicubic_conv3(kx,ky);
383  }
384  if(mode == 2 || mode == 3 || mode == 4){
385  /*Ableitung in y ausgeben*/
386  for(int i = 0; i < 6; i++){
387  kx.weights[i] = Cubic_conv_kernel3(dx+2-i);
388  ky.weights[i] = Cubic_conv_kernel3_d(dy+2-i);
389  }
390  dhdy = Bicubic_conv3(kx,ky);
391  }
392 
393  if(hc != NULL) //Hoehe mitgeben
394  *((float*)hc) = h;
395  if(mode == 0) //Hoehe ausgeben
396  return h;
397  if(mode == 1) //x-Ableitung ausgeben
398  return dhdx;
399  if(mode == 2) //y-Ableitung ausgeben
400  return dhdy;
401  if(mode == 3) //Gradientenrichtung ausgeben
402  return atan2(double(dhdx),double(dhdy*cos(lat/180*PI)));
403  if(mode == 4){ //Vertrauensindex mit- u. Hoehe ausgeben
404  //((double*)hc)[0] = abs(((double*)hc)[2]*dhdy-((double*)hc)[1]*dhdx/cos(lat/180*PI));
405  double hn[4];
406  double dlat = 0.00007, dlon = dlat/cos(lat/180.*PI);
407  hn[0] = height_t(lon-2*dlon*((double*)hc)[1],lat+2*dlat*((double*)hc)[2],NULL,0);
408  hn[1] = height_t(lon-1*dlon*((double*)hc)[1],lat+1*dlat*((double*)hc)[2],NULL,0);
409  hn[2] = height_t(lon+1*dlon*((double*)hc)[1],lat-1*dlat*((double*)hc)[2],NULL,0);
410  hn[3] = height_t(lon+2*dlon*((double*)hc)[1],lat-2*dlat*((double*)hc)[2],NULL,0);
411  ((double*)hc)[0]= (2*hn[0]-hn[1]-2*h-hn[2]+2*hn[3])/14.*
412  qAbs(((double*)hc)[2]*dhdy-((double*)hc)[1]*dhdx/cos(lat/180*PI));
413  return h;
414  }
415  else
416  return 0;
417 }
418 double HeightData::height(double lon,double lat,void* hc,short mode){
419  set_rect(lat,lat,lon,lon);
420  nReserve = 3;
421  create_grid();
422  double result = height_t(lon,lat,hc,mode);
423  free_grid();
424  return result;
425 }
426 void HeightData::paek_pro(double lon, double lat, qint64** c,float* Ahigher,
427  int cc, int ac,float d,float mperlon, float mperlat){
428  int j,i,x,y,ymin,ymax,xmin,xmax, Ac = 0, AcF = 0;
429  float x0 = lon*1200 - GridRect.minlon; //Gipfellage
430  float y0 = lat*1200 - GridRect.minlat;
431  float sx = mperlon/1000; //Skalierung km je x
432  float sy = mperlat/1000; //Saklierung km je y
433  float h0 = height_t(lon,lat,NULL,0);
434 
435  /*Punktzaehler anlegen*/
436  int **wc = new int*[ac];
437  for(i = 0; i < ac; i++)
438  wc[i] = new int[cc];
439 
440  /*Höhenzaehler und Punktzaehler initialisieren*/
441  for(i = 0; i < ac; i++)
442  for(j = 0; j < cc; j++){
443  c[i][j] = 0; //Zählt alle Höhen
444  wc[i][j] = 0; //Anzahl enthaltener Punkte
445  }
446 
447  xmin = qMax(0 ,int(floor(x0-(d*cc)/sx)));
448  xmax = qMin(this->sx-1,int (ceil(x0+(d*cc)/sx)));
449  for(x = xmin; x <= xmax; x++){
450  ymin = qMax(0 ,int(floor(y0-sqrt((d*cc*d*cc)-((x-x0)*sx*(x-x0)*sx)+1)/sx)));
451  ymax = qMin(this->sy-1,int(ceil (y0+sqrt((d*cc*d*cc)-((x-x0)*sx*(x-x0)*sx)+1)/sx)));
452  for(y = ymin; y <= ymax; y++){
453  i = int(floor((atan2(y-y0,x-x0)/PI+1)/2*ac));
454  j = int(floor(sqrt(((y-y0)*sy*(y-y0)*sy)+((x-x0)*sx*(x-x0)*sx))/d));
455  if(i < cc && grid[y][x] < 9000
456  && i < ac && j < cc){
457  c[i][j] = c[i][j] + long(grid[y][x]);
458  wc[i][j]++;
459  if(grid[y][x] > h0) Ac++;
460  AcF++;
461  }
462  }
463  }
464  /*Durchschnitt bilden*/
465  for(i = 0; i < ac; i++)
466  for(j = 0; j < cc; j++)
467  if(wc[i][j] > 0)
468  c[i][j] = c[i][j] / wc[i][j];
469  else
470  c[i][j] = 12000;
471 
472  /*Punktzaehler löschen*/
473  for(i = 0; i < ac; i++)
474  delete(wc[i]);
475  delete(wc);
476 
477  *Ahigher = float(1.*Ac/AcF);
478 
479 }
480 void HeightData::interpolate_linear(int width,int height,double alpha,
481  double** h, double** dhdx, double** dhdy){
482  nReserve = 3;
483  create_grid();
484 
485  int ix,iy;
486  for(int x = 0; x < width; x++)
487  for(int y = 0; y < height; y++){
488  ix = ltx(DEMRect.minlon + (DEMRect.maxlon-DEMRect.minlon)*(1.*x/(width -1)));
489  iy = lty(DEMRect.minlat + (DEMRect.maxlat-DEMRect.minlat)*(1.*y/(height-1)));
490  h[x][height-1-y] = 1.*( + grid[iy ][ix] + grid[iy ][ix+1]
491  + grid[iy+1][ix] + grid[iy+1][ix+1])/4.;
492  dhdx[x][height-1-y] = 1.*( - grid[iy ][ix] + grid[iy ][ix+1]
493  - grid[iy+1][ix] + grid[iy+1][ix+1])/2.;
494  dhdy[x][height-1-y] =-1.*( - grid[iy ][ix] - grid[iy ][ix+1]
495  + grid[iy+1][ix] + grid[iy+1][ix+1])/2.;
496  }
497 }
498 void HeightData::interpolate_linear_mercator(int width,int height,double alpha,
499  double** h, double** dhdx, double** dhdy){
500  nReserve = 3;
501  create_grid();
502  double mmax = lattomercator(DEMRect.maxlat);
503  double mmin = lattomercator(DEMRect.minlat);
504  int ix,iy;
505  for(int y = 0; y < height; y++){
506  iy = lty(mercatortolat(mmin + (mmax-mmin)*((0.5+y)/height)));
507  for(int x = 0; x < width; x++){
508  ix = ltx(DEMRect.minlon + (DEMRect.maxlon-DEMRect.minlon)*((0.5+x)/width));
509  h[x][height-1-y] = 1.*( + grid[iy ][ix] + grid[iy ][ix+1]
510  + grid[iy+1][ix] + grid[iy+1][ix+1])/4.;
511  dhdx[x][height-1-y] = 1.*( - grid[iy ][ix] + grid[iy ][ix+1]
512  - grid[iy+1][ix] + grid[iy+1][ix+1])/2.;
513  dhdy[x][height-1-y] =-1.*( - grid[iy ][ix] - grid[iy ][ix+1]
514  + grid[iy+1][ix] + grid[iy+1][ix+1])/2.;
515  }
516  }
517 }
518 
519 void HeightData::get_sorted_grid_data(short** heights, double** lons, double** lats,
520  long* count){
521  /*extracts all grid values with KOO into three list, sorted by height
522  beginning with the highest*/
523  *count = sy*sx;
524  qsh = *heights = new short[*count];
525  qslon = *lons = new double[*count];
526  qslat= *lats = new double[*count];
527  /*Extracting*/
528  for(int ii = 0; ii < sy; ii++)
529  for( int jj = 0; jj < sx; jj++){
530  (*heights)[ii*sx+jj] = grid[ii][jj];
531  (*lons)[ii*sx+jj] = 1.*(GridRect.minlon+jj)/1200.;
532  (*lats)[ii*sx+jj] = 1.*(GridRect.minlat+ii)/1200.;
533  }
534  /*QuickSort*/
535  sort_GridList(0,*count-1);
536 }
537 void HeightData::sort_GridList(long s, long e){
538  short sb;
539  long l = s, r = e;
540  double p = qsh[(s+e)/2], db;
541  do{
542  while(qsh[l] < p) l++;
543  while(qsh[r] > p) r--;
544  if(l<=r){
545  sb = qsh[l]; qsh[l] = qsh[r]; qsh[r] = sb;
546  db = qslon[l]; qslon[l] = qslon[r]; qslon[r] = db;
547  db = qslat[l]; qslat[l] = qslat[r]; qslat[r] = db;
548  l++;
549  r--;
550  }
551  }while(l <= r);
552  if(s < r) sort_GridList(s,r);
553  if(l < e) sort_GridList(l,e);
554 }
555 recti HeightData::get_GridRect(){
556  return GridRect;
557 }
558 int HeightData::correct_hgt(char* filename){
559  char fn[100];
560  short badheight = -1000;
561  int lnl = strlen(filename);
562  strcpy(&(fn[0]),&(filename[strlen(filename)-11]));
563  int lon = -190,lat = -190;
564  fn[7]='\0';
565  if(fn[3]=='E') lon = atoi(&(fn[4]));
566  if(fn[3]=='W') lon = - atoi(&(fn[4]));
567  fn[3]='\0';
568  if(fn[0]=='N') lat = atoi(&(fn[1]));
569  if(fn[0]=='S') lat = - atoi(&(fn[1]));
570  bool corrected = false;
571  if(lon >= -180 && lon < 180 && lat >= -90 && lat < 90){
572  set_rect(lat,lat+1,lon,lon+1);
573  nReserve = 100;
574  create_grid();
575  for(int x = 99; x < 1300; x++)
576  for(int y = 99; y < 1300; y++)
577  if(grid[y][x] < badheight){
578  corrected = true;
579  int ym = y, yp = y, xm = x, xp = x;
580  double hy = badheight-100., hx = badheight-100.;
581  while(grid[ym][x] < badheight && ym > 0) ym--;
582  while(grid[yp][x] < badheight && yp < sy-1) yp++;
583  while(grid[y][xm] < badheight && xm > 0) xm--;
584  while(grid[y][xp] < badheight && xp < sx-1) xp++;
585  short y0 = grid[y][x];
586  short y1 = grid[yp][x];
587  short y2 = grid[ym][x];
588  short y3 = grid[y][xp];
589  short y4 = grid[y][xm];
590  if(grid[ym][x] > badheight && grid[yp][x] > badheight)
591  hy = ((yp-y)*grid[ym][x]+(y-ym)*grid[yp][x])/(yp-ym);
592  if(grid[ym][x] > badheight && grid[yp][x] > badheight)
593  hx = ((xp-x)*grid[y][xm]+(x-xm)*grid[y][xp])/(xp-xm);
594  if(hy > badheight && hx > badheight)
595  grid[y][x] = (hx+hy)/2;
596  else{
597  if(hy > badheight)
598  grid[y][x] = hy;
599  else
600  grid[y][x] = hx;
601  }
602  }
603  if(corrected){
604  /*Datei existiert*/
605  FILE* f = fopen(filename,"wb");
606  for(int y = 1300; y >= 99; y--){
607  for(int x = 99; x < 1300; x++)
608  grid[y][x] = (((grid[y][x])>>8)&0xff)+(((grid[y][x]) << 8)&0xff00);
609  fwrite(&(grid[y][100]),2,1201,f);
610  }
611  fclose(f);
612  }
613  return corrected;
614  }
615  else
616  return 0;
617 }
HeightData::grid
short ** grid
Definition: HeightData.h:110
HeightData::get_sorted_grid_data
void get_sorted_grid_data(short **heights, double **lons, double **lats, long *count)
returns sorted raw data heights within the specified area
Definition: HeightData.cpp:519
HeightData::DEMRect
rect DEMRect
Definition: HeightData.h:105
HeightData::SRTM3folder
char SRTM3folder[1000]
Definition: HeightData.h:107
HeightData::height_t
double height_t(double lon, double lat, void *hc, short mode)
returns height data operating at the precalculated data grid
Definition: HeightData.cpp:343
HeightData::interpolate_linear
void interpolate_linear(int width, int height, double alpha, double **h, double **dhdx, double **dhdy)
returns heigth data grid h with size (width,height)
Definition: HeightData.cpp:480
HeightData::height
double height(double lon, double lat, void *hc, short mode)
returns height data
Definition: HeightData.cpp:418
HeightData::interpolate_bicubic_conv3
void interpolate_bicubic_conv3(int width, int height, double alpha, double **h, double **dhdx, double **dhdy)
returns heigth data grid h with size (width,height)
Definition: HeightData.cpp:106
HeightData::Cubic_conv_kernel3
double Cubic_conv_kernel3(double x)
calculates weights of the bicubic convolution for the absolute height
Definition: HeightData.cpp:291
HeightData::ltx
int ltx(double lon)
Definition: HeightData.cpp:97
CubicKernel3
holds weights for cubic interpolation
Definition: HeightData.h:13
HeightData::free_grid
void free_grid()
Definition: HeightData.cpp:85
recti
holds data grid information with angles multiplied by 1200
Definition: HeightData.h:11
HeightData::GridRect
recti GridRect
Definition: HeightData.h:106
HeightData::create_grid
void create_grid()
Definition: HeightData.cpp:22
HeightData::paek_pro
void paek_pro(double lon, double lat, qint64 **c, float *Ahigher, int cc, int ac, float d, float mperlon, float mperlat)
evaluates prominence of summits
Definition: HeightData.cpp:426
HeightData::Bicubic_conv3
double Bicubic_conv3(CubicKernel3 kx, CubicKernel3 ky)
executes the bicubic convolution with the kernels kx and ky
Definition: HeightData.cpp:316
HeightData::set_rect
void set_rect(double minlat, double maxlat, double minlon, double maxlon)
specifies area of the precalculated data grid
Definition: HeightData.cpp:11
HeightData::reset_t
void reset_t(double minlat, double maxlat, double minlon, double maxlon)
generates the data grid being
Definition: HeightData.cpp:327
HeightData::interpolate_bicubic_conv3_mercator
void interpolate_bicubic_conv3_mercator(int width, int height, double alpha, double **h, double **dhdx, double **dhdy)
returns heigth data grid h with size (width,height)
Definition: HeightData.cpp:194
HeightData::ltdx
double ltdx(double lon)
Definition: HeightData.cpp:103
HeightData::nReserve
int nReserve
Definition: HeightData.h:104
HeightData::set_SRTM3_folder
void set_SRTM3_folder(char *folder)
specifies the folder where hgt-files of SRTM3 data are searched
Definition: HeightData.cpp:324
HeightData::ltdy
double ltdy(double lat)
Definition: HeightData.cpp:100
HeightData::Cubic_conv_kernel3_d
double Cubic_conv_kernel3_d(double x)
calculates weights of the bicubic convolution for the gradient
Definition: HeightData.cpp:299
HeightData::interpolate_linear_mercator
void interpolate_linear_mercator(int width, int height, double alpha, double **h, double **dhdx, double **dhdy)
returns heigth data grid h with size (width,height)
Definition: HeightData.cpp:498
HeightData::max_t
int max_t()
maximum height of the precalculated data grid
Definition: HeightData.cpp:332
HeightData::lty
int lty(double lat)
Definition: HeightData.cpp:94
HeightData::sort_GridList
void sort_GridList(long s, long e)
quick sort used by get_sorted_grid_data()
Definition: HeightData.cpp:537
Generated on Sun Mar 2 2014 12:09:53 for BTP by   doxygen 1.8.5