GIS空间分析实现源码

GIS空间分析实现源码

http://blog.sina.com.cn/s/blog_70214d730100nw3e.html

GIS空间分析实现源码

#include
#include
#include
#include
#include
#include
#include "geometryc.h"

int angle_contains_ray_2d ( float x1, float y1, float x2, float y2,
  float x3, float y3, float x, float y )


{   float a1;
  float a2;
  a1 = angle_deg_2d ( x1, y1, x2, y2, x, y );
  a2 = angle_deg_2d ( x1, y1, x2, y2, x3, y3 );
  if ( a1 <= a2 ) {
    return TRUE;
  } else {
    return FALSE;
  }
}

float angle_deg_2d ( float x1, float y1, float x2, float y2, float x3,
  float y3 )


{   float angle_rad_2d;
  float value;
  float x;
  float y;
  x = ( x1 - x2 ) * ( x3 - x2 ) + ( y1 - y2 ) * ( y3 - y2 );
  y = ( x1 - x2 ) * ( y3 - y2 ) - ( y1 - y2 ) * ( x3 - x2 );
  if ( x == 0.0 && y == 0.0 ) {
    value = 0.0;
  }
  else {
    angle_rad_2d = atan2 ( y, x );
    if ( angle_rad_2d < 0.0 ) {
      angle_rad_2d = angle_rad_2d + 2.0 * PI;
    }
    value = angle_rad_2d * RAD_TO_DEG;
  }
  return value;
}
float angle_rad_2d ( float x1, float y1, float x2, float y2,
  float x3, float y3 )


{   float value;
  float x;
  float y;
  x = ( x1 - x2 ) * ( x3 - x2 ) + ( y1 - y2 ) * ( y3 - y2 );
  y = ( x1 - x2 ) * ( y3 - y2 ) - ( y1 - y2 ) * ( x3 - x2 );
  if ( x == 0.0 && y == 0.0 ) {
    value = 0.0;
  }
  else {
    value = atan2 ( y, x );
    if ( value < 0.0 ) {
      value = value + 2.0 * PI;
    }
  }
  return value;
}
float angle_rad_3d ( float x1, float y1, float z1, float x2, float y2,
  float z2, float x3, float y3, float z3 )


{   float dot;
  float v1norm;
  float v2norm;
  float value;
  dot = dot0_3d ( x2, y2, z2, x1, y1, z1, x3, y3, z3 );
  v1norm = enorm0_3d ( x1, y1, z1, x2, y2, z2 );
  v2norm = enorm0_3d ( x3, y3, z3, x2, y2, z2 );
  if ( v1norm == 0.0 || v2norm == 0.0 ) {
    value = 0.0;
  }
  else {
    value = acos ( dot / ( v1norm * v2norm ) );
  }
  return value;
}
float angle_rad_nd ( int n, float vec1[], float vec2[] )


{   float dot;
  float v1norm;
  float v2norm;
  float value;
  dot = dot_nd ( n, vec1, vec2 );
  v1norm = enorm_nd ( n, vec1 );
  v2norm = enorm_nd ( n, vec2 );
  if ( v1norm == 0.0 || v2norm == 0.0 ) {
    value = 0.0;
  }
  else {
    value = acos ( dot / ( v1norm * v2norm ) );
  }
  return value;
}
float anglei_deg_2d ( float x1, float y1, float x2, float y2, float x3,
  float y3 )


{   float value;
  float x;
  float y;
  x = ( x1 - x2 ) * ( x3 - x2 ) + ( y1 - y2 ) * ( y3 - y2 );
  y = ( x1 - x2 ) * ( y3 - y2 ) - ( y1 - y2 ) * ( x3 - x2 );
  if ( x == 0.0 && y == 0.0 ) {
    value = 0.0;
  }
  else {
    value = atan2 ( y, x );
    if ( value < 0.0 ) {
      value = value + 2.0 * PI;
    }
    value = value * RAD_TO_DEG;
    if ( value > 180.0 ) {
      value = 360.0 - value;
    }
  }
  return value;
}
float anglei_rad_2d ( float x1, float y1, float x2, float y2, float x3,
  float y3 )


