[GRASS-dev] Re: [GRASS GIS] #1575: r.horizon bugfix and speedup

GRASS GIS trac at osgeo.org
Tue Feb 14 11:23:46 EST 2012


#1575: r.horizon bugfix and speedup
-------------------------+--------------------------------------------------
 Reporter:  sprice       |       Owner:  grass-dev@…              
     Type:  enhancement  |      Status:  new                      
 Priority:  normal       |   Milestone:  7.0.0                    
Component:  Default      |     Version:  svn-trunk                
 Keywords:               |    Platform:  Unspecified              
      Cpu:  OSX/Intel    |  
-------------------------+--------------------------------------------------

Comment(by sprice):

 Well let's try that patch again...
 {{{
 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;
         }

 }}}

-- 
Ticket URL: <http://trac.osgeo.org/grass/ticket/1575#comment:1>
GRASS GIS <http://grass.osgeo.org>



More information about the grass-dev mailing list