00001
#include <volume_io/internal_volume_io.h>
00002
#include <bicpl/geom.h>
00003
00004
#undef USE_FAST_APPROX
00005 #define USE_FAST_APPROX
00006
00007
#ifdef USE_FAST_APPROX
00008 #define FAST_SQRT fast_approx_sqrt
00009
#else
00010
#define FAST_SQRT sqrt
00011
#endif
00012
00013 #define MIN_SQRT_ARG 1.0
00014 #define MAX_SQRT_ARG 1000.0
00015 #define N_SQRT_PRECOMPUTE 30000
00016
00017 public double fast_approx_sqrt(
00018
double y )
00019 {
00020 Real x;
00021
static BOOLEAN first =
TRUE;
00022
static struct {
double scale;
double trans; }
00023 precomp_sqrt[
N_SQRT_PRECOMPUTE], *
lookup;
00024
00025
if( y <= MIN_SQRT_ARG || y >=
MAX_SQRT_ARG )
00026 {
00027
return( sqrt(
y ) );
00028 }
00029
00030
if( first )
00031 {
00032
int i;
00033 Real next_sqrt, this_sqrt;
00034
00035 first =
FALSE;
00036
00037 next_sqrt = 0.0;
00038 for_less( i, 0,
N_SQRT_PRECOMPUTE )
00039 {
00040 this_sqrt = next_sqrt;
00041 next_sqrt = sqrt( (Real) (i+1) / (Real)
N_SQRT_PRECOMPUTE *
00042
MAX_SQRT_ARG );
00043
00044 precomp_sqrt[i].scale = next_sqrt - this_sqrt;
00045 precomp_sqrt[i].trans = this_sqrt -
00046 precomp_sqrt[i].scale * (Real) i;
00047 }
00048 }
00049
00050 x = ((Real)
N_SQRT_PRECOMPUTE /
MAX_SQRT_ARG) *
y;
00051
00052
lookup = &precomp_sqrt[(
int) x];
00053
00054
return(
lookup->trans +
lookup->scale * x );
00055 }
00056
00057 public Real
fast_approx_distance_between_points(
00058 Point *p1,
00059 Point *p2 )
00060 {
00061 Real dx, dy, dz;
00062
00063 dx = (Real) Point_x(*p2) - (Real) Point_x(*p1);
00064 dy = (Real) Point_y(*p2) - (Real) Point_y(*p1);
00065 dz = (Real) Point_z(*p2) - (Real) Point_z(*p1);
00066
00067
return(
FAST_SQRT( dx * dx + dy * dy + dz * dz ) );
00068 }