{   float value;
  float x;
  float y;
  x = ( x1 - x2 ) * ( x3 - x2 ) + ( y1 - y2 ) * ( y3 - y2 );
  y = ( x1 - x2 ) * ( y3 - y2 ) - ( y1 - y2 ) * ( x3 - x2 );
  if ( x == 0.0 && y == 0.0 ) {
    value = 0.0;
  }
  else {
    value = atan2 ( y, x );
    if ( value < 0.0 ) {
      value = value + 2.0 * PI;
    }
    if ( value > PI ) {
      value = 2.0 * PI - value;
    }
  }
  return value;
}
float arc_cosine ( float c )


{   if ( c <= -1.0 ) {
    return PI;
  } else if ( c >= 1.0 ) {
    return 0.0;
  } else {
    return acos ( c );
  }
}

float atan4 ( float y, float x )


{

  if ( x == 0.0E+00 ) {
    if ( y > 0.0E+00 ) {
      return ( PI / 2.0E+00 );
    } else if ( y < 0.0E+00 ) {
      return ( 3.0E+00 * PI / 2.0E+00 );
    } else if ( y == 0.0E+00 ) {
      return ( 0.0E+00 );
    }
  } else if ( y == 0.0E+00 ) {
    if ( x > 0.0E+00 ) {
      return 0.0E+00;
    } else if ( x < 0.0E+00 ) {
      return PI;
    }
  }

  if        ( x > 0.0E+00 && y > 0.0E+00 ) {
    return                  atan2 (  y,  x );
  } else if ( x < 0.0E+00 && y > 0.0E+00 ) {
    return (           PI - atan2 (  y, -x ) );
  } else if ( x < 0.0E+00 && y < 0.0E+00 ) {
    return (           PI + atan2 ( -y, -x ) );
  } else if ( x > 0.0E+00 && y < 0.0E+00 ) {
    return ( 2.0E+00 * PI - atan2 ( -y,  x ) );
  }
  return 0.0;
}
int box_contains_point_3d ( float x1, float y1, float z1, float x2,
  float y2, float z2, float x3, float y3, float z3, float x4, float y4,
  float z4, float x, float y, float z )


{   float dot;
  float normsq;
  dot = dot0_3d ( x1, y1, z1, x2, y2, z2, x, y, z );
  if ( dot < 0.0 ) {
    return FALSE;
  }
  normsq = enormsq0_3d ( x1, y1, z1, x2, y2, z2 );
  if ( normsq < dot ) {
    return FALSE;
  }
  dot = dot0_3d ( x1, y1, z1, x3, y3, z3, x, y, z );
  if ( dot < 0.0 ) {
    return FALSE;
  }
  normsq = enormsq0_3d ( x1, y1, z1, x3, y3, z3 );
  if ( normsq < dot ) {
    return FALSE;
  }
  dot = dot0_3d ( x1, y1, z1, x4, y4, z4, x, y, z );
  if ( dot < 0.0 ) {
    return FALSE;
  }
  normsq = enormsq0_3d ( x1, y1, z1, x4, y4, z4 );
  if ( normsq < dot ) {
    return FALSE;
  }
  return TRUE;
}
float box_point_dist_3d ( float x1, float y1, float z1, float x2,
  float y2, float z2, float x3, float y3, float z3, float x4, float y4,
  float z4, float x, float y, float z )


