From 693913ffdbb28014bef1b9fa82fcceb04ae1e52b Mon Sep 17 00:00:00 2001 From: Homer Reid Date: Tue, 10 Jul 2018 18:28:31 -0400 Subject: [PATCH 1/8] update point_in_prism to avoid double-counting of polygon edge intersections when plumb lines pass through vertices --- utils/geom.c | 12 ++++++++---- 1 file changed, 8 insertions(+), 4 deletions(-) diff --git a/utils/geom.c b/utils/geom.c index 80e9bf0..a2619f8 100644 --- a/utils/geom.c +++ b/utils/geom.c @@ -2149,10 +2149,14 @@ boolean intersect_line_with_segment(double px, double py, double dx, double dy, double t = (M00*RHSy-M10*RHSx)/DetM; if (s) *s = (M11*RHSx-M01*RHSy)/DetM; - if (t<0.0 || t>1.0) // intersection of lines does not lie between vertices - return 0; - - return 1; + // the plumb line intersects the segment if 0<=t<=1, with t==0,1 + // corresponding to the endpoints; for our purposes we count + // the intersection if the plumb line runs through the t==0 vertex, but + // NOT the t==1 vertex, to avoid double-counting for complete polygons. +#define DELTAT 1.0e-6 +#define TMIN (0.0) +#define TMAX (1.0-DELTAT) + return ( ( tTMAX ) ? 0 : 1 ); } // like the previous routine, but only count intersections if s>=0 From 893208c4b5eecb1ef2f41202de67ef8201f7169d Mon Sep 17 00:00:00 2001 From: Homer Reid Date: Thu, 12 Jul 2018 08:04:47 -0400 Subject: [PATCH 2/8] revised normal_to_object implementation for prisms to account for in-plane distance from point to prism face when determining closest face --- utils/geom.c | 138 +++++++++++++++++++++++++++++++-------------- utils/test-prism.c | 20 ++++++- 2 files changed, 116 insertions(+), 42 deletions(-) diff --git a/utils/geom.c b/utils/geom.c index a2619f8..8e8f168 100644 --- a/utils/geom.c +++ b/utils/geom.c @@ -2301,27 +2301,85 @@ double intersect_line_segment_with_prism(prism *prsm, vector3 p, vector3 d, doub } /***************************************************************/ -/* compute the minimum distance from p to the plane */ -/* containing o (origin) and spanned by basis vectors v1,v2 */ -/* algorithm: solve the 3x3 system p-s*v3 = o + t*v1 + u*v2 */ -/* where v3 = v1 x v2 and s,t,u are unknowns. */ -/* return value is unit normal vector to plane and *s is such */ -/* that p-(*s)*v3 lies in the plane */ +/* compute the minimum distance from a 3D point p to the */ +/* line segment with endpoints v1,v2. */ +/* algorithm: let pLine = v1 + d*(v2-v1) be the point on the */ +/* line closest to p; d is defined by minimizing |p-pLine|^2. */ +/* --> |p-v1|^2 + d^2 |v2-v1|^2 - 2*d*dot(p-v1,v2-v1) = min */ +/* --> 2d |v2-v1|^2 - 2*dot(p-v1,v2-v1) = 0 */ +/* --> d = dot(p-v1,v2-v1) / |v2-v1|^2 */ /***************************************************************/ -vector3 normal_to_plane(vector3 o, vector3 v1, vector3 v2, vector3 p, double *s) +double min_distance_to_line_segment(vector3 p, vector3 v1, vector3 v2) +{ vector3 v2mv1 = vector3_minus(v2,v1); + vector3 pmv1 = vector3_minus(p,v1); + double d = vector3_dot(v2mv1,pmv1) / vector3_dot(v2mv1,v2mv1); + if (d<0.0) d=0.0; // if pProj lies outside the line segment, + if (d>1.0) d=1.0; // displace it to whichever vertex is closer + vector3 pLine = vector3_plus(v1, vector3_scale(d,v2mv1)); + return vector3_norm(vector3_minus(p,pLine)); +} + +/***************************************************************/ +/* compute the minimum distance from a 3D point p to the */ +/* planar quadrilateral with vertices {o, o+l1, o+l2, o+l1+l2} */ +/* (if quadrilateral==1) */ +/* or the planar triangle with vertices {o, o+l1, o+l2} */ +/* (if quadrilateral==0). */ +/* algorithm: */ +/* (A) solve a 3x3 system to compute the projection of p into */ +/* the plane of the triangle/quadrilateral (call it pPlane)*/ +/* pPlane = p-s*l3 = o + t*l1 + u*l2 */ +/* where l3 is the normal to the surface and s,t,u are */ +/* unknowns. */ +/* (B) if pPlane lies within the triangle/quadrilateral, the */ +/* minimum distance is simply s. */ +/* (C) otherwise compute the minimum distance d from pProj */ +/* to any edge of the triangle/quadrilateral and set */ +/* minimum_distance = sqrt(s^2 + d^2). */ +/* return value is normal to plane. */ +/***************************************************************/ +vector3 normal_to_plane(vector3 o, vector3 l1, vector3 l2, vector3 p, + int quadrilateral, double *min_distance) { - vector3 v3 = unit_vector3(vector3_cross(v1, v2)); - CHECK( (vector3_norm(v3)>1.0e-6), "degenerate plane in normal_to_plane" ); + vector3 l3 = unit_vector3(vector3_cross(l1, l2)); + CHECK( (vector3_norm(l3)>1.0e-6), "degenerate plane in normal_to_plane" ); matrix3x3 M; - M.c0 = v1; - M.c1 = v2; - M.c2 = v3; + M.c0 = l1; + M.c1 = l2; + M.c2 = l3; vector3 RHS = vector3_minus(p,o); vector3 tus = matrix3x3_vector3_mult(matrix3x3_inverse(M),RHS); // "t, u, s" - *s = tus.z; - return v3; + double t=tus.x, u=tus.y, s = tus.z; + int inside = ( ( 0.0<=t && t<=1.0 && 0.0<=u && u<=1.0 ) ? 1 : 0 ); + if (!quadrilateral && (t+u)>1.0) inside=0; // need (t+u)<1 for triangle interior + double d=0.0; + if(inside==0) + { vector3 pPlane = vector3_minus(p, vector3_scale(s,l3) ); + vector3 v01 = vector3_plus(o,l1); + vector3 v10 = vector3_plus(o,l2); + d=min_distance_to_line_segment(pPlane, o, v01); + d=fmin(d,min_distance_to_line_segment(pPlane, o, v10)); + if (quadrilateral) + { vector3 v11 = vector3_plus(v01,l2); + d=fmin(d,min_distance_to_line_segment(pPlane, v01, v11)); + d=fmin(d,min_distance_to_line_segment(pPlane, v11, v10)); + } + else + d=fmin(d,min_distance_to_line_segment(pPlane, v01, v10)); + } + *min_distance=sqrt(s*s+d*d); + return l3; } +vector3 normal_to_quadrilateral(vector3 o, vector3 v1, vector3 v2, vector3 p, + double *min_distance) +{ return normal_to_plane(o,v1,v2,p,1,min_distance); } + +vector3 normal_to_triangle(vector3 o, vector3 v1, vector3 v2, vector3 p, + double *min_distance) +{ return normal_to_plane(o,v1,v2,p,0,min_distance); } + + /***************************************************************/ /* find the face of the prism for which the normal distance */ /* from x to the plane of that face is the shortest, then */ @@ -2335,43 +2393,41 @@ vector3 normal_to_prism(prism *prsm, vector3 xc) int num_vertices = prsm->vertices.num_items; vector3 xp = prism_coordinate_c2p(prsm, xc); - vector3 axis = {0,0,0}; axis.z=height; + vector3 zhat = {0,0,1.0}, axis=vector3_scale(height, zhat); + if (height==0.0) + return prism_vector_p2c(prsm, zhat); vector3 retval; double min_distance=HUGE_VAL; + int nv; + for(nv=0; nv0.0) - { int nv; - for(nv=0; nv=argc) usage("too few arguments to --line"); + vector3 v1,v2; + sscanf(argv[narg+1],"%le",&(v1.x)); + sscanf(argv[narg+2],"%le",&(v1.y)); + sscanf(argv[narg+3],"%le",&(v1.z)); + sscanf(argv[narg+4],"%le",&(v2.x)); + sscanf(argv[narg+5],"%le",&(v2.y)); + sscanf(argv[narg+6],"%le",&(v2.z)); + printf("Min distance=%e\n",min_distance_to_line_segment(test_point,v1,v2)); + narg+=6; + } else if (!strcmp(argv[narg],"--dir")) { if (narg+5>=argc) usage("too few arguments to --lineseg"); sscanf(argv[narg+1],"%le",&(test_dir.x)); From 7c2368f94c46cd31d8b5ecf8b5a006162e137955 Mon Sep 17 00:00:00 2001 From: Homer Reid Date: Thu, 12 Jul 2018 12:10:57 -0400 Subject: [PATCH 3/8] re-updated normal_to_prism to behave correctly when the projection of the point into the plane does not lie within the prism face in question --- utils/geom.c | 134 ++++++++++++++++++++++++--------------------- utils/test-prism.c | 3 + 2 files changed, 75 insertions(+), 62 deletions(-) diff --git a/utils/geom.c b/utils/geom.c index 8e8f168..1d5f411 100644 --- a/utils/geom.c +++ b/utils/geom.c @@ -2320,64 +2320,77 @@ double min_distance_to_line_segment(vector3 p, vector3 v1, vector3 v2) } /***************************************************************/ -/* compute the minimum distance from a 3D point p to the */ -/* planar quadrilateral with vertices {o, o+l1, o+l2, o+l1+l2} */ -/* (if quadrilateral==1) */ -/* or the planar triangle with vertices {o, o+l1, o+l2} */ -/* (if quadrilateral==0). */ -/* algorithm: */ -/* (A) solve a 3x3 system to compute the projection of p into */ -/* the plane of the triangle/quadrilateral (call it pPlane)*/ -/* pPlane = p-s*l3 = o + t*l1 + u*l2 */ -/* where l3 is the normal to the surface and s,t,u are */ -/* unknowns. */ -/* (B) if pPlane lies within the triangle/quadrilateral, the */ -/* minimum distance is simply s. */ -/* (C) otherwise compute the minimum distance d from pProj */ -/* to any edge of the triangle/quadrilateral and set */ -/* minimum_distance = sqrt(s^2 + d^2). */ -/* return value is normal to plane. */ +/* compute the projection of a 3D point p into the plane */ +/* containing the three points {o, o+v1, o+v2}. */ +/* algorithm: solve a 3x3 system to compute the projection of */ +/* p into the plane (call it pPlane) */ +/* pPlane = p-s*v3 = o + t*v1 + u*v2 */ +/* where v3 is the normal to the plane and s,t,u */ +/* are unknowns. */ +/* on return, s is the normal distance. the return value is */ +/* the normal vector to the plane. */ +/* if in_quadrilateral is non-null it is set to 1 or 0 */ +/* according as pPlane does or does not lie in the */ +/* quadrilateral with vertices (o, o+v1, o+v2, o+v1+v2). */ /***************************************************************/ -vector3 normal_to_plane(vector3 o, vector3 l1, vector3 l2, vector3 p, - int quadrilateral, double *min_distance) +vector3 normal_to_plane(vector3 p, vector3 o, vector3 v1, vector3 v2, + double *s, int *in_quadrilateral) { - vector3 l3 = unit_vector3(vector3_cross(l1, l2)); - CHECK( (vector3_norm(l3)>1.0e-6), "degenerate plane in normal_to_plane" ); + vector3 v3 = unit_vector3(vector3_cross(v1, v2)); + CHECK( (vector3_norm(v3)>1.0e-6), "degenerate plane in project_point_into_plane" ); matrix3x3 M; - M.c0 = l1; - M.c1 = l2; - M.c2 = l3; + M.c0 = v1; + M.c1 = v2; + M.c2 = v3; vector3 RHS = vector3_minus(p,o); vector3 tus = matrix3x3_vector3_mult(matrix3x3_inverse(M),RHS); // "t, u, s" - double t=tus.x, u=tus.y, s = tus.z; - int inside = ( ( 0.0<=t && t<=1.0 && 0.0<=u && u<=1.0 ) ? 1 : 0 ); - if (!quadrilateral && (t+u)>1.0) inside=0; // need (t+u)<1 for triangle interior - double d=0.0; + double t=tus.x, u=tus.y; + *s = tus.z; + if (in_quadrilateral) + *in_quadrilateral = ( ( 0.0<=t && t<=1.0 && 0.0<=u && u<=1.0 ) ? 1 : 0 ); + return v3; +} + +vector3 normal_to_quadrilateral(vector3 p, vector3 o, vector3 v1, vector3 v2, + double *min_distance) +{ + double s, d=0.0; + int inside; + vector3 v3=normal_to_plane(p, o, v1, v2, &s, &inside); if(inside==0) - { vector3 pPlane = vector3_minus(p, vector3_scale(s,l3) ); - vector3 v01 = vector3_plus(o,l1); - vector3 v10 = vector3_plus(o,l2); - d=min_distance_to_line_segment(pPlane, o, v01); - d=fmin(d,min_distance_to_line_segment(pPlane, o, v10)); - if (quadrilateral) - { vector3 v11 = vector3_plus(v01,l2); - d=fmin(d,min_distance_to_line_segment(pPlane, v01, v11)); - d=fmin(d,min_distance_to_line_segment(pPlane, v11, v10)); - } - else - d=fmin(d,min_distance_to_line_segment(pPlane, v01, v10)); + { vector3 pPlane = vector3_minus(p, vector3_scale(s,v3) ); + vector3 p01 = vector3_plus(o,v1); + vector3 p10 = vector3_plus(o,v2); + vector3 p11 = vector3_plus(p01,v2); + d= min_distance_to_line_segment(pPlane, o, p01); + d=fmin(d,min_distance_to_line_segment(pPlane, o, p10)); + d=fmin(d,min_distance_to_line_segment(pPlane, p01, p11)); + d=fmin(d,min_distance_to_line_segment(pPlane, p11, p10)); } *min_distance=sqrt(s*s+d*d); - return l3; + return v3; } -vector3 normal_to_quadrilateral(vector3 o, vector3 v1, vector3 v2, vector3 p, - double *min_distance) -{ return normal_to_plane(o,v1,v2,p,1,min_distance); } +vector3 normal_to_prism_roof_or_ceiling(vector3 p, vector3 *vertices, int num_vertices, + double height, double *min_distance) +{ + vector3 o = {0.0,0.0,0.0}; o.z=height; + double s; + vector3 v3 = normal_to_plane(p, o, vertices[0], vertices[1], &s, 0); + vector3 pPlane = vector3_minus(p, vector3_scale(s,v3) ); -vector3 normal_to_triangle(vector3 o, vector3 v1, vector3 v2, vector3 p, - double *min_distance) -{ return normal_to_plane(o,v1,v2,p,0,min_distance); } + double d = 0.0; + if (point_in_polygon(pPlane.x,pPlane.y,vertices,num_vertices)==0) + { int nv; + d=min_distance_to_line_segment(pPlane,vertices[0],vertices[1] ); + for(nv=0; nv Date: Thu, 12 Jul 2018 14:34:50 -0400 Subject: [PATCH 4/8] updates --- utils/geom.c | 7 ++----- 1 file changed, 2 insertions(+), 5 deletions(-) diff --git a/utils/geom.c b/utils/geom.c index 1d5f411..cc597ee 100644 --- a/utils/geom.c +++ b/utils/geom.c @@ -2153,10 +2153,7 @@ boolean intersect_line_with_segment(double px, double py, double dx, double dy, // corresponding to the endpoints; for our purposes we count // the intersection if the plumb line runs through the t==0 vertex, but // NOT the t==1 vertex, to avoid double-counting for complete polygons. -#define DELTAT 1.0e-6 -#define TMIN (0.0) -#define TMAX (1.0-DELTAT) - return ( ( tTMAX ) ? 0 : 1 ); + return t < 0 || ((float)t) >= 1 ? 0 : 1; } // like the previous routine, but only count intersections if s>=0 @@ -2344,7 +2341,7 @@ vector3 normal_to_plane(vector3 p, vector3 o, vector3 v1, vector3 v2, M.c2 = v3; vector3 RHS = vector3_minus(p,o); vector3 tus = matrix3x3_vector3_mult(matrix3x3_inverse(M),RHS); // "t, u, s" - double t=tus.x, u=tus.y; + float t=tus.x, u=tus.y; *s = tus.z; if (in_quadrilateral) *in_quadrilateral = ( ( 0.0<=t && t<=1.0 && 0.0<=u && u<=1.0 ) ? 1 : 0 ); From d42b418ec54f2e8c85ca8e03b5cc0bce037d0115 Mon Sep 17 00:00:00 2001 From: Homer Reid Date: Sun, 22 Jul 2018 23:58:39 -0400 Subject: [PATCH 5/8] revise algorithm for intersecting line segment with prism --- utils/geom.c | 168 ++++++++++++++++++++++++++++----------------- utils/test-prism.c | 45 +++++++++++- 2 files changed, 149 insertions(+), 64 deletions(-) diff --git a/utils/geom.c b/utils/geom.c index cc597ee..b59597a 100644 --- a/utils/geom.c +++ b/utils/geom.c @@ -2204,19 +2204,19 @@ boolean point_in_prism(prism *prsm, vector3 xc) /***************************************************************/ /* insert s into slist, which is sorted in increasing order */ -/* has maximum length 2 */ +/* of |s-s_0| and has maximum length 2 */ /***************************************************************/ -int insert_s_in_list(double s, double slist[2], int length) +int insert_s_in_list(double s, double s0, double slist[2], int length) { if (length==0) { slist[0]=s; return 1; } - if ( s < slist[0] ) + if ( fabs(s-s0) < fabs(slist[0]-s0) ) { slist[1]=slist[0]; slist[0]=s; } - else if (length==1 || s < slist[1] ) + else if (length==1 || (fabs(s-s0) < fabs(slist[1]-s0)) ) slist[1]=s; return 2; @@ -2225,9 +2225,10 @@ int insert_s_in_list(double s, double slist[2], int length) /***************************************************************/ /* compute all values of s at which the line p+s*d intersects */ /* a prism face, retaining the first two entries in a list */ -/* of s values sorted in increasing order m */ +/* of s values sorted in increasing order of distance to s0. */ /***************************************************************/ -int intersect_line_with_prism(prism *prsm, vector3 p, vector3 d, double slist[2]) +int intersect_line_with_prism(prism *prsm, vector3 p, vector3 d, double slist[2], + double smin, double smax, double s0) { int num_vertices = prsm->vertices.num_items; vector3 *vertices = prsm->vertices.items; @@ -2257,6 +2258,8 @@ int intersect_line_with_prism(prism *prsm, vector3 p, vector3 d, double slist[2] vertices[nv], vertices[nvp1], &s) ) continue; + if ( ssmax ) continue; + // OK, we know the XY-plane projection of the line intersects the polygon edge; // now go back to 3D, ask for the z-coordinate at which the line intersects // the z-axis extrusion of the edge, and determine whether this point lies @@ -2268,7 +2271,7 @@ int intersect_line_with_prism(prism *prsm, vector3 p, vector3 d, double slist[2] if ( (z_intersectz_max) ) continue; - num_intersections=insert_s_in_list(s, slist, num_intersections); + num_intersections=insert_s_in_list(s, s0, slist, num_intersections); } // identify intersections with prism ceiling and floor faces @@ -2278,9 +2281,10 @@ int intersect_line_with_prism(prism *prsm, vector3 p, vector3 d, double slist[2] for(LowerUpper=0; LowerUpper<2; LowerUpper++) { double z0 = LowerUpper ? height : 0.0; double s = (z0 - p_prsm.z)/dz; + if ( ssmax) continue; if (!point_in_polygon(p_prsm.x + s*d_prsm.x, p_prsm.y+s*d_prsm.y, vertices, num_vertices)) continue; - num_intersections=insert_s_in_list(s,slist,num_intersections); + num_intersections=insert_s_in_list(s,s0,slist,num_intersections); } return num_intersections; } @@ -2290,11 +2294,25 @@ int intersect_line_with_prism(prism *prsm, vector3 p, vector3 d, double slist[2] /***************************************************************/ double intersect_line_segment_with_prism(prism *prsm, vector3 p, vector3 d, double a, double b) { - double slist[2]={0.0, 0.0}; - if (2!=intersect_line_with_prism(prsm, p, d, slist)) - return 0.0; - double ds = fmin(slist[1],b) - fmax(slist[0],a); - return ds > 0.0 ? ds : 0.0; + boolean a_inside = point_in_prism(prsm, vector3_plus(p, vector3_scale(a,d))); + boolean b_inside = point_in_prism(prsm, vector3_plus(p, vector3_scale(b,d))); + + // case with segment endpoints both inside or both outside is easy + if (a_inside==b_inside) + return a_inside ? fabs(b-a) : 0.0; + + double slist[2]; + if (a_inside==1 && b_inside==0) + { int ns=intersect_line_with_prism(prsm, p, d, slist, a, b, a); + CHECK(ns>0, "internal error 1 in ilsw_prism") + return slist[0]-a; + } + else // if (a_inside==0 && b_inside==1) + { int ns=intersect_line_with_prism(prsm, p, d, slist, a, b, b); + CHECK(ns>0, "internal error 2 in ilsw_prism") + return b-slist[0]; + } + return 0.0; // never get here } /***************************************************************/ @@ -2318,22 +2336,23 @@ double min_distance_to_line_segment(vector3 p, vector3 v1, vector3 v2) /***************************************************************/ /* compute the projection of a 3D point p into the plane */ -/* containing the three points {o, o+v1, o+v2}. */ +/* that contains the three points {o, o+v1, o+v2} and has */ +/* normal vector v3. */ /* algorithm: solve a 3x3 system to compute the projection of */ /* p into the plane (call it pPlane) */ /* pPlane = p-s*v3 = o + t*v1 + u*v2 */ /* where v3 is the normal to the plane and s,t,u */ /* are unknowns. */ -/* on return, s is the normal distance. the return value is */ -/* the normal vector to the plane. */ -/* if in_quadrilateral is non-null it is set to 1 or 0 */ -/* according as pPlane does or does not lie in the */ +/* the return value is the value of s (where pPlane = p-s*v3), */ +/* i.e. the minimum distance from p to the plane. */ +/* if in_quadrilateral is non-null it is set to 0 */ +/* or 1 according as pPlane does or does not lie in the */ /* quadrilateral with vertices (o, o+v1, o+v2, o+v1+v2). */ /***************************************************************/ -vector3 normal_to_plane(vector3 p, vector3 o, vector3 v1, vector3 v2, - double *s, int *in_quadrilateral) +double normal_distance_to_plane(vector3 p, + vector3 o, vector3 v1, vector3 v2, vector3 v3, + int *in_quadrilateral) { - vector3 v3 = unit_vector3(vector3_cross(v1, v2)); CHECK( (vector3_norm(v3)>1.0e-6), "degenerate plane in project_point_into_plane" ); matrix3x3 M; M.c0 = v1; @@ -2341,54 +2360,78 @@ vector3 normal_to_plane(vector3 p, vector3 o, vector3 v1, vector3 v2, M.c2 = v3; vector3 RHS = vector3_minus(p,o); vector3 tus = matrix3x3_vector3_mult(matrix3x3_inverse(M),RHS); // "t, u, s" - float t=tus.x, u=tus.y; - *s = tus.z; + float t=tus.x, u=tus.y, s=tus.z; if (in_quadrilateral) *in_quadrilateral = ( ( 0.0<=t && t<=1.0 && 0.0<=u && u<=1.0 ) ? 1 : 0 ); - return v3; + return s; } -vector3 normal_to_quadrilateral(vector3 p, vector3 o, vector3 v1, vector3 v2, - double *min_distance) +// like normal_distance_to_plane, but if pPlane (projection of point into plane) +// lies outside the quadrilateral {o,o+v1,o+v2,o+v1+v2} then take into account +// the in-plane distance from pPlane to the quadrilateral +double min_distance_to_quadrilateral(vector3 p, + vector3 o, vector3 v1, vector3 v2, vector3 v3) { - double s, d=0.0; int inside; - vector3 v3=normal_to_plane(p, o, v1, v2, &s, &inside); - if(inside==0) - { vector3 pPlane = vector3_minus(p, vector3_scale(s,v3) ); - vector3 p01 = vector3_plus(o,v1); - vector3 p10 = vector3_plus(o,v2); - vector3 p11 = vector3_plus(p01,v2); - d= min_distance_to_line_segment(pPlane, o, p01); - d=fmin(d,min_distance_to_line_segment(pPlane, o, p10)); - d=fmin(d,min_distance_to_line_segment(pPlane, p01, p11)); - d=fmin(d,min_distance_to_line_segment(pPlane, p11, p10)); - } - *min_distance=sqrt(s*s+d*d); - return v3; + double s=normal_distance_to_plane(p, o, v1, v2, v3, &inside); + if(inside==1) + return s; + vector3 pPlane = vector3_minus(p, vector3_scale(s,v3) ); + vector3 p01 = vector3_plus(o,v1); + vector3 p10 = vector3_plus(o,v2); + vector3 p11 = vector3_plus(p01,v2); + double d = min_distance_to_line_segment(pPlane, o, p01) ; + d = fmin(d, min_distance_to_line_segment(pPlane, o, p10) ); + d = fmin(d, min_distance_to_line_segment(pPlane, p01, p11) ); + d = fmin(d, min_distance_to_line_segment(pPlane, p11, p10) ); + return sqrt(s*s+d*d); +} + +double min_distance_to_prism_roof_or_ceiling(vector3 p, + vector3 *vertices, int num_vertices, + double height) +{ + vector3 o = {0.0,0.0,0.0}; o.z=height; + vector3 v3 = {0,0,1.0}; + double s = normal_distance_to_plane(p,o,vertices[0],vertices[1],v3,0); + vector3 pPlane = vector3_minus(p, vector3_scale(s,v3) ); + if (point_in_polygon(pPlane.x,pPlane.y,vertices,num_vertices)==1) + return s; + + int nv; + double d=min_distance_to_line_segment(pPlane,vertices[0],vertices[1] ); + for(nv=1; nvvertices.items; + int num_vertices = prsm->vertices.num_items; + double height = prsm->height; + + box->low = box->high = prism_coordinate_p2c(prsm, vertices[0]); + int nv, fc; + for(nv=0; nvlow.x = fmin(box->low.x, vc.x); + box->low.y = fmin(box->low.y, vc.y); + box->low.z = fmin(box->low.z, vc.z); + + box->high.x = fmax(box->high.x, vc.x); + box->high.y = fmax(box->high.y, vc.y); + box->high.z = fmax(box->high.z, vc.z); + } +} + static vector3 make_vector3(double x, double y, double z) { vector3 v; @@ -426,7 +453,7 @@ int main(int argc, char *argv[]) narg+=6; } else if (!strcmp(argv[narg],"--dir")) - { if (narg+5>=argc) usage("too few arguments to --lineseg"); + { if (narg+3>=argc) usage("too few arguments to --lineseg"); sscanf(argv[narg+1],"%le",&(test_dir.x)); sscanf(argv[narg+2],"%le",&(test_dir.y)); sscanf(argv[narg+3],"%le",&(test_dir.z)); @@ -475,6 +502,22 @@ int main(int argc, char *argv[]) printf("Wrote prism description to GNUPLOT file test-prism.gp.\n"); printf("Wrote prism description to GMSH file test-prism.geo.\n"); + geom_box prism_box; + my_get_prism_bounding_box(prsm, &prism_box); + f=fopen("test-prism-bb.gp","w"); + fprintf(f,"%e %e %e\n",prism_box.low.x, prism_box.low.y, prism_box.low.z); + fprintf(f,"%e %e %e\n",prism_box.high.x, prism_box.low.y, prism_box.low.z); + fprintf(f,"%e %e %e\n",prism_box.high.x, prism_box.high.y, prism_box.low.z); + fprintf(f,"%e %e %e\n",prism_box.low.x, prism_box.high.y, prism_box.low.z); + fprintf(f,"%e %e %e\n\n\n",prism_box.low.x, prism_box.low.y, prism_box.low.z); + + fprintf(f,"%e %e %e\n",prism_box.low.x, prism_box.low.y, prism_box.high.z); + fprintf(f,"%e %e %e\n",prism_box.high.x, prism_box.low.y, prism_box.high.z); + fprintf(f,"%e %e %e\n",prism_box.high.x, prism_box.high.y, prism_box.high.z); + fprintf(f,"%e %e %e\n",prism_box.low.x, prism_box.high.y, prism_box.high.z); + fprintf(f,"%e %e %e\n\n\n",prism_box.low.x, prism_box.low.y, prism_box.high.z); + printf("Wrote bounding box to GNUPLOT file test-prism-bb.gp.\n"); + /***************************************************************/ /* test point inclusion, normal to object, and line-segment */ /* intersection with specified data */ From ec96dea85b535556b0513a85bdb3729fa99b2401 Mon Sep 17 00:00:00 2001 From: Homer Reid Date: Mon, 23 Jul 2018 02:50:15 -0400 Subject: [PATCH 6/8] overhauled intersect_line_segment_with_prism to do the full exact calculation with no approximations or assumptions about proximity --- utils/geom.c | 65 ++++++++++++++++++++++++++++++++-------------------- 1 file changed, 40 insertions(+), 25 deletions(-) diff --git a/utils/geom.c b/utils/geom.c index b59597a..089a61e 100644 --- a/utils/geom.c +++ b/utils/geom.c @@ -69,7 +69,7 @@ using namespace ctlio; static vector3 prism_vector_p2c(prism *prsm, vector3 vp); static vector3 prism_vector_c2p(prism *prsm, vector3 vc); static void get_prism_bounding_box(prism *prsm, geom_box *box); -static double intersect_line_segment_with_prism(prism *prsm, vector3 p, vector3 d, double a, double b); +double intersect_line_segment_with_prism(prism *prsm, vector3 p, vector3 d, double a, double b); static boolean point_in_polygon(double px, double py, vector3 *vertices, int num_vertices); static boolean point_in_prism(prism *prsm, vector3 p); static vector3 normal_to_prism(prism *prsm, vector3 p); @@ -2222,13 +2222,21 @@ int insert_s_in_list(double s, double s0, double slist[2], int length) return 2; } +// comparator for qsort +static int dcmp(const void *pd1, const void *pd2) +{ double d1=*((double *)pd1), d2=*((double *)pd2); + return (d1d2) ? 1.0 : 0.0; +} + /***************************************************************/ /* compute all values of s at which the line p+s*d intersects */ -/* a prism face, retaining the first two entries in a list */ -/* of s values sorted in increasing order of distance to s0. */ +/* a prism face. */ +/* slist is a caller-allocated buffer with enough room for */ +/* at least num_vertices+2 doubles. on return it contains */ +/* the intersection s-values sorted in ascending order. */ +/* the return value is the number of intersections. */ /***************************************************************/ -int intersect_line_with_prism(prism *prsm, vector3 p, vector3 d, double slist[2], - double smin, double smax, double s0) +int intersect_line_with_prism(prism *prsm, vector3 p, vector3 d, double *slist) { int num_vertices = prsm->vertices.num_items; vector3 *vertices = prsm->vertices.items; @@ -2258,8 +2266,6 @@ int intersect_line_with_prism(prism *prsm, vector3 p, vector3 d, double slist[2] vertices[nv], vertices[nvp1], &s) ) continue; - if ( ssmax ) continue; - // OK, we know the XY-plane projection of the line intersects the polygon edge; // now go back to 3D, ask for the z-coordinate at which the line intersects // the z-axis extrusion of the edge, and determine whether this point lies @@ -2271,7 +2277,7 @@ int intersect_line_with_prism(prism *prsm, vector3 p, vector3 d, double slist[2] if ( (z_intersectz_max) ) continue; - num_intersections=insert_s_in_list(s, s0, slist, num_intersections); + slist[num_intersections++]=s; } // identify intersections with prism ceiling and floor faces @@ -2281,11 +2287,12 @@ int intersect_line_with_prism(prism *prsm, vector3 p, vector3 d, double slist[2] for(LowerUpper=0; LowerUpper<2; LowerUpper++) { double z0 = LowerUpper ? height : 0.0; double s = (z0 - p_prsm.z)/dz; - if ( ssmax) continue; if (!point_in_polygon(p_prsm.x + s*d_prsm.x, p_prsm.y+s*d_prsm.y, vertices, num_vertices)) continue; - num_intersections=insert_s_in_list(s,s0,slist,num_intersections); + slist[num_intersections++]=s; } + + qsort((void *)slist,num_intersections,sizeof(double),dcmp); return num_intersections; } @@ -2294,25 +2301,33 @@ int intersect_line_with_prism(prism *prsm, vector3 p, vector3 d, double slist[2] /***************************************************************/ double intersect_line_segment_with_prism(prism *prsm, vector3 p, vector3 d, double a, double b) { - boolean a_inside = point_in_prism(prsm, vector3_plus(p, vector3_scale(a,d))); - boolean b_inside = point_in_prism(prsm, vector3_plus(p, vector3_scale(b,d))); + double *slist = (double *)malloc( (2+prsm->vertices.num_items)*sizeof(double) ); + int num_intersections=intersect_line_with_prism(prsm, p, d, slist); - // case with segment endpoints both inside or both outside is easy - if (a_inside==b_inside) - return a_inside ? fabs(b-a) : 0.0; + // na=smallest index such that slist[na] > a + int na=-1; + int ns; + for(ns=0; na==-1 && nsa) + na=ns; - double slist[2]; - if (a_inside==1 && b_inside==0) - { int ns=intersect_line_with_prism(prsm, p, d, slist, a, b, a); - CHECK(ns>0, "internal error 1 in ilsw_prism") - return slist[0]-a; + if (na==-1) + { free(slist); + return 0.0; } - else // if (a_inside==0 && b_inside==1) - { int ns=intersect_line_with_prism(prsm, p, d, slist, a, b, b); - CHECK(ns>0, "internal error 2 in ilsw_prism") - return b-slist[0]; + + int inside = ( (na%2)==0 ? 0 : 1); + double last_s=a; + double ds=0.0; + for(ns=na; ns Date: Mon, 23 Jul 2018 02:55:33 -0400 Subject: [PATCH 7/8] update acquisition of array slices for DFT fields to behave properly in the presence of symmetries (2) --- utils/geom.c | 20 -------------------- 1 file changed, 20 deletions(-) diff --git a/utils/geom.c b/utils/geom.c index 089a61e..bed7e27 100644 --- a/utils/geom.c +++ b/utils/geom.c @@ -2202,26 +2202,6 @@ boolean point_in_prism(prism *prsm, vector3 xc) return point_in_polygon(xp.x, xp.y, vertices, num_vertices); } -/***************************************************************/ -/* insert s into slist, which is sorted in increasing order */ -/* of |s-s_0| and has maximum length 2 */ -/***************************************************************/ -int insert_s_in_list(double s, double s0, double slist[2], int length) -{ - if (length==0) - { slist[0]=s; - return 1; - } - if ( fabs(s-s0) < fabs(slist[0]-s0) ) - { slist[1]=slist[0]; - slist[0]=s; - } - else if (length==1 || (fabs(s-s0) < fabs(slist[1]-s0)) ) - slist[1]=s; - - return 2; -} - // comparator for qsort static int dcmp(const void *pd1, const void *pd2) { double d1=*((double *)pd1), d2=*((double *)pd2); From 7f18509a96a318e9c623afe8402c61ce2b508f3b Mon Sep 17 00:00:00 2001 From: Homer Reid Date: Mon, 23 Jul 2018 03:43:53 -0400 Subject: [PATCH 8/8] update to intersect_line_segment_with_prism to fix failing unit test --- utils/geom.c | 4 ++-- utils/test-prism.c | 13 ++++++++++++- 2 files changed, 14 insertions(+), 3 deletions(-) diff --git a/utils/geom.c b/utils/geom.c index bed7e27..be2bb1f 100644 --- a/utils/geom.c +++ b/utils/geom.c @@ -69,7 +69,7 @@ using namespace ctlio; static vector3 prism_vector_p2c(prism *prsm, vector3 vp); static vector3 prism_vector_c2p(prism *prsm, vector3 vc); static void get_prism_bounding_box(prism *prsm, geom_box *box); -double intersect_line_segment_with_prism(prism *prsm, vector3 p, vector3 d, double a, double b); +static double intersect_line_segment_with_prism(prism *prsm, vector3 p, vector3 d, double a, double b); static boolean point_in_polygon(double px, double py, vector3 *vertices, int num_vertices); static boolean point_in_prism(prism *prsm, vector3 p); static vector3 normal_to_prism(prism *prsm, vector3 p); @@ -2307,7 +2307,7 @@ double intersect_line_segment_with_prism(prism *prsm, vector3 p, vector3 d, doub last_s = this_s; } free(slist); - return ds; + return ds > 0.0 ? ds : 0.0; } /***************************************************************/ diff --git a/utils/test-prism.c b/utils/test-prism.c index b61df17..c67b6ea 100644 --- a/utils/test-prism.c +++ b/utils/test-prism.c @@ -321,7 +321,18 @@ int test_line_segment_intersection(geometric_object the_block, geometric_object num_failed++; if (f) - fprintf(f," %e %e %s\n",sblock,sprism,fabs(sblock-sprism) > 1.0e-6*fmax(fabs(sblock),fabs(sprism)) ? "fail" : "success"); + { + int success = fabs(sblock-sprism) <= 1.0e-6*fmax(fabs(sblock),fabs(sprism)); + fprintf(f," %e %e %s\n",sblock,sprism,success ? "success" : "fail"); + if (success==0) + { fprintf(f,"#%e %e %e %e %e %e %e %e\n",p.x,p.y,p.z,d.x,d.y,d.z,a,b); + fprintf(f,"%e %e %e\n%e %e %e\n%e %e %e\n", + p.x,p.y,p.z, + p.x+a*d.x,p.y+a*d.y,p.z+a*d.z, + p.x+b*d.x,p.y+b*d.y,p.z+b*d.z); + } + fprintf(f,"\n"); + } } if (f) fclose(f);