Menu
▾
▴
[Meshdb-cvs] src/geo/alt main.c,1.3.2.6,1.3.2.7
|
From: <le...@pr...> - 2005-01-23 01:52:42
|
Update of /cvsroot/meshdb/src/geo/alt In directory sc8-pr-cvs1.sourceforge.net:/tmp/cvs-serv24905 Modified Files: Tag: leonard-dev main.c Log Message: replace stupid neighbour algorithm with a more realstic projection algorithm Index: main.c =================================================================== RCS file: /cvsroot/meshdb/src/geo/alt/main.c,v retrieving revision 1.3.2.6 retrieving revision 1.3.2.7 diff -u -d -r1.3.2.6 -r1.3.2.7 --- main.c 22 Jan 2005 09:51:37 -0000 1.3.2.6 +++ main.c 23 Jan 2005 01:52:23 -0000 1.3.2.7 @@ -264,7 +264,7 @@ hm = (float *)malloc((w+2) * (h+2) * sizeof (float)); /* Shift hm so that H[-h2-1,-w2-1] == h[0] (step of w+1) */ -#define MAP(m,x,y) (m)[(x) + (w+2)*(y) - (w+2)*(h+2)/2] +#define MAP(m,x,y) (m)[((x) + (w2+1)) + (w+2)*((y) + (h2+1))] #define H(x,y) MAP(hm,x,y) /* Compute the height map */ @@ -272,61 +272,149 @@ n = n0 - y / scale; for (x = -w2 -1; x <= w2; x++) { e = e0 + x / scale; - H(x,y) = alt_fastat(&alt, e, n); +#if 1 + H(x,y) = alt_fastat(&alt, e, n); +#else + /* Synthetic test features */ + { + float d = sqrtf((float)x * x + (float)y * y); + float h = 0; + float th = atan2f(y,x); + + if ((0.5 < th && th < 0.7) || + (1.8 < th && th < 2.2)) { + if (75 < d && d <= 80) + h = 4 * (d-75); + else if (80 < d && d <= 85) + h = 4 * (85-d); + } + H(x,y) = h; + } +#endif } } if (lflag) { + /* + * The line-of-sight algorithm works like this: + * Imaging a vertical wall surrounding the + * grid. On this wall we mark the 'shadow' of + * processed cells. + * We process cells in (ractangular) rings. + * Each ring's cells are processed like this: + * Step 1. compute that cell's 'rise over run' + * r = (h(x,y)-h0)/sqrt(x^2+y^2) + * Step 2. Project the cell onto the wall to compare + * its r value against the previous ring's shadow. + * If r is bigger, then the cell is visible + * + * After processing all the cells in the ring, + * the next step is to re-process the wall shadow. + * Step 3. Project each wall cell onto the ring. + * Compute the projected point's rise/run (perhaps + * interpolating from the height map) + * Step 4. Update the wall shadow if the rise/run + * has increased. + * + * Unlike previous algorithms, this one retains + * directional precision. Also of benefit is that + * the running time is still O(n^2). + */ + int n, nmax = w2 > h2 ? w2 : h2; - azm = (float *)malloc((w+2) * (h+2) * sizeof (float)); - azmaxm = (float *)malloc((w+2) * (h+2) * sizeof (float)); + /* Initialise the wall to a r/r of ZENITH (-Inf) */ +#define ZENITH (-Inf) + float *wall_e = (float*)malloc((h+2) * sizeof (float)); + float *wall_w = (float*)malloc((h+2) * sizeof (float)); + float *wall_n = (float*)malloc((h+2) * sizeof (float)); + float *wall_s = (float*)malloc((h+2) * sizeof (float)); + float *az_e = (float *)malloc((h+2) * sizeof (float)); + float *az_w = (float *)malloc((h+2) * sizeof (float)); + float *az_n = (float *)malloc((w+2) * sizeof (float)); + float *az_s = (float *)malloc((w+2) * sizeof (float)); + + for (x = -w2-1; x < w2+1; x++) + wall_n[x+w2+1] = wall_s[x+w2+1] = ZENITH; + for (y = -h2-1; y < h2+1; y++) + wall_e[y+h2+1] = wall_w[y+h2+1] = ZENITH; + visiblem = (char *)malloc((w+2) * (h+2) * sizeof (char)); -#define AZ(x,y) MAP(azm,x,y) -#define AZMAX(x,y) MAP(azmaxm,x,y) #define VISIBLE(x,y) MAP(visiblem,x,y) -#define ZENITH (-Inf) - AZ(0,0) = AZMAX(0,0) = ZENITH; VISIBLE(0,0) = 1; - for (n = 1; n < nmax; n++) { int p; - float d; -#define LOS(x,y,xp,yp,xq,yq) \ - if (-w2<(x) && (x)<w2 && -h2<(y) && (y)<h2) { \ - float az, azprev, azmaxprev; \ - az = (H(x,y) - h0) / d; \ - if (isnanf(az)) az = ZENITH; \ - azprev = AZ(xp,yp)*(p/n) + \ - AZ(xq,yq)*(1-p/n); \ - azmaxprev = AZMAX(xp,yp)*(p/n) + \ - AZMAX(xq,yq)*(1-p/n); \ - if (az < azmaxprev) { \ - AZMAX(x,y) = azmaxprev; \ - VISIBLE(x,y) = 0; \ - } else { \ - AZMAX(x,y) = az; \ - VISIBLE(x,y) = 1; \ - } \ - AZ(x,y) = az; \ - } + for (y = -n; y <= n; y++) { + float d = sqrt(n*n + y*y); + float aze = (H(n,y)-h0) / d; + float azw = (H(-n,y)-h0) / d; + float yp = y * h2 / (float)n; + int ypI = floor(yp); + float ypM = yp - ypI; + float shadow_e = wall_e[ypI+h2+1] * (1-ypM) + + wall_e[ypI+h2+1+1] * ypM; + float shadow_w = wall_w[ypI+h2+1] * (1-ypM) + + wall_w[ypI+h2+1+1] * ypM; + + if (isnanf(aze)) aze = ZENITH; + if (isnanf(azw)) azw = ZENITH; + az_e[y+h2+1] = aze; + az_w[y+h2+1] = azw; + VISIBLE(n,y) = (aze >= shadow_e); + VISIBLE(-n,y) = (azw >= shadow_w); + } + + for (y = -h2; y < h2; y++) { + float yr = y * n / (float)h2; + int yrI = floor(yr); + float yrM = yr - yrI; + float aze = az_e[yrI+h2+1] * (1-yrM) + + az_e[yrI+h2+1+1] * yrM; + float azw = az_w[yrI+h2+1] * (1-yrM) + + az_w[yrI+h2+1+1] * yrM; + if (wall_e[y+h2+1] < aze) + wall_e[y+h2+1] = aze; + if (wall_w[y+h2+1] < azw) + wall_w[y+h2+1] = azw; + } + + for (x = -n; x <= n; x++) { + float d = sqrt(n*n + x*x); + float azn = (H(x,n)-h0) / d; + float azs = (H(x,-n)-h0) / d; + float xp = x * w2 / (float)n; + int xpI = floor(xp); + float xpM = xp - xpI; + float shadow_n = wall_n[xpI+w2+1] * (1-xpM) + + wall_n[xpI+w2+1+1] * xpM; + float shadow_s = wall_s[xpI+w2+1] * (1-xpM) + + wall_s[xpI+w2+1+1] * xpM; + + if (isnanf(azn)) azn = ZENITH; + if (isnanf(azs)) azs = ZENITH; + az_n[x+w2+1] = azn; + az_s[x+w2+1] = azs; + VISIBLE(x,n) = (azn >= shadow_n); + VISIBLE(x,-n) = (azs >= shadow_s); + } + + for (x = -w2; x < w2; x++) { + float xr = x * n / (float)w2; + int xrI = floor(xr); + float xrM = xr - xrI; + float azn = az_n[xrI+w2+1] * (1-xrM) + + az_n[xrI+w2+1+1] * xrM; + float azs = az_s[xrI+w2+1] * (1-xrM) + + az_s[xrI+w2+1+1] * xrM; + if (wall_n[x+w2+1] < azn) + wall_n[x+w2+1] = azn; + if (wall_s[x+w2+1] < azs) + wall_s[x+w2+1] = azs; + } + - for (p = 1; p <= n; p++) { - d = sqrt(n*n + p*p); - LOS(n,p,n-1,p-1,n-1,p) /* ESE */ - LOS(-p,n,-(p-1),n-1,-p,n-1) /* SSW */ - LOS(-n,-p,-(n-1),-(p-1),-(n-1),-p) /* WNW */ - LOS(p,-n,p-1,-(n-1),p,-(n-1)) /* NNE */ - } - for (p = 0; p < n; p++) { - d = sqrt(n*n + p*p); - LOS(n,-p,n-1,-(p-1),n-1,-p) /* ENE */ - LOS(p,n,p-1,n-1,p,n-1) /* SSE */ - LOS(-n,p,-(n-1),p-1,-(n-1),p) /* WSW */ - LOS(-p,-n,-(p-1),-(n-1),-p,-(n-1)) /* NNW */ - } } } @@ -355,7 +443,7 @@ dy = y - cross[i].y; if ((dx == dy || dx == -dy) && (-CROSS2 < dx && dx < CROSS2)) - iscross = 1; + iscross = 2; break; case 'x': @@ -378,8 +466,14 @@ /* Crosses (and circles in red) */ if (iscross) { - rgb[0] = 255; - rgb[1] = rgb[2] = 0; + if (iscross == 1) { + rgb[0] = 255; + rgb[1] = rgb[2] = 0; + } else { + rgb[0] = 255; + rgb[1] = 0; + rgb[2] = 128; + } goto plot; } |