{   float dis;
  float dist;
  float x5;
  float x6;
  float x7;
  float x8;
  float y5;
  float y6;
  float y7;
  float y8;
  float z5;
  float z6;
  float z7;
  float z8;

  x5 = x2 + x3 - x1;
  y5 = y2 + y3 - y1;
  z5 = z2 + z3 - z1;
  x6 = x2 + x4 - x1;
  y6 = y2 + y4 - y1;
  z6 = z2 + z4 - z1;
  x7 = x3 + x4 - x1;
  y7 = y3 + y4 - y1;
  z7 = z3 + z4 - z1;
  x8 = x2 + x3 + x4 - 2.0 * x1;
  y8 = y2 + y3 + y4 - 2.0 * y1;
  z8 = z2 + z3 + z4 - 2.0 * z1;

  dis = para_point_dist_3d ( x1, y1, z1, x2, y2, z2, x3, y3, z3, x, y, z );
  dist = dis;
  dis = para_point_dist_3d ( x1, y1, z1, x2, y2, z2, x4, y4, z4, x, y, z );
  if ( dis < dist ) {
    dist = dis;
  }
  dis = para_point_dist_3d ( x1, y1, z1, x3, y3, z3, x4, y4, z4, x, y, z );
  if ( dis < dist ) {
    dist = dis;
  }
  dis = para_point_dist_3d ( x8, y8, z8, x5, y5, z5, x6, y6, z6, x, y, z );
  if ( dis < dist ) {
    dist = dis;
  }
  dis = para_point_dist_3d ( x8, y8, z8, x5, y5, z5, x7, y7, z7, x, y, z );
  if ( dis < dist ) {
    dist = dis;
  }
  dis = para_point_dist_3d ( x8, y8, z8, x6, y6, z6, x7, y7, z7, x, y, z );
  if ( dis < dist ) {
    dist = dis;
  }
  return dist;
}
void circle_dia2imp_2d ( float x1, float y1, float x2, float y2, float *r,
  float *xc, float *yc )


{   *r = 0.5 * enorm0_2d ( x1, y1, x2, y2 );
  *xc = 0.5 * ( x1 + x2 );
  *yc = 0.5 * ( y1 + y2 );
  return;
}
int circle_exp_contains_point_2d ( float x1, float y1, float x2, float y2,
  float x3, float y3, float x, float y )


{   float a[4][4];
  float det;
  int inside;

  if ( x1 == x2 && y1 == y2 ) {
    if ( x1 == x3 && y1 == y3 ) {
      if ( x1 == x && y1 == y ) {
        inside = 6;
      }
      else {
        inside = 7;
      }
    }
    else {
      det = ( x1 - x3 ) * ( y - y3 ) - ( x - x3 ) * ( y1 - y3 );
      if ( det == 0.0 ) {
        inside = 4;
      }
      else {
        inside = 5;
      }
    }
    return inside;
  }

  if ( x1 == x3 && y1 == y3 ) {
    det = ( x1 - x2 ) * ( y - y2 ) - ( x - x2 ) * ( y1 - y2 );
    if ( det == 0.0 ) {
      inside = 4;
    }
    else {
      inside = 5;
    }
    return inside;
  }

  det = ( x1 - x2 ) * ( y3 - y2 ) - ( x3 - x2 ) * ( y1 - y2 );
  if ( det == 0.0 ) {
    det = ( x1 - x2 ) * ( y - y2 ) - ( x - x2 ) * ( y1 - y2 );
    if ( det == 0.0 ) {
      inside = 2;
    }
    else {
      inside = 3;
    }
    return inside;
  }

  a[0][0] = x1;
  a[0][1] = y1;
  a[0][2] = x1 * x1 + y1 * y1;
  a[0][3] = 1.0;
  a[1][0] = x2;
  a[1][1] = y2;
  a[1][2] = x2 * x2 + y2 * y2;
  a[1][3] = 1.0;
  a[2][0] = x3;
  a[2][1] = y3;
  a[2][2] = x3 * x3 + y3 * y3;
  a[2][3] = 1.0;
  a[3][0] = x;
  a[3][1] = y;
  a[3][2] = x * x + y * y;
  a[3][3] = 1.0;
  det = rmat4_det ( a );
  if ( det < 0.0 ) {
    inside = 1;
  }
  else if ( det == 0.0 ) {
    inside = 0;
  }
  else {
    inside = -1;
  }
  return inside;
}

