[GRASS-dev] Re: [GRASS GIS] #1575: r.horizon bugfix and speedup
GRASS GIS
trac at osgeo.org
Tue Feb 14 18:02:22 EST 2012
#1575: r.horizon bugfix and speedup
-------------------------+--------------------------------------------------
Reporter: sprice | Owner: grass-dev@…
Type: enhancement | Status: new
Priority: normal | Milestone: 7.0.0
Component: Raster | Version: svn-trunk
Keywords: r.horizon | Platform: All
Cpu: OSX/Intel |
-------------------------+--------------------------------------------------
Changes (by hamish):
* keywords: => r.horizon
* platform: Unspecified => All
* component: Default => Raster
Old description:
> I mentioned this on the mailing list, but didn't get a response. So
> here's a more official report:
>
> I made some tweaks to r.horizon/main.c to almost half its runtime and
> remove a bug. By my tests, it's output is identical to before, but a 219
> minute run has been reduced to 110 minutes. Mostly due to removal of
> unneeded floor() calls. I could use some people to test and make sure I'm
> not introducing any bugs, though.
> Thanks,
> Seth
>
> --- main.orig 2011-08-14 06:32:49.000000000 -0600
> +++ main.c 2012-02-12 19:20:15.000000000 -0700
> @@ -144,12 +144,15 @@
> /* use G_distance() instead ??!?! */
> double distance(double x1, double x2, double y1, double y2)
> {
> + double diffX = x1 - x2;
> + double diffY = y1 - y2;
> +
> if (ll_correction) {
> - return DEGREEINMETERS * sqrt(coslatsq * (x1 - x2) * (x1 - x2)
> - + (y1 - y2) * (y1 - y2));
> + return DEGREEINMETERS * sqrt(coslatsq * diffX * diffX
> + + diffY * diffY);
> }
> else {
> - return sqrt((x1 - x2) * (x1 - x2) + (y1 - y2) * (y1 - y2));
> + return sqrt(diffX * diffX + diffY * diffY);
> }
> }
>
> @@ -841,36 +844,28 @@
>
> int new_point()
> {
> - int iold, jold;
> - int succes = 1, succes2 = 1;
> - double sx, sy;
> - double dx, dy;
> -
> - iold = ip;
> - jold = jp;
> + int iold = ip;
> + int jold = jp;
>
> - while (succes) {
> + while (1) {
> yy0 += stepsinangle;
> xx0 += stepcosangle;
>
> /* offset 0.5 cell size to get the right cell i, j */
> - sx = xx0 * invstepx + offsetx;
> - sy = yy0 * invstepy + offsety;
> - ip = (int)sx;
> - jp = (int)sy;
> + ip = (int) (xx0 * invstepx + offsetx);
> + jp = (int) (yy0 * invstepy + offsety);
>
> /* test outside of raster */
> if ((ip < 0) || (ip >= n) || (jp < 0) || (jp >= m))
> return (3);
>
> if ((ip != iold) || (jp != jold)) {
> - dx = (double)ip *stepx;
> - dy = (double)jp *stepy;
> + double dx = (double)ip *stepx;
> + double dy = (double)jp *stepy;
>
> length = distance(xg0, dx, yg0, dy); /* dist from orig. grid
> point
> to the current grid point */
> - succes2 = test_low_res();
> - if (succes2 == 1) {
> + if (test_low_res() == 1) {
> zp = z[jp][ip];
> return (1);
> }
> @@ -882,54 +877,44 @@
>
> int test_low_res()
> {
> - int iold100, jold100;
> - double sx, sy;
> - int delx, dely, mindel;
> - double zp100, z2, curvature_diff;
> -
> - iold100 = ip100;
> - jold100 = jp100;
> - ip100 = floor(ip / 100.);
> - jp100 = floor(jp / 100.);
> + int iold100 = ip100;
> + int jold100 = jp100;
> + ip100 = ip * .01;
> + jp100 = jp * .01;
> /*test the new position with low resolution */
> if ((ip100 != iold100) || (jp100 != jold100)) {
> /* G_debug(3,"%d %d %d %d\n",ip,jp, iold100,jold100); */
> /* replace with approximate version
> curvature_diff = EARTHRADIUS*(1.-cos(length/EARTHRADIUS));
> - */
> +
> + Code folded into the 'if' statement:
> curvature_diff = 0.5 * length * length * invEarth;
> z2 = z_orig + curvature_diff + length * tanh0;
> zp100 = z100[jp100][ip100];
> + if (zp100 <= z2) {...
> + */
> /*G_debug(3,"%d %d %lf %lf \n",ip,jp,z2,zp100); */
>
> - if (zp100 <= z2)
> + if (z100[jp100][ip100] <= z_orig + 0.5 * length * length *
> invEarth +
> length * tanh0)
> /*skip to the next lowres cell */
> {
> - delx = 32000;
> - dely = 32000;
> + int mindel;
> + int delx = 32000;
> + int dely = 32000;
> + double sx = (xx0 * invstepx + offsetx) * .01;
> + double sy = (yy0 * invstepy + offsety) * .01;
> +
> if (cosangle > 0.) {
> - sx = xx0 * invstepx + offsetx;
> - delx =
> - floor(fabs
> - ((ceil(sx / 100.) - (sx / 100.)) *
> distcosangle));
> + delx = fabs((ceil(sx) - sx) * distcosangle);
> }
> - if (cosangle < 0.) {
> - sx = xx0 * invstepx + offsetx;
> - delx =
> - floor(fabs
> - ((floor(sx / 100.) - (sx / 100.)) *
> distcosangle));
> + else if (cosangle < 0.) {
> + delx = fabs((floor(sx) - sx) * distcosangle);
> }
> if (sinangle > 0.) {
> - sy = yy0 * invstepy + offsety;
> - dely =
> - floor(fabs
> - ((ceil(sy / 100.) - (sy / 100.)) *
> distsinangle));
> + dely = fabs((ceil(sy) - sy) * distsinangle);
> }
> else if (sinangle < 0.) {
> - sy = yy0 * invstepy + offsety;
> - dely =
> - floor(fabs
> - ((floor(jp / 100.) - (sy / 100.)) *
> distsinangle));
> + dely = fabs((floor(sy) - sy) * distsinangle);
> }
>
> mindel = min(delx, dely);
> @@ -953,17 +938,14 @@
>
> double searching()
> {
> - double z2;
> - double curvature_diff;
> - int succes = 1;
> -
> if (zp == UNDEFZ)
> return 0;
>
> while (1) {
> - succes = new_point();
> + double z2;
> + double curvature_diff;
>
> - if (succes != 1) {
> + if (new_point() != 1) {
> break;
> }
New description:
I mentioned this on the mailing list, but didn't get a response. So here's
a more official report:
I made some tweaks to r.horizon/main.c to almost half its runtime and
remove a bug. By my tests, it's output is identical to before, but a 219
minute run has been reduced to 110 minutes. Mostly due to removal of
unneeded floor() calls. I could use some people to test and make sure I'm
not introducing any bugs, though.
Thanks,
Seth
--
--
Ticket URL: <https://trac.osgeo.org/grass/ticket/1575#comment:2>
GRASS GIS <http://grass.osgeo.org>
More information about the grass-dev
mailing list