From fec6a053be28c539f5868acffac2d53f5a00ef53 Mon Sep 17 00:00:00 2001 From: cvs-fast-export Date: Tue, 21 Dec 2004 07:53:42 +0000 Subject: [PATCH] Synthetic commit for incomplete tag rel-7-0-1 svn:revision:23684 cvs:account:cvs-fast-export --- librt/const.c | 52 -- librt/hist.c | 170 ---- librt/plane.c | 2379 ----------------------------------------------- librt/polylib.c | 551 ----------- librt/polyno.h | 43 - librt/snoise.c | 1797 ----------------------------------- 6 files changed, 4992 deletions(-) delete mode 100644 librt/const.c delete mode 100644 librt/hist.c delete mode 100644 librt/plane.c delete mode 100644 librt/polylib.c delete mode 100644 librt/polyno.h delete mode 100644 librt/snoise.c diff --git a/librt/const.c b/librt/const.c deleted file mode 100644 index e50afe91e0a..00000000000 --- a/librt/const.c +++ /dev/null @@ -1,52 +0,0 @@ -/* - * C O N S T . C - * - * Constants used by the ray tracing library. - * - * Source - - * SECAD/VLD Computing Consortium, Bldg 394 - * The U. S. Army Ballistic Research Laboratory - * Aberdeen Proving Ground, Maryland 21005 - * - * Copyright Notice - - * This software is Copyright (C) 1989 by the United States Army. - * All rights reserved. - */ -#ifndef lint -static char RCSmat[] = "@(#)$Header$ (BRL)"; -#endif - -#include "conf.h" - -#include -#include "machine.h" -#include "vmath.h" -#include "raytrace.h" - -CONST double rt_pi = 3.14159265358979323846; /* pi */ -CONST double rt_twopi = 6.28318530717958647692; /* pi*2 */ -CONST double rt_halfpi = 1.57079632679489661923; /* pi/2 */ -CONST double rt_quarterpi=0.78539816339744830961; /* pi/4 */ -CONST double rt_invpi = 0.318309886183790671538; /* 1/pi */ -CONST double rt_inv2pi = 0.159154943091895335769; /* 1/(pi*2) */ -CONST double rt_inv4pi = 0.07957747154594766788; /* 1/(pi*4) */ - -CONST double rt_inv255 = 1.0/255.0; - -CONST double rt_degtorad = 0.0174532925199433; /* (pi*2)/360 */ -CONST double rt_radtodeg = 57.29577951308230698802; /* 360/(pi*2) */ - -CONST mat_t rt_identity = { - 1, 0, 0, 0, - 0, 1, 0, 0, - 0, 0, 1, 0, - 0, 0, 0, 1 -}; - -#ifdef ipsc860 -/* Hack to work around the broken iPSC/860 loader. It doesn't - * see the definition in global.c, probably because there is no - * .text or .data there that it needs. - */ -struct rt_g rt_g; -#endif diff --git a/librt/hist.c b/librt/hist.c deleted file mode 100644 index cc77fc1d170..00000000000 --- a/librt/hist.c +++ /dev/null @@ -1,170 +0,0 @@ -/* - * H I S T . C - * - * General purpose histogram handling routines - * - * The macro RT_HISTOGRAM_TALLY is used to record items that - * live in a single "bin", while the subroutine rt_hist_range() - * is used to record items that may extend across multiple "bin"s. - * - * Author - - * Michael John Muuss - * - * Source - - * SECAD/VLD Computing Consortium, Bldg 394 - * The U. S. Army Ballistic Research Laboratory - * Aberdeen Proving Ground, Maryland 21005-5066 - * - * Copyright Notice - - * This software is Copyright (C) 1990 by the United States Army. - * All rights reserved. - */ -#ifndef lint -static char RCShist[] = "@(#)$Header$ (BRL)"; -#endif - -#include "conf.h" - -#include -#include -#include "machine.h" -#include "vmath.h" -#include "raytrace.h" -#include "./debug.h" - -/* - * R T _ H I S T _ F R E E - */ -void -rt_hist_free( histp ) -struct histogram *histp; -{ - /* The 'clean' command can be called before the first rt_prep(), - * resulting in zero magic numbers in rtip->rti_hist_cellsize. - */ - if( histp && histp->magic == 0 ) return; - RT_CK_HISTOGRAM(histp); - if( histp->hg_bins ) - rt_free( (char *)histp->hg_bins, "old histogram bins"); - histp->hg_bins = (long *)0; - histp->hg_nbins = 0; - histp->magic = -1; /* sanity */ -} - -/* - * R T _ H I S T _ I N I T - * - * Initialize a histogram structure. - * It is expected that the structure is junk upon entry. - */ -void -rt_hist_init( histp, min, max, nbins ) -struct histogram *histp; -fastf_t min, max; -int nbins; -{ - - if( max <= min ) max = min+1; - if( nbins < 2 ) { - nbins = 2; - } else if( nbins > 10000 ) { - nbins = 10000; - } - - histp->hg_min = floor(min); - histp->hg_max = ceil(max); - histp->hg_nbins = nbins; - - /* When max-min <= nbins, clumpsize should be 1 */ - histp->hg_clumpsize = ((max-min)/nbins)+1; - if( histp->hg_clumpsize <= 0 ) histp->hg_clumpsize = 1; - - histp->hg_nsamples = 0L; - histp->hg_bins = (long *)rt_calloc( nbins+1, sizeof(long), "histogram bins"); - histp->magic = RT_HISTOGRAM_MAGIC; -} - -/* - * R T _ H I S T _ R A N G E - */ -void -rt_hist_range( hp, low, high ) -register struct histogram *hp; -fastf_t low; -fastf_t high; -{ - long a; - long b; - register int i; - - RT_CK_HISTOGRAM(hp); - if( low <= hp->hg_min ) - a = 0; - else - a = (low - hp->hg_min) / hp->hg_clumpsize; - if( high >= hp->hg_max ) - b = hp->hg_nbins-1; - else - b = (high - hp->hg_min) / hp->hg_clumpsize; - if( a < 0 || b >= hp->hg_nbins ) rt_bomb("rt_hist_range() out of range\n"); - for( i=a; i <= b; i++ ) { - hp->hg_bins[i]++; - } - hp->hg_nsamples++; -} - -/* - * R T _ H I S T _ P R - */ -void -rt_hist_pr( histp, title ) -register struct histogram *histp; -CONST char *title; -{ - register int i; - long maxcount; - static CONST char marks[] = "################################################################"; -#define NMARKS 50 - char buf[256]; - int percent; - int mark_count; - int val; - int nbins; - - RT_CK_HISTOGRAM(histp); - - /* Find entry with highest count */ - maxcount = 0L; - for( i=0; ihg_nbins; i++ ) { - if( histp->hg_bins[i] > maxcount ) - maxcount = histp->hg_bins[i]; - } - if( maxcount <= 0 ) maxcount = 1; - - /* Supress trailing bins with zero counts. nbins s/b >= 1 */ - for( nbins = histp->hg_nbins-1; nbins >= 1; nbins-- ) - if(histp->hg_bins[nbins] > 0) break; - - /* 12345678 12345678 123 .... */ - rt_log("\nHistogram of %s\nmin=%g, max=%g, nbins=%d, clumpsize=%g\n%d samples collected, highest count was %d\n\n Value Count Rel%%| Bar Graph\n", - title, - histp->hg_min, histp->hg_max, - histp->hg_nbins, histp->hg_clumpsize, - histp->hg_nsamples, maxcount ); - - /* Print each bin. Final bins with zero counts are supressed. */ - for( i=0; i <= nbins; i++ ) { - percent = (int)(((double)histp->hg_bins[i])*100.0/maxcount); - mark_count = percent*NMARKS/100; - if( mark_count <= 0 ) { - buf[0] = '\0'; - } else { - bcopy( marks, buf, mark_count ); - buf[mark_count] = '\0'; - } - val = histp->hg_min + i*histp->hg_clumpsize; - rt_log("%8d %8d %3d |%s\n", - val, - histp->hg_bins[i], percent, buf ); - } -} diff --git a/librt/plane.c b/librt/plane.c deleted file mode 100644 index 79860afc35b..00000000000 --- a/librt/plane.c +++ /dev/null @@ -1,2379 +0,0 @@ -/* - * P L A N E . C - * - * Some useful routines for dealing with planes and lines. - * - * Author - - * Michael John Muuss - * - * Source - - * SECAD/VLD Computing Consortium, Bldg 394 - * The U. S. Army Ballistic Research Laboratory - * Aberdeen Proving Ground, Maryland 21005 - * - * Distribution Status - - * Public Domain, Distribution Unlimited - */ -#ifndef lint -static char RCSplane[] = "@(#)$Header$ (BRL)"; -#endif - -#include "conf.h" - -#include -#include - -#include "machine.h" -#include "vmath.h" -#include "raytrace.h" -#include "./debug.h" - -/* - * R T _ D I S T _ P T 3 _ P T 3 - * - * Returns distance between two points. - */ -double -rt_dist_pt3_pt3( a, b ) -CONST point_t a; -CONST point_t b; -{ - vect_t diff; - - VSUB2( diff, a, b ); - return MAGNITUDE( diff ); -} - -/* - * R T _ P T 3 _ P T 3 _ E Q U A L - * - * Returns - - * 1 if the two points are equal, within the tolerance - * 0 if the two points are not "the same" - */ -int -rt_pt3_pt3_equal( a, b, tol ) -CONST point_t a; -CONST point_t b; -CONST struct rt_tol *tol; -{ - vect_t diff; - - RT_CK_TOL(tol); - VSUB2( diff, b, a ); - if( MAGSQ( diff ) < tol->dist_sq ) return 1; - return 0; -} - -/* - * R T _ P T 2 _ P T 2 _ E Q U A L - * - * Returns - - * 1 if the two points are equal, within the tolerance - * 0 if the two points are not "the same" - */ -int -rt_pt2_pt2_equal( a, b, tol ) -CONST point_t a; -CONST point_t b; -CONST struct rt_tol *tol; -{ - vect_t diff; - - RT_CK_TOL(tol); - VSUB2_2D( diff, b, a ); - if( MAGSQ_2D( diff ) < tol->dist_sq ) return 1; - return 0; -} - -/* - * R T _ 3 P T S _ C O L L I N E A R - * - * Check to see if three points are collinear. - * - * The algorithm is designed to work properly regardless of the - * order in which the points are provided. - * - * Returns (boolean) - - * 1 If 3 points are collinear - * 0 If they are not - */ -int -rt_3pts_collinear(a, b, c, tol) -point_t a, b, c; -struct rt_tol *tol; -{ - fastf_t mag_ab, mag_bc, mag_ca, max_len, dist_sq; - fastf_t cos_a, cos_b, cos_c; - vect_t ab, bc, ca; - int max_edge_no; - - VSUB2(ab, b, a); - VSUB2(bc, c, b); - VSUB2(ca, a, c); - mag_ab = MAGNITUDE(ab); - mag_bc = MAGNITUDE(bc); - mag_ca = MAGNITUDE(ca); - - /* find longest edge */ - max_len = mag_ab; - max_edge_no = 1; - - if( mag_bc > max_len ) - { - max_len = mag_bc; - max_edge_no = 2; - } - - if( mag_ca > max_len ) - { - max_len = mag_ca; - max_edge_no = 3; - } - - switch( max_edge_no ) - { - default: - case 1: - cos_b = (-VDOT( ab , bc ))/(mag_ab * mag_bc ); - dist_sq = mag_bc*mag_bc*( 1.0 - cos_b*cos_b); - break; - case 2: - cos_c = (-VDOT( bc , ca ))/(mag_bc * mag_ca ); - dist_sq = mag_ca*mag_ca*(1.0 - cos_c*cos_c); - break; - case 3: - cos_a = (-VDOT( ca , ab ))/(mag_ca * mag_ab ); - dist_sq = mag_ab*mag_ab*(1.0 - cos_a*cos_a); - break; - } - - if( dist_sq <= tol->dist_sq ) - return( 1 ); - else - return( 0 ); -} - - -/* - * R T _ 3 P T S _ D I S T I N C T - * - * Check to see if three points are all distinct, i.e., - * ensure that there is at least sqrt(dist_tol_sq) distance - * between every pair of points. - * - * Returns (boolean) - - * 1 If all three points are distinct - * 0 If two or more points are closer together than dist_tol_sq - */ -int -rt_3pts_distinct( a, b, c, tol ) -CONST point_t a, b, c; -CONST struct rt_tol *tol; -{ - vect_t B_A; - vect_t C_A; - vect_t C_B; - - RT_CK_TOL(tol); - VSUB2( B_A, b, a ); - if( MAGSQ( B_A ) <= tol->dist_sq ) return(0); - VSUB2( C_A, c, a ); - if( MAGSQ( C_A ) <= tol->dist_sq ) return(0); - VSUB2( C_B, c, b ); - if( MAGSQ( C_B ) <= tol->dist_sq ) return(0); - return(1); -} - -/* - * R T _ M K _ P L A N E _ 3 P T S - * - * Find the equation of a plane that contains three points. - * Note that normal vector created is expected to point out (see vmath.h), - * so the vector from A to C had better be counter-clockwise - * (about the point A) from the vector from A to B. - * This follows the BRL-CAD outward-pointing normal convention, and the - * right-hand rule for cross products. - * - * - * C - * * - * |\ - * | \ - * ^ N | \ - * | \ | \ - * | \ | \ - * |C-A \ | \ - * | \ | \ - * | \ | \ - * \| \ - * *---------* - * A B - * -----> - * B-A - * - * If the points are given in the order A B C (eg, *counter*-clockwise), - * then the outward pointing surface normal N = (B-A) x (C-A). - * - * Explicit Return - - * 0 OK - * -1 Failure. At least two of the points were not distinct, - * or all three were colinear. - * - * Implicit Return - - * plane The plane equation is stored here. - */ -int -rt_mk_plane_3pts( plane, a, b, c, tol ) -plane_t plane; -CONST point_t a, b, c; -CONST struct rt_tol *tol; -{ - vect_t B_A; - vect_t C_A; - vect_t C_B; - register fastf_t mag; - - RT_CK_TOL(tol); - - VSUB2( B_A, b, a ); - if( MAGSQ( B_A ) <= tol->dist_sq ) return(-1); - VSUB2( C_A, c, a ); - if( MAGSQ( C_A ) <= tol->dist_sq ) return(-1); - VSUB2( C_B, c, b ); - if( MAGSQ( C_B ) <= tol->dist_sq ) return(-1); - - VCROSS( plane, B_A, C_A ); - - /* Ensure unit length normal */ - if( (mag = MAGNITUDE(plane)) <= SMALL_FASTF ) - return(-1); /* FAIL */ - mag = 1/mag; - VSCALE( plane, plane, mag ); - - /* Find distance from the origin to the plane */ - /* XXX Should do with pt that has smallest magnitude (closest to origin) */ - plane[3] = VDOT( plane, a ); - -#if 0 - /* Check to be sure that angle between A-Origin and N vect < 89 degrees */ - /* XXX Could complain for pts on axis-aligned plane, far from origin */ - mag = MAGSQ(a); - if( mag > tol->dist_sq ) { - /* cos(89 degrees) = 0.017452406, reciprocal is 57.29 */ - if( plane[3]/sqrt(mag) < 0.017452406 ) { - rt_log("rt_mk_plane_3pts() WARNING: plane[3] value is suspect\n"); - } - } -#endif - return(0); /* OK */ -} - -/* - * R T _ M K P O I N T _ 3 P L A N E S - * - * Given the description of three planes, compute the point of intersection, - * if any. - * - * Find the solution to a system of three equations in three unknowns: - * - * Px * Ax + Py * Ay + Pz * Az = -A3; - * Px * Bx + Py * By + Pz * Bz = -B3; - * Px * Cx + Py * Cy + Pz * Cz = -C3; - * - * or - * - * [ Ax Ay Az ] [ Px ] [ -A3 ] - * [ Bx By Bz ] * [ Py ] = [ -B3 ] - * [ Cx Cy Cz ] [ Pz ] [ -C3 ] - * - * - * Explitic Return - - * 0 OK - * -1 Failure. Intersection is a line or plane. - * - * Implicit Return - - * pt The point of intersection is stored here. - */ -int -rt_mkpoint_3planes( pt, a, b, c ) -point_t pt; -CONST plane_t a, b, c; -{ - vect_t v1, v2, v3; - register fastf_t det; - - /* Find a vector perpendicular to both planes B and C */ - VCROSS( v1, b, c ); - - /* If that vector is perpendicular to A, - * then A is parallel to either B or C, and no intersection exists. - * This dot&cross product is the determinant of the matrix M. - * (I suspect there is some deep significance to this!) - */ - det = VDOT( a, v1 ); - if( NEAR_ZERO( det, SMALL_FASTF ) ) return(-1); - - VCROSS( v2, a, c ); - VCROSS( v3, a, b ); - - det = 1/det; - pt[X] = det*(a[3]*v1[X] - b[3]*v2[X] + c[3]*v3[X]); - pt[Y] = det*(a[3]*v1[Y] - b[3]*v2[Y] + c[3]*v3[Y]); - pt[Z] = det*(a[3]*v1[Z] - b[3]*v2[Z] + c[3]*v3[Z]); - return(0); -} - -/* - * R T _ 2 L I N E 3 _ C O L I N E A R - * - * Returns non-zero if the 3 lines are colinear to within tol->dist - * over the given distance range. - * - * Range should be at least one model diameter for most applications. - * 1e5 might be OK for a default for "vehicle sized" models. - * - * The direction vectors do not need to be unit length. - */ -int -rt_2line3_colinear( p1, d1, p2, d2, range, tol ) -CONST point_t p1; -CONST vect_t d1; -CONST point_t p2; -CONST vect_t d2; -double range; -CONST struct rt_tol *tol; -{ - fastf_t mag1; - fastf_t mag2; - point_t tail; - - RT_CK_TOL(tol); - - if( (mag1 = MAGNITUDE(d1)) < SMALL_FASTF ) rt_bomb("rt_2line3_colinear() mag1 zero\n"); - if( (mag2 = MAGNITUDE(d2)) < SMALL_FASTF ) rt_bomb("rt_2line3_colinear() mag2 zero\n"); - - /* Impose a general angular tolerance to reject "obviously" non-parallel lines */ - /* tol->para and RT_DOT_TOL are too tight a tolerance. 0.1 is 5 degrees */ - if( fabs(VDOT(d1, d2)) < 0.9 * mag1 * mag2 ) goto fail; - - /* See if start points are within tolerance of other line */ - if( rt_distsq_line3_pt3( p1, d1, p2 ) > tol->dist_sq ) goto fail; - if( rt_distsq_line3_pt3( p2, d2, p1 ) > tol->dist_sq ) goto fail; - - VJOIN1( tail, p1, range/mag1, d1 ); - if( rt_distsq_line3_pt3( p2, d2, tail ) > tol->dist_sq ) goto fail; - - VJOIN1( tail, p2, range/mag2, d2 ); - if( rt_distsq_line3_pt3( p1, d1, tail ) > tol->dist_sq ) goto fail; - - if( rt_g.debug & DEBUG_MATH ) { - rt_log("rt_2line3colinear(range=%g) ret=1\n",range); - } - return 1; -fail: - if( rt_g.debug & DEBUG_MATH ) { - rt_log("rt_2line3colinear(range=%g) ret=0\n",range); - } - return 0; -} - -/* - * R T _ I S E C T _ L I N E 3 _ P L A N E - * - * Intersect an infinite line (specified in point and direction vector form) - * with a plane that has an outward pointing normal. - * The direction vector need not have unit length. - * - * Explicit Return - - * -2 missed (ray is outside halfspace) - * -1 missed (ray is inside) - * 0 line lies on plane - * 1 hit (ray is entering halfspace) - * 2 hit (ray is leaving) - * - * Implicit Return - - * The value at *dist is set to the parametric distance of the intercept - */ -int -rt_isect_line3_plane( dist, pt, dir, plane, tol ) -fastf_t *dist; -CONST point_t pt; -CONST vect_t dir; -CONST plane_t plane; -CONST struct rt_tol *tol; -{ - register fastf_t slant_factor; - register fastf_t norm_dist; - - RT_CK_TOL(tol); - - norm_dist = plane[3] - VDOT( plane, pt ); - slant_factor = VDOT( plane, dir ); - - if( slant_factor < -SMALL_FASTF ) { - *dist = norm_dist/slant_factor; - return 1; /* HIT, entering */ - } else if( slant_factor > SMALL_FASTF ) { - *dist = norm_dist/slant_factor; - return 2; /* HIT, leaving */ - } - - /* - * Ray is parallel to plane when dir.N == 0. - */ - *dist = 0; /* sanity */ - if( norm_dist < -tol->dist ) - return -2; /* missed, outside */ - if( norm_dist > tol->dist ) - return -1; /* missed, inside */ - return 0; /* Ray lies in the plane */ -} - -/* - * R T _ I S E C T _ 2 P L A N E S - * - * Given two planes, find the line of intersection between them, - * if one exists. - * The line of intersection is returned in parametric line - * (point & direction vector) form. - * - * In order that all the geometry under consideration be in "front" - * of the ray, it is necessary to pass the minimum point of the model - * RPP. If this convention is unnecessary, just pass (0,0,0) as rpp_min. - * - * Explicit Return - - * 0 OK, line of intersection stored in `pt' and `dir'. - * -1 FAIL, planes are identical (co-planar) - * -2 FAIL, planes are parallel and distinct - * -3 FAIL, unable to find line of intersection - * - * Implicit Returns - - * pt Starting point of line of intersection - * dir Direction vector of line of intersection (unit length) - */ -int -rt_isect_2planes( pt, dir, a, b, rpp_min, tol ) -point_t pt; -vect_t dir; -CONST plane_t a; -CONST plane_t b; -CONST vect_t rpp_min; -CONST struct rt_tol *tol; -{ - LOCAL vect_t abs_dir; - LOCAL plane_t pl; - int i; - - if( (i = rt_coplanar( a, b, tol )) != 0 ) { - if( i > 0 ) - return(-1); /* FAIL -- coplanar */ - return(-2); /* FAIL -- parallel & distinct */ - } - - /* Direction vector for ray is perpendicular to both plane normals */ - VCROSS( dir, a, b ); - VUNITIZE( dir ); /* safety */ - - /* - * Select an axis-aligned plane which has it's normal pointing - * along the same axis as the largest magnitude component of - * the direction vector. - * If the largest magnitude component is negative, reverse the - * direction vector, so that model is "in front" of start point. - */ - abs_dir[X] = (dir[X] >= 0) ? dir[X] : (-dir[X]); - abs_dir[Y] = (dir[Y] >= 0) ? dir[Y] : (-dir[Y]); - abs_dir[Z] = (dir[Z] >= 0) ? dir[Z] : (-dir[Z]); - - if( abs_dir[X] >= abs_dir[Y] ) { - if( abs_dir[X] >= abs_dir[Z] ) { - VSET( pl, 1, 0, 0 ); /* X */ - pl[3] = rpp_min[X]; - if( dir[X] < 0 ) { - VREVERSE( dir, dir ); - } - } else { - VSET( pl, 0, 0, 1 ); /* Z */ - pl[3] = rpp_min[Z]; - if( dir[Z] < 0 ) { - VREVERSE( dir, dir ); - } - } - } else { - if( abs_dir[Y] >= abs_dir[Z] ) { - VSET( pl, 0, 1, 0 ); /* Y */ - pl[3] = rpp_min[Y]; - if( dir[Y] < 0 ) { - VREVERSE( dir, dir ); - } - } else { - VSET( pl, 0, 0, 1 ); /* Z */ - pl[3] = rpp_min[Z]; - if( dir[Z] < 0 ) { - VREVERSE( dir, dir ); - } - } - } - - /* Intersection of the 3 planes defines ray start point */ - if( rt_mkpoint_3planes( pt, pl, a, b ) < 0 ) - return(-3); /* FAIL -- no intersection */ - - return(0); /* OK */ -} - -/* - * R T _ I S E C T _ L I N E 2 _ L I N E 2 - * - * Intersect two lines, each in given in parametric form: - * - * X = P + t * D - * and - * X = A + u * C - * - * While the parametric form is usually used to denote a ray - * (ie, positive values of the parameter only), in this case - * the full line is considered. - * - * The direction vectors C and D need not have unit length. - * - * Explicit Return - - * -1 no intersection, lines are parallel. - * 0 lines are co-linear - * dist[0] gives distance from P to A, - * dist[1] gives distance from P to (A+C) [not same as below] - * 1 intersection found (t and u returned) - * dist[0] gives distance from P to isect, - * dist[1] gives distance from A to isect. - * - * Implicit Returns - - * When explicit return > 0, dist[0] and dist[1] are the - * line parameters of the intersection point on the 2 rays. - * The actual intersection coordinates can be found by - * substituting either of these into the original ray equations. - * - * Note that for lines which are very nearly parallel, but not - * quite parallel enough to have the determinant go to "zero", - * the intersection can turn up in surprising places. - * (e.g. when det=1e-15 and det1=5.5e-17, t=0.5) - */ -int -rt_isect_line2_line2( dist, p, d, a, c, tol ) -fastf_t *dist; /* dist[2] */ -CONST point_t p; -CONST vect_t d; -CONST point_t a; -CONST vect_t c; -CONST struct rt_tol *tol; -{ - fastf_t hx, hy; /* A - P */ - register fastf_t det; - register fastf_t det1; - vect_t unit_d; - vect_t unit_c; - vect_t unit_h; - fastf_t dot; - int parallel; - int parallel1; - - RT_CK_TOL(tol); - if( rt_g.debug & DEBUG_MATH ) { - rt_log("rt_isect_line2_line2() p=(%g,%g), d=(%g,%g)\n\t\t\ta=(%g,%g), c=(%g,%g)\n", - V2ARGS(p), V2ARGS(d), V2ARGS(a), V2ARGS(c) ); - } - - /* - * From the two components q and r, form a system - * of 2 equations in 2 unknowns. - * Solve for t and u in the system: - * - * Px + t * Dx = Ax + u * Cx - * Py + t * Dy = Ay + u * Cy - * or - * t * Dx - u * Cx = Ax - Px - * t * Dy - u * Cy = Ay - Py - * - * Let H = A - P, resulting in: - * - * t * Dx - u * Cx = Hx - * t * Dy - u * Cy = Hy - * - * or - * - * [ Dx -Cx ] [ t ] [ Hx ] - * [ ] * [ ] = [ ] - * [ Dy -Cy ] [ u ] [ Hy ] - * - * This system can be solved by direct substitution, or by - * finding the determinants by Cramers rule: - * - * [ Dx -Cx ] - * det(M) = det [ ] = -Dx * Cy + Cx * Dy - * [ Dy -Cy ] - * - * If det(M) is zero, then the lines are parallel (perhaps colinear). - * Otherwise, exactly one solution exists. - */ - det = c[X] * d[Y] - d[X] * c[Y]; - - /* - * det(M) is non-zero, so there is exactly one solution. - * Using Cramer's rule, det1(M) replaces the first column - * of M with the constant column vector, in this case H. - * Similarly, det2(M) replaces the second column. - * Computation of the determinant is done as before. - * - * Now, - * - * [ Hx -Cx ] - * det [ ] - * det1(M) [ Hy -Cy ] -Hx * Cy + Cx * Hy - * t = ------- = --------------- = ------------------ - * det(M) det(M) -Dx * Cy + Cx * Dy - * - * and - * - * [ Dx Hx ] - * det [ ] - * det2(M) [ Dy Hy ] Dx * Hy - Hx * Dy - * u = ------- = --------------- = ------------------ - * det(M) det(M) -Dx * Cy + Cx * Dy - */ - hx = a[X] - p[X]; - hy = a[Y] - p[Y]; - det1 = (c[X] * hy - hx * c[Y]); - - unit_d[0] = d[0]; - unit_d[1] = d[1]; - unit_d[2] = 0.0; - VUNITIZE( unit_d ); - unit_c[0] = c[0]; - unit_c[1] = c[1]; - unit_c[2] = 0.0; - VUNITIZE( unit_c ); - unit_h[0] = hx; - unit_h[1] = hy; - unit_h[2] = 0.0; - VUNITIZE( unit_h ); - - if( fabs( VDOT( unit_d, unit_c ) ) >= tol->para ) - parallel = 1; - else - parallel = 0; - - if( fabs( VDOT( unit_h, unit_c ) ) >= tol->para ) - parallel1 = 1; - else - parallel1 = 0; - - /* XXX This zero tolerance here should actually be - * XXX determined by something like - * XXX max(c[X], c[Y], d[X], d[Y]) / MAX_FASTF_DYNAMIC_RANGE - * XXX In any case, nothing smaller than 1e-16 - */ -#define DETERMINANT_TOL 1.0e-14 /* XXX caution on non-IEEE machines */ - if( parallel || NEAR_ZERO( det, DETERMINANT_TOL ) ) { - /* Lines are parallel */ - if( !parallel1 && !NEAR_ZERO( det1, DETERMINANT_TOL ) ) { - /* Lines are NOT co-linear, just parallel */ - if( rt_g.debug & DEBUG_MATH ) { - rt_log("\tparallel, not co-linear. det=%e, det1=%g\n", det, det1); - } - return -1; /* parallel, no intersection */ - } - - /* - * Lines are co-linear. - * Determine t as distance from P to A. - * Determine u as distance from P to (A+C). [special!] - * Use largest direction component, for numeric stability - * (and avoiding division by zero). - */ - if( fabs(d[X]) >= fabs(d[Y]) ) { - dist[0] = hx/d[X]; - dist[1] = (hx + c[X]) / d[X]; - } else { - dist[0] = hy/d[Y]; - dist[1] = (hy + c[Y]) / d[Y]; - } - if( rt_g.debug & DEBUG_MATH ) { - rt_log("\tcolinear, t = %g, u = %g\n", dist[0], dist[1] ); - } - return 0; /* Lines co-linear */ - } - if( rt_g.debug & DEBUG_MATH ) { - /* XXX This print is temporary */ -rt_log("\thx=%g, hy=%g, det=%g, det1=%g, det2=%g\n", hx, hy, det, det1, (d[X] * hy - hx * d[Y]) ); - } - det = 1/det; - dist[0] = det * det1; - dist[1] = det * (d[X] * hy - hx * d[Y]); - if( rt_g.debug & DEBUG_MATH ) { - rt_log("\tintersection, t = %g, u = %g\n", dist[0], dist[1] ); - } - -#if 0 - /* XXX This isn't any good. - * 1) Sometimes, dist[0] is very large. Only caller can tell whether - * that is useful to him or not. - * 2) Sometimes, the difference between the two hit points is - * not much more than tol->dist. Either hit point is perfectly - * good; the caller just needs to be careful and not use *both*. - */ - { - point_t hit1, hit2; - vect_t diff; - fastf_t dist_sq; - - VJOIN1_2D( hit1, p, dist[0], d ); - VJOIN1_2D( hit2, a, dist[1], c ); - VSUB2_2D( diff, hit1, hit2 ); - dist_sq = MAGSQ_2D( diff ); - if( dist_sq >= tol->dist_sq ) { - if( rt_g.debug & DEBUG_MATH || dist_sq < 100*tol->dist_sq ) { - rt_log("rt_isect_line2_line2(): dist=%g >%g, inconsistent solution, hit1=(%g,%g), hit2=(%g,%g)\n", - sqrt(dist_sq), tol->dist, - hit1[X], hit1[Y], hit2[X], hit2[Y]); - } - return -2; /* s/b -1? */ - } - } -#endif - - return 1; /* Intersection found */ -} - -/* - * R T _ I S E C T _ L I N E 2 _ L S E G 2 - * - * Intersect a line in parametric form: - * - * X = P + t * D - * - * with a line segment defined by two distinct points A and B=(A+C). - * - * XXX probably should take point B, not vector C. Sigh. - * - * Explicit Return - - * -4 A and B are not distinct points - * -3 Lines do not intersect - * -2 Intersection exists, but outside segemnt, < A - * -1 Intersection exists, but outside segment, > B - * 0 Lines are co-linear (special meaning of dist[1]) - * 1 Intersection at vertex A - * 2 Intersection at vertex B (A+C) - * 3 Intersection between A and B - * - * Implicit Returns - - * t When explicit return >= 0, t is the parameter that describes - * the intersection of the line and the line segment. - * The actual intersection coordinates can be found by - * solving P + t * D. However, note that for return codes - * 1 and 2 (intersection exactly at a vertex), it is - * strongly recommended that the original values passed in - * A or B are used instead of solving P + t * D, to prevent - * numeric error from creeping into the position of - * the endpoints. - */ -int -rt_isect_line2_lseg2( dist, p, d, a, c, tol ) -fastf_t *dist; /* dist[2] */ -CONST point_t p; -CONST vect_t d; -CONST point_t a; -CONST vect_t c; -CONST struct rt_tol *tol; -{ - register fastf_t f; - fastf_t ctol; - int ret; - point_t b; - - RT_CK_TOL(tol); - if( rt_g.debug & DEBUG_MATH ) { - rt_log("rt_isect_line2_lseg2() p=(%g,%g), pdir=(%g,%g)\n\t\t\ta=(%g,%g), adir=(%g,%g)\n", - V2ARGS(p), V2ARGS(d), V2ARGS(a), V2ARGS(c) ); - } - - /* - * To keep the values of u between 0 and 1, - * C should NOT be scaled to have unit length. - * However, it is a good idea to make sure that - * C is a non-zero vector, (ie, that A and B are distinct). - */ - if( (ctol = MAGSQ_2D(c)) <= tol->dist_sq ) { - ret = -4; /* points A and B are not distinct */ - goto out; - } - - /* - * Detecting colinearity is difficult, and very very important. - * As a first step, check to see if both points A and B lie - * within tolerance of the line. If so, then the line segment AC - * is ON the line. - */ - VADD2_2D( b, a, c ); - if( rt_distsq_line2_point2( p, d, a ) <= tol->dist_sq && - (ctol=rt_distsq_line2_point2( p, d, b )) <= tol->dist_sq ) { - if( rt_g.debug & DEBUG_MATH ) { -rt_log("b=(%g, %g), b_dist_sq=%g\n", V2ARGS(b), ctol); - rt_log("rt_isect_line2_lseg2() pts A and B within tol of line\n"); - } - /* Find the parametric distance along the ray */ - dist[0] = rt_dist_pt2_along_line2( p, d, a ); - dist[1] = rt_dist_pt2_along_line2( p, d, b ); - ret = 0; /* Colinear */ - goto out; - } - - if( (ret = rt_isect_line2_line2( dist, p, d, a, c, tol )) < 0 ) { - /* Lines are parallel, non-colinear */ - ret = -3; /* No intersection found */ - goto out; - } - if( ret == 0 ) { - fastf_t dtol; - /* Lines are colinear */ - /* If P within tol of either endpoint (0, 1), make exact. */ - dtol = tol->dist / sqrt(MAGSQ_2D(d)); - if( rt_g.debug & DEBUG_MATH ) { - rt_log("rt_isect_line2_lseg2() dtol=%g, dist[0]=%g, dist[1]=%g\n", - dtol, dist[0], dist[1]); - } - if( dist[0] > -dtol && dist[0] < dtol ) dist[0] = 0; - else if( dist[0] > 1-dtol && dist[0] < 1+dtol ) dist[0] = 1; - - if( dist[1] > -dtol && dist[1] < dtol ) dist[1] = 0; - else if( dist[1] > 1-dtol && dist[1] < 1+dtol ) dist[1] = 1; - ret = 0; /* Colinear */ - goto out; - } - - /* - * The two lines are claimed to intersect at a point. - * First, validate that hit point represented by dist[0] - * is in fact on and between A--B. - * (Nearly parallel lines can result in odd situations here). - * The performance hit of doing this is vastly preferable - * to returning wrong answers. Know a faster algorithm? - */ - { - fastf_t ab_dist = 0; - point_t hit_pt; - point_t hit2; - - VJOIN1_2D( hit_pt, p, dist[0], d ); - VJOIN1_2D( hit2, a, dist[1], c ); - /* Check both hit point value calculations */ - if( rt_pt2_pt2_equal( a, hit_pt, tol ) || - rt_pt2_pt2_equal( a, hit2, tol ) ) { - dist[1] = 0; - ret = 1; /* Intersect is at A */ - } - if( rt_pt2_pt2_equal( b, hit_pt, tol ) || - rt_pt2_pt2_equal( b, hit_pt, tol ) ) { - dist[1] = 1; - ret = 2; /* Intersect is at B */ - } - - ret = rt_isect_pt2_lseg2( &ab_dist, a, b, hit_pt, tol ); - if( rt_g.debug & DEBUG_MATH ) { - /* XXX This is temporary */ - V2PRINT("a", a); - V2PRINT("hit", hit_pt); - V2PRINT("b", b); -rt_log("rt_isect_pt2_lseg2() hit2d=(%g,%g) ab_dist=%g, ret=%d\n", hit_pt[X], hit_pt[Y], ab_dist, ret); -rt_log("\tother hit2d=(%g,%g)\n", hit2[X], hit2[Y] ); - } - if( ret <= 0 ) { - if( ab_dist < 0 ) { - ret = -2; /* Intersection < A */ - } else { - ret = -1; /* Intersection >B */ - } - goto out; - } - if( ret == 1 ) { - dist[1] = 0; - ret = 1; /* Intersect is at A */ - goto out; - } - if( ret == 2 ) { - dist[1] = 1; - ret = 2; /* Intersect is at B */ - goto out; - } - /* ret == 3, hit_pt is between A and B */ - - if( !rt_between( a[X], hit_pt[X], b[X], tol ) || - !rt_between( a[Y], hit_pt[Y], b[Y], tol ) ) { - rt_bomb("rt_isect_line2_lseg2() hit_pt not between A and B!\n"); - } - } - - /* - * If the dist[1] parameter is outside the range (0..1), - * reject the intersection, because it falls outside - * the line segment A--B. - * - * Convert the tol->dist into allowable deviation in terms of - * (0..1) range of the parameters. - */ - ctol = tol->dist / sqrt(ctol); - if( rt_g.debug & DEBUG_MATH ) { - rt_log("rt_isect_line2_lseg2() ctol=%g, dist[1]=%g\n", ctol, dist[1]); - } - if( dist[1] < -ctol ) { - ret = -2; /* Intersection < A */ - goto out; - } - if( (f=(dist[1]-1)) > ctol ) { - ret = -1; /* Intersection > B */ - goto out; - } - - /* Check for ctoly intersection with one of the verticies */ - if( dist[1] < ctol ) { - dist[1] = 0; - ret = 1; /* Intersection at A */ - goto out; - } - if( f >= -ctol ) { - dist[1] = 1; - ret = 2; /* Intersection at B */ - goto out; - } - ret = 3; /* Intersection between A and B */ -out: - if( rt_g.debug & DEBUG_MATH ) { - rt_log("rt_isect_line2_lseg2() dist[0]=%g, dist[1]=%g, ret=%d\n", - dist[0], dist[1], ret); - } - return ret; -} - -/* - * R T _ I S E C T _ L S E G 2 _ L S E G 2 - * - * Intersect two 2D line segments, defined by two points and two vectors. - * The vectors are unlikely to be unit length. - * - * Explicit Return - - * -2 missed (line segments are parallel) - * -1 missed (colinear and non-overlapping) - * 0 hit (line segments colinear and overlapping) - * 1 hit (normal intersection) - * - * Implicit Return - - * The value at dist[] is set to the parametric distance of the intercept - * dist[0] is parameter along p, range 0 to 1, to intercept. - * dist[1] is parameter along q, range 0 to 1, to intercept. - * If within distance tolerance of the endpoints, these will be - * exactly 0.0 or 1.0, to ease the job of caller. - * - * Special note: when return code is "0" for co-linearity, dist[1] has - * an alternate interpretation: it's the parameter along p (not q) - * which takes you from point p to the point (q + qdir), i.e., it's - * the endpoint of the q linesegment, since in this case there may be - * *two* intersections, if q is contained within span p to (p + pdir). - * And either may be -10 if the point is outside the span. - */ -int -rt_isect_lseg2_lseg2( dist, p, pdir, q, qdir, tol ) -fastf_t *dist; -CONST point_t p; -CONST vect_t pdir; -CONST point_t q; -CONST vect_t qdir; -struct rt_tol *tol; -{ - fastf_t ptol, qtol; /* length in parameter space == tol->dist */ - int status; - - RT_CK_TOL(tol); - if( rt_g.debug & DEBUG_MATH ) { - rt_log("rt_isect_lseg2_lseg2() p=(%g,%g), pdir=(%g,%g)\n\t\tq=(%g,%g), qdir=(%g,%g)\n", - V2ARGS(p), V2ARGS(pdir), V2ARGS(q), V2ARGS(qdir) ); - } - - status = rt_isect_line2_line2( dist, p, pdir, q, qdir, tol ); - if( status < 0 ) { - /* Lines are parallel, non-colinear */ - return -1; /* No intersection */ - } - if( status == 0 ) { - int nogood = 0; - /* Lines are colinear */ - /* If P within tol of either endpoint (0, 1), make exact. */ - ptol = tol->dist / sqrt( MAGSQ_2D(pdir) ); - if( rt_g.debug & DEBUG_MATH ) { - rt_log("ptol=%g\n", ptol); - } - if( dist[0] > -ptol && dist[0] < ptol ) dist[0] = 0; - else if( dist[0] > 1-ptol && dist[0] < 1+ptol ) dist[0] = 1; - - if( dist[1] > -ptol && dist[1] < ptol ) dist[1] = 0; - else if( dist[1] > 1-ptol && dist[1] < 1+ptol ) dist[1] = 1; - - if( dist[1] < 0 || dist[1] > 1 ) nogood = 1; - if( dist[0] < 0 || dist[0] > 1 ) nogood++; - if( nogood >= 2 ) - return -1; /* colinear, but not overlapping */ - if( rt_g.debug & DEBUG_MATH ) { - rt_log(" HIT colinear!\n"); - } - return 0; /* colinear and overlapping */ - } - /* Lines intersect */ - /* If within tolerance of an endpoint (0, 1), make exact. */ - ptol = tol->dist / sqrt( MAGSQ_2D(pdir) ); - if( dist[0] > -ptol && dist[0] < ptol ) dist[0] = 0; - else if( dist[0] > 1-ptol && dist[0] < 1+ptol ) dist[0] = 1; - - qtol = tol->dist / sqrt( MAGSQ_2D(qdir) ); - if( dist[1] > -qtol && dist[1] < qtol ) dist[1] = 0; - else if( dist[1] > 1-qtol && dist[1] < 1+qtol ) dist[1] = 1; - - if( rt_g.debug & DEBUG_MATH ) { - rt_log("ptol=%g, qtol=%g\n", ptol, qtol); - } - if( dist[0] < 0 || dist[0] > 1 || dist[1] < 0 || dist[1] > 1 ) { - if( rt_g.debug & DEBUG_MATH ) { - rt_log(" MISS\n"); - } - return -1; /* missed */ - } - if( rt_g.debug & DEBUG_MATH ) { - rt_log(" HIT!\n"); - } - return 1; /* hit, normal intersection */ -} - -/* - * R T _ I S E C T _ L S E G 3 _ L S E G 3 - * - * Intersect two 3D line segments, defined by two points and two vectors. - * The vectors are unlikely to be unit length. - * - * Explicit Return - - * -2 missed (line segments are parallel) - * -1 missed (colinear and non-overlapping) - * 0 hit (line segments colinear and overlapping) - * 1 hit (normal intersection) - * - * Implicit Return - - * The value at dist[] is set to the parametric distance of the intercept - * dist[0] is parameter along p, range 0 to 1, to intercept. - * dist[1] is parameter along q, range 0 to 1, to intercept. - * If within distance tolerance of the endpoints, these will be - * exactly 0.0 or 1.0, to ease the job of caller. - * - * Special note: when return code is "0" for co-linearity, dist[1] has - * an alternate interpretation: it's the parameter along p (not q) - * which takes you from point p to the point (q + qdir), i.e., it's - * the endpoint of the q linesegment, since in this case there may be - * *two* intersections, if q is contained within span p to (p + pdir). - * And either may be -10 if the point is outside the span. - */ -int -rt_isect_lseg3_lseg3( dist, p, pdir, q, qdir, tol ) -fastf_t *dist; -CONST point_t p; -CONST vect_t pdir; -CONST point_t q; -CONST vect_t qdir; -struct rt_tol *tol; -{ - fastf_t ptol, qtol; /* length in parameter space == tol->dist */ - fastf_t pmag, qmag; - int status; - - RT_CK_TOL(tol); - if( rt_g.debug & DEBUG_MATH ) { - rt_log("rt_isect_lseg3_lseg3() p=(%g,%g), pdir=(%g,%g)\n\t\tq=(%g,%g), qdir=(%g,%g)\n", - V2ARGS(p), V2ARGS(pdir), V2ARGS(q), V2ARGS(qdir) ); - } - - status = rt_isect_line3_line3( &dist[0], &dist[1], p, pdir, q, qdir, tol ); - if( status < 0 ) { - /* Lines are parallel, non-colinear */ - return -1; /* No intersection */ - } - pmag = MAGNITUDE(pdir); - if( pmag < SMALL_FASTF ) - rt_bomb("rt_isect_lseg3_lseg3: |p|=0\n"); - if( status == 0 ) { - int nogood = 0; - /* Lines are colinear */ - /* If P within tol of either endpoint (0, 1), make exact. */ - ptol = tol->dist / pmag; - if( rt_g.debug & DEBUG_MATH ) { - rt_log("ptol=%g\n", ptol); - } - if( dist[0] > -ptol && dist[0] < ptol ) dist[0] = 0; - else if( dist[0] > 1-ptol && dist[0] < 1+ptol ) dist[0] = 1; - - if( dist[1] > -ptol && dist[1] < ptol ) dist[1] = 0; - else if( dist[1] > 1-ptol && dist[1] < 1+ptol ) dist[1] = 1; - - if( dist[1] < 0 || dist[1] > 1 ) nogood = 1; - if( dist[0] < 0 || dist[0] > 1 ) nogood++; - if( nogood >= 2 ) - return -1; /* colinear, but not overlapping */ - if( rt_g.debug & DEBUG_MATH ) { - rt_log(" HIT colinear!\n"); - } - return 0; /* colinear and overlapping */ - } - /* Lines intersect */ - /* If within tolerance of an endpoint (0, 1), make exact. */ - ptol = tol->dist / pmag; - if( dist[0] > -ptol && dist[0] < ptol ) dist[0] = 0; - else if( dist[0] > 1-ptol && dist[0] < 1+ptol ) dist[0] = 1; - - qmag = MAGNITUDE(qdir); - if( qmag < SMALL_FASTF ) - rt_bomb("rt_isect_lseg3_lseg3: |q|=0\n"); - qtol = tol->dist / qmag; - if( dist[1] > -qtol && dist[1] < qtol ) dist[1] = 0; - else if( dist[1] > 1-qtol && dist[1] < 1+qtol ) dist[1] = 1; - - if( rt_g.debug & DEBUG_MATH ) { - rt_log("ptol=%g, qtol=%g\n", ptol, qtol); - } - if( dist[0] < 0 || dist[0] > 1 || dist[1] < 0 || dist[1] > 1 ) { - if( rt_g.debug & DEBUG_MATH ) { - rt_log(" MISS\n"); - } - return -1; /* missed */ - } - if( rt_g.debug & DEBUG_MATH ) { - rt_log(" HIT!\n"); - } - return 1; /* hit, normal intersection */ -} - -/* - * R T _ I S E C T _ L I N E 3 _ L I N E 3 - * - * Intersect two lines, each in given in parametric form: - * - * X = P + t * D - * and - * X = A + u * C - * - * While the parametric form is usually used to denote a ray - * (ie, positive values of the parameter only), in this case - * the full line is considered. - * - * The direction vectors C and D need not have unit length. - * - * Explicit Return - - * -2 no intersection, lines are parallel. - * -1 no intersection - * 0 lines are co-linear (t returned for u=0 to give distance to A) - * 1 intersection found (t and u returned) - * - * Implicit Returns - - * - * t,u When explicit return >= 0, t and u are the - * line parameters of the intersection point on the 2 rays. - * The actual intersection coordinates can be found by - * substituting either of these into the original ray equations. - * - * XXX It would be sensible to change the t,u pair to dist[2]. - */ -int -rt_isect_line3_line3( t, u, p, d, a, c, tol ) -fastf_t *t; -fastf_t *u; -CONST point_t p; -CONST vect_t d; -CONST point_t a; -CONST vect_t c; -CONST struct rt_tol *tol; -{ - LOCAL vect_t n; - LOCAL vect_t abs_n; - LOCAL vect_t h; - register fastf_t det; - register fastf_t det1; - register short int q,r,s; - - RT_CK_TOL(tol); - - /* - * Any intersection will occur in the plane with surface - * normal D cross C, which may not have unit length. - * The plane containing the two lines will be a constant - * distance from a plane with the same normal that contains - * the origin. Therfore, the projection of any point on the - * plane along N has the same length. - * Verify that this holds for P and A. - * If N dot P != N dot A, there is no intersection, because - * P and A must lie on parallel planes that are different - * distances from the origin. - */ - VCROSS( n, d, c ); - det = VDOT( n, p ) - VDOT( n, a ); - if( !NEAR_ZERO( det, tol->dist ) ) { - return(-1); /* No intersection */ - } - - if( NEAR_ZERO( MAGSQ( n ) , SMALL_FASTF ) ) - { - vect_t a_to_p; - - /* lines are parallel, must find another way to get normal vector */ - VSUB2( a_to_p , p , a ); - VCROSS( n , a_to_p , d ); - - /* if normal still has zero length, then lines are parallel and collinear - * and the following code will work OK */ - } - - /* - * Solve for t and u in the system: - * - * Px + t * Dx = Ax + u * Cx - * Py + t * Dy = Ay + u * Cy - * Pz + t * Dz = Az + u * Cz - * - * This system is over-determined, having 3 equations in 2 unknowns. - * However, the intersection problem is really only a 2-dimensional - * problem, being located in the surface of a plane. - * Therefore, the "least important" of these equations can - * be initially ignored, leaving a system of 2 equations in - * 2 unknowns. - * - * Find the component of N with the largest magnitude. - * This component will have the least effect on the parameters - * in the system, being most nearly perpendicular to the plane. - * Denote the two remaining components by the - * subscripts q and r, rather than x,y,z. - * Subscript s is the smallest component, used for checking later. - */ - abs_n[X] = (n[X] >= 0) ? n[X] : (-n[X]); - abs_n[Y] = (n[Y] >= 0) ? n[Y] : (-n[Y]); - abs_n[Z] = (n[Z] >= 0) ? n[Z] : (-n[Z]); - if( abs_n[X] >= abs_n[Y] ) { - if( abs_n[X] >= abs_n[Z] ) { - /* X is largest in magnitude */ - q = Y; - r = Z; - s = X; - } else { - /* Z is largest in magnitude */ - q = X; - r = Y; - s = Z; - } - } else { - if( abs_n[Y] >= abs_n[Z] ) { - /* Y is largest in magnitude */ - q = X; - r = Z; - s = Y; - } else { - /* Z is largest in magnitude */ - q = X; - r = Y; - s = Z; - } - } - -#if 0 - /* XXX Use rt_isect_line2_line2() here */ - /* move the 2d vectors around */ - rt_isect_line2_line2( &dist, p, d, a, c, tol ); -#endif - - /* - * From the two components q and r, form a system - * of 2 equations in 2 unknowns: - * - * Pq + t * Dq = Aq + u * Cq - * Pr + t * Dr = Ar + u * Cr - * or - * t * Dq - u * Cq = Aq - Pq - * t * Dr - u * Cr = Ar - Pr - * - * Let H = A - P, resulting in: - * - * t * Dq - u * Cq = Hq - * t * Dr - u * Cr = Hr - * - * or - * - * [ Dq -Cq ] [ t ] [ Hq ] - * [ ] * [ ] = [ ] - * [ Dr -Cr ] [ u ] [ Hr ] - * - * This system can be solved by direct substitution, or by - * finding the determinants by Cramers rule: - * - * [ Dq -Cq ] - * det(M) = det [ ] = -Dq * Cr + Cq * Dr - * [ Dr -Cr ] - * - * If det(M) is zero, then the lines are parallel (perhaps colinear). - * Otherwise, exactly one solution exists. - */ - VSUB2( h, a, p ); - det = c[q] * d[r] - d[q] * c[r]; - det1 = (c[q] * h[r] - h[q] * c[r]); /* see below */ - /* XXX This should be no smaller than 1e-16. See rt_isect_line2_line2 for details */ - if( NEAR_ZERO( det, DETERMINANT_TOL ) ) { - /* Lines are parallel */ - if( !NEAR_ZERO( det1, DETERMINANT_TOL ) ) { - /* Lines are NOT co-linear, just parallel */ - return -2; /* parallel, no intersection */ - } - - /* Lines are co-linear */ - /* Compute t for u=0 as a convenience to caller */ - *u = 0; - /* Use largest direction component */ - if( fabs(d[q]) >= fabs(d[r]) ) { - *t = h[q]/d[q]; - } else { - *t = h[r]/d[r]; - } - return(0); /* Lines co-linear */ - } - - /* - * det(M) is non-zero, so there is exactly one solution. - * Using Cramer's rule, det1(M) replaces the first column - * of M with the constant column vector, in this case H. - * Similarly, det2(M) replaces the second column. - * Computation of the determinant is done as before. - * - * Now, - * - * [ Hq -Cq ] - * det [ ] - * det1(M) [ Hr -Cr ] -Hq * Cr + Cq * Hr - * t = ------- = --------------- = ------------------ - * det(M) det(M) -Dq * Cr + Cq * Dr - * - * and - * - * [ Dq Hq ] - * det [ ] - * det2(M) [ Dr Hr ] Dq * Hr - Hq * Dr - * u = ------- = --------------- = ------------------ - * det(M) det(M) -Dq * Cr + Cq * Dr - */ - det = 1/det; - *t = det * det1; - *u = det * (d[q] * h[r] - h[q] * d[r]); - - /* - * Check that these values of t and u satisfy the 3rd equation - * as well! - * XXX It isn't clear that "det" is exactly a model-space distance. - */ - det = *t * d[s] - *u * c[s] - h[s]; - if( !NEAR_ZERO( det, tol->dist ) ) { - /* XXX This tolerance needs to be much less loose than - * XXX SQRT_SMALL_FASTF. What about DETERMINANT_TOL? - */ - /* Inconsistent solution, lines miss each other */ - return(-1); - } - - /* To prevent errors, check the answer. - * Not returning bogus results to our caller is worth the extra time. - */ - { - point_t hit1, hit2; - - VJOIN1( hit1, p, *t, d ); - VJOIN1( hit2, a, *u, c ); - if( !rt_pt3_pt3_equal( hit1, hit2, tol ) ) { -/* rt_log("rt_isect_line3_line3(): BOGUS RESULT, hit1=(%g,%g,%g), hit2=(%g,%g,%g)\n", - hit1[X], hit1[Y], hit1[Z], hit2[X], hit2[Y], hit2[Z]); */ - return -1; - } - } - - return(1); /* Intersection found */ -} - -/* - * R T _ I S E C T _ L I N E _ L S E G - * - * Intersect a line in parametric form: - * - * X = P + t * D - * - * with a line segment defined by two distinct points A and B. - * - * Explicit Return - - * -4 A and B are not distinct points - * -3 Intersection exists, < A (t is returned) - * -2 Intersection exists, > B (t is returned) - * -1 Lines do not intersect - * 0 Lines are co-linear (t for A is returned) - * 1 Intersection at vertex A - * 2 Intersection at vertex B - * 3 Intersection between A and B - * - * Implicit Returns - - * t When explicit return >= 0, t is the parameter that describes - * the intersection of the line and the line segment. - * The actual intersection coordinates can be found by - * solving P + t * D. However, note that for return codes - * 1 and 2 (intersection exactly at a vertex), it is - * strongly recommended that the original values passed in - * A or B are used instead of solving P + t * D, to prevent - * numeric error from creeping into the position of - * the endpoints. - */ -/* XXX should probably be called rt_isect_line3_lseg3() */ -/* XXX should probably be changed to return dist[2] */ -int -rt_isect_line_lseg( t, p, d, a, b, tol ) -fastf_t *t; -CONST point_t p; -CONST vect_t d; -CONST point_t a; -CONST point_t b; -CONST struct rt_tol *tol; -{ - LOCAL vect_t c; /* Direction vector from A to B */ - auto fastf_t u; /* As in, A + u * C = X */ - register fastf_t f; - register int ret; - fastf_t fuzz; - - RT_CK_TOL(tol); - - VSUB2( c, b, a ); - /* - * To keep the values of u between 0 and 1, - * C should NOT be scaled to have unit length. - * However, it is a good idea to make sure that - * C is a non-zero vector, (ie, that A and B are distinct). - */ - if( (fuzz = MAGSQ(c)) < tol->dist_sq ) { - return(-4); /* points A and B are not distinct */ - } - - /* - * Detecting colinearity is difficult, and very very important. - * As a first step, check to see if both points A and B lie - * within tolerance of the line. If so, then the line segment AC - * is ON the line. - */ - if( rt_distsq_line3_pt3( p, d, a ) <= tol->dist_sq && - rt_distsq_line3_pt3( p, d, b ) <= tol->dist_sq ) { - if( rt_g.debug & DEBUG_MATH ) { - rt_log("rt_isect_line3_lseg3() pts A and B within tol of line\n"); - } - /* Find the parametric distance along the ray */ - *t = rt_dist_pt3_along_line3( p, d, a ); - /*** dist[1] = rt_dist_pt3_along_line3( p, d, b ); ***/ - return 0; /* Colinear */ - } - - if( (ret = rt_isect_line3_line3( t, &u, p, d, a, c, tol )) < 0 ) { - /* No intersection found */ - return( -1 ); - } - if( ret == 0 ) { - /* co-linear (t was computed for point A, u=0) */ - return( 0 ); - } - - /* - * The two lines intersect at a point. - * If the u parameter is outside the range (0..1), - * reject the intersection, because it falls outside - * the line segment A--B. - * - * Convert the tol->dist into allowable deviation in terms of - * (0..1) range of the parameters. - */ - fuzz = tol->dist / sqrt(fuzz); - if( u < -fuzz ) - return(-3); /* Intersection < A */ - if( (f=(u-1)) > fuzz ) - return(-2); /* Intersection > B */ - - /* Check for fuzzy intersection with one of the verticies */ - if( u < fuzz ) - return( 1 ); /* Intersection at A */ - if( f >= -fuzz ) - return( 2 ); /* Intersection at B */ - - return(3); /* Intersection between A and B */ -} - -/* - * R T _ D I S T _ L I N E 3_ P T 3 - * - * Given a parametric line defined by PT + t * DIR and a point A, - * return the closest distance between the line and the point. - * - * 'dir' need not have unit length. - * - * Find parameter for PCA along line with unitized DIR: - * d = VDOT(f, dir) / MAGNITUDE(dir); - * Find distance g from PCA to A using Pythagoras: - * g = sqrt( MAGSQ(f) - d**2 ) - * - * Return - - * Distance - */ -double -rt_dist_line3_pt3( pt, dir, a ) -CONST point_t pt; -CONST vect_t dir; -CONST point_t a; -{ - LOCAL vect_t f; - register fastf_t FdotD; - - if( (FdotD = MAGNITUDE(dir)) <= SMALL_FASTF ) { - FdotD = 0.0; - goto out; - } - VSUB2( f, a, pt ); - FdotD = VDOT( f, dir ) / FdotD; - FdotD = MAGSQ( f ) - FdotD * FdotD; - if( FdotD <= SMALL_FASTF ) { - FdotD = 0.0; - goto out; - } - FdotD = sqrt(FdotD); -out: - if( rt_g.debug & DEBUG_MATH ) { - rt_log("rt_dist_line3_pt3() ret=%g\n", FdotD); - } - return FdotD; -} - -/* - * R T _ D I S T S Q _ L I N E 3 _ P T 3 - * - * Given a parametric line defined by PT + t * DIR and a point A, - * return the square of the closest distance between the line and the point. - * - * 'dir' need not have unit length. - * - * Return - - * Distance squared - */ -double -rt_distsq_line3_pt3( pt, dir, a ) -CONST point_t pt; -CONST vect_t dir; -CONST point_t a; -{ - LOCAL vect_t f; - register fastf_t FdotD; - - VSUB2( f, pt, a ); - if( (FdotD = MAGNITUDE(dir)) <= SMALL_FASTF ) { - FdotD = 0.0; - goto out; - } - FdotD = VDOT( f, dir ) / FdotD; - if( (FdotD = VDOT( f, f ) - FdotD * FdotD ) <= SMALL_FASTF ) { - FdotD = 0.0; - } -out: - if( rt_g.debug & DEBUG_MATH ) { - rt_log("rt_distsq_line3_pt3() ret=%g\n", FdotD); - } - return FdotD; -} - -/* - * R T _ D I S T _ L I N E _ O R I G I N - * - * Given a parametric line defined by PT + t * DIR, - * return the closest distance between the line and the origin. - * - * 'dir' need not have unit length. - * - * Return - - * Distance - */ -double -rt_dist_line_origin( pt, dir ) -CONST point_t pt; -CONST vect_t dir; -{ - register fastf_t PTdotD; - - if( (PTdotD = MAGNITUDE(dir)) <= SMALL_FASTF ) - return 0.0; - PTdotD = VDOT( pt, dir ) / PTdotD; - if( (PTdotD = VDOT( pt, pt ) - PTdotD * PTdotD ) <= SMALL_FASTF ) - return(0.0); - return( sqrt(PTdotD) ); -} -/* - * R T _ D I S T _ L I N E 2 _ P O I N T 2 - * - * Given a parametric line defined by PT + t * DIR and a point A, - * return the closest distance between the line and the point. - * - * 'dir' need not have unit length. - * - * Return - - * Distance - */ -double -rt_dist_line2_point2( pt, dir, a ) -CONST point_t pt; -CONST vect_t dir; -CONST point_t a; -{ - LOCAL vect_t f; - register fastf_t FdotD; - - VSUB2_2D( f, pt, a ); - if( (FdotD = sqrt(MAGSQ_2D(dir))) <= SMALL_FASTF ) - return 0.0; - FdotD = VDOT_2D( f, dir ) / FdotD; - if( (FdotD = VDOT_2D( f, f ) - FdotD * FdotD ) <= SMALL_FASTF ) - return(0.0); - return( sqrt(FdotD) ); -} - -/* - * R T _ D I S T S Q _ L I N E 2 _ P O I N T 2 - * - * Given a parametric line defined by PT + t * DIR and a point A, - * return the closest distance between the line and the point, squared. - * - * 'dir' need not have unit length. - * - * Return - - * Distance squared - */ -double -rt_distsq_line2_point2( pt, dir, a ) -CONST point_t pt; -CONST vect_t dir; -CONST point_t a; -{ - LOCAL vect_t f; - register fastf_t FdotD; - - VSUB2_2D( f, pt, a ); - if( (FdotD = sqrt(MAGSQ_2D(dir))) <= SMALL_FASTF ) - return 0.0; - FdotD = VDOT_2D( f, dir ) / FdotD; - if( (FdotD = VDOT_2D( f, f ) - FdotD * FdotD ) <= SMALL_FASTF ) - return(0.0); - return( FdotD ); -} - -/* - * R T _ A R E A _ O F _ T R I A N G L E - * - * Returns the area of a triangle. - * Algorithm by Jon Leech 3/24/89. - */ -double -rt_area_of_triangle( a, b, c ) -register CONST point_t a, b, c; -{ - register double t; - register double area; - - t = a[Y] * (b[Z] - c[Z]) - - b[Y] * (a[Z] - c[Z]) + - c[Y] * (a[Z] - b[Z]); - area = t*t; - t = a[Z] * (b[X] - c[X]) - - b[Z] * (a[X] - c[X]) + - c[Z] * (a[X] - b[X]); - area += t*t; - t = a[X] * (b[Y] - c[Y]) - - b[X] * (a[Y] - c[Y]) + - c[X] * (a[Y] - b[Y]); - area += t*t; - - return( 0.5 * sqrt(area) ); -} - - -/* - * R T _ I S E C T _ P T _ L S E G - * - * Intersect a point P with the line segment defined by two distinct - * points A and B. - * - * Explicit Return - * -2 P on line AB but outside range of AB, - * dist = distance from A to P on line. - * -1 P not on line of AB within tolerance - * 1 P is at A - * 2 P is at B - * 3 P is on AB, dist = distance from A to P on line. - * - * B * - * | - * P'*-tol-*P - * | / _ - * dist / /| - * | / / - * | / / AtoP - * |/ / - * A * / - * - * tol = distance limit from line to pt P; - * dist = distance from A to P' - */ -int rt_isect_pt_lseg(dist, a, b, p, tol) -fastf_t *dist; /* distance along line from A to P */ -CONST point_t a, b, p; /* points for line and intersect */ -CONST struct rt_tol *tol; -{ - vect_t AtoP, - BtoP, - AtoB, - ABunit; /* unit vector from A to B */ - fastf_t APprABunit; /* Mag of projection of AtoP onto ABunit */ - fastf_t distsq; - - RT_CK_TOL(tol); - - VSUB2(AtoP, p, a); - if (MAGSQ(AtoP) < tol->dist_sq) - return(1); /* P at A */ - - VSUB2(BtoP, p, b); - if (MAGSQ(BtoP) < tol->dist_sq) - return(2); /* P at B */ - - VSUB2(AtoB, b, a); - VMOVE(ABunit, AtoB); - distsq = MAGSQ(ABunit); - if( distsq < tol->dist_sq ) - return -1; /* A equals B, and P isn't there */ - distsq = 1/sqrt(distsq); - VSCALE( ABunit, ABunit, distsq ); - - /* Similar to rt_dist_line_pt, except we - * never actually have to do the sqrt that the other routine does. - */ - - /* find dist as a function of ABunit, actually the projection - * of AtoP onto ABunit - */ - APprABunit = VDOT(AtoP, ABunit); - - /* because of pythgorean theorem ... */ - distsq = MAGSQ(AtoP) - APprABunit * APprABunit; - if (distsq > tol->dist_sq) - return(-1); /* dist pt to line too large */ - - /* Distance from the point to the line is within tolerance. */ - *dist = VDOT(AtoP, AtoB) / MAGSQ(AtoB); - - if (*dist > 1.0 || *dist < 0.0) /* P outside AtoB */ - return(-2); - - return(3); /* P on AtoB */ -} - -/* - * R T _ I S E C T _ P T 2 _ L S E G 2 - * - * Intersect a point P with the line segment defined by two distinct - * points A and B. - * - * Explicit Return - * -2 P on line AB but outside range of AB, - * dist = distance from A to P on line. - * -1 P not on line of AB within tolerance - * 1 P is at A - * 2 P is at B - * 3 P is on AB, dist = distance from A to P on line. - * - * B * - * | - * P'*-tol-*P - * | / _ - * dist / /| - * | / / - * | / / AtoP - * |/ / - * A * / - * - * tol = distance limit from line to pt P; - * dist = distance from A to P' - */ -int -rt_isect_pt2_lseg2(dist, a, b, p, tol) -fastf_t *dist; /* distance along line from A to P */ -CONST point_t a, b, p; /* points for line and intersect */ -CONST struct rt_tol *tol; -{ - vect_t AtoP, - BtoP, - AtoB, - ABunit; /* unit vector from A to B */ - fastf_t APprABunit; /* Mag of projection of AtoP onto ABunit */ - fastf_t distsq; - - RT_CK_TOL(tol); - - VSUB2_2D(AtoP, p, a); - if (MAGSQ_2D(AtoP) < tol->dist_sq) - return(1); /* P at A */ - - VSUB2_2D(BtoP, p, b); - if (MAGSQ_2D(BtoP) < tol->dist_sq) - return(2); /* P at B */ - - VSUB2_2D(AtoB, b, a); - VMOVE_2D(ABunit, AtoB); - distsq = MAGSQ_2D(ABunit); - if( distsq < tol->dist_sq ) { - if( rt_g.debug & DEBUG_MATH ) { - rt_log("distsq A=%g\n", distsq); - } - return -1; /* A equals B, and P isn't there */ - } - distsq = 1/sqrt(distsq); - VSCALE_2D( ABunit, ABunit, distsq ); - - /* Similar to rt_dist_line_pt, except we - * never actually have to do the sqrt that the other routine does. - */ - - /* find dist as a function of ABunit, actually the projection - * of AtoP onto ABunit - */ - APprABunit = VDOT_2D(AtoP, ABunit); - - /* because of pythgorean theorem ... */ - distsq = MAGSQ_2D(AtoP) - APprABunit * APprABunit; - if (distsq > tol->dist_sq) { - if( rt_g.debug & DEBUG_MATH ) { - V2PRINT("ABunit", ABunit); - rt_log("distsq B=%g\n", distsq); - } - return(-1); /* dist pt to line too large */ - } - - /* Distance from the point to the line is within tolerance. */ - *dist = VDOT_2D(AtoP, AtoB) / MAGSQ_2D(AtoB); - - if (*dist > 1.0 || *dist < 0.0) /* P outside AtoB */ - return(-2); - - return(3); /* P on AtoB */ -} - -/* - * R T _ D I S T _ P T 3 _ L S E G 3 - * - * Find the distance from a point P to a line segment described - * by the two endpoints A and B, and the point of closest approach (PCA). - * - * P - * * - * /. - * / . - * / . - * / . (dist) - * / . - * / . - * *------*--------* - * A PCA B - * - * There are six distinct cases, with these return codes - - * 0 P is within tolerance of lseg AB. *dist isn't 0: (SPECIAL!!!) - * *dist = parametric dist = |PCA-A| / |B-A|. pca=computed. - * 1 P is within tolerance of point A. *dist = 0, pca=A. - * 2 P is within tolerance of point B. *dist = 0, pca=B. - * 3 P is to the "left" of point A. *dist=|P-A|, pca=A. - * 4 P is to the "right" of point B. *dist=|P-B|, pca=B. - * 5 P is "above/below" lseg AB. *dist=|PCA-P|, pca=computed. - * - * This routine was formerly called rt_dist_pt_lseg(). - * - * XXX For efficiency, a version of this routine that provides the - * XXX distance squared would be faster. - */ -int -rt_dist_pt3_lseg3( dist, pca, a, b, p, tol ) -fastf_t *dist; -point_t pca; -CONST point_t a, b, p; -CONST struct rt_tol *tol; -{ - vect_t PtoA; /* P-A */ - vect_t PtoB; /* P-B */ - vect_t AtoB; /* B-A */ - fastf_t P_A_sq; /* |P-A|**2 */ - fastf_t P_B_sq; /* |P-B|**2 */ - fastf_t B_A; /* |B-A| */ - fastf_t t; /* distance along ray of projection of P */ - - RT_CK_TOL(tol); - - if( rt_g.debug & DEBUG_MATH ) { - rt_log("rt_dist_pt3_lseg3() a=(%g,%g,%g) b=(%g,%g,%g)\n\tp=(%g,%g,%g), tol->dist=%g sq=%g\n", - V3ARGS(a), - V3ARGS(b), - V3ARGS(p), - tol->dist, tol->dist_sq ); - } - - /* Check proximity to endpoint A */ - VSUB2(PtoA, p, a); - if( (P_A_sq = MAGSQ(PtoA)) < tol->dist_sq ) { - /* P is within the tol->dist radius circle around A */ - VMOVE( pca, a ); - if( rt_g.debug & DEBUG_MATH ) rt_log(" at A\n"); - *dist = 0.0; - return 1; - } - - /* Check proximity to endpoint B */ - VSUB2(PtoB, p, b); - if( (P_B_sq = MAGSQ(PtoB)) < tol->dist_sq ) { - /* P is within the tol->dist radius circle around B */ - VMOVE( pca, b ); - if( rt_g.debug & DEBUG_MATH ) rt_log(" at B\n"); - *dist = 0.0; - return 2; - } - - VSUB2(AtoB, b, a); - B_A = sqrt( MAGSQ(AtoB) ); - - /* compute distance (in actual units) along line to PROJECTION of - * point p onto the line: point pca - */ - t = VDOT(PtoA, AtoB) / B_A; - if( rt_g.debug & DEBUG_MATH ) { - rt_log("rt_dist_pt3_lseg3() B_A=%g, t=%g\n", - B_A, t ); - } - - if( t <= 0 ) { - /* P is "left" of A */ - if( rt_g.debug & DEBUG_MATH ) rt_log(" left of A\n"); - VMOVE( pca, a ); - *dist = sqrt(P_A_sq); - return 3; - } - if( t < B_A ) { - /* PCA falls between A and B */ - register fastf_t dsq; - fastf_t param_dist; /* parametric dist */ - - /* Find PCA */ - param_dist = t / B_A; /* Range 0..1 */ - VJOIN1(pca, a, param_dist, AtoB); - - /* Find distance from PCA to line segment (Pythagorus) */ - if( (dsq = P_A_sq - t * t ) <= tol->dist_sq ) { - if( rt_g.debug & DEBUG_MATH ) rt_log(" ON lseg\n"); - /* Distance from PCA to lseg is zero, give param instead */ - *dist = param_dist; /* special! */ - return 0; - } - if( rt_g.debug & DEBUG_MATH ) rt_log(" closest to lseg\n"); - *dist = sqrt(dsq); - return 5; - } - /* P is "right" of B */ - if( rt_g.debug & DEBUG_MATH ) rt_log(" right of B\n"); - VMOVE(pca, b); - *dist = sqrt(P_B_sq); - return 4; -} - -/* - * R T _ D I S T _ P T 2 _ L S E G 2 - * - * Find the distance from a point P to a line segment described - * by the two endpoints A and B, and the point of closest approach (PCA). - * - * P - * * - * /. - * / . - * / . - * / . (dist) - * / . - * / . - * *------*--------* - * A PCA B - * - * There are six distinct cases, with these return codes - - * 0 P is within tolerance of lseg AB. *dist isn't 0: (SPECIAL!!!) - * *dist = parametric dist = |PCA-A| / |B-A|. pca=computed. - * 1 P is within tolerance of point A. *dist = 0, pca=A. - * 2 P is within tolerance of point B. *dist = 0, pca=B. - * 3 P is to the "left" of point A. *dist=|P-A|**2, pca=A. - * 4 P is to the "right" of point B. *dist=|P-B|**2, pca=B. - * 5 P is "above/below" lseg AB. *dist=|PCA-P|**2, pca=computed. - * - * - * Patterned after rt_dist_pt3_lseg3(). - */ -int -rt_dist_pt2_lseg2( dist_sq, pca, a, b, p, tol ) -fastf_t *dist_sq; -fastf_t pca[2]; -CONST point_t a, b, p; -CONST struct rt_tol *tol; -{ - vect_t PtoA; /* P-A */ - vect_t PtoB; /* P-B */ - vect_t AtoB; /* B-A */ - fastf_t P_A_sq; /* |P-A|**2 */ - fastf_t P_B_sq; /* |P-B|**2 */ - fastf_t B_A; /* |B-A| */ - fastf_t t; /* distance along ray of projection of P */ - - RT_CK_TOL(tol); - - if( rt_g.debug & DEBUG_MATH ) { - rt_log("rt_dist_pt3_lseg3() a=(%g,%g,%g) b=(%g,%g,%g)\n\tp=(%g,%g,%g), tol->dist=%g sq=%g\n", - V3ARGS(a), - V3ARGS(b), - V3ARGS(p), - tol->dist, tol->dist_sq ); - } - - - /* Check proximity to endpoint A */ - VSUB2_2D(PtoA, p, a); - if( (P_A_sq = MAGSQ_2D(PtoA)) < tol->dist_sq ) { - /* P is within the tol->dist radius circle around A */ - V2MOVE( pca, a ); - if( rt_g.debug & DEBUG_MATH ) rt_log(" at A\n"); - *dist_sq = 0.0; - return 1; - } - - /* Check proximity to endpoint B */ - VSUB2_2D(PtoB, p, b); - if( (P_B_sq = MAGSQ_2D(PtoB)) < tol->dist_sq ) { - /* P is within the tol->dist radius circle around B */ - V2MOVE( pca, b ); - if( rt_g.debug & DEBUG_MATH ) rt_log(" at B\n"); - *dist_sq = 0.0; - return 2; - } - - VSUB2_2D(AtoB, b, a); - B_A = sqrt( MAGSQ_2D(AtoB) ); - - /* compute distance (in actual units) along line to PROJECTION of - * point p onto the line: point pca - */ - t = VDOT_2D(PtoA, AtoB) / B_A; - if( rt_g.debug & DEBUG_MATH ) { - rt_log("rt_dist_pt3_lseg3() B_A=%g, t=%g\n", - B_A, t ); - } - - if( t <= 0 ) { - /* P is "left" of A */ - if( rt_g.debug & DEBUG_MATH ) rt_log(" left of A\n"); - V2MOVE( pca, a ); - *dist_sq = P_A_sq; - return 3; - } - if( t < B_A ) { - /* PCA falls between A and B */ - register fastf_t dsq; - fastf_t param_dist; /* parametric dist */ - - /* Find PCA */ - param_dist = t / B_A; /* Range 0..1 */ - V2JOIN1(pca, a, param_dist, AtoB); - - /* Find distance from PCA to line segment (Pythagorus) */ - if( (dsq = P_A_sq - t * t ) <= tol->dist_sq ) { - if( rt_g.debug & DEBUG_MATH ) rt_log(" ON lseg\n"); - /* Distance from PCA to lseg is zero, give param instead */ - *dist_sq = param_dist; /* special! Not squared. */ - return 0; - } - if( rt_g.debug & DEBUG_MATH ) rt_log(" closest to lseg\n"); - *dist_sq = dsq; - return 5; - } - /* P is "right" of B */ - if( rt_g.debug & DEBUG_MATH ) rt_log(" right of B\n"); - V2MOVE(pca, b); - *dist_sq = P_B_sq; - return 4; -} - -/* - * R T _ R O T A T E _ B B O X - * - * Transform a bounding box (RPP) by the given 4x4 matrix. - * There are 8 corners to the bounding RPP. - * Each one needs to be transformed and min/max'ed. - * This is not minimal, but does fully contain any internal object, - * using an axis-aligned RPP. - */ -void -rt_rotate_bbox( omin, omax, mat, imin, imax ) -point_t omin; -point_t omax; -CONST mat_t mat; -CONST point_t imin; -CONST point_t imax; -{ - point_t local; /* vertex point in local coordinates */ - point_t model; /* vertex point in model coordinates */ - -#define ROT_VERT( a, b, c ) \ - VSET( local, a[X], b[Y], c[Z] ); \ - MAT4X3PNT( model, mat, local ); \ - VMINMAX( omin, omax, model ) \ - - ROT_VERT( imin, imin, imin ); - ROT_VERT( imin, imin, imax ); - ROT_VERT( imin, imax, imin ); - ROT_VERT( imin, imax, imax ); - ROT_VERT( imax, imin, imin ); - ROT_VERT( imax, imin, imax ); - ROT_VERT( imax, imax, imin ); - ROT_VERT( imax, imax, imax ); -#undef ROT_VERT -} - -/* - * R T _ R O T A T E _ P L A N E - * - * Transform a plane equation by the given 4x4 matrix. - */ -void -rt_rotate_plane( oplane, mat, iplane ) -plane_t oplane; -CONST mat_t mat; -CONST plane_t iplane; -{ - point_t orig_pt; - point_t new_pt; - - /* First, pick a point that lies on the original halfspace */ - VSCALE( orig_pt, iplane, iplane[3] ); - - /* Transform the surface normal */ - MAT4X3VEC( oplane, mat, iplane ); - - /* Transform the point from original to new halfspace */ - MAT4X3PNT( new_pt, mat, orig_pt ); - - /* - * The transformed normal is all that is required. - * The new distance is found from the transformed point on the plane. - */ - oplane[3] = VDOT( new_pt, oplane ); -} - -/* - * R T _ C O P L A N A R - * - * Test if two planes are identical. If so, their dot products will be - * either +1 or -1, with the distance from the origin equal in magnitude. - * - * Returns - - * -1 not coplanar, parallel but distinct - * 0 not coplanar, not parallel. Planes intersect. - * +1 coplanar, same normal direction - * +2 coplanar, opposite normal direction - */ -int -rt_coplanar( a, b, tol ) -CONST plane_t a; -CONST plane_t b; -CONST struct rt_tol *tol; -{ - register fastf_t f; - register fastf_t dot; - - RT_CK_TOL(tol); - - /* Check to see if the planes are parallel */ - dot = VDOT( a, b ); - if( dot >= 0 ) { - /* Normals head in generally the same directions */ - if( dot < tol->para ) - return(0); /* Planes intersect */ - - /* Planes have "exactly" the same normal vector */ - f = a[3] - b[3]; - if( NEAR_ZERO( f, tol->dist ) ) { - return(1); /* Coplanar, same direction */ - } - return(-1); /* Parallel but distinct */ - } - /* Normals head in generally opposite directions */ - if( -dot < tol->para ) - return(0); /* Planes intersect */ - - /* Planes have "exactly" opposite normal vectors */ - f = a[3] + b[3]; - if( NEAR_ZERO( f, tol->dist ) ) { - return(2); /* Coplanar, opposite directions */ - } - return(-1); /* Parallel but distinct */ -} - -/* - * R T _ A N G L E _ M E A S U R E - * - * Using two perpendicular vectors (x_dir and y_dir) which lie - * in the same plane as 'vec', return the angle (in radians) of 'vec' - * from x_dir, going CCW around the perpendicular x_dir CROSS y_dir. - * - * Trig note - - * - * theta = atan2(x,y) returns an angle in the range -pi to +pi. - * Here, we need an angle in the range of 0 to 2pi. - * This could be implemented by adding 2pi to theta when theta is negative, - * but this could have nasty numeric ambiguity right in the vicinity - * of theta = +pi, which is a very critical angle for the applications using - * this routine. - * So, an alternative formulation is to compute gamma = atan2(-x,-y), - * and then theta = gamma + pi. Now, any error will occur in the - * vicinity of theta = 0, which can be handled much more readily. - * - * If theta is negative, or greater than two pi, - * wrap it around. - * These conditions only occur if there are problems in atan2(). - * - * Returns - - * vec == x_dir returns 0, - * vec == y_dir returns pi/2, - * vec == -x_dir returns pi, - * vec == -y_dir returns 3*pi/2. - * - * In all cases, the returned value is between 0 and rt_twopi. - */ -double -rt_angle_measure( vec, x_dir, y_dir ) -vect_t vec; -CONST vect_t x_dir; -CONST vect_t y_dir; -{ - fastf_t xproj, yproj; - fastf_t gamma; - fastf_t ang; - - xproj = -VDOT( vec, x_dir ); - yproj = -VDOT( vec, y_dir ); - gamma = atan2( yproj, xproj ); /* -pi..+pi */ - ang = rt_pi + gamma; /* 0..+2pi */ - if( ang < 0 ) { - do { - ang += rt_twopi; - } while( ang < 0 ); - } else if( ang > rt_twopi ) { - do { - ang -= rt_twopi; - } while( ang > rt_twopi ); - } - if( ang < 0 || ang > rt_twopi ) rt_bomb("rt_angle_measure() angle out of range\n"); - return ang; -} - -/* - * R T _ D I S T _ P T 3 _ A L O N G _ L I N E 3 - * - * Return the parametric distance t of a point X along a line defined - * as a ray, i.e. solve X = P + t * D. - * If the point X does not lie on the line, then t is the distance of - * the perpendicular projection of point X onto the line. - */ -double -rt_dist_pt3_along_line3( p, d, x ) -CONST point_t p; -CONST vect_t d; -CONST point_t x; -{ - vect_t x_p; - - VSUB2( x_p, x, p ); - return VDOT( x_p, d ); -} - - -/* - * R T _ D I S T _ P T 2 _ A L O N G _ L I N E 2 - * - * Return the parametric distance t of a point X along a line defined - * as a ray, i.e. solve X = P + t * D. - * If the point X does not lie on the line, then t is the distance of - * the perpendicular projection of point X onto the line. - */ -double -rt_dist_pt2_along_line2( p, d, x ) -CONST point_t p; -CONST vect_t d; -CONST point_t x; -{ - vect_t x_p; - double ret; - - VSUB2_2D( x_p, x, p ); - ret = VDOT_2D( x_p, d ); - if( rt_g.debug & DEBUG_MATH ) { - rt_log("rt_dist_pt2_along_line2() p=(%g, %g), d=(%g, %g), x=(%g, %g) ret=%g\n", - V2ARGS(p), - V2ARGS(d), - V2ARGS(x), - ret ); - } - return ret; -} - -/* - * Returns - - * 1 if left <= mid <= right - * 0 if mid is not in the range. - */ -int -rt_between( left, mid, right, tol ) -double left; -double mid; -double right; -CONST struct rt_tol *tol; -{ - RT_CK_TOL(tol); - - if( left < right ) { - if( NEAR_ZERO(left-right, tol->dist*0.1) ) { - left -= tol->dist*0.1; - right += tol->dist*0.1; - } - if( mid < left || mid > right ) goto fail; - return 1; - } - /* The 'right' value is lowest */ - if( NEAR_ZERO(left-right, tol->dist*0.1) ) { - right -= tol->dist*0.1; - left += tol->dist*0.1; - } - if( mid < right || mid > left ) goto fail; - return 1; -fail: - if( rt_g.debug & DEBUG_MATH ) { - rt_log("rt_between( %.17e, %.17e, %.17e ) ret=0 FAIL\n", - left, mid, right); - } - return 0; -} diff --git a/librt/polylib.c b/librt/polylib.c deleted file mode 100644 index 063f06231b9..00000000000 --- a/librt/polylib.c +++ /dev/null @@ -1,551 +0,0 @@ -/* - * P O L Y L I B . C - * - * Library for dealing with polynomials. - * - * Functions - - * rt_poly_mul Multiply two polynomials - * rt_poly_scale Scale a polynomial - * rt_poly_add Add two polynomials - * rt_poly_sub Subtract two polynomials - * rt_poly_synthetic_division Divide 1 poly into another using Synthetic Division - * quadratic Solve quadratic formula - * cubic Solve cubic forumla - * rt_pr_poly Print a polynomial - * - * Author - - * Jeff Hanes - * - * Source - - * SECAD/VLD Computing Consortium, Bldg 394 - * The U. S. Army Ballistic Research Laboratory - * Aberdeen Proving Ground, Maryland 21005 - * - * Copyright Notice - - * This software is Copyright (C) 1985 by the United States Army. - * All rights reserved. - */ -#ifndef lint -static char RCSpolylib[] = "@(#)$Header$ (BRL)"; -#endif - -#include "conf.h" - -#include -#include -#include -#include "machine.h" -#include "externs.h" -#include "vmath.h" -#include "raytrace.h" -#include "rtstring.h" -#include "./polyno.h" -#include "./complex.h" - -extern void rt_pr_poly(); - -#ifndef M_PI -#define M_PI 3.14159265358979323846 -#endif -#define PI_DIV_3 (M_PI/3.0) -static poly Zpoly = { 0, 0.0 }; - -/* - * rt_poly_mul -- multiply two polynomials - */ -poly * -rt_poly_mul(m1,m2,product) -register poly *m1, *m2, *product; -{ - if( m1->dgr == 1 && m2->dgr == 1 ) { - product->dgr = 2; - product->cf[0] = m1->cf[0] * m2->cf[0]; - product->cf[1] = m1->cf[0] * m2->cf[1] + - m1->cf[1] * m2->cf[0]; - product->cf[2] = m1->cf[1] * m2->cf[1]; - return(product); - } - if( m1->dgr == 2 && m2->dgr == 2 ) { - product->dgr = 4; - product->cf[0] = m1->cf[0] * m2->cf[0]; - product->cf[1] = m1->cf[0] * m2->cf[1] + - m1->cf[1] * m2->cf[0]; - product->cf[2] = m1->cf[0] * m2->cf[2] + - m1->cf[1] * m2->cf[1] + - m1->cf[2] * m2->cf[0]; - product->cf[3] = m1->cf[1] * m2->cf[2] + - m1->cf[2] * m2->cf[1]; - product->cf[4] = m1->cf[2] * m2->cf[2]; - return(product); - } - - /* Not one of the common (or easy) cases. */ - { - register int ct1, ct2; - - *product = Zpoly; - - /* If the degree of the product will be larger than the - * maximum size allowed in "polyno.h", then return a null - * pointer to indicate failure. - */ - if ( (product->dgr = m1->dgr + m2->dgr) > MAXP ) - return PM_NULL; - - for ( ct1=0; ct1 <= m1->dgr; ++ct1 ){ - for ( ct2=0; ct2 <= m2->dgr; ++ct2 ){ - product->cf[ct1+ct2] += - m1->cf[ct1] * m2->cf[ct2]; - } - } - } - return product; -} - - -/* - * rt_poly_scale -- scale a polynomial - */ -poly * -rt_poly_scale(eqn,factor) -register poly *eqn; -double factor; -{ - register int cnt; - - for ( cnt=0; cnt <= eqn->dgr; ++cnt ){ - eqn->cf[cnt] *= factor; - } - return eqn; -} - - -/* - * rt_poly_add -- add two polynomials - */ -poly * -rt_poly_add(poly1,poly2,sum) -register poly *poly1, *poly2, *sum; -{ - LOCAL poly tmp; - register int i, offset; - - offset = Abs(poly1->dgr - poly2->dgr); - - tmp = Zpoly; - - if ( poly1->dgr >= poly2->dgr ){ - *sum = *poly1; - for ( i=0; i <= poly2->dgr; ++i ){ - tmp.cf[i+offset] = poly2->cf[i]; - } - } else { - *sum = *poly2; - for ( i=0; i <= poly1->dgr; ++i ){ - tmp.cf[i+offset] = poly1->cf[i]; - } - } - - for ( i=0; i <= sum->dgr; ++i ){ - sum->cf[i] += tmp.cf[i]; - } - return sum; -} - - -/* - * rt_poly_sub -- subtract two polynomials - */ -poly * -rt_poly_sub(poly1,poly2,diff) -register poly *poly1, *poly2, *diff; -{ - LOCAL poly tmp; - register int i, offset; - - offset = Abs(poly1->dgr - poly2->dgr); - - *diff = Zpoly; - tmp = Zpoly; - - if ( poly1->dgr >= poly2->dgr ){ - *diff = *poly1; - for ( i=0; i <= poly2->dgr; ++i ){ - tmp.cf[i+offset] = poly2->cf[i]; - } - } else { - diff->dgr = poly2->dgr; - for ( i=0; i <= poly1->dgr; ++i ){ - diff->cf[i+offset] = poly1->cf[i]; - } - tmp = *poly2; - } - - for ( i=0; i <= diff->dgr; ++i ){ - diff->cf[i] -= tmp.cf[i]; - } - return diff; -} - - -/* >>> s y n D i v ( ) <<< - * Divides any polynomial into any other polynomial using synthetic - * division. Both polynomials must have real coefficients. - */ -void -rt_poly_synthetic_division(dvdend,dvsor,quo,rem) -register poly *dvdend, *dvsor, *quo, *rem; -{ - register int div; - register int n; - - *quo = *dvdend; - *rem = Zpoly; - - if ((quo->dgr = dvdend->dgr - dvsor->dgr) < 0) - quo->dgr = -1; - if ((rem->dgr = dvsor->dgr - 1) > dvdend->dgr) - rem->dgr = dvdend->dgr; - - for ( n=0; n <= quo->dgr; ++n){ - quo->cf[n] /= dvsor->cf[0]; - for ( div=1; div <= dvsor->dgr; ++div){ - quo->cf[n+div] -= quo->cf[n] * dvsor->cf[div]; - } - } - for ( n=1; n<=(rem->dgr+1); ++n){ - rem->cf[n-1] = quo->cf[quo->dgr+n]; - quo->cf[quo->dgr+n] = 0; - } -} - - -/* >>> q u a d r a t i c ( ) <<< - * - * Uses the quadratic formula to find the roots (in `complex' form) - * of any quadratic equation with real coefficients. - */ -void -rt_poly_quadratic_roots( quadrat, root ) -register poly *quadrat; -register complex root[]; -{ - LOCAL fastf_t discrim, denom, rad; - - if( NEAR_ZERO( quadrat->cf[0], SMALL ) ) { - /* root = -cf[2] / cf[1] */ - if( NEAR_ZERO( quadrat->cf[1], SMALL ) ) { - /* No solution. Now what? */ - rt_log("rt_poly_quadratic_roots(): ERROR, no solution\n"); - return; - } - /* Fake it as a repeated root. */ - root[0].re = root[1].re = -quadrat->cf[2]/quadrat->cf[1]; - root[0].im = root[1].im = 0.0; - return; - } - /* What to do if cf[1] > SQRT_MAX_FASTF ? */ - - discrim = quadrat->cf[1]*quadrat->cf[1] - 4.0* quadrat->cf[0]*quadrat->cf[2]; - denom = 0.5 / quadrat->cf[0]; - if ( discrim >= 0.0 ){ - rad = sqrt( discrim ); - root[0].re = ( -quadrat->cf[1] + rad ) * denom; - root[1].re = ( -quadrat->cf[1] - rad ) * denom; - root[1].im = root[0].im = 0.0; - } else { - root[1].re = root[0].re = -quadrat->cf[1] * denom; - root[1].im = -(root[0].im = sqrt( -discrim ) * denom); - } -} - - -#define SQRT3 1.732050808 -#define THIRD 0.333333333333333333333333333 -#define INV_TWENTYSEVEN 0.037037037037037037037037037 -#define CUBEROOT( a ) (( (a) >= 0.0 ) ? pow( a, THIRD ) : -pow( -(a), THIRD )) - -/* >>> c u b i c ( ) <<< - * - * Uses the cubic formula to find the roots ( in `complex' form ) - * of any cubic equation with real coefficients. - * - * to solve a polynomial of the form: - * - * X**3 + c1*X**2 + c2*X + c3 = 0, - * - * first reduce it to the form: - * - * Y**3 + a*Y + b = 0, - * - * where - * Y = X + c1/3, - * and - * a = c2 - c1**2/3, - * b = ( 2*c1**3 - 9*c1*c2 + 27*c3 )/27. - * - * Then we define the value delta, D = b**2/4 + a**3/27. - * - * If D > 0, there will be one real root and two conjugate - * complex roots. - * If D = 0, there will be three real roots at least two of - * which are equal. - * If D < 0, there will be three unequal real roots. - * - * Returns 1 for success, 0 for fail. - */ -static int expecting_fpe = 0; -static jmp_buf abort_buf; -HIDDEN void rt_catch_FPE(sig) -int sig; -{ - if( !expecting_fpe ) - rt_bomb("unexpected SIGFPE!"); - if( !rt_g.rtg_parallel ) - (void)signal(SIGFPE, rt_catch_FPE); /* Renew handler */ - longjmp(abort_buf, 1); /* return error code */ -} - -int -rt_poly_cubic_roots( eqn, root ) -register poly *eqn; -register complex root[]; -{ - LOCAL fastf_t a, b, c1, c1_3rd, delta; - register int i; - static int first_time = 1; - - if( !rt_g.rtg_parallel ) { - /* abort_buf is NOT parallel! */ - if( first_time ) { - first_time = 0; - (void)signal(SIGFPE, rt_catch_FPE); - } - expecting_fpe = 1; - if( setjmp( abort_buf ) ) { - (void)signal(SIGFPE, rt_catch_FPE); - rt_log("rt: rt_poly_cubic_roots() Floating Point Error\n"); - return(0); /* FAIL */ - } - } - - c1 = eqn->cf[1]; - if( Abs(c1) > SQRT_MAX_FASTF ) return(0); /* FAIL */ - c1_3rd = c1 * THIRD; - a = eqn->cf[2] - c1*c1_3rd; - if( Abs(a) > SQRT_MAX_FASTF ) return(0); /* FAIL */ - b = (2.0*c1*c1*c1 - 9.0*c1*eqn->cf[2] + 27.0*eqn->cf[3])*INV_TWENTYSEVEN; - if( Abs(b) > SQRT_MAX_FASTF ) return(0); /* FAIL */ - - if( (delta = a*a) > SQRT_MAX_FASTF ) return(0); /* FAIL */ - delta = b*b*0.25 + delta*a*INV_TWENTYSEVEN; - - if ( delta > 0.0 ){ - LOCAL fastf_t r_delta, A, B; - - r_delta = sqrt( delta ); - A = B = -0.5 * b; - A += r_delta; - B -= r_delta; - - A = CUBEROOT( A ); - B = CUBEROOT( B ); - - root[2].re = root[1].re = -0.5 * ( root[0].re = A + B ); - - root[0].im = 0.0; - root[2].im = -( root[1].im = (A - B)*SQRT3*0.5 ); - } else if ( delta == 0.0 ){ - LOCAL fastf_t b_2; - b_2 = -0.5 * b; - - root[0].re = 2.0* CUBEROOT( b_2 ); - root[2].re = root[1].re = -0.5 * root[0].re; - root[2].im = root[1].im = root[0].im = 0.0; - } else { - LOCAL fastf_t phi, fact; - LOCAL fastf_t cs_phi, sn_phi_s3; - - if( a >= 0.0 ) { - fact = 0.0; - phi = 0.0; - cs_phi = 1.0; /* cos( phi ); */ - sn_phi_s3 = 0.0; /* sin( phi ) * SQRT3; */ - } else { - FAST fastf_t f; - a *= -THIRD; - fact = sqrt( a ); - if( (f = b * (-0.5) / (a*fact)) >= 1.0 ) { - phi = 0.0; - cs_phi = 1.0; /* cos( phi ); */ - sn_phi_s3 = 0.0; /* sin( phi ) * SQRT3; */ - } else if( f <= -1.0 ) { - phi = PI_DIV_3; - cs_phi = cos( phi ); - sn_phi_s3 = sin( phi ) * SQRT3; - } else { - phi = acos( f ) * THIRD; - cs_phi = cos( phi ); - sn_phi_s3 = sin( phi ) * SQRT3; - } - } - - root[0].re = 2.0*fact*cs_phi; - root[1].re = fact*( sn_phi_s3 - cs_phi); - root[2].re = fact*( -sn_phi_s3 - cs_phi); - root[2].im = root[1].im = root[0].im = 0.0; - } - for ( i=0; i < 3; ++i ) - root[i].re -= c1_3rd; - expecting_fpe = 0; - return(1); /* OK */ -} - - -/* >>> q u a r t i c ( ) <<< - * - * Uses the quartic formula to find the roots ( in `complex' form ) - * of any quartic equation with real coefficients. - * - * Returns 1 for success, 0 for fail. - */ -int -rt_poly_quartic_roots( eqn, root ) -register poly *eqn; -register complex root[]; -{ - LOCAL poly cube, quad1, quad2; - LOCAL complex u[3]; - LOCAL fastf_t U, p, q, q1, q2; -#define Max3(a,b,c) ((c)>((a)>(b)?(a):(b)) ? (c) : ((a)>(b)?(a):(b))) - - cube.dgr = 3; - cube.cf[0] = 1.0; - cube.cf[1] = -eqn->cf[2]; - cube.cf[2] = eqn->cf[3]*eqn->cf[1] - - 4*eqn->cf[4]; - cube.cf[3] = -eqn->cf[3]*eqn->cf[3] - - eqn->cf[4]*eqn->cf[1]*eqn->cf[1] - + 4*eqn->cf[4]*eqn->cf[2]; - - if( !rt_poly_cubic_roots( &cube, u ) ) { - return( 0 ); /* FAIL */ - } - if ( u[1].im != 0.0 ){ - U = u[0].re; - } else { - U = Max3( u[0].re, u[1].re, u[2].re ); - } - - p = eqn->cf[1]*eqn->cf[1]*0.25 + U - eqn->cf[2]; - U *= 0.5; - q = U*U - eqn->cf[4]; - if( p < 0 ) { - if( p < -1.0e-8 ) { - return(0); /* FAIL */ - } - p = 0; - } else { - p = sqrt( p ); - } - if( q < 0 ) { - if( q < -1.0e-8 ) { - return(0); /* FAIL */ - } - q = 0; - } else { - q = sqrt( q ); - } - - quad1.dgr = quad2.dgr = 2; - quad1.cf[0] = quad2.cf[0] = 1.0; - quad1.cf[1] = eqn->cf[1]*0.5; - quad2.cf[1] = quad1.cf[1] + p; - quad1.cf[1] -= p; - - q1 = U - q; - q2 = U + q; - - p = quad1.cf[1]*q2 + quad2.cf[1]*q1 - eqn->cf[3]; - if( Abs( p ) < 1.0e-8){ - quad1.cf[2] = q1; - quad2.cf[2] = q2; - } else { - q = quad1.cf[1]*q1 + quad2.cf[1]*q2 - eqn->cf[3]; - if( Abs( q ) < 1.0e-8 ){ - quad1.cf[2] = q2; - quad2.cf[2] = q1; - } else { - return(0); /* FAIL */ - } - } - - rt_poly_quadratic_roots( &quad1, root ); - rt_poly_quadratic_roots( &quad2, &root[2] ); - return(1); /* SUCCESS */ -} - - -/* - * R T _ P R _ P O L Y - */ -void -rt_pr_poly(title, eqn) -char *title; -register poly *eqn; -{ - register int n; - register int exp; - struct rt_vls str; - char buf[48]; - - rt_vls_init( &str ); - rt_vls_extend( &str, 196 ); - rt_vls_strcat( &str, title ); - sprintf(buf, " polynomial, degree = %d\n", eqn->dgr); - rt_vls_strcat( &str, buf ); - - exp = eqn->dgr; - for ( n=0; n<=eqn->dgr; n++,exp-- ) { - register double coeff = eqn->cf[n]; - if( n > 0 ) { - if( coeff < 0 ) { - rt_vls_strcat( &str, " - " ); - coeff = -coeff; - } else { - rt_vls_strcat( &str, " + " ); - } - } - sprintf( buf, "%g", coeff ); - rt_vls_strcat( &str, buf ); - if( exp > 1 ) { - sprintf( buf, " *X^%d", exp ); - rt_vls_strcat( &str, buf ); - } else if( exp == 1 ) { - - rt_vls_strcat( &str, " *X" ); - } else { - /* For constant term, add nothing */ - } - } - rt_vls_strcat( &str, "\n" ); - rt_log( "%s", rt_vls_addr(&str) ); - rt_vls_free( &str ); -} - -/* - * R T _ P R _ R O O T S - */ -void -rt_pr_roots( title, roots, n ) -char *title; -complex roots[]; -int n; -{ - register int i; - - rt_log("%s: %d roots:\n", title, n ); - for( i=0; i= 0 ? (a) : -(a)) -#define Max( a, b ) ((a) > (b) ? (a) : (b)) - -/* error return value for 'polyMul' */ -#define PM_NULL ((poly *)0) - -/* polynomial data type */ -typedef struct { - int dgr; - double cf[MAXP+1]; -} poly; - - -/* library functions in polylib.c */ -extern poly *polyAdd(), *polySub(), *polyMul(), *polyScal(); -extern void quadratic(), synDiv(), prntpoly(), pr_poly(); -extern int quartic(), cubic(); - -#endif /* POLYNO_H */ diff --git a/librt/snoise.c b/librt/snoise.c deleted file mode 100644 index fd773331c4e..00000000000 --- a/librt/snoise.c +++ /dev/null @@ -1,1797 +0,0 @@ -/* S N O I S E . C - * - * Signed noise functions: - * - * Assembled from Darwin Peachey's examples in - * "Texturing and Modeling: A Procedural Approach" - * and from examples provided by F. Kenton Musgrave. - * - * noise_g Gradient Noise - * noise_v Value Noise - * noise_gv Gradient Value Noise - * noise_vc Value Convolution Noise - * noise_sc Sparse Convolution Noise - * noise_perlin Robert Skinner's Perlin-style "Noise" function - * noise_vec Vector-valued noise - * - * Spectral Noise functions - * noise_fbm fractional Brownian motion. Based on noise_perlin - * noise_turb turbulence. Based on noise_perlin - * - * Author - Lee A. Butler - * - * Source - - * The U. S. Army Research Laboratory - * Aberdeen Proving Ground, Maryland 21005-5068 USA - * - * Distribution Status - - * Public Domain, Distribution Unlimited. - */ -#ifndef lint -static char RCSid[] = "@(#)$Header$ (ARL)"; -#endif - -#include "conf.h" - -#include -#include -#include "machine.h" -#include "externs.h" -#include "vmath.h" -#include "bu.h" -#include "raytrace.h" - - -#define TABSIZE 256 /* must be power of 2 */ -#define TABMASK (TABSIZE-1) -#define PERM(x) perm[ (x)&TABMASK ] -#define INDEX(ix,iy,iz) PERM( (ix)+PERM( (iy)+PERM(iz) ) ) -#define LERP(t, x0, x1) ( (x0) + (t) * ((x1) - (x0)) ) -#define SMOOTHSTEP(x) ( (x) * (x) * (3 - 2*(x)) ) - -#define CLAMP(_x,_a,_b) (_x < _a ? _a : (_x > _b ? _b : _x)) -#define FLOOR(x) ( (int)(x) - ( (x) < 0 && (x) != (int)(x) ) ) -#define CEIL(x) ( (int)(x) + ( (x) > 0 && (x) != (int)(x) ) ) - - -/* - * The noise function has a finite domain, which can be exceeded when - * using fractal textures with very high frequencies. This routine is - * designed to extend the effective domain of the function, albeit by - * introducing periodicity. -FKM 4/93 - */ -#define MAXVAL 2147483647. /* max value for noise integers */ -#define TWICE_MAXVAL 4294967294. -#define FILTER_ARGS( src) {\ - register int i; \ - point_t dst; \ - \ - for (i=0 ; i < 3 ; i++) { \ - /* assure values are positive */ \ - if (src[i] < 0) dst[i] = -src[i]; \ - else dst[i] = src[i]; \ - \ - /* fold space */ \ - while (dst[i] > MAXVAL || dst[i]<0) \ - if (dst[i] > MAXVAL) \ - dst[i] = TWICE_MAXVAL - dst[i]; \ - else \ - dst[i] = -dst[i]; \ - } \ -\ - x = dst[0]; ix = FLOOR(x); fx = x - ix; \ - y = dst[1]; iy = FLOOR(y); fy = y - iy; \ - z = dst[2]; iz = FLOOR(z); fz = z - iz; \ -} - -/* Permutation table */ -static unsigned char perm[TABSIZE] = { - 101, 215, 90, 114, 81, 226, 3, 149, 40, 98, 176, 15, 230, 41, - 40, 136, 154, 149, 69, 99, 75, 190, 76, 19, 103, 219, 241, 169, - 216, 0, 118, 136, 201, 67, 251, 78, 153, 155, 54, 226, 77, 38, - 86, 99, 164, 192, 154, 136, 117, 167, 83, 242, 94, 241, 1, 132, - 69, 6, 151, 52, 224, 15, 66, 77, 228, 127, 181, 73, 221, 172, - 117, 245, 198, 96, 58, 90, 76, 171, 184, 144, 211, 99, 209, 216, - 46, 241, 108, 133, 16, 233, 225, 194, 102, 176, 194, 103, 32, 124, - 56, 223, 135, 0, 220, 4, 208, 62, 80, 247, 239, 207, 125, 56, - 147, 73, 82, 66, 44, 0, 11, 61, 106, 179, 56, 129, 17, 100, - 122, 55, 56, 234, 89, 49, 54, 162, 13, 200, 7, 113, 45, 238, - 232, 121, 223, 178, 238, 116, 102, 228, 177, 214, 189, 166, 173, 147, - 69, 239, 169, 12, 95, 158, 38, 96, 165, 6, 215, 19, 190, 65, - 230, 96, 81, 54, 166, 64, 58, 64, 241, 34, 69, 140, 82, 221, - 75, 174, 213, 224, 166, 18, 230, 64, 156, 215, 212, 95, 193, 27, - 217, 143, 219, 87, 177, 88, 228, 245, 31, 251, 14, 157, 9, 96, - 134, 72, 143, 155, 209, 114, 6, 120, 72, 74, 50, 4, 212, 146, - 26, 187, 30, 57, 242, 189, 210, 211, 64, 65, 86, 99, 134, 68, - 102, 222, 11, 75, 100, 143, 79, 210, 121, 23, 67, 234, 250, 85, - 230, 61, 95, 192 -}; - -/* Coefficients of (Catmull spline) basis matrix. */ -#define CR00 -0.5 -#define CR01 1.5 -#define CR02 -1.5 -#define CR03 0.5 - -#define CR10 1.0 -#define CR11 -2.5 -#define CR12 2.0 -#define CR13 -0.5 - -#define CR20 -0.5 -#define CR21 0.0 -#define CR22 0.5 -#define CR23 0.0 - -#define CR30 0.0 -#define CR31 1.0 -#define CR32 0.0 -#define CR33 0.0 - - -static double -spline(x, nknots, knot) -double x; -int nknots; -double *knot; -{ - int span; - int nspans = nknots - 3; - double c0, c1, c2, c3; /* coefficients of the cubic.*/ - - if (nspans < 1) { /* illegal */ - /* fprintf(stderr, "Spline has too few knots.\n"); */ - abort(); - return 0; - } - - /* Find the appropriate 4-point span of the spline. */ - x = CLAMP(x, 0, 1) * nspans; - span = (int) x; - if (span >= nknots - 3) - span = nknots - 3; - x -= span; - knot += span; - - /* Evaluate the span cubic at x using Horner's rule. - * Note: Coefficients which are 0 are not calculated - */ - c3 = CR00*knot[0] + CR01*knot[1] + CR02*knot[2] + CR03*knot[3]; - c2 = CR10*knot[0] + CR11*knot[1] + CR12*knot[2] + CR13*knot[3]; - c1 = CR20*knot[0] /* +CR21*knot[1]*/ + CR22*knot[2]/*+ CR23*knot[3]*/; - c0 =/*CR30*knot[0]+*/ CR31*knot[1] /*+ CR32*knot[2] + CR33*knot[3] */; - - return ((c3*x + c2)*x + c1)*x + c0; -} - -static double -catrom2(d) -double d; -{ -#define SAMPRATE 100 -#define NENTRIES (4*SAMPRATE+1) - double x; - int i; - static double table[NENTRIES]; - static int initialized=0; - - if (d >= 4.) return 0; - - if (!initialized) { - initialized = 1; - for (i=0 ; i < NENTRIES ; i++) { - x = sqrt( (double)i / (double)SAMPRATE ); - if (x < 1) - table[i] = 0.5 * (2.+x*x*(-5+x*3)); - else - table[i] = 0.5 * (4+x*(-8+x*(5-x))); - } - } - - d = d*SAMPRATE + 0.5; - i = FLOOR(d); - if (i >= NENTRIES) return 0.; - return table[i]; -} - -/* - * This is our table of random numbers. Rather than calling drand48() or - * random() or rand() we just pick numbers out of this table. This table - * has 4096 entries. - */ -CONST float rt_rand_table[RT_RAND_TABSIZE] = { - 0.39646477, 0.84048537, 0.35333610, 0.44658343, 0.31869277, 0.88642843, - 0.01558285, 0.58409022, 0.15936863, 0.38371587, 0.69100437, 0.05885891, - 0.89985431, 0.16354595, 0.15907150, 0.53306471, 0.60414419, 0.58269902, - 0.26997112, 0.39047820, 0.29340057, 0.74237741, 0.29852561, 0.07553808, - 0.40498263, 0.85737794, 0.94196832, 0.66283066, 0.84647578, 0.00275508, - 0.46237925, 0.53259602, 0.78787662, 0.26561223, 0.98275226, 0.30678513, - 0.60085514, 0.60871565, 0.21243880, 0.88589513, 0.30465710, 0.15185986, - 0.33766190, 0.38747695, 0.64360983, 0.75355328, 0.60361610, 0.53162825, - 0.45936032, 0.65248845, 0.32718116, 0.94637049, 0.36803987, 0.94389034, - 0.00742826, 0.51659995, 0.27277095, 0.02429916, 0.59195450, 0.20496351, - 0.87769335, 0.05936869, 0.26084255, 0.30282918, 0.89149522, 0.49819806, - 0.71002558, 0.28641399, 0.86492358, 0.67554067, 0.45848997, 0.95963556, - 0.77467541, 0.37655128, 0.22863912, 0.35453388, 0.30031825, 0.66976583, - 0.71896657, 0.56595451, 0.82446531, 0.39061191, 0.81876631, 0.84400846, - 0.18046777, 0.94339589, 0.42488677, 0.52066578, 0.06564375, 0.91350817, - 0.88258457, 0.76136413, 0.39892255, 0.68825684, 0.76154830, 0.40500880, - 0.12525114, 0.48463390, 0.22246255, 0.87312117, 0.52882190, 0.00141396, - 0.86051381, 0.01869740, 0.81489194, 0.24288402, 0.31457184, 0.96573252, - 0.93556011, 0.80943170, 0.49210915, 0.22013551, 0.57635374, 0.28902906, - 0.32106698, 0.26132334, 0.17398786, 0.00181730, 0.04478414, 0.24117455, - 0.41545119, 0.70162465, 0.22184569, 0.50391034, 0.06703021, 0.39306344, - 0.47947653, 0.21814221, 0.21951128, 0.91620319, 0.35022175, 0.19269394, - 0.21123498, 0.63368163, 0.05356539, 0.78341141, 0.03067392, 0.44409660, - 0.17641289, 0.93218022, 0.90964827, 0.47284483, 0.87169546, 0.69556736, - 0.93018962, 0.45509056, 0.39856678, 0.89330400, 0.69354675, 0.83873438, - 0.73969914, 0.65126955, 0.67815424, 0.57721232, 0.27311828, 0.93538805, - 0.66193889, 0.04741251, 0.37303839, 0.61819478, 0.14862799, 0.37730558, - 0.64459140, 0.02568782, 0.84138051, 0.07701881, 0.74263392, 0.25598614, - 0.90184423, 0.37764759, 0.31952992, 0.21143499, 0.64853250, 0.25131508, - 0.22864432, 0.25092218, 0.94322116, 0.13671508, 0.27006077, 0.54870719, - 0.32402145, 0.86508002, 0.29670335, 0.68005934, 0.83314725, 0.87630751, - 0.64966456, 0.07311548, 0.89854697, 0.25358126, 0.61131863, 0.84188993, - 0.83201930, 0.37258709, 0.75704822, 0.10921491, 0.85081198, 0.55932933, - 0.85758046, 0.34309698, 0.69176960, 0.34519729, 0.89358465, 0.95870788, - 0.12173520, 0.98154436, 0.05502479, 0.61477628, 0.03812668, 0.37647260, - 0.52560486, 0.28192396, 0.56053796, 0.60724297, 0.81643958, 0.44653499, - 0.02711168, 0.47185401, 0.28489831, 0.29272368, 0.19568014, 0.01757942, - 0.82959923, 0.57316544, 0.10478060, 0.73297648, 0.11909752, 0.22388840, - 0.94725286, 0.73905162, 0.82135698, 0.82613156, 0.25097997, 0.25649573, - 0.33819320, 0.38822167, 0.52731848, 0.26640913, 0.40122053, 0.87089457, - 0.04589179, 0.29457382, 0.39424656, 0.56000239, 0.31082114, 0.82257687, - 0.47548149, 0.09114359, 0.26240109, 0.91704409, 0.97843590, 0.33209335, - 0.90207403, 0.24063487, 0.37316437, 0.75228603, 0.45782739, 0.90062788, - 0.42289451, 0.56592348, 0.97046216, 0.35387088, 0.43150651, 0.17909136, - 0.21531911, 0.33717677, 0.45436775, 0.04472570, 0.68388245, 0.06214914, - 0.55137528, 0.31617124, 0.26821790, 0.50076599, 0.06253696, 0.96482654, - 0.72957579, 0.80631556, 0.15197441, 0.70501796, 0.72995745, 0.84256584, - 0.61792858, 0.95544778, 0.47938813, 0.52708549, 0.02938659, 0.60761211, - 0.04756937, 0.59276965, 0.40264693, 0.90439440, 0.33013107, 0.31083745, - 0.85494474, 0.73414953, 0.09862653, 0.34045859, 0.40173949, 0.98516706, - 0.62959478, 0.83196239, 0.43711973, 0.79122669, 0.05946186, 0.13513490, - 0.54382390, 0.93045165, 0.88489902, 0.38989459, 0.41455752, 0.85686265, - 0.28675858, 0.10775672, 0.57667600, 0.77760218, 0.14225093, 0.80867719, - 0.74809422, 0.41002486, 0.47511577, 0.34445954, 0.41679885, 0.02364775, - 0.52581101, 0.71789154, 0.23602268, 0.69486508, 0.77995061, 0.76275041, - 0.12156272, 0.37051796, 0.92825143, 0.76600634, 0.28957953, 0.98366071, - 0.42810692, 0.03781309, 0.59832645, 0.01559561, 0.75672384, 0.13067900, - 0.37456716, 0.21749617, 0.57779700, 0.20566339, 0.78702094, 0.74836919, - 0.45533640, 0.28526833, 0.85022008, 0.45000752, 0.83501192, 0.43024603, - 0.33099577, 0.95087211, 0.24887282, 0.35679405, 0.41787068, 0.49705314, - 0.88529295, 0.56316535, 0.68550024, 0.96775937, 0.18236876, 0.56780538, - 0.02966464, 0.20396187, 0.15776493, 0.54759438, 0.20399063, 0.21128662, - 0.06729404, 0.19309738, 0.90152153, 0.78204947, 0.30258985, 0.42200777, - 0.98952543, 0.79712060, 0.38747161, 0.18194028, 0.22473139, 0.73738293, - 0.53350707, 0.66473599, 0.20183098, 0.46882384, 0.56649724, 0.09746379, - 0.27702369, 0.01731167, 0.56388407, 0.88693968, 0.07019007, 0.06202093, - 0.25330677, 0.30706457, 0.89663986, 0.13407325, 0.94770981, 0.42945465, - 0.07812638, 0.78366722, 0.44367027, 0.42118211, 0.70188645, 0.18982070, - 0.12622022, 0.24805934, 0.65792792, 0.13750542, 0.32277563, 0.82251428, - 0.01384098, 0.92370735, 0.37635347, 0.58206852, 0.26137310, 0.41370918, - 0.68769666, 0.75372087, 0.17132775, 0.12920425, 0.15899024, 0.29275826, - 0.85759328, 0.27887709, 0.68842366, 0.83888960, 0.81497662, 0.04861429, - 0.07294221, 0.43097879, 0.74982575, 0.15195442, 0.41226738, 0.10251594, - 0.76112259, 0.64178617, 0.73490959, 0.54741372, 0.09988188, 0.73232412, - 0.52882975, 0.51947453, 0.28649967, 0.94732048, 0.42613666, 0.49429513, - 0.70247499, 0.25822254, 0.81136807, 0.88043637, 0.00112771, 0.10709923, - 0.27325164, 0.98398617, 0.36112317, 0.81679947, 0.85526165, 0.28941299, - 0.99948609, 0.00875904, 0.06378548, 0.96769399, 0.68316238, 0.30155620, - 0.72542307, 0.92358475, 0.64901375, 0.91576760, 0.53612168, 0.97667749, - 0.23347556, 0.41855917, 0.58187885, 0.78029750, 0.20017003, 0.90105181, - 0.12913710, 0.92866606, 0.76486375, 0.47974114, 0.55637733, 0.01208822, - 0.38373625, 0.37766509, 0.72855173, 0.54380673, 0.67755567, 0.58976557, - 0.28276397, 0.10782662, 0.29873343, 0.52877053, 0.90501287, 0.43730920, - 0.95117220, 0.23234421, 0.51999770, 0.61653656, 0.18827655, 0.97925561, - 0.50319736, 0.16005923, 0.53953297, 0.19256178, 0.24906409, 0.33018760, - 0.15088122, 0.62258310, 0.21172293, 0.94382294, 0.62201213, 0.21367579, - 0.47367458, 0.84272338, 0.44118411, 0.11902201, 0.92041789, 0.67426564, - 0.16472098, 0.78002358, 0.70952318, 0.18482036, 0.77673340, 0.21519147, - 0.25450596, 0.66286094, 0.96061541, 0.94208271, 0.80800541, 0.79734687, - 0.98645123, 0.69988945, 0.66452100, 0.11786686, 0.68969624, 0.15539494, - 0.32768682, 0.45121898, 0.41156526, 0.50746722, 0.84627410, 0.82489765, - 0.56879108, 0.69197801, 0.14044106, 0.84732625, 0.70125160, 0.20754978, - 0.06695332, 0.43922141, 0.84497649, 0.30101358, 0.33906069, 0.53742714, - 0.03292710, 0.87869117, 0.18144436, 0.96476526, 0.72751967, 0.44501378, - 0.79927650, 0.49539015, 0.33107384, 0.15359297, 0.73721792, 0.44891119, - 0.35028905, 0.61833266, 0.94150374, 0.13788360, 0.75867440, 0.04480321, - 0.58663764, 0.65751747, 0.00990359, 0.96680331, 0.33626699, 0.84459331, - 0.61178133, 0.01235191, 0.85230017, 0.51120210, 0.23446091, 0.67044493, - 0.80213760, 0.06027408, 0.76044009, 0.61856870, 0.92402669, 0.36177601, - 0.20350460, 0.32559878, 0.55615096, 0.87441762, 0.94812956, 0.03036690, - 0.40538316, 0.76681572, 0.66196149, 0.86851554, 0.79405332, 0.56904565, - 0.58707982, 0.57556754, 0.57311427, 0.28071644, 0.74909981, 0.41933082, - 0.80136088, 0.00111404, 0.16827220, 0.53737595, 0.96463389, 0.92216821, - 0.09812967, 0.79352987, 0.32609668, 0.84324004, 0.04262168, 0.37991233, - 0.16672271, 0.40084932, 0.15288743, 0.32741607, 0.48020341, 0.64494132, - 0.12378731, 0.63773727, 0.74284151, 0.89023276, 0.89511602, 0.09476246, - 0.87901346, 0.39630206, 0.01325157, 0.77360169, 0.50697619, 0.11875299, - 0.07917336, 0.55373911, 0.07678998, 0.77695809, 0.15475466, 0.41119012, - 0.02868890, 0.81818189, 0.63057161, 0.15637346, 0.61123570, 0.55551929, - 0.88786545, 0.48341520, 0.33705686, 0.11591437, 0.08057005, 0.47721493, - 0.13099915, 0.66802435, 0.68941569, 0.00186024, 0.28043875, 0.28696567, - 0.78111478, 0.97353869, 0.14759043, 0.34874942, 0.46671974, 0.34925594, - 0.89451620, 0.39183381, 0.17273889, 0.57633823, 0.92675197, 0.31266868, - 0.49794719, 0.22441621, 0.43281549, 0.71798237, 0.56260682, 0.92951971, - 0.91421196, 0.39609148, 0.50271111, 0.62784723, 0.42192285, 0.61555041, - 0.71070714, 0.24412578, 0.87688152, 0.25047744, 0.36988002, 0.56628297, - 0.21786419, 0.27199749, 0.01926812, 0.41508006, 0.63293902, 0.65415704, - 0.26981691, 0.47568941, 0.06062798, 0.20574399, 0.27379098, 0.25985808, - 0.84287801, 0.75766222, 0.23469406, 0.71481863, 0.55056790, 0.77546769, - 0.69366582, 0.48165395, 0.83847807, 0.11059789, 0.80322752, 0.86518332, - 0.80342971, 0.79956996, 0.51615567, 0.23402665, 0.87779053, 0.63146171, - 0.99512074, 0.98836098, 0.28000438, 0.45942571, 0.45960435, 0.96967357, - 0.53161393, 0.77654627, 0.70570180, 0.43782638, 0.43191766, 0.02038338, - 0.19004420, 0.59289301, 0.13056690, 0.04003430, 0.59268429, 0.30004982, - 0.85682375, 0.54608257, 0.34994998, 0.80346149, 0.33174901, 0.57462389, - 0.66467700, 0.04838776, 0.61427891, 0.70450838, 0.79870563, 0.91906587, - 0.90011267, 0.16049770, 0.07501304, 0.53333152, 0.71782816, 0.96830335, - 0.46142079, 0.39409409, 0.74205247, 0.97017616, 0.07466423, 0.00396880, - 0.58934667, 0.99720926, 0.58605443, 0.71450257, 0.07735063, 0.57777176, - 0.68445195, 0.45805182, 0.43843468, 0.75711266, 0.26837314, 0.34164253, - 0.39050042, 0.44338766, 0.40946091, 0.25185489, 0.11267285, 0.97801954, - 0.82009196, 0.90234104, 0.46356207, 0.66688069, 0.23750772, 0.56604339, - 0.95380290, 0.88238627, 0.92761603, 0.87819739, 0.36008849, 0.14993526, - 0.23871129, 0.02006514, 0.59040581, 0.72558514, 0.57712148, 0.63396868, - 0.70329103, 0.43137116, 0.31465760, 0.45511417, 0.39940817, 0.58188493, - 0.40628451, 0.74308326, 0.97149811, 0.36783513, 0.03334644, 0.66450604, - 0.57675455, 0.98217792, 0.48378778, 0.21773421, 0.82009278, 0.52050219, - 0.03005330, 0.64034290, 0.37184202, 0.46378676, 0.51480345, 0.43166945, - 0.63149750, 0.57813972, 0.15443478, 0.72663540, 0.06082055, 0.82803596, - 0.12397288, 0.23570890, 0.59340780, 0.39786342, 0.40065108, 0.92158799, - 0.66062077, 0.91787441, 0.17697285, 0.00918329, 0.61035890, 0.27714125, - 0.66997465, 0.45554272, 0.21857858, 0.02568675, 0.06757884, 0.16490334, - 0.26484245, 0.88407773, 0.58429941, 0.68611608, 0.75178376, 0.62053331, - 0.51920603, 0.93093194, 0.86157311, 0.86773765, 0.89415988, 0.29835242, - 0.55578762, 0.94440944, 0.72822905, 0.06934604, 0.31213773, 0.54624471, - 0.50631544, 0.77525638, 0.65406240, 0.16871710, 0.61095184, 0.56033928, - 0.24714263, 0.22761744, 0.79364983, 0.11285600, 0.61761327, 0.22499263, - 0.50691925, 0.17943005, 0.27910714, 0.65678093, 0.64408139, 0.79398722, - 0.84840117, 0.12803671, 0.93833231, 0.32620397, 0.78173470, 0.80191970, - 0.65537508, 0.10627944, 0.69074281, 0.59494370, 0.92864953, 0.42493397, - 0.89982518, 0.69431363, 0.98919228, 0.99236097, 0.41161203, 0.03012915, - 0.94302675, 0.71058484, 0.75897544, 0.31802531, 0.66040998, 0.75409221, - 0.17377417, 0.49430617, 0.72906845, 0.23812349, 0.00534743, 0.60928373, - 0.45668616, 0.09767355, 0.61494335, 0.10446236, 0.20843376, 0.50872015, - 0.87316633, 0.47013591, 0.25844337, 0.20276389, 0.54166009, 0.42224112, - 0.09816068, 0.19333972, 0.95124378, 0.32691462, 0.26088813, 0.72628485, - 0.54834457, 0.77973233, 0.12514232, 0.05238664, 0.69970434, 0.17129662, - 0.48737869, 0.91037489, 0.07998871, 0.16452908, 0.53800937, 0.29948063, - 0.55422577, 0.21020295, 0.02739935, 0.25399627, 0.57120802, 0.24000881, - 0.11774320, 0.33810484, 0.14748690, 0.22051884, 0.63259469, 0.52898741, - 0.81488187, 0.31214854, 0.04166442, 0.01178180, 0.77614946, 0.45558760, - 0.66076230, 0.27389164, 0.64036896, 0.57078570, 0.76501193, 0.25159684, - 0.74658299, 0.48265951, 0.87039814, 0.28268290, 0.31429221, 0.73108738, - 0.55014797, 0.27930268, 0.64367509, 0.64395311, 0.23733366, 0.02054776, - 0.81614523, 0.28009259, 0.39549466, 0.32990267, 0.88445160, 0.92871285, - 0.79019956, 0.34517222, 0.28648062, 0.89354353, 0.89462838, 0.19452477, - 0.14186463, 0.91962389, 0.44289582, 0.38952085, 0.49719762, 0.22982353, - 0.99052801, 0.75256173, 0.35398875, 0.13417099, 0.93994345, 0.97088491, - 0.15145498, 0.23067090, 0.68079389, 0.06146777, 0.06999866, 0.52000537, - 0.32453850, 0.43101508, 0.18499239, 0.26243931, 0.47672436, 0.50326530, - 0.82668450, 0.08310974, 0.78932407, 0.38344360, 0.01115167, 0.91702712, - 0.82245722, 0.94873009, 0.09789368, 0.62852837, 0.37634165, 0.44663951, - 0.82864701, 0.05819955, 0.73995364, 0.89338660, 0.11417093, 0.37565563, - 0.75531821, 0.67543002, 0.69414089, 0.74687905, 0.68142244, 0.04385994, - 0.51100119, 0.77355292, 0.22839082, 0.06815969, 0.05094052, 0.00808088, - 0.94285300, 0.88879444, 0.33144257, 0.18649060, 0.18456602, 0.22334871, - 0.80063730, 0.01252259, 0.02064522, 0.37546320, 0.38256972, 0.24588788, - 0.45865628, 0.94715004, 0.80626959, 0.13753732, 0.97789551, 0.93135690, - 0.32503039, 0.57205053, 0.43713130, 0.37169472, 0.92359696, 0.54477280, - 0.22771783, 0.71945435, 0.16115116, 0.25758243, 0.29726726, 0.93445346, - 0.27133177, 0.36625920, 0.13804928, 0.19654112, 0.97941075, 0.97407387, - 0.46617237, 0.39003165, 0.47168203, 0.32578068, 0.51571010, 0.75585681, - 0.91173252, 0.81464588, 0.24271236, 0.87835430, 0.01931690, 0.63908696, - 0.33705804, 0.47851562, 0.26493030, 0.01715658, 0.53805493, 0.95698606, - 0.76987435, 0.33433295, 0.05839757, 0.79255870, 0.69877446, 0.52635188, - 0.85949752, 0.42037089, 0.77838817, 0.71323587, 0.46137467, 0.14416121, - 0.14052842, 0.94796973, 0.98743793, 0.67506477, 0.36351923, 0.53379089, - 0.26690813, 0.62243055, 0.84943780, 0.82378338, 0.27679276, 0.41757404, - 0.32523395, 0.63698809, 0.68577877, 0.84739694, 0.67453815, 0.31379589, - 0.15186396, 0.37020105, 0.70158059, 0.16282464, 0.66481828, 0.33021097, - 0.08527924, 0.31027457, 0.15759895, 0.87033761, 0.31447165, 0.20500829, - 0.83470948, 0.94522718, 0.99684706, 0.72804066, 0.88038525, 0.17993617, - 0.62488891, 0.73909916, 0.06873046, 0.11556673, 0.76069763, 0.21640605, - 0.15796581, 0.48442746, 0.09468491, 0.92543900, 0.12134504, 0.80691614, - 0.94569913, 0.82037552, 0.26075648, 0.68414236, 0.33516296, 0.43606808, - 0.03136254, 0.44873334, 0.69641195, 0.18220942, 0.14579246, 0.24633508, - 0.81187530, 0.28393200, 0.84646005, 0.45769613, 0.52241915, 0.71916830, - 0.45972586, 0.58219356, 0.77058682, 0.74303172, 0.75418072, 0.12514247, - 0.34625224, 0.83923767, 0.69874176, 0.64076205, 0.35110990, 0.33977209, - 0.88909188, 0.24841571, 0.51734057, 0.66385747, 0.82095555, 0.67049021, - 0.56372325, 0.56861555, 0.92660721, 0.03544243, 0.28554481, 0.91515749, - 0.65079961, 0.38916163, 0.91714585, 0.93924922, 0.07104736, 0.49253405, - 0.92021835, 0.59492656, 0.61510267, 0.12424425, 0.00754476, 0.19782523, - 0.19883383, 0.71731799, 0.25937293, 0.83351979, 0.70701056, 0.83824376, - 0.30481323, 0.08143559, 0.72728114, 0.26186776, 0.97840780, 0.81913488, - 0.08659232, 0.33743060, 0.70536365, 0.76283466, 0.24527282, 0.00551574, - 0.27899429, 0.03136477, 0.39592097, 0.37734966, 0.15993748, 0.22152365, - 0.06417676, 0.78249634, 0.01609090, 0.09766092, 0.69161364, 0.82481149, - 0.83044214, 0.13766314, 0.62567611, 0.34648124, 0.50469635, 0.46731761, - 0.87926085, 0.51360233, 0.78158268, 0.54688146, 0.13776192, 0.12223265, - 0.37783116, 0.75481451, 0.67187032, 0.93993968, 0.17106400, 0.41689789, - 0.90095430, 0.59401434, 0.90337106, 0.04283122, 0.38839579, 0.48373914, - 0.26639679, 0.51098156, 0.53388900, 0.11106455, 0.10973833, 0.65497409, - 0.35861673, 0.41304702, 0.38684881, 0.38861834, 0.71072661, 0.35234868, - 0.91908994, 0.90923215, 0.30708016, 0.53153287, 0.04584024, 0.16035961, - 0.08241219, 0.68597392, 0.24152430, 0.14540544, 0.85151099, 0.96272446, - 0.15209020, 0.66167826, 0.76963298, 0.18306187, 0.64496103, 0.88322602, - 0.17004280, 0.22970732, 0.11707908, 0.03975910, 0.18739250, 0.10153946, - 0.78273129, 0.31261781, 0.96681610, 0.48529984, 0.11189916, 0.35750204, - 0.98943970, 0.24222239, 0.23454868, 0.44351865, 0.09693648, 0.47525380, - 0.51832390, 0.71061832, 0.60536174, 0.69772938, 0.79641372, 0.14991409, - 0.59948054, 0.57581969, 0.82693393, 0.52655141, 0.18396973, 0.16417968, - 0.54113774, 0.03120916, 0.43118240, 0.03705477, 0.84939584, 0.16805586, - 0.75367692, 0.16629029, 0.16355482, 0.48379094, 0.89362949, 0.38454127, - 0.13249740, 0.69961885, 0.06388892, 0.25840339, 0.22312307, 0.62661127, - 0.60324095, 0.78472088, 0.88830085, 0.69151093, 0.88099625, 0.50713854, - 0.90149956, 0.75684538, 0.86599794, 0.06284109, 0.58778359, 0.23410381, - 0.47837951, 0.33577834, 0.80380595, 0.82515916, 0.48783997, 0.97754757, - 0.60459354, 0.91419468, 0.68190788, 0.11616024, 0.14968469, 0.01484386, - 0.37683694, 0.30244626, 0.94105415, 0.36542821, 0.11073204, 0.73076121, - 0.42387834, 0.88505394, 0.77469184, 0.24249763, 0.12031183, 0.14650880, - 0.18774278, 0.69582623, 0.18594764, 0.69561026, 0.59282427, 0.72189474, - 0.78652647, 0.12374727, 0.08488700, 0.33580533, 0.27487649, 0.46930717, - 0.12223051, 0.12011581, 0.41926361, 0.23858772, 0.01560299, 0.93002825, - 0.08669509, 0.06738898, 0.33651073, 0.10413094, 0.83875428, 0.52611200, - 0.18995741, 0.71445975, 0.95188495, 0.05483756, 0.11389985, 0.30728189, - 0.85544678, 0.51712048, 0.30763636, 0.11454252, 0.89421276, 0.41433937, - 0.96683330, 0.73914275, 0.50028227, 0.70423969, 0.67976363, 0.60954897, - 0.21377238, 0.23127718, 0.28652452, 0.97570832, 0.76138373, 0.58160878, - 0.94258919, 0.18452705, 0.40439755, 0.31681769, 0.68522652, 0.37554639, - 0.64835593, 0.87804841, 0.85441157, 0.41283934, 0.88924031, 0.36214337, - 0.75029458, 0.79508875, 0.70715172, 0.99521715, 0.85166791, 0.83045884, - 0.93819323, 0.62410306, 0.41186321, 0.53502520, 0.81646770, 0.33511208, - 0.88094574, 0.41259867, 0.91696616, 0.33103261, 0.16989513, 0.53980369, - 0.21089029, 0.15537184, 0.52074117, 0.64391736, 0.71001110, 0.96722336, - 0.62655582, 0.19352873, 0.19652090, 0.82619206, 0.56010807, 0.88545505, - 0.75812174, 0.34549266, 0.12631225, 0.40808958, 0.95156690, 0.87926432, - 0.01091339, 0.48462036, 0.65641666, 0.90953759, 0.81884974, 0.08313831, - 0.13739252, 0.03700178, 0.99623409, 0.34896134, 0.66954850, 0.02122909, - 0.12884641, 0.43257745, 0.78478221, 0.48941593, 0.07950125, 0.51011847, - 0.83836371, 0.08273074, 0.44183966, 0.11850539, 0.06567619, 0.04636041, - 0.37731752, 0.54127132, 0.57587074, 0.21930263, 0.79686471, 0.61110745, - 0.11239614, 0.43445976, 0.23344420, 0.55955468, 0.78025702, 0.90587759, - 0.51424044, 0.70176886, 0.00438665, 0.93089201, 0.32227454, 0.94433412, - 0.52882808, 0.35897344, 0.98082307, 0.63672559, 0.40675592, 0.73029611, - 0.79718033, 0.83289500, 0.35912491, 0.73523154, 0.26524352, 0.90985766, - 0.11540680, 0.92642634, 0.50243902, 0.91701275, 0.99210843, 0.17884858, - 0.84208972, 0.35576873, 0.94030752, 0.75677768, 0.32701883, 0.05703666, - 0.10718289, 0.54239482, 0.46669026, 0.81808486, 0.94145871, 0.80625755, - 0.49432001, 0.75453334, 0.43996266, 0.16898902, 0.16110126, 0.33294432, - 0.02137160, 0.19508086, 0.95404022, 0.58019988, 0.04106973, 0.46852629, - 0.34709349, 0.46393951, 0.27337984, 0.38641706, 0.07065264, 0.24067118, - 0.00788229, 0.14663229, 0.39338238, 0.15231094, 0.28686524, 0.14091922, - 0.48755575, 0.08600964, 0.38103706, 0.85224310, 0.23203833, 0.20498248, - 0.81848563, 0.31811682, 0.38677373, 0.50907972, 0.18362144, 0.69818721, - 0.24296072, 0.33306600, 0.85275074, 0.67875562, 0.36115710, 0.86867365, - 0.77337390, 0.57887241, 0.30549375, 0.47475260, 0.09579983, 0.57751925, - 0.76935336, 0.28016905, 0.77684416, 0.46015672, 0.77270791, 0.98904684, - 0.37299274, 0.45216906, 0.75338760, 0.86066275, 0.87934727, 0.84870953, - 0.13656602, 0.19058923, 0.23393366, 0.31299008, 0.06451366, 0.78993980, - 0.65567868, 0.46164525, 0.32113817, 0.61085142, 0.30780485, 0.01553258, - 0.46506638, 0.85429998, 0.44394019, 0.71193856, 0.86622522, 0.31727764, - 0.59739528, 0.95999038, 0.75236607, 0.97717895, 0.98845944, 0.00110685, - 0.74519566, 0.59948487, 0.95196711, 0.00215843, 0.38529653, 0.60866461, - 0.20766827, 0.96787199, 0.19853626, 0.34042337, 0.41711054, 0.78376394, - 0.28499823, 0.74542367, 0.97893854, 0.77999021, 0.75856188, 0.43043826, - 0.02399550, 0.39406678, 0.62818864, 0.78321739, 0.27475616, 0.27616528, - 0.59092547, 0.65346187, 0.85093651, 0.25997639, 0.35494196, 0.36609947, - 0.33324756, 0.26593693, 0.38501432, 0.57710566, 0.86694471, 0.42853425, - 0.67233089, 0.88955319, 0.80095512, 0.57766286, 0.28801410, 0.11220959, - 0.14994893, 0.95346215, 0.05115330, 0.33923894, 0.94689118, 0.50818931, - 0.79636759, 0.57244585, 0.56162556, 0.92657841, 0.08503698, 0.70268735, - 0.11428299, 0.17253175, 0.09354997, 0.26086290, 0.85628473, 0.80523556, - 0.10953611, 0.27818728, 0.55005375, 0.87786380, 0.05452842, 0.83779496, - 0.50301807, 0.07853461, 0.66340314, 0.82328198, 0.63196771, 0.05333895, - 0.02496850, 0.98335094, 0.20832115, 0.98792252, 0.71504034, 0.63470091, - 0.29950013, 0.45574646, 0.52719557, 0.15773614, 0.24066173, 0.97188136, - 0.07420548, 0.33221847, 0.18687574, 0.34072899, 0.24511444, 0.58394877, - 0.17897235, 0.73976862, 0.56236628, 0.33084542, 0.55535200, 0.04026661, - 0.99673609, 0.92993273, 0.63380399, 0.39647042, 0.32466895, 0.72892237, - 0.62670658, 0.91935778, 0.91565539, 0.55105699, 0.82288511, 0.83738029, - 0.13230596, 0.61409905, 0.33102089, 0.99830501, 0.73176906, 0.56883094, - 0.62208830, 0.71231805, 0.64613903, 0.59484706, 0.79923468, 0.62173911, - 0.86671859, 0.50805419, 0.85831456, 0.61735511, 0.39075307, 0.70206796, - 0.15457926, 0.43149269, 0.04098807, 0.69213851, 0.16256298, 0.81511340, - 0.62268765, 0.67352228, 0.28687992, 0.55765997, 0.22728199, 0.10512828, - 0.65946148, 0.70695750, 0.36741255, 0.03321968, 0.32985075, 0.38065681, - 0.07875382, 0.96898541, 0.13379212, 0.48460052, 0.64668432, 0.30857389, - 0.50388801, 0.43682097, 0.41085104, 0.85599317, 0.08505214, 0.87399606, - 0.32124118, 0.19004074, 0.83550096, 0.04106973, 0.92198789, 0.75867300, - 0.25784983, 0.65099019, 0.03336388, 0.43112107, 0.01940191, 0.98235807, - 0.10095235, 0.17339695, 0.81894218, 0.65605885, 0.44965030, 0.81102496, - 0.35040995, 0.21115750, 0.18789907, 0.91891778, 0.34097893, 0.36455768, - 0.37642212, 0.53795720, 0.21970575, 0.78290199, 0.13097548, 0.10732775, - 0.14852556, 0.67926283, 0.34850206, 0.05789852, 0.14360772, 0.44732418, - 0.71207344, 0.41249598, 0.38369556, 0.50358289, 0.97667044, 0.35192068, - 0.67292223, 0.18083142, 0.40289416, 0.53982370, 0.35673004, 0.82667366, - 0.67993037, 0.45448807, 0.81705398, 0.64754801, 0.05645953, 0.59119537, - 0.36112078, 0.50578894, 0.42578084, 0.10586197, 0.80944494, 0.04437628, - 0.23939812, 0.84029307, 0.15338027, 0.66655224, 0.75501772, 0.16363826, - 0.13954332, 0.03640591, 0.86677965, 0.28598036, 0.72432050, 0.65843566, - 0.67715250, 0.64161301, 0.78138684, 0.52122744, 0.81396682, 0.23525096, - 0.82170979, 0.80058409, 0.98835937, 0.01858467, 0.46001722, 0.77418237, - 0.34889885, 0.08227679, 0.92343489, 0.73222398, 0.18593124, 0.76536543, - 0.61432600, 0.12130894, 0.55341636, 0.39637870, 0.00580329, 0.04224794, - 0.43154582, 0.87518874, 0.43856396, 0.66271913, 0.38368254, 0.40928257, - 0.50337870, 0.99416772, 0.95021297, 0.55086359, 0.01463005, 0.72636909, - 0.65745074, 0.00081230, 0.31475020, 0.58547970, 0.84259510, 0.10144061, - 0.08338780, 0.67066174, 0.04987421, 0.32221045, 0.57289357, 0.88458883, - 0.53546253, 0.53181886, 0.44564044, 0.37404670, 0.12974876, 0.26756614, - 0.16929104, 0.57378808, 0.99836600, 0.27063637, 0.33202868, 0.25075950, - 0.57452442, 0.10636918, 0.69773359, 0.89339249, 0.41417913, 0.97495142, - 0.33599781, 0.91406857, 0.96511787, 0.72399812, 0.56673759, 0.72997817, - 0.20698347, 0.37463587, 0.66511363, 0.67187469, 0.98050220, 0.08198928, - 0.51876553, 0.42077819, 0.20283958, 0.67506557, 0.53338976, 0.88423472, - 0.38724614, 0.70622265, 0.68994838, 0.95663912, 0.22730750, 0.12360727, - 0.49092116, 0.83835762, 0.95062040, 0.23869974, 0.64576479, 0.65754310, - 0.62815059, 0.96356518, 0.36700222, 0.20384731, 0.28158357, 0.34737755, - 0.53614047, 0.06755230, 0.05247672, 0.14866479, 0.28903988, 0.31944371, - 0.67790589, 0.57856175, 0.62768840, 0.77013400, 0.30464148, 0.87181519, - 0.93195563, 0.53223975, 0.21788411, 0.05765129, 0.99131898, 0.06017602, - 0.26455795, 0.17380561, 0.75140265, 0.56132862, 0.72154615, 0.29656591, - 0.05020634, 0.92560398, 0.65562838, 0.64075559, 0.45900735, 0.39006217, - 0.70693959, 0.74997666, 0.02610552, 0.58936327, 0.91733292, 0.47355250, - 0.40897745, 0.50457804, 0.13804328, 0.17749904, 0.06127244, 0.94163021, - 0.55958132, 0.01817415, 0.88280129, 0.39097941, 0.58603669, 0.31529723, - 0.93874684, 0.55246857, 0.15425351, 0.18270817, 0.48239620, 0.07952807, - 0.40710687, 0.83808682, 0.87551600, 0.59674891, 0.55392539, 0.03504402, - 0.97194628, 0.18143794, 0.67766843, 0.10131283, 0.82108480, 0.94600016, - 0.01859663, 0.27462766, 0.94480291, 0.30233267, 0.92181684, 0.69331861, - 0.05823558, 0.88596393, 0.67830721, 0.20233253, 0.63420421, 0.93615781, - 0.11603123, 0.39861027, 0.67326915, 0.40397403, 0.16285015, 0.27441518, - 0.03520218, 0.94431420, 0.06258755, 0.77586065, 0.63977594, 0.75688511, - 0.28935800, 0.31968745, 0.26006937, 0.50700603, 0.98583888, 0.18808550, - 0.84299415, 0.07812566, 0.41436091, 0.02345137, 0.71058959, 0.99727363, - 0.12446444, 0.78746386, 0.36106969, 0.11155264, 0.42898452, 0.33740698, - 0.21802004, 0.57290743, 0.16798465, 0.65356366, 0.71504037, 0.77013605, - 0.40427630, 0.06165743, 0.77514496, 0.40741624, 0.61728445, 0.01127253, - 0.60112461, 0.10910375, 0.57821900, 0.96309956, 0.86472382, 0.38844948, - 0.53658283, 0.21352301, 0.82690294, 0.47602811, 0.71681049, 0.44390239, - 0.67476991, 0.80054953, 0.70823144, 0.96394017, 0.78723809, 0.05443785, - 0.87865947, 0.95444005, 0.68347773, 0.97883715, 0.90192083, 0.07582916, - 0.83628749, 0.75288972, 0.26586664, 0.05695956, 0.19525090, 0.17120912, - 0.82631786, 0.91645865, 0.33871985, 0.40803659, 0.15353169, 0.28942814, - 0.41416447, 0.87885413, 0.34385331, 0.56938583, 0.01625497, 0.56807972, - 0.85604628, 0.75341863, 0.55232899, 0.75836933, 0.69832404, 0.66260711, - 0.44490425, 0.70209484, 0.31040385, 0.65555568, 0.28800599, 0.76348664, - 0.83016930, 0.75660831, 0.20379669, 0.54107301, 0.13507721, 0.86895963, - 0.61195760, 0.58039434, 0.33572608, 0.16392659, 0.22823332, 0.39811898, - 0.95758536, 0.94928874, 0.73039579, 0.56786612, 0.26446383, 0.08062597, - 0.37762072, 0.73810111, 0.39469104, 0.93939481, 0.22902952, 0.35188639, - 0.74214395, 0.95012693, 0.09319190, 0.97114877, 0.31624107, 0.72727100, - 0.14509092, 0.94418697, 0.35804516, 0.05112000, 0.29093107, 0.42974436, - 0.92698174, 0.68388500, 0.94285437, 0.75508718, 0.36084068, 0.40383196, - 0.31187609, 0.75575897, 0.76025188, 0.80065454, 0.64304524, 0.15325564, - 0.97040890, 0.22270266, 0.83498072, 0.47163203, 0.89327637, 0.57073643, - 0.51766784, 0.28848874, 0.53180422, 0.30856932, 0.14697621, 0.33899699, - 0.56585343, 0.41396431, 0.82916494, 0.73832739, 0.81281029, 0.96974427, - 0.05917982, 0.21840658, 0.91254925, 0.53218439, 0.06645393, 0.16531006, - 0.47794297, 0.28891123, 0.06295557, 0.53541282, 0.70308362, 0.69990253, - 0.01026776, 0.57888682, 0.65946586, 0.81155424, 0.02354035, 0.13110837, - 0.82641089, 0.64456998, 0.74636133, 0.98449562, 0.58992637, 0.02167233, - 0.82620632, 0.22933757, 0.75107065, 0.35149935, 0.43752698, 0.17202851, - 0.46091576, 0.29182233, 0.90923845, 0.48813696, 0.34770286, 0.84031617, - 0.26399543, 0.14689222, 0.02920177, 0.34707090, 0.54496359, 0.61956665, - 0.26919258, 0.70301586, 0.44668159, 0.23909403, 0.00407794, 0.66622638, - 0.06307568, 0.85994542, 0.54972977, 0.96810940, 0.67886697, 0.23926025, - 0.91871877, 0.65875971, 0.94975913, 0.65160416, 0.20613966, 0.66100585, - 0.60929384, 0.00522157, 0.01352552, 0.90211290, 0.90162043, 0.25762367, - 0.40315991, 0.76769302, 0.24327867, 0.85877565, 0.64089239, 0.20163279, - 0.46519584, 0.67262439, 0.17891084, 0.17694250, 0.25833530, 0.74946432, - 0.83851853, 0.74415649, 0.85081077, 0.11552736, 0.81762537, 0.57567877, - 0.26693121, 0.60067890, 0.01233623, 0.15631428, 0.87382427, 0.59430703, - 0.70377557, 0.20281627, 0.00138292, 0.47313408, 0.12032132, 0.84990472, - 0.74679998, 0.51810834, 0.89505329, 0.82878266, 0.24452955, 0.18172719, - 0.35489688, 0.63923570, 0.62290247, 0.09870198, 0.71123217, 0.32111017, - 0.08470608, 0.78402612, 0.08433680, 0.24911436, 0.37382601, 0.33110304, - 0.35175922, 0.68089127, 0.73467109, 0.92210611, 0.73252090, 0.99720778, - 0.07482037, 0.67228038, 0.92699402, 0.49055883, 0.72152609, 0.34418167, - 0.35158891, 0.63078540, 0.06488613, 0.14723501, 0.31633401, 0.30036503, - 0.72829681, 0.05845699, 0.87348552, 0.02539326, 0.25196998, 0.50332778, - 0.95103492, 0.59909741, 0.74772927, 0.47884426, 0.41408444, 0.21179913, - 0.54099793, 0.20658793, 0.96035529, 0.66398440, 0.49498658, 0.61067629, - 0.54752438, 0.98758676, 0.87684259, 0.52783592, 0.11672464, 0.76540677, - 0.21831488, 0.79426615, 0.16858546, 0.19106641, 0.41482809, 0.70180519, - 0.66895505, 0.91394555, 0.14538004, 0.77514315, 0.41737956, 0.56770803, - 0.75813436, 0.15904930, 0.46208621, 0.58046103, 0.80785918, 0.90202113, - 0.15332375, 0.04309245, 0.03447005, 0.12827243, 0.70827435, 0.59431141, - 0.39800417, 0.06703147, 0.63431108, 0.13141487, 0.82364136, 0.48582269, - 0.39520588, 0.73048386, 0.12827687, 0.59747988, 0.38218825, 0.15486439, - 0.59315178, 0.58943926, 0.59683647, 0.33332472, 0.28878887, 0.68436264, - 0.69709523, 0.92674187, 0.35072549, 0.97494610, 0.88582730, 0.69271971, - 0.32716111, 0.09010421, 0.92091697, 0.55996787, 0.68065513, 0.21743633, - 0.63786980, 0.78798902, 0.81155176, 0.81645532, 0.34404285, 0.10267134, - 0.19513239, 0.96259308, 0.26516098, 0.44542620, 0.57998340, 0.13321638, - 0.59982457, 0.77383876, 0.70222766, 0.11914411, 0.31961169, 0.91732398, - 0.72773347, 0.86942841, 0.49756895, 0.67454255, 0.06894273, 0.47004901, - 0.96594027, 0.38349396, 0.33943464, 0.15089975, 0.50759678, 0.98405657, - 0.84076861, 0.50786957, 0.54062452, 0.90070138, 0.03366471, 0.83141347, - 0.53084483, 0.81520570, 0.16843795, 0.12887521, 0.78475530, 0.88110020, - 0.45518844, 0.70968712, 0.16068297, 0.09958544, 0.48806982, 0.93078700, - 0.86715515, 0.26766775, 0.82364627, 0.25013288, 0.37606389, 0.10082662, - 0.70266187, 0.25066555, 0.49869588, 0.19280217, 0.52736817, 0.93269276, - 0.49281815, 0.41286925, 0.32089329, 0.75229927, 0.27956918, 0.05566025, - 0.53712326, 0.54550001, 0.24385260, 0.59024878, 0.12345935, 0.32655143, - 0.91833481, 0.63447686, 0.10699877, 0.70384469, 0.05920560, 0.87868694, - 0.44766627, 0.30347571, 0.11489109, 0.06202175, 0.82823492, 0.62989753, - 0.97335763, 0.86681851, 0.10777320, 0.49161164, 0.70908275, 0.23551309, - 0.60672474, 0.75987321, 0.09989137, 0.26367431, 0.30157532, 0.35772294, - 0.16248163, 0.50672150, 0.54602107, 0.05879029, 0.28254053, 0.42157995, - 0.59298071, 0.73934319, 0.56294588, 0.30366231, 0.03507543, 0.01771209, - 0.35686198, 0.67041865, 0.20584056, 0.84284487, 0.11988366, 0.67218765, - 0.56356658, 0.18432884, 0.43821122, 0.00674587, 0.41605092, 0.92134950, - 0.55271774, 0.59395440, 0.91755903, 0.73496854, 0.78939640, 0.19003227, - 0.78695134, 0.94964802, 0.86586766, 0.33665778, 0.44436511, 0.03459996, - 0.11533141, 0.83835059, 0.35730399, 0.81740206, 0.75063630, 0.02308454, - 0.01786634, 0.35111835, 0.58220386, 0.68711139, 0.66164455, 0.04724864, - 0.25770826, 0.06745381, 0.18418331, 0.89118043, 0.46511939, 0.77699125, - 0.63955381, 0.81103190, 0.44353764, 0.12067766, 0.72246568, 0.10708233, - 0.09087837, 0.97098644, 0.07340694, 0.06523760, 0.40972294, 0.69058746, - 0.51670676, 0.78765022, 0.66148728, 0.86556383, 0.18341522, 0.95580709, - 0.52653754, 0.65883644, 0.20187022, 0.79084036, 0.61920148, 0.85717020, - 0.40428252, 0.25285955, 0.66010414, 0.50310245, 0.52124889, 0.40172148, - 0.44534206, 0.85289989, 0.16955704, 0.46745253, 0.90727202, 0.15137941, - 0.71185881, 0.20673931, 0.17282442, 0.87290987, 0.73776642, 0.47536244, - 0.68307623, 0.04455988, 0.52080416, 0.06363940, 0.43628161, 0.58221055, - 0.70868204, 0.92696831, 0.00272832, 0.16417507, 0.64793179, 0.21615101, - 0.58386089, 0.61074942, 0.89856647, 0.45767367, 0.45340821, 0.49407289, - 0.50872672, 0.85351140, 0.23121707, 0.37324385, 0.67001538, 0.72990276, - 0.66069740, 0.62797860, 0.21628282, 0.44713658, 0.65092174, 0.98400849, - 0.64501364, 0.63555922, 0.02031386, 0.00929679, 0.59231579, 0.01022071, - 0.55945264, 0.77706162, 0.96908244, 0.56693662, 0.08787933, 0.59517717, - 0.83456816, 0.02077345, 0.18912276, 0.88808394, 0.93990270, 0.52442191, - 0.29777014, 0.57378119, 0.19158909, 0.12458072, 0.01375421, 0.33234240, - 0.86555754, 0.97093450, 0.61667268, 0.13185021, 0.62575445, 0.89180105, - 0.13459692, 0.02148055, 0.47554430, 0.76589505, 0.32828543, 0.16909003, - 0.98152539, 0.83938627, 0.12693798, 0.20211634, 0.07196474, 0.21561980, - 0.39267931, 0.28655092, 0.61832650, 0.65124178, 0.09804221, 0.87047680, - 0.87880966, 0.26238775, 0.72714869, 0.33809599, 0.81605767, 0.02901823, - 0.81351039, 0.81983860, 0.72401183, 0.10102948, 0.31495573, 0.32721513, - 0.51613604, 0.88865318, 0.01183643, 0.11133231, 0.27194933, 0.05973819, - 0.26419079, 0.45585551, 0.12209796, 0.68125743, 0.44589051, 0.92370725, - 0.51333894, 0.41090484, 0.81104227, 0.68760397, 0.15117297, 0.46166782, - 0.76525005, 0.93394293, 0.11496696, 0.88813101, 0.79008934, 0.27447053, - 0.50187519, 0.95396444, 0.08474622, 0.60813218, 0.53346299, 0.07717461, - 0.27880135, 0.31468634, 0.98254009, 0.27047696, 0.97228850, 0.04084992, - 0.50532963, 0.52143510, 0.96424822, 0.55834341, 0.53186723, 0.12928126, - 0.63297547, 0.27126104, 0.63969912, 0.33668652, 0.27709607, 0.35836771, - 0.25491284, 0.38491428, 0.26669517, 0.32550927, 0.36079575, 0.36511807, - 0.49128306, 0.87623104, 0.99660848, 0.43801121, 0.93286371, 0.58746863, - 0.37384043, 0.41534371, 0.79734039, 0.39497126, 0.67976503, 0.62179610, - 0.76864614, 0.48372765, 0.16114198, 0.93289364, 0.82314392, 0.86650463, - 0.04565417, 0.97683336, 0.07594888, 0.88371585, 0.75437729, 0.52480218, - 0.19795215, 0.62587778, 0.75264199, 0.92017403, 0.61766826, 0.06217101, - 0.17457026, 0.12225818, 0.53390645, 0.54202126, 0.01198264, 0.19100763, - 0.44512524, 0.54868460, 0.48303749, 0.20688546, 0.63374014, 0.74273044, - 0.06976231, 0.11215656, 0.71631144, 0.69785226, 0.34457628, 0.44907262, - 0.82426273, 0.45101351, 0.88011689, 0.64522688, 0.25851709, 0.35359473, - 0.36842835, 0.31087286, 0.06225657, 0.40061083, 0.25722948, 0.93695573, - 0.07910673, 0.55250415, 0.26306064, 0.65134612, 0.86274307, 0.76040336, - 0.64224656, 0.05527222, 0.63465387, 0.81360946, 0.01756218, 0.13881754, - 0.15914602, 0.49059240, 0.37613461, 0.39581518, 0.43134116, 0.19574371, - 0.90546394, 0.32237752, 0.92070005, 0.57794676, 0.03063173, 0.80312781, - 0.44847771, 0.96873535, 0.09604634, 0.54405053, 0.72452231, 0.65307756, - 0.91683153, 0.79547570, 0.95218087, 0.44093434, 0.05370701, 0.48341108, - 0.50808599, 0.07961491, 0.36456962, 0.40821134, 0.49771478, 0.21312319, - 0.70771900, 0.01802712, 0.41240443, 0.72332481, 0.74566548, 0.07783475, - 0.30667350, 0.88765206, 0.58271763, 0.25596251, 0.96705511, 0.25861868, - 0.50020562, 0.67737653, 0.74914168, 0.05167603, 0.69451201, 0.05067172, - 0.86524277, 0.05057197, 0.61743668, 0.28345598, 0.36116131, 0.17466187, - 0.03632483, 0.04397176, 0.23042168, 0.36248230, 0.12146970, 0.50634381, - 0.68427804, 0.91287052, 0.90732209, 0.29937753, 0.58556695, 0.88216188, - 0.11688003, 0.24658516, 0.22663618, 0.19731068, 0.16830221, 0.88593867, - 0.70200537, 0.33457263, 0.86041968, 0.55521554, 0.01403669, 0.39669796, - 0.74867152, 0.35181729, 0.12678652, 0.25407686, 0.14244217, 0.86932867, - 0.52456770, 0.79324301, 0.37434139, 0.88719778, 0.46589506, 0.82649027, - 0.98527239, 0.27073185, 0.28186258, 0.58663122, 0.53793964, 0.11772120, - 0.43517569, 0.07928799, 0.29563707, 0.57265089, 0.01973573, 0.65034563, - 0.19073691, 0.89248855, 0.28166232, 0.59772870, 0.64440408, 0.97792355, - 0.12489173, 0.30452369, 0.49994647, 0.13477852, 0.42507341, 0.37532567, - 0.78926164, 0.77184657, 0.01459860, 0.44082754, 0.28556373, 0.80897727, - 0.07255588, 0.36900776, 0.90460144, 0.01381822, 0.21266146, 0.43734677, - 0.42766950, 0.32381884, 0.11933000, 0.82247758, 0.73775231, 0.47414474, - 0.14972808, 0.43534136, 0.61812880, 0.23988686, 0.86370833, 0.44568770, - 0.07184413, 0.12585575, 0.13290306, 0.71467184, 0.31987516, 0.97842943, - 0.93703174, 0.09402039, 0.94296405, 0.34609999, 0.50484272, 0.10105256, - 0.69247764, 0.33889079, 0.37738624, 0.67766999, 0.03290055, 0.60766965, - 0.74394835, 0.56580320, 0.54444604, 0.02416489, 0.96053432, 0.64989130, - 0.53865297, 0.84842780, 0.14497033, 0.43756262, 0.08073748, 0.04542752, - 0.74231497, 0.57219133, 0.60601763, 0.89489031, 0.76328121, 0.93546812, - 0.32818512, 0.83250531, 0.26383686, 0.93520592, 0.10694382, 0.28593984, - 0.50627131, 0.05459252, 0.64207800, 0.11109143, 0.02410949, 0.62591678, - 0.13898630, 0.59362647, 0.61176530, 0.18697022, 0.62387837, 0.12471922, - 0.63432733, 0.17502842, 0.13006353, 0.87023567, 0.26605380, 0.24533266, - 0.11307102, 0.57725687, 0.53384239, 0.19273866, 0.84286351, 0.16428121, - 0.78307944, 0.47171256, 0.65539683, 0.75078757, 0.26993092, 0.20680054, - 0.89062224, 0.14280701, 0.59151188, 0.55936556, 0.99523038, 0.61409118, - 0.12178031, 0.45368713, 0.59963240, 0.73335856, 0.25789168, 0.19399655, - 0.69161961, 0.37854504, 0.28369241, 0.04183542, 0.18372278, 0.35621091, - 0.91171838, 0.62906761, 0.91551318, 0.88111322, 0.91332010, 0.30535886, - 0.23208163, 0.75591723, 0.73230334, 0.45000298, 0.84995876, 0.37819480, - 0.09782391, 0.09728269, 0.45418018, 0.66096794, 0.60740309, 0.25339599, - 0.88556381, 0.81945357, 0.19415600, 0.23327006, 0.47418517, 0.50132744, - 0.68009447, 0.56986153, 0.98026763, 0.43735312, 0.05924920, 0.28997095, - 0.22213214, 0.19584179, 0.26304806, 0.99296719, 0.13548181, 0.53597256, - 0.40170888, 0.48654007, 0.13026393, 0.58498728, 0.51830166, 0.14171475, - 0.92991362, 0.69178004, 0.85867589, 0.43353399, 0.91032197, 0.42915646, - 0.37583094, 0.78033317, 0.44789492, 0.18939256, 0.42984912, 0.12031170, - 0.16121025, 0.01865476, 0.39565974, 0.14668304, 0.74079686, 0.40695221, - 0.58522660, 0.63022693, 0.43633514, 0.34880725, 0.66594773, 0.29936163, - 0.16811860, 0.03860058, 0.70394403, 0.72231471, 0.68605676, 0.50205188, - 0.37343846, 0.08363673, 0.57909586, 0.49563442, 0.49866058, 0.59883851, - 0.18307059, 0.34862191, 0.38284623, 0.13762649, 0.31618320, 0.56327763, - 0.73326133, 0.88997858, 0.13717060, 0.92839688, 0.01820541, 0.21268153, - 0.44852148, 0.40660574, 0.85121682, 0.13010594, 0.21072828, 0.90872810, - 0.25190127, 0.96571483, 0.47651001, 0.31701091, 0.09775980, 0.77450085, - 0.05735449, 0.32652099, 0.67604650, 0.29626359, 0.98786046, 0.16969690, - 0.34111008, 0.11103827, 0.82359404, 0.41758347, 0.95606322, 0.24735670, - 0.39101760, 0.42304520, 0.07629842, 0.25827725, 0.27713570, 0.71222586, - 0.94777521, 0.46680045, 0.90358975, 0.63229476, 0.02971087, 0.80955489, - 0.85490303, 0.73332454, 0.29701079, 0.23428380, 0.43143945, 0.87282965, - 0.37329514, 0.58209996, 0.46052435, 0.57153801, 0.13466181, 0.09974419, - 0.83834289, 0.17247041, 0.80444093, 0.36084503, 0.28306307, 0.33209075, - 0.43033955, 0.24972348, 0.37663435, 0.03751623, 0.65086670, 0.48853011, - 0.25752207, 0.94851652, 0.67421394, 0.82098282, 0.72600968, 0.42388709, - 0.71941515, 0.64685489, 0.04251911, 0.14097982, 0.81931728, 0.51443430, - 0.17067073, 0.05256610, 0.82348062, 0.19785087, 0.55166256, 0.36027010, - 0.83553959, 0.02748641, 0.93425043, 0.79764408, 0.79586139, 0.01030830, - 0.14254004, 0.40329066, 0.80173747, 0.01918919, 0.58746605, 0.98562010, - 0.32215059, 0.78300744, 0.54086497, 0.01183996, 0.19382398, 0.55010409, - 0.74520593, 0.83621285, 0.11434937, 0.49466741, 0.45720908, 0.17371025, - 0.91440231, 0.19761280, 0.97776279, 0.94778974, 0.78520668, 0.75993995, - 0.63988166, 0.03360466, 0.65345133, 0.33741810, 0.97518463, 0.68462187, - 0.28085992, 0.55143958, 0.79989913, 0.29110049, 0.36976620, 0.40327895, - 0.52422195, 0.55852474, 0.88267141, 0.95971580, 0.44262137, 0.71386312, - 0.03930301, 0.28956130, 0.78410790, 0.21404605, 0.59872571, 0.40623665, - 0.14551409, 0.92654144, 0.32927120, 0.44497467, 0.68041561, 0.42421355, - 0.23127308, 0.99886755, 0.23262914, 0.94925216, 0.52576316, 0.62090428, - 0.97549509, 0.58295736, 0.52847499, 0.52499294, 0.11377772, 0.05428723, - 0.86790567, 0.68403694, 0.71949956, 0.54067301, 0.60712753, 0.91434256, - 0.78000421, 0.38565474, 0.42296304, 0.38051836, 0.33776915, 0.27099260, - 0.26760921, 0.87174820, 0.07685093, 0.44329070, 0.77451723, 0.44766852, - 0.56834743, 0.77968343, 0.91311962, 0.77947077, 0.41303724, 0.97225516, - 0.08979550, 0.00932069, 0.82766448, 0.87069481, 0.53234612, 0.89940173, - 0.50455909, 0.01695322, 0.58241143, 0.76213893, 0.35292745, 0.86972532, - 0.60801566, 0.36735676, 0.60284786, 0.90139751, 0.99323060, 0.06280776, - 0.64651800, 0.17228444, 0.87335941, 0.46484393, 0.41051032, 0.67875582, - 0.47476726, 0.65977698, 0.33070318, 0.53107402, 0.55433807, 0.17810236, - 0.54930412, 0.15690358, 0.18346389, 0.47828350, 0.51022477, 0.32661650, - 0.10754613, 0.10987161, 0.01238629, 0.95142795, 0.95940083, 0.56366110, - 0.13583935, 0.67540852, 0.49203346, 0.39626615, 0.50650484, 0.14350648, - 0.70442828, 0.01548833, 0.66363946, 0.12148005, 0.69785327, 0.39929313, - 0.96968647, 0.02731004, 0.45740848, 0.50391355, 0.51299419, 0.83009065, - 0.59496543, 0.59496592, 0.83610618, 0.12005625, 0.70920114, 0.74953527, - 0.44020479, 0.65639624, 0.16319464, 0.25762984, 0.81220737, 0.07756619, - 0.20022046, 0.16950506, 0.96547942, 0.19727100, 0.81726336, 0.25303056, - 0.02677381, 0.40915369, 0.29983745, 0.74505943, 0.41758834, 0.95107107, - 0.85048932, 0.49581990, 0.42253783, 0.91790322, 0.07959654, 0.04947580, - 0.99356913, 0.00797764, 0.98655835, 0.40160741, 0.16804166, 0.60586867, - 0.36116068, 0.29048722, 0.47010892, 0.27265925, 0.65085293, 0.96933078, - 0.68585172, 0.37845229, 0.34509481, 0.16747920, 0.84873888, 0.08134411, - 0.85113768, 0.96663149, 0.61685880, 0.25819913, 0.90804227, 0.32133128, - 0.26873363, 0.21699529, 0.49197305, 0.05908679, 0.67027919, 0.97546938, - 0.29354478, 0.16355479, 0.97961621, 0.74913276, 0.53804889, 0.62356279, - 0.84510548, 0.74192709, 0.71570122, 0.91141771, 0.02942715, 0.82894818, - 0.68869987, 0.41661766, 0.95113782, 0.73251084, 0.34505623, 0.72703619, - 0.59542382, 0.55570975, 0.12918494, 0.23079788, 0.70635969, 0.30252786, - 0.25255492, 0.96095156, 0.56063415, 0.29705582, 0.48388354, 0.16483471, - 0.94230491, 0.24874036, 0.15009486, 0.63450074, 0.02269665, 0.39362799, - 0.31238144, 0.15816203, 0.61019541, 0.48678874, 0.18515253, 0.46032698, - 0.91067635, 0.00759546, 0.36148925, 0.99584060, 0.34801208, 0.36202345, - 0.62215612, 0.89666678, 0.78105422, 0.98497424, 0.31455745, 0.07688824, - 0.34229393, 0.43679342, 0.97472210, 0.95858828, 0.16584540, 0.99324651, - 0.96402143, 0.85717907, 0.09834350, 0.00446583, 0.04056849, 0.76814498, - 0.55380931, 0.86863656, 0.57393094, 0.02748927, 0.51864933, 0.76892838, - 0.33486268, 0.11544810, 0.19997272, 0.06150180, 0.69628830, 0.16793087, - 0.38171904, 0.09805189, 0.61795730, 0.62394989, 0.47919435, 0.30364850, - 0.85834017, 0.83192898, 0.91214951, 0.01056750, 0.56991278, 0.46576893, - 0.42625775, 0.13798368, 0.34753798, 0.75606213, 0.63764341, 0.58803950, - 0.17690054, 0.95667807, 0.43798673, 0.74876615, 0.78670237, 0.16649930, - 0.13981391, 0.81709112, 0.06352629, 0.77541428, 0.85213562, 0.60982110, - 0.31399744, 0.95843141, 0.62300280, 0.96176742, 0.83030776, 0.99825887, - 0.26331429, 0.82843437, 0.96064077, 0.86180896, 0.44896924, 0.44189585, - 0.43698718, 0.61737576, 0.98796578, 0.70424403, 0.25384554, 0.57486171, - 0.91519003, 0.85053158, 0.81005336, 0.72018429, 0.04842224, 0.94529842, - 0.16331144, 0.24629028, 0.59064443, 0.18205138, 0.14679423, 0.04708226, - 0.41032253, 0.18755142, 0.25665135, 0.61768362, 0.46533965, 0.56218639, - 0.83839556, 0.21612406, 0.97818276, 0.94941727, 0.33246099, 0.59963034, - 0.01846613, 0.77226043, 0.54104703, 0.50579748, 0.03066894, 0.91960885, - 0.85930916, 0.04815602, 0.06495406, 0.64212750, 0.68671226, 0.72083225, - 0.36012848, 0.52433477, 0.26474655, 0.03253265, 0.90133587, 0.94977825, - 0.64416448, 0.33101214, 0.37019532, 0.04903235, 0.49050120, 0.94117034, - 0.83920204, 0.81031207, 0.81608501, 0.76925160, 0.78079121, 0.20016505, - 0.25834919, 0.35921566, 0.68817023, 0.66294782, 0.18946707, 0.50273705, - 0.39160343, 0.88421442, 0.17571952, 0.76964486, 0.13017095, 0.87197292, - 0.24426241, 0.24972187, 0.56312704, 0.15384830, 0.15091202, 0.77905437, - 0.21298530, 0.64600698, 0.66915420, 0.50460329, 0.68861779, 0.72534858, - 0.58487818, 0.75001514, 0.38335367, 0.44968277, 0.03353412, 0.35519846, - 0.95080913, 0.20673317, 0.45983426, 0.43377117, 0.57860654, 0.33282387, - 0.42120774, 0.42867583, 0.08400530, 0.19260875, 0.75454705, 0.44022695, - 0.90065178, 0.63442488, 0.04404682, 0.68638750, 0.69516945, 0.67964213, - 0.58195254, 0.15658237, 0.30728546, 0.33221676, 0.02275314, 0.70288666, - 0.93688437, 0.83563966, 0.06683834, 0.03163784, 0.57317598, 0.90332618, - 0.84709326, 0.01070375, 0.97057062, 0.01557599, 0.25964063, 0.97970641, - 0.49650424, 0.93751463, 0.50940926, 0.06999394, 0.71010489, 0.05951450, - 0.31697736, 0.02623365, 0.95712525, 0.90170935, 0.94290571, 0.43095741, - 0.79398933, 0.33021382, 0.77035889, 0.21664025, 0.96824961, 0.39301477, - 0.19015828, 0.69697318, 0.59848021, 0.65231613, 0.95158850, 0.24056626, - 0.30602577, 0.93711994, 0.67748509, 0.96767628, 0.37459222, 0.30036449, - 0.54668466, 0.18561001, 0.65003720, 0.88999350, 0.07726681, 0.54374686, - 0.26865245, 0.44949814, 0.80972786, 0.32781876, 0.42261100, 0.28298487, - 0.50299699, 0.80896938, 0.06908590, 0.05052620, 0.18385489, 0.89060223, - 0.67596765, 0.51492647, 0.20078962, 0.37015859, 0.01137236, 0.92018180, - 0.94869365, 0.25940426, 0.44362413, 0.63955852, 0.47626578, 0.73183833, - 0.86862197, 0.94771009, 0.04185220, 0.22638184, 0.94609569, 0.62257461, - 0.32718704, 0.01009828, 0.85697497, 0.71054544, 0.38664571, 0.45750474, - 0.92082412, 0.34981261, 0.13989817, 0.44884271, 0.59303862, 0.46468432, - 0.89906563, 0.52390336, 0.13239346, 0.35659210, 0.85158740, 0.79168434, - 0.45483232, 0.42946229, 0.84406538, 0.51995600, 0.81147661, 0.21792629, - 0.18056740, 0.88402572, 0.21821975, 0.52179478, 0.36455743, 0.50252126, - 0.36958935, 0.29031253, 0.88322581, 0.23240925, 0.76756442, 0.21154312, - 0.81608690, 0.72888508, 0.24161096, 0.02803250, 0.50898931, 0.63200632, - 0.01235317, 0.96751872, 0.00596693, 0.94540599, 0.11916351, 0.41036952, - 0.79359038, 0.38815488, 0.88078982, 0.84273991, 0.55773169, 0.78450643, - 0.15532463, 0.11087690, 0.65612160, 0.10414411, 0.40686884, 0.10790753, - 0.40008297, 0.66457177, 0.94478144, 0.91580812, 0.74680723, 0.48498351, - 0.47919701, 0.48366682, 0.81407669, 0.45548292, 0.21341618, 0.55604577, - 0.63011741, 0.11568323, 0.77783864, 0.44507494, 0.56311602, 0.57192246, - 0.18744672, 0.55534316, 0.98179553, 0.71261114, 0.78971778, 0.27579468, - 0.37738804, 0.08725811, 0.55462815, 0.91741883, 0.62856459, 0.91778153, - 0.56361874, 0.03488687, 0.46770121, 0.92272177, 0.77865432, 0.47542075, - 0.86572382, 0.30846909, 0.04258394, 0.36567526, 0.80325130, 0.15976596, - 0.28278103, 0.50950362, 0.25054190, 0.92833120, 0.39451597, 0.12400924, - 0.60366358, 0.04789403, 0.67982555, 0.80173938, 0.98829883, 0.31126026, - 0.75104860, 0.13235162, 0.86772958, 0.33592419, 0.79307959, 0.64230727, - 0.10302327, 0.66547720, 0.16025926, 0.97646982, 0.78957686, 0.14156171, - 0.68630613, 0.42688665, 0.81641807, 0.61291622, 0.36055321, 0.65311069, - 0.59758071, 0.69231761, 0.30087287, 0.39240225, 0.94616151, 0.63773967, - 0.19105764, 0.04989065, 0.69214686, 0.92022614, 0.58420835, 0.72674799, - 0.67871745, 0.18786155, 0.06812167, 0.69273440, 0.95556472, 0.07347767, - 0.18614898, 0.75784219, 0.77198267, 0.19464658, 0.97802581, 0.93154697, - 0.76045086, 0.21999738, 0.89837883, 0.17834361, 0.09954517, 0.62849139, - 0.86276277, 0.75140560, 0.51287292, 0.76116437, 0.01704602, 0.29458402, - 0.46102228, 0.59178293, 0.10875866, 0.69440919, 0.53206421, 0.84827982, - 0.71977298, 0.10174997, 0.75189310, 0.36576877, 0.08766211, 0.95642505, - 0.69070985, 0.90005695, 0.36566039, 0.58207272, 0.91143468, 0.77416264, - 0.34383601, 0.58332389, 0.51810539, 0.25631948, 0.35051741, 0.99876187, - 0.54881699, 0.20384774, 0.07494492, 0.83289562, 0.40305926, 0.37125948, - 0.77943371, 0.11564495, 0.78510687, 0.63182952, 0.01003871, 0.19149036, - 0.12423076, 0.84901943, 0.75494166, 0.21095160, 0.64287662, 0.04762861, - 0.65379082, 0.23853011, 0.18266543, 0.05502760, 0.52139959, 0.47340771, - 0.26530790, 0.73159343, 0.13605749, 0.53738871, 0.31461027, 0.99842950, - 0.05486034, 0.86978128, 0.35312733, 0.95215166, 0.58268196, 0.99773876, - 0.12226111, 0.67412280, 0.94257148, 0.85524844, 0.47040517, 0.77880797, - 0.00647563, 0.70231593, 0.97591539, 0.07731658, 0.55738706, 0.36142466, - 0.28919391, 0.87273538, 0.20656569, 0.02015289, 0.21307323, 0.65920031, - 0.59975298, 0.34163260, 0.41817890, 0.31309590, 0.05980695, 0.93868608, - 0.18835050, 0.07922173, 0.05356133, 0.77895470, 0.19195842, 0.74476566, - 0.88971302, 0.93874078, 0.91263775, 0.72763697, 0.48205163, 0.46984523, - 0.68526601, 0.39192863, 0.30637009, 0.40925093, 0.86075896, 0.64141035, - 0.37606968, 0.08455302, 0.26039296, 0.77204206, 0.02505414, 0.55267207, - 0.02395482, 0.11720525, 0.95937772, 0.97533250, 0.15626459, 0.26596589, - 0.69178053, 0.96739722, 0.43683515, 0.53712165, 0.45434231, 0.95439768, - 0.98406868, 0.05258748, 0.67077224, 0.61008697, 0.21801553, 0.66599979, - 0.31620906, 0.79579965, 0.50703031, 0.30006193, 0.08184154, 0.63867687, - 0.62947374, 0.95352375, 0.35530642, 0.52888537, 0.75087233, 0.92082748, - 0.11493688, 0.71592649, 0.78173144, 0.18002254, 0.98390321, 0.68572155, - 0.31669596, 0.64384449, 0.79461531, 0.71372631, 0.39028739, 0.12178050, - 0.04807762, 0.54084435, 0.74946165, 0.63504544, 0.91895815, 0.39997999, - 0.42329829, 0.77983289, 0.83483829, 0.03892205, 0.95294535, 0.40577275, - 0.85545545, 0.61273127, 0.57132416, 0.11563644, 0.88087046, 0.54439432, - 0.93599067, 0.28901987, 0.24064290, 0.63030622, 0.21702494, 0.16766900, - 0.15166367, 0.76687850, 0.62507295, 0.33458056 -}; - -/* - * Value Noise - * - */ -static double valueTab[TABSIZE]; -static int valueTabInitialized = 0; -#define vlattice(ix, iy, iz) (valueTab[INDEX(ix, iy, iz)]) - - -void -valueTabInit(seed) -long seed; -{ - register double *table = valueTab; - register int i; - int rndtabi = RT_RAND_TABSIZE - 1; - - if (valueTabInitialized) return; - valueTabInitialized = 1; - - RT_RANDSEED(rndtabi, seed); - - for (i=0 ; i < TABSIZE ; i++ ) - *table++ = 1. - 2.* RT_RANDOM(rndtabi); -} - -double -noise_v(point) -point_t point; -{ - int i, j, k; - int ix, iy, iz; /* lower integer lattice point */ - double x, y, z; /* corrected point */ - double fx, fy, fz; /* distance above integer lattice point */ - double xknots[4], yknots[4], zknots[4]; - static int initialized = 0; - double pt[3]; - - if (!valueTabInitialized) valueTabInit(665L); - - FILTER_ARGS( point); /* sets x,y,z, ix,iy,iz, fx,fy,fz */ - - for (k = -1; k <= 2; k++) { - for (j = -1; j <= 2; j++) { - for (i = -1; i <= 2; i++) - xknots[i+1] = vlattice(ix+i,iy+j,iz+k); - yknots[j+1] = spline(fx, 4, xknots); - } - zknots[k+1] = spline(fy, 4, yknots); - } - return spline(fz, 4, zknots); -} - - /* - * Lattice Convolution Noise - * - * - */ -double -noise_vc(point) -point_t point; -{ - int ix, iy, iz; /* lower integer lattice point */ - double x, y, z; /* corrected point */ - double fx, fy, fz; /* distance above integer lattice point */ - int i, j, k; - double dx, dy, dz; - double sum = 0; - - - if (!valueTabInitialized) valueTabInit(665L); - - FILTER_ARGS( point); /* sets x,y,z, ix,iy,iz, fx,fy,fz */ - - for (k= -1 ; k <= 2 ; k++) { - dz = k - fz; - dz = dz*dz; - for (j= -1; j <= 2 ; j++) { - dy = j - fy; - dy = dy*dy; - for (i=-1 ; i <= 2 ; i++) { - dx = i - fx; - dx = dx*dx; - sum += vlattice(ix+i,iy+j,iz+k) * - catrom2(dx+dy+dz); - } - } - } - return sum; -} - - - -/* Gradient Noise - * - * - * - */ -static double gradientTab[TABSIZE*3]; -static int gradientTabInitialized; - -static void -gradientTabInit(seed) -long seed; -{ - double *table = gradientTab; - double r, z, theta; - int i; - int rndtabi = RT_RAND_TABSIZE - 1; - - if (gradientTabInitialized) return; - gradientTabInitialized = 1; - - RT_RANDSEED(rndtabi, seed); - - for (i=0 ; i < TABSIZE ; i++) { - - /* z is [-1:1] */ - z = 1. - 2. * RT_RANDOM(rndtabi); - - - /* r is radius of x,y circle */ - r = sqrt(1 - z*z); - - /* theta is agnle in (x,y) */ - theta = 2 * rt_pi * RT_RANDOM(rndtabi); - - *table++ = r * cos(theta); - *table++ = r * sin(theta); - *table++ = z; - } -} - - -#define glattice(ix, iy, iz, fx, fy, fz) (\ - gradientTab[INDEX(ix, iy, iz)*3] * fx + \ - gradientTab[INDEX(ix, iy, iz)*3+1] * fy + \ - gradientTab[INDEX(ix, iy, iz)*3+2] * fz ) - - -double -noise_g(point) -point_t point; -{ - int ix, iy, iz; /* lower integer lattice point */ - double x, y, z; /* corrected point */ - double fx, fy, fz; /* distance above integer lattice point */ - double fx1, fy1, fz1; - double wx, wy, wz; - double vx0, vx1, vy0, vy1, vz0, vz1; - static int initialized = 0; - - if (!gradientTabInitialized) gradientTabInit(665L); - - FILTER_ARGS( point); /* sets x,y,z, ix,iy,iz, fx,fy,fz */ - - fx1 = fx - 1; - wx = SMOOTHSTEP(fx); - - fy1 = fy - 1; - wy = SMOOTHSTEP(fy); - - fz1 = fz - 1; - wz = SMOOTHSTEP(fz); - - vx0 = glattice(ix,iy,iz,fx,fy,fz); - vx1 = glattice(ix+1,iy,iz,fx1,fy,fz); - vy0 = LERP(wx, vx0, vx1); - vx0 = glattice(ix,iy+1,iz,fx,fy1,fz); - vx1 = glattice(ix+1,iy+1,iz,fx1,fy1,fz); - vy1 = LERP(wx, vx0, vx1); - vz0 = LERP(wy, vy0, vy1); - - vx0 = glattice(ix,iy,iz+1,fx,fy,fz1); - vx1 = glattice(ix+1,iy,iz+1,fx1,fy,fz1); - vy0 = LERP(wx, vx0, vx1); - vx0 = glattice(ix,iy+1,iz+1,fx,fy1,fz1); - vx1 = glattice(ix+1,iy+1,iz+1,fx1,fy1,fz1); - vy1 = LERP(wx, vx0, vx1); - vz1 = LERP(wy, vy0, vy1); - - return LERP(wz, vz0, vz1); -} - -/* - * Gradient Value Noise - * - */ -double -noise_gv(pt) -point_t pt; -{ - return 0.3 * noise_v(pt) + 0.7 * noise_g(pt); -} - - -/* Sparse Convolution Noise - * - * - */ - -#define NEXT(h) (((h)+1) & TABMASK) -#define NIMPULSES 3 - -static double impulseTab[TABSIZE*4]; -static int impulseTabInitialized; - -static void -impulseTabInit(seed) -long seed; -{ - int i; - double *f = impulseTab; - int rndtabi = RT_RAND_TABSIZE - 1; - - if (impulseTabInitialized) return; - impulseTabInitialized = 1; - - RT_RANDSEED(rndtabi, seed); - - for (i = 0; i < TABSIZE; i++) { - *f++ = RT_RANDOM(rndtabi); - *f++ = RT_RANDOM(rndtabi); - *f++ = RT_RANDOM(rndtabi); - *f++ = 1. - 2.*RT_RANDOM(rndtabi); - } - -} - -double -noise_sc(point) -point_t point; -{ - static int initialized; - double *fp; - int i, j, k, h, n; - int ix, iy, iz; /* lower integer lattice point */ - double x, y, z; /* corrected point */ - double fx, fy, fz; /* distance above integer lattice point */ - double sum = 0; - double dx, dy, dz, distsq; - - /* Initialize the random impulse table if necessary. */ - if (!impulseTabInitialized) impulseTabInit(665L); - - FILTER_ARGS(point); /* sets x,y,z, ix,iy,iz, fx,fy,fz */ - - /* Perform the sparse convolution. */ - for (i = -2; i <= 2; i++) { - for (j = -2; j <= 2; j++) { - for (k = -2; k <= 2; k++) { - /* Compute voxel hash code. */ - h = INDEX(ix+i,iy+j,iz+k); - - for (n = NIMPULSES; n > 0; n--, h = NEXT(h)) { - /* Convolve filter and impulse. */ - fp = &impulseTab[h*4]; - dx = fx - (i + *fp++); - dy = fy - (j + *fp++); - dz = fz - (k + *fp++); - distsq = dx*dx + dy*dy + dz*dz; - sum += catrom2(distsq) * *fp; - } - } - } - } - - return sum / NIMPULSES; -} - -void -init_noise_funcs() -{ - static int done=0; - RES_ACQUIRE(&rt_g.res_model); - - if (done) { - RES_RELEASE(&rt_g.res_model); - return; - } - done = 1; - - valueTabInit(665L); - gradientTabInit(665L); - impulseTabInit(665L); - - - RES_RELEASE(&rt_g.res_model); -} - -/**************************************************************** - * * - * F. Kenton Musgrave's noise functions * - * * - ****************************************************************/ - -/* Big/little endian is one heck of a pain in the butt. BYTE_SWAP should - be defined for big-endian machines--most of them. PC's, DECStations(?), - and other random architectures must not have it defined */ -#ifndef __linux__ -#define BYTE_SWAP -#endif - -#define MINX -10000 -#define MINY MINX -#define MINZ MINX - -#define REALSCALE ( 2.0 / 65536.0 ) -#define NREALSCALE ( 2.0 / 4096.0 ) -#define Hash3d(a,b,c) \ - ht.hashTable[ ht.hashTable[ ht.hashTable[(a) & 0xfff] ^ ((b) & 0xfff) ] ^ ((c) & 0xfff) ] - -#define Hash(a,b,c) (xtab[(xtab[(xtab[(a) & 0xff] ^ (b)) & 0xff] ^ (c)) & 0xff] & 0xff) - -#define INCRSUM(m,s,x,y,z) ((s)*(RTable[m]*0.5 \ - + RTable[m+1]*(x) \ - + RTable[m+2]*(y) \ - + RTable[m+3]*(z))) - -#define MAXSIZE 267 - -double RTable[MAXSIZE]; - -#define MAGIC_STRHT1 1771561 -#define MAGIC_STRHT2 1651771 -#define MAGIC_TAB1 9823 -#define MAGIC_TAB2 784642 -#define CK_HT() { \ - RT_CKMAG(&ht.magic, MAGIC_STRHT1, "struct str_ht ht 1"); \ - RT_CKMAG(&ht.magic_end, MAGIC_STRHT2, "struct str_ht ht 2"); \ - RT_CKMAG(ht.hashTableMagic1, MAGIC_TAB1, "hashTable Magic 1"); \ - RT_CKMAG(ht.hashTableMagic2, MAGIC_TAB2, "hashTable Magic 2"); \ - if (ht.hashTable != (short *)&ht.hashTableMagic1[1] ) \ - rt_bomb("ht.hashTable changed rel ht.hashTableMagic1"); \ - if (ht.hashTableMagic2 != (long *)&ht.hashTable[4096] ) \ - rt_bomb("ht.hashTable changed rel ht.hashTableMagic2"); \ -} - -struct str_ht { - long magic; - char hashTableValid; - long *hashTableMagic1; - short *hashTable; - long *hashTableMagic2; - long magic_end; -}; - -static struct str_ht ht; - -static unsigned short xtab[256] = -{ - 0x0000, 0xc0c1, 0xc181, 0x0140, 0xc301, 0x03c0, 0x0280, 0xc241, - 0xc601, 0x06c0, 0x0780, 0xc741, 0x0500, 0xc5c1, 0xc481, 0x0440, - 0xcc01, 0x0cc0, 0x0d80, 0xcd41, 0x0f00, 0xcfc1, 0xce81, 0x0e40, - 0x0a00, 0xcac1, 0xcb81, 0x0b40, 0xc901, 0x09c0, 0x0880, 0xc841, - 0xd801, 0x18c0, 0x1980, 0xd941, 0x1b00, 0xdbc1, 0xda81, 0x1a40, - 0x1e00, 0xdec1, 0xdf81, 0x1f40, 0xdd01, 0x1dc0, 0x1c80, 0xdc41, - 0x1400, 0xd4c1, 0xd581, 0x1540, 0xd701, 0x17c0, 0x1680, 0xd641, - 0xd201, 0x12c0, 0x1380, 0xd341, 0x1100, 0xd1c1, 0xd081, 0x1040, - 0xf001, 0x30c0, 0x3180, 0xf141, 0x3300, 0xf3c1, 0xf281, 0x3240, - 0x3600, 0xf6c1, 0xf781, 0x3740, 0xf501, 0x35c0, 0x3480, 0xf441, - 0x3c00, 0xfcc1, 0xfd81, 0x3d40, 0xff01, 0x3fc0, 0x3e80, 0xfe41, - 0xfa01, 0x3ac0, 0x3b80, 0xfb41, 0x3900, 0xf9c1, 0xf881, 0x3840, - 0x2800, 0xe8c1, 0xe981, 0x2940, 0xeb01, 0x2bc0, 0x2a80, 0xea41, - 0xee01, 0x2ec0, 0x2f80, 0xef41, 0x2d00, 0xedc1, 0xec81, 0x2c40, - 0xe401, 0x24c0, 0x2580, 0xe541, 0x2700, 0xe7c1, 0xe681, 0x2640, - 0x2200, 0xe2c1, 0xe381, 0x2340, 0xe101, 0x21c0, 0x2080, 0xe041, - 0xa001, 0x60c0, 0x6180, 0xa141, 0x6300, 0xa3c1, 0xa281, 0x6240, - 0x6600, 0xa6c1, 0xa781, 0x6740, 0xa501, 0x65c0, 0x6480, 0xa441, - 0x6c00, 0xacc1, 0xad81, 0x6d40, 0xaf01, 0x6fc0, 0x6e80, 0xae41, - 0xaa01, 0x6ac0, 0x6b80, 0xab41, 0x6900, 0xa9c1, 0xa881, 0x6840, - 0x7800, 0xb8c1, 0xb981, 0x7940, 0xbb01, 0x7bc0, 0x7a80, 0xba41, - 0xbe01, 0x7ec0, 0x7f80, 0xbf41, 0x7d00, 0xbdc1, 0xbc81, 0x7c40, - 0xb401, 0x74c0, 0x7580, 0xb541, 0x7700, 0xb7c1, 0xb681, 0x7640, - 0x7200, 0xb2c1, 0xb381, 0x7340, 0xb101, 0x71c0, 0x7080, 0xb041, - 0x5000, 0x90c1, 0x9181, 0x5140, 0x9301, 0x53c0, 0x5280, 0x9241, - 0x9601, 0x56c0, 0x5780, 0x9741, 0x5500, 0x95c1, 0x9481, 0x5440, - 0x9c01, 0x5cc0, 0x5d80, 0x9d41, 0x5f00, 0x9fc1, 0x9e81, 0x5e40, - 0x5a00, 0x9ac1, 0x9b81, 0x5b40, 0x9901, 0x59c0, 0x5880, 0x9841, - 0x8801, 0x48c0, 0x4980, 0x8941, 0x4b00, 0x8bc1, 0x8a81, 0x4a40, - 0x4e00, 0x8ec1, 0x8f81, 0x4f40, 0x8d01, 0x4dc0, 0x4c80, 0x8c41, - 0x4400, 0x84c1, 0x8581, 0x4540, 0x8701, 0x47c0, 0x4680, 0x8641, - 0x8201, 0x42c0, 0x4380, 0x8341, 0x4100, 0x81c1, 0x8081, 0x4040 -}; - -/* - * Note that passing a double to Crc16 and interpreting it as - * an array of chars means that machines with different floating-point - * representation schemes will evaluate Noise(point) differently. - */ -static int -Crc16(buf, count) -register char *buf; -register size_t count; -{ - register unsigned int crc = 0; - - while (count--) -#ifdef BYTE_SWAP - crc = (crc >> 8) ^ xtab[ (unsigned char)(crc ^ *(buf+count))]; -#else - crc = (crc >> 8) ^ xtab[ (unsigned char) (crc ^ *buf++) ]; -#endif - - return crc; -} - -static int -Randomizer(v) -point_t v; -{ - point_t t; - - t[0] = v[0] * .12345; - t[1] = v[1] * .12345; - t[2] = v[2] * .12345; - - return Crc16((char *)t, sizeof(t) ); -} - - -void -noise_init() -{ - point_t rp; - int i, j, k, temp; - int rndtabi = RT_RAND_TABSIZE - 1; - - RES_ACQUIRE(&rt_g.res_model); - - if (ht.hashTableValid) { - RES_RELEASE(&rt_g.res_model); - return; - } - - RT_RANDSEED(rndtabi, (RT_RAND_TABSIZE-1) ); - ht.hashTableMagic1 = (long *) rt_malloc( - 2*sizeof(long) + 4096*sizeof(short int), - "noise hashTable"); - ht.hashTable = (short *)&ht.hashTableMagic1[1]; - ht.hashTableMagic2 = (long *)&ht.hashTable[4096]; - - *ht.hashTableMagic1 = MAGIC_TAB1; - *ht.hashTableMagic2 = MAGIC_TAB2; - - ht.magic_end = MAGIC_STRHT2; - ht.magic = MAGIC_STRHT1; - - for (i = 0; i < 4096; i++) - ht.hashTable[i] = i; - - /* scramble the hash table */ - for (k = 0, i = 4095; i > 0; i--, k++) { - j = (int)(RT_RANDOM(rndtabi) * 4096.0); - - temp = ht.hashTable[i]; - ht.hashTable[i] = ht.hashTable[j]; - ht.hashTable[j] = temp; - } - - for (i = 0; i < MAXSIZE; i++) { - rp[0] = rp[1] = rp[2] = (double)i; - RTable[i] = Randomizer(rp) * REALSCALE - 1.0; - } - - ht.hashTableValid = 1; - - RES_RELEASE(&rt_g.res_model); - - CK_HT(); -} - - - -/* - * Robert Skinner's Perlin-style "Noise" function - */ -double -noise_perlin(point) -point_t point; -{ - register int jx, jy, jz; - int ix, iy, iz; /* lower integer lattice point */ - double x, y, z; /* corrected point */ - double fx, fy, fz; /* distance above integer lattice point */ - double sx, sy, sz, tx, ty, tz; - double sum; - short m; - point_t pt; - - if (!ht.hashTableValid) noise_init(); - else { - CK_HT(); - } - FILTER_ARGS( point); /* sets x,y,z, ix,iy,iz, fx,fy,fz */ - - jx = ix + 1; - jy = iy + 1; - jz = iz + 1; - - sx = SMOOTHSTEP(fx); - sy = SMOOTHSTEP(fy); - sz = SMOOTHSTEP(fz); - - /* the complement values of sx,sy,sz */ - tx = 1.0 - sx; - ty = 1.0 - sy; - tz = 1.0 - sz; - - /* - * interpolate! - */ - m = Hash3d( ix, iy, iz ) & 0xFF; - sum = INCRSUM(m,(tx*ty*tz),(x-ix),(y-iy),(z-iz)); - - m = Hash3d( jx, iy, iz ) & 0xFF; - sum += INCRSUM(m,(sx*ty*tz),(x-jx),(y-iy),(z-iz)); - - m = Hash3d( ix, jy, iz ) & 0xFF; - sum += INCRSUM(m,(tx*sy*tz),(x-ix),(y-jy),(z-iz)); - - m = Hash3d( jx, jy, iz ) & 0xFF; - sum += INCRSUM(m,(sx*sy*tz),(x-jx),(y-jy),(z-iz)); - - m = Hash3d( ix, iy, jz ) & 0xFF; - sum += INCRSUM(m,(tx*ty*sz),(x-ix),(y-iy),(z-jz)); - - m = Hash3d( jx, iy, jz ) & 0xFF; - sum += INCRSUM(m,(sx*ty*sz),(x-jx),(y-iy),(z-jz)); - - m = Hash3d( ix, jy, jz ) & 0xFF; - sum += INCRSUM(m,(tx*sy*sz),(x-ix),(y-jy),(z-jz)); - - m = Hash3d( jx, jy, jz ) & 0xFF; - sum += INCRSUM(m,(sx*sy*sz),(x-jx),(y-jy),(z-jz)); - - return sum; - -} - -/* - * Vector-valued "Noise" - */ -void -noise_vec(point, result) -point_t point; -point_t result; -{ - register int jx, jy, jz; - int ix, iy, iz; /* lower integer lattice point */ - double x, y, z; /* corrected point */ - double fx, fy, fz; /* distance above integer lattice point */ - double px, py, pz, s; - double sx, sy, sz, tx, ty, tz; - short m; - point_t pt; - - - if ( ! ht.hashTableValid ) noise_init(); - - - FILTER_ARGS( point); /* sets x,y,z, ix,iy,iz, fx,fy,fz */ - - jx = ix+1; jy = iy + 1; jz = iz + 1; - - sx = SMOOTHSTEP(x - ix); - sy = SMOOTHSTEP(y - iy); - sz = SMOOTHSTEP(z - iz); - - /* the complement values of sx,sy,sz */ - tx = 1.0 - sx; - ty = 1.0 - sy; - tz = 1.0 - sz; - - /* - * interpolate! - */ - m = Hash3d( ix, iy, iz ) & 0xFF; - px = x-ix; - py = y-iy; - pz = z-iz; - s = tx*ty*tz; - result[0] = INCRSUM(m,s,px,py,pz); - result[1] = INCRSUM(m+4,s,px,py,pz); - result[2] = INCRSUM(m+8,s,px,py,pz); - - m = Hash3d( jx, iy, iz ) & 0xFF; - px = x-jx; - s = sx*ty*tz; - result[0] += INCRSUM(m,s,px,py,pz); - result[1] += INCRSUM(m+4,s,px,py,pz); - result[2] += INCRSUM(m+8,s,px,py,pz); - - m = Hash3d( jx, jy, iz ) & 0xFF; - py = y-jy; - s = sx*sy*tz; - result[0] += INCRSUM(m,s,px,py,pz); - result[1] += INCRSUM(m+4,s,px,py,pz); - result[2] += INCRSUM(m+8,s,px,py,pz); - - m = Hash3d( ix, jy, iz ) & 0xFF; - px = x-ix; - s = tx*sy*tz; - result[0] += INCRSUM(m,s,px,py,pz); - result[1] += INCRSUM(m+4,s,px,py,pz); - result[2] += INCRSUM(m+8,s,px,py,pz); - - m = Hash3d( ix, jy, jz ) & 0xFF; - pz = z-jz; - s = tx*sy*sz; - result[0] += INCRSUM(m,s,px,py,pz); - result[1] += INCRSUM(m+4,s,px,py,pz); - result[2] += INCRSUM(m+8,s,px,py,pz); - - m = Hash3d( jx, jy, jz ) & 0xFF; - px = x-jx; - s = sx*sy*sz; - result[0] += INCRSUM(m,s,px,py,pz); - result[1] += INCRSUM(m+4,s,px,py,pz); - result[2] += INCRSUM(m+8,s,px,py,pz); - - m = Hash3d( jx, iy, jz ) & 0xFF; - py = y-iy; - s = sx*ty*sz; - result[0] += INCRSUM(m,s,px,py,pz); - result[1] += INCRSUM(m+4,s,px,py,pz); - result[2] += INCRSUM(m+8,s,px,py,pz); - - m = Hash3d( ix, iy, jz ) & 0xFF; - px = x-ix; - s = tx*ty*sz; - result[0] += INCRSUM(m,s,px,py,pz); - result[1] += INCRSUM(m+4,s,px,py,pz); - result[2] += INCRSUM(m+8,s,px,py,pz); -} -/************************************************************* - * - * Spectral Noise functions - * - ************************************************************* - * - * The Spectral Noise functions cache the values of the - * term: - * (-h_val) - * freq - * - * Which on some systems is rather expensive to compute. - * - *************************************************************/ -struct fbm_spec { - long magic; - double octaves; - double lacunarity; - double h_val; - double remainder; - double *spec_wgts; -}; -#define MAGIC_fbm_spec_wgt 0x837592 - -static struct fbm_spec *etbl = (struct fbm_spec *)NULL; -static int etbl_next = 0; -static int etbl_size = 0; - -#define PSCALE(_p, _s) _p[0] *= _s; _p[1] *= _s; _p[2] *= _s -#define PCOPY(_d, _s) _d[0] = _s[0]; _d[1] = _s[1]; _d[2] = _s[2] - - - -static struct fbm_spec * -build_spec_tbl(h_val, lacunarity, octaves) -double h_val, lacunarity, octaves; -{ - struct fbm_spec *ep; - double *spec_wgts; - double frequency; - int i; - - /* The right spectral weights table for these parameters has not been - * pre-computed. As a result, we compute the table now and save it - * with the knowledge that we'll likely want it again later. - */ - - /* allocate storage for new tables if needed */ - if (etbl_next >= etbl_size) { - if (etbl_size) { - etbl_size *= 2; - etbl = (struct fbm_spec *)rt_realloc((char *)etbl, - etbl_size*sizeof(struct fbm_spec), - "spectral weights table"); - } else - etbl = (struct fbm_spec *)rt_calloc(etbl_size = 10, - sizeof(struct fbm_spec), - "spectral weights table"); - - if (!etbl) abort(); - } - - /* set up the next available table */ - ep = &etbl[etbl_next]; - ep->magic = MAGIC_fbm_spec_wgt; ep->octaves = octaves; - ep->h_val = h_val; ep->lacunarity = lacunarity; - spec_wgts = ep->spec_wgts = - (double *)rt_malloc( ((int)(octaves+1)) * sizeof(double), - "spectral weights" ); - - /* precompute and store spectral weights table */ - for (frequency = 1.0, i=0 ; i < octaves ; i++) { - /* compute weight for each frequency */ - spec_wgts[i] = pow(frequency, -h_val); - frequency *= lacunarity; - } - - etbl_next++; /* saved for last in case we're running multi-threaded */ - return ep; -} - -/* The first order of business is to see if we have pre-computed - * the spectral weights table for these parameters in a previous - * invocation. If not, the we compute them and save them for - * possible future use - */ - -struct fbm_spec * -find_spec_wgt(h, l, o) -double h; -double l; -double o; -{ - struct fbm_spec *ep; - int i; - - - for (ep = etbl, i=0 ; i < etbl_next ; i++, ep++) { - if (ep->magic != MAGIC_fbm_spec_wgt) rt_bomb("find_spec_wgt"); - if (ep->lacunarity == l && ep->h_val == h && - ep->octaves >= o ) - return ep; - } - - /* we didn't find the table we wanted so we've got to semaphore on - * the list to wait our turn to add what we want to the table. - */ - - RES_ACQUIRE(&rt_g.res_model); - - /* We search the list one more time in case the last process to - * hold the semaphore just created the table we were about to add - */ - for (ep = etbl, i=0 ; i < etbl_next ; i++, ep++) { - if (ep->magic != MAGIC_fbm_spec_wgt) rt_bomb("find_spec_wgt"); - if (ep->lacunarity == l && ep->h_val == h && - ep->octaves >= o ) - break; - } - - if (i >= etbl_next) ep = build_spec_tbl(h, l, o); - - RES_RELEASE(&rt_g.res_model); - - return (ep); -} -/* - * Procedural fBm evaluated at "point"; returns value stored in "value". - * - * Parameters: - * ``h_val'' fractal increment parameter - * ``lacunarity'' gap between successive frequencies - * ``octaves'' number of frequencies in the fBm - * - * The function call pow() is relatively expensive. Therfore, this function - * pre-computes and saves the spectral weights in a table for re-use in - * successive invocations. - */ -double -noise_fbm(point, h_val, lacunarity, octaves) -point_t point; -double h_val; -double lacunarity; -double octaves; -{ - struct fbm_spec *ep; - double value, frequency, remainder, *spec_wgts; - point_t pt; - double Nosie3(); - int i, oct; - - /* The first order of business is to see if we have pre-computed - * the spectral weights table for these parameters in a previous - * invocation. If not, the we compute them and save them for - * possible future use - */ - - ep = find_spec_wgt(h_val, lacunarity, octaves); - - /* now we're ready to compute the fBm value */ - - value = 0.0; /* initialize vars to proper values */ - frequency = 1.0; - /* copy the point so we don't corrupt the caller's version */ - PCOPY(pt, point); - - spec_wgts = ep->spec_wgts; - - /* inner loop of spectral construction */ - oct=(int)octaves; /* save repeating double->int cast */ - for (i=0 ; i < oct ; i++) { - value += noise_perlin( pt ) * spec_wgts[i]; - PSCALE(pt, lacunarity); - } - - remainder = octaves - (int)octaves; - if ( remainder ) { - /* add in ``octaves'' remainder - * ``i'' and spatial freq. are preset in loop above - */ - value += remainder * noise_perlin( pt ) * spec_wgts[i]; - } - - return( value ); - -} /* noise_fbm() */ - - -/* - * Procedural turbulence evaluated at "point"; - * - * returns value stored in "value". - * - * Parameters: - * ``h_val'' fractal increment parameter - * ``lacunarity'' gap between successive frequencies - * ``octaves'' number of frequencies in the fBm - * - * The function call pow() is relatively expensive. Therfore, this function - * pre-computes and saves the spectral weights in a table for re-use in - * successive invocations. - */ -double -noise_turb( point, h_val, lacunarity, octaves) -point_t point; -double h_val; -double lacunarity; -double octaves; -{ - struct fbm_spec *ep; - double value, frequency, remainder, *spec_wgts; - point_t pt; - double Nosie3(); - int i, oct; - - - /* The first order of business is to see if we have pre-computed - * the spectral weights table for these parameters in a previous - * invocation. If not, the we compute them and save them for - * possible future use - */ - -#define CACHE_SPECTRAL_WGTS 1 -#ifdef CACHE_SPECTRAL_WGTS - - ep = find_spec_wgt(h_val, lacunarity, octaves); - - /* now we're ready to compute the fBm value */ - - value = 0.0; /* initialize vars to proper values */ - frequency = 1.0; - /* copy the point so we don't corrupt - * the caller's copy of the variable - */ - PCOPY(pt, point); - spec_wgts = ep->spec_wgts; - - /* inner loop of spectral construction */ - oct=(int)octaves; /* save repeating double->int cast */ - for (i=0 ; i < oct ; i++) { - value += fabs(noise_perlin( pt )) * spec_wgts[i]; - PSCALE(pt, lacunarity); - } - - remainder = octaves - (int)octaves; - if ( remainder ) { - /* add in ``octaves'' remainder - * ``i'' and spatial freq. are preset in loop above - */ - value += remainder * noise_perlin( pt ) * spec_wgts[i]; - } -#else - PCOPY(pt, point); - - value = 0.0; /* initialize vars to proper values */ - frequency = 1.0; - - oct=(int)octaves; /* save repeating double->int cast */ - for (i=0 ; i < oct ; i++) { - value += fabs(noise_perlin( pt )) * pow(frequency, -h_val); - frequency *= lacunarity; - PSCALE(pt, lacunarity); - } - - remainder = octaves - (int)octaves; - if ( remainder ) { - /* add in ``octaves'' remainder - * ``i'' and spatial freq. are preset in loop above - */ - value += remainder * noise_perlin( pt ) * pow(frequency, -h_val); - } -#endif - return( value ); - -} /* noise_turb() */