int circle_imp_contains_point_2d ( float radius, float xc, float yc,
  float x, float y )


{   if ( enorm0_2d ( x, y, xc, yc ) <= radius ) {
    return TRUE;
  }
  else {
    return FALSE;
  }
}
void circle_imp_line_par_int_2d ( float radius, float xc, float yc,
  float x0, float y0, float f, float g, int *nint, float x[], float y[] )


{   float root;
  float t;
  root = radius * radius * ( f * f + g * g )
    - ( f * ( yc - y0 ) - g * ( xc - x0 ) )
    * ( f * ( yc - y0 ) - g * ( xc - x0 ) );
  if ( root < 0.0 ) {
    *nint = 0;
  }
  else if ( root == 0.0 ) {
    *nint = 1;
    t = ( f * ( xc - x0 ) + g * ( yc - y0 ) ) / ( f * f + g * g );
    x[0] = x0 + f * t;
    y[0] = y0 + g * t;
  }
  else if ( root > 0.0 ) {
    *nint = 2;
    t = ( ( f * ( xc - x0 ) + g * ( yc - y0 ) ) - sqrt ( root ) )
      / ( f * f + g * g );
    x[0] = x0 + f * t;
    y[0] = y0 + g * t;
    t = ( ( f * ( xc - x0 ) + g * ( yc - y0 ) ) + sqrt ( root ) )
      / ( f * f + g * g );
    x[1] = x0 + f * t;
    y[1] = y0 + g * t;
  }
  return;
}

void circle_points_2d ( int n, float x[], float y[] )


{   float angle;
  int i;
  if ( n <= 0 ) {
    printf ( "\n" );
    printf ( "CIRCLE_POINTS_2D - Fatal error\n" );
    printf ( "  The input value of N is %d\n", n );
    printf ( "  N must be at least 1.\n" );
    exit ( EXIT_FAILURE );
  }
  for ( i = 0; i < n; i++ ) {
    angle = ( 2.0 * PI * ( float ) i ) / ( float ) n ;
    x[i] = cos ( angle );
    y[i] = sin ( angle );
  }
  return;
}
float cone_area_3d ( float h, float r ) {


  float area;
  area = PI * r * sqrt ( h * h + r * r );
  return area;
}
float cone_volume_3d ( float h, float r )


{   float volume;
  volume = PI * r * r * h / 3.0;
  return volume;
}
float cot_rad ( float angle )


{   float value;
  value = cos ( angle ) / sin ( angle );
  return value;
}
float cot_deg ( float angle )


{   float angle_rad;
  float value;
  angle_rad = angle * DEG_TO_RAD;
  value  = cos ( angle_rad ) / sin ( angle_rad );
  return value;
}
float cross_2d ( float x1, float y1, float x2, float y2 )


{   float value;
  value = x1 * y2 - y1 * x2;
  return value;
}
void cross_3d ( float x1, float y1, float z1, float x2, float y2, float z2,
  float *x3, float *y3, float *z3 )


{   *x3 = y1 * z2 - z1 * y2;
  *y3 = z1 * x2 - x1 * z2;
  *z3 = x1 * y2 - y1 * x2;
  return;
}
float cross0_2d ( float x0, float y0, float x1, float y1, float x2,
  float y2 )


{   float value;
  value =
      ( x1 - x0 ) * ( y2 - y0 )
    - ( y1 - y0 ) * ( x2 - x0 );
  return value;
}
void cross0_3d ( float x0, float y0, float z0, float x1, float y1, float z1,
  float x2, float y2, float z2, float *x3, float *y3, float *z3 )


{   *x3 =
      ( y1 - y0 ) * ( z2 - z0 )
    - ( z1 - z0 ) * ( y2 - y0 );
  *y3 =
      ( z1 - z0 ) * ( x2 - x0 )
    - ( x1 - x0 ) * ( z2 - z0 );
  *z3 =
      ( x1 - x0 ) * ( y2 - y0 )
    - ( y1 - y0 ) * ( x2 - x0 );
  return;
}
void direction_pert_3d ( float sigma, int *iseed, float vbase[3],
  float vran[3] )


{   float dphi;
  int i;
  float phi;
  float psi;
  float theta;
  float vdot;
  float x;
  float x2;
  float y;
  float y2;
  float z;
  float z2;

  if ( sigma >= 1.0 ) {
    for ( i = 0; i < 3; i++ ) {
      vran[i] = vbase[i];
    }
  }
  else if ( sigma <= 0.0 ) {
    vdot = 2.0 * uniform_01_sample ( iseed ) - 1.0;
    phi = acos ( vdot );
    theta = 2.0 * PI * uniform_01_sample ( iseed );
    vran[0] = cos ( theta ) * sin ( phi );
    vran[1] = sin ( theta ) * sin ( phi );
    vran[2] = cos ( phi );
  }
  else {
    phi = acos ( vbase[2] );
    theta = atan2 ( vbase[1], vbase[0] );

    x = uniform_01_sample ( iseed );
    x = exp ( ( 1.0 - sigma ) * log ( x ) );
    vdot = 2.0 * x - 1.0;
    dphi = acos ( vdot );

    x = cos ( theta ) * sin ( phi + dphi );
    y = sin ( theta ) * sin ( phi + dphi );
    z = cos ( phi + dphi );

    psi = 2.0 * PI * uniform_01_sample ( iseed );
    vector_rotate_3d ( x, y, z, &x2, &y2, &z2, vbase[0], vbase[1], vbase[2], psi
);
    vran[0] = x2;
    vran[1] = y2;
    vran[2] = z2;
  }
  return;
}
void direction_random_3d ( int *iseed, float vran[3] )


{   float phi;
  float theta;
  float vdot;

  vdot = 2.0 * uniform_01_sample ( iseed ) - 1.0;
  phi = acos ( vdot );

  theta = 2.0 * PI * uniform_01_sample ( iseed );
  vran[0] = cos ( theta ) * sin ( phi );
  vran[1] = sin ( theta ) * sin ( phi );
  vran[2] = cos ( phi );
  return;
}
void direction_random_nd ( int n, int *iseed, float w[] )


{   int i;
  float norm;

  for ( i = 0; i < n; i++ ) {
    w[i] = normal_01_sample ( iseed );
  }

  norm = 0.0;
  for ( i = 0; i < n; i++ ) {
    norm = norm + w[i] * w[i];
  }
  norm = sqrt ( norm );

  for ( i = 0; i < n; i++ ) {
    w[i] = w[i] / norm;
  }
  return;
}
float dot_2d ( float x1, float y1, float x2, float y2 )


{   float value;
  value = x1 * x2 + y1 * y2;
  return value;
}
float dot_3d ( float x1, float y1, float z1, float x2, float y2, float z2 )


{   float value;
  value =
    x1 * x2 +
    y1 * y2 +
    z1 * z2;
  return value;
}
float dot_nd ( int n, float vec1[], float vec2[] )


{   int i;
  float value;
  value = 0.0;
  for ( i = 0; i < n; i++ ) {
    value = value + vec1[i] * vec2[i];
  }
  return value;
}
float dot0_2d ( float x0, float y0, float x1, float y1, float x2,
  float y2 )


{   float value;
  value =
    ( x1 - x0 ) * ( x2 - x0 ) +
    ( y1 - y0 ) * ( y2 - y0 );
  return value;
}
float dot0_3d ( float x0, float y0, float z0, float x1, float y1, float z1,
  float x2, float y2, float z2 )


{   float value;
  value =
    ( x1 - x0 ) * ( x2 - x0 ) +
    ( y1 - y0 ) * ( y2 - y0 ) +
    ( z1 - z0 ) * ( z2 - z0 );
  return value;
}
float enorm_2d ( float x1, float y1 )


{   float value;
  value = sqrt ( x1 * x1 + y1 * y1 );
  return value;
}
float enorm_3d ( float x1, float y1, float z1 )


{   float value;
  value = sqrt (
    x1 * x1 +
    y1 * y1 +
    z1 * z1 );
  return value;
}
float enorm_nd ( int n, float x[] )


{   int i;
  float value;
  value = 0.0;
  for ( i = 0; i < n; i++ ) {
    value = value + x[i] * x[i];
  }
  value = sqrt ( value );
  return value;
}
float enorm0_2d ( float x0, float y0, float x1, float y1 )


{   float value;
  value = sqrt (
    ( x1 - x0 ) * ( x1 - x0 ) +
    ( y1 - y0 ) * ( y1 - y0 ) );
  return value;
}
float enorm0_3d ( float x0, float y0, float z0, float x1, float y1, float z1 )


{   float value;
  value = sqrt (
    ( x1 - x0 ) * ( x1 - x0 ) +
    ( y1 - y0 ) * ( y1 - y0 ) +
    ( z1 - z0 ) * ( z1 - z0 ) );
  return value;
}
float enorm0_nd ( int n, float x[], float y[] )


{   int i;
  float value;
  value = 0.0;
  for ( i = 0; i < n; i++ ) {
    value = value + ( x[i] - y[i] ) * ( x[i] - y[i] );
  }
  value = sqrt ( value );
  return value;
}
float enormsq0_2d ( float x0, float y0, float x1, float y1 )


{   float value;
  value =
    ( x1 - x0 ) * ( x1 - x0 ) +
    ( y1 - y0 ) * ( y1 - y0 );
  return value;
}
float enormsq0_3d ( float x0, float y0, float z0, float x1, float y1,
  float z1 )


{   float value;
  value =
    ( x1 - x0 ) * ( x1 - x0 ) +
    ( y1 - y0 ) * ( y1 - y0 ) +
    ( z1 - z0 ) * ( z1 - z0 );
  return value;
}
float enormsq0_nd ( int n, float x0[], float x1[] )


{   float value;
  int i;
  value = 0.0;
  for ( i = 0; i < n; i++ ) {
    value = value + ( x1[i] - x0[i] ) * ( x1[i] - x0[i] );
  }
  return value;
}
int get_seed ( void )


{   time_t clock;
  int i;
  int ihour;
  int imin;
  int isec;
  int iseed;
  static int iseed_internal = 0;
  struct tm *lt;
  time_t tloc;

  if ( iseed_internal == 0 ) {
    clock = time ( &tloc );
    lt = localtime ( &clock );

    ihour = lt->tm_hour;
    if ( ihour > 12 ) {
      ihour = ihour - 12;
    }

    ihour = ihour - 1;
    imin = lt->tm_min;
    isec = lt->tm_sec;
    iseed_internal = isec + 60 * ( imin + 60 * ihour );

    iseed_internal = iseed_internal + 1;

    iseed_internal = ( int ) (
      ( float ) iseed_internal * ( float ) I_MAX / ( 60.0 * 60.0 * 12.0 )
    );
  }

  if ( iseed_internal == 0 ) {
    iseed_internal = 1;
  }
  if ( iseed_internal == I_MAX ) {
    iseed_internal = I_MAX - 1;
  }

  for ( i = 0; i < 10; i++ ) {
    uniform_01_sample ( &iseed_internal );
  }
  iseed = iseed_internal;
  return iseed;
}
int halfspace_imp_triangle_int_3d ( float x1, float y1, float z1, float x2,
  float y2, float z2, float x3, float y3, float z3, float a, float b,
  float c, float d, float x[], float y[], float z[] )


{   float dist1;
  float dist2;
  float dist3;
  int num_int;

  dist1 = a * x1 + b * y1 + c * z1 + d;
  dist2 = a * x2 + b * y2 + c * z2 + d;
  dist3 = a * x3 + b * y3 + c * z3 + d;

  num_int = halfspace_triangle_int_3d ( x1, y1, z1, x2, y2,
    z2, x3, y3, z3, dist1, dist2, dist3, x, y, z );
  return num_int;
}
int halfspace_norm_triangle_int_3d ( float x1, float y1, float z1, float x2,
  float y2, float z2, float x3, float y3, float z3, float xp, float yp,
  float zp, float xn, float yn, float zn, float x[], float y[], float z[] )


{   float d;
  float dist1;
  float dist2;
  float dist3;
  int num_int;

  d = - xn * xp - yn * yp - zn * zp;
  dist1 = xn * x1 + yn * y1 + zn * z1 + d;
  dist2 = xn * x2 + yn * y2 + zn * z2 + d;
  dist3 = xn * x3 + yn * y3 + zn * z3 + d;

  num_int = halfspace_triangle_int_3d ( x1, y1, z1, x2, y2,
    z2, x3, y3, z3, dist1, dist2, dist3, x, y, z );
  return num_int;
}
int halfspace_triangle_int_3d ( float x1, float y1, float z1,
  float x2, float y2, float z2, float x3, float y3, float z3, float dist1,
  float dist2, float dist3, float x[], float y[], float z[] )


{   int num_int;

  num_int = 0;
  if ( dist1 <= 0.0 ) {
    x[num_int] = x1;
    y[num_int] = y1;
    z[num_int] = z1;
    num_int = num_int + 1;
  }
  if ( dist1 * dist2 < 0.0 ) {
    x[num_int] = ( dist1 * x2 - dist2 * x1 ) / ( dist1 - dist2 );
    y[num_int] = ( dist1 * y2 - dist2 * y1 ) / ( dist1 - dist2 );
    z[num_int] = ( dist1 * z2 - dist2 * z1 ) / ( dist1 - dist2 );
    num_int = num_int + 1;
  }
  if ( dist2 <= 0.0 ) {
    x[num_int] = x2;
    y[num_int] = y2;
    z[num_int] = z2;
    num_int = num_int + 1;
  }
  if ( dist2 * dist3 < 0.0 ) {
    x[num_int] = ( dist2 * x3 - dist3 * x2 ) / ( dist2 - dist3 );
    y[num_int] = ( dist2 * y3 - dist3 * y2 ) / ( dist2 - dist3 );
    z[num_int] = ( dist2 * z3 - dist3 * z2 ) / ( dist2 - dist3 );
    num_int = num_int + 1;
  }
  if ( dist3 <= 0.0 ) {
    x[num_int] = x3;
    y[num_int] = y3;
    z[num_int] = z3;
    num_int = num_int + 1;
  }
  if ( dist3 * dist1 < 0.0 ) {
    x[num_int] = ( dist3 * x1 - dist1 * x3 ) / ( dist3 - dist1 );
    y[num_int] = ( dist3 * y1 - dist1 * y3 ) / ( dist3 - dist1 );
    z[num_int] = ( dist3 * z1 - dist1 * z3 ) / ( dist3 - dist1 );
    num_int = num_int + 1;
  }
  return num_int;
}
int i_random ( int ilo, int ihi, int *iseed )


{   int i;
  float r;
  float rhi;
  float rlo;

  r = uniform_01_sample ( iseed );

  rlo = ( ( float ) ilo ) - 0.5;
  rhi = ( ( float ) ihi ) + 0.5;

  r = ( ( 1.0 - r ) * rlo + r * rhi );
  if ( r < 0.0 ) {
    r = r - 0.5;
  }
  else {
    r = r + 0.5;
  }
  i = ( int ) r;

  if ( i < ilo ) {
    i = ilo;
  }
  if ( i > ihi ) {
    i = ihi;
  }
  return i;
}
hello! I am a gis fans!!
2003-12-14 21:17:00举报帖子
使用道具
    gisjian

等级：侠之大者
威望：2
文章：306
积分：766
门派：无门无派
注册：2002年12月8日第 3 楼       

float line_exp_point_dist_2d ( float x1, float y1, float x2, float y2,
  float x, float y )


{   float bot;
  float dist;
  float dot;
  float t;
  float xn;
  float yn;
  bot = enormsq0_2d ( x1, y1, x2, y2 );
  if ( bot == 0.0 ) {
    xn = x1;
    yn = y1;
  }

  else {
    dot =
        ( x - x1 ) * ( x2 - x1 )
      + ( y - y1 ) * ( y2 - y1 );
    t = dot / bot;
    xn = x1 + t * ( x2 - x1 );
    yn = y1 + t * ( y2 - y1 );
  }
  dist = enorm0_2d ( xn, yn, x, y );
  return dist;
}
float line_exp_point_dist_3d ( float x1, float y1, float z1, float x2,
  float y2, float z2, float x, float y, float z )


{   float bot;
  float dist;
  float t;
  float xn;
  float yn;
  float zn;
  bot = enormsq0_3d ( x1, y1, z1, x2, y2, z2 );
  if ( bot == 0.0 ) {
    xn = x1;
    yn = y1;
    zn = z1;
  }

  else {
    t = (
      ( x - x1 ) * ( x2 - x1 ) +
      ( y - y1 ) * ( y2 - y1 ) +
      ( z - z1 ) * ( z2 - z1 ) ) / bot;
    xn = x1 + t * ( x2 - x1 );
    yn = y1 + t * ( y2 - y1 );
    zn = z1 + t * ( z2 - z1 );
  }

  dist = enorm0_3d ( x, y, z, xn, yn, zn );
  return dist;
}
float line_exp_point_dist_signed_2d ( float x1, float y1, float x2,
  float y2, float x, float y )


{   float a;
  float b;
  float c;
  float dist_signed;

  line_exp2imp_2d ( x1, y1, x2, y2, &a, &b, &c );

  dist_signed = ( a * x + b * y + c ) / sqrt ( a * a + b * b );
  return dist_signed;
}
void line_exp2imp_2d ( float x1, float y1, float x2, float y2, float *a,
  float *b, float *c )


{
  if ( x1 == x2 && y1 == y2 ) {
    printf ( "\n" );
    printf ( "LINE_EXP2IMP_2D - Fatal error!\n" );
    printf ( "  (X1,Y1) = (X2,Y2)\n" );
    printf ( "  (X1,Y1) = %f %f\n", x1, y1 );
    printf ( "  (X2,Y2) = %f %f\n", x2, y2 );
    exit ( EXIT_FAILURE );
  }
  *a = y2 - y1;
  *b = x1 - x2;
  *c = x2 * y1 - x1 * y2;
  return;
}
void line_exp2par_3d ( float x1, float y1, float z1, float x2, float y2, 
  float z2, float *f, float *g, float *h, float *x0, float *y0, float *z0 )


{   float norm;

  if ( x1 == x2 && y1 == y2 && z1 == z2 ) {
    *f = 0.0;
    *g = 0.0;
    *h = 0.0;
    *x0 = x1;
    *y0 = y1;
    *z0 = z1;
  }
  else {
    norm = sqrt (
        ( x2 - x1 ) * ( x2 - x1 )
      + ( y2 - y1 ) * ( y2 - y1 )
      + ( z2 - z1 ) * ( z2 - z1 ) );
    *f = ( x2 - x1 ) / norm;
    *g = ( y2 - y1 ) / norm;
    *h = ( z2 - z1 ) / norm;
    *x0 = x1;
    *y0 = y1;
    *z0 = z1;
  }
  return;
}

分享