star-2d

commit ae991ae98d757135bf8584f7f25d283ec1edb643
parent 5ce3cf61ba1c35e1585a913d8d494b307718b842
Author: Christophe Coustet <christophe.coustet@meso-star.com>
Date:   Thu, 28 Jan 2021 09:36:54 +0100

Improve accuracy of closest point requests for huge segments

Diffstat:
M
src/s2d_scene_view_closest_point.c
 | 
59
+++++++++++++++++++++++++++++++++--------------------------
M
src/test_s2d_closest_point.c
 | 
2
+-

2 files changed, 34 insertions(+), 27 deletions(-)
diff --git a/src/s2d_scene_view_closest_point.c b/src/s2d_scene_view_closest_point.c
@@ -33,6 +33,7 @@
 #include "s2d_scene_view_c.h"
 
 #include <rsys/float2.h>
+#include <rsys/double2.h>
 
 struct point_query_context {
   struct RTCPointQueryContext rtc;
@@ -45,40 +46,48 @@ struct point_query_context {
  ******************************************************************************/
 static INLINE float*
 closest_point_segment
-  (const float p[2], /* Position */
-   const float v0[2], /* 1st segment vertex */
-   const float v1[2], /* 2nd segment vertex */
-   const float E[2], /* v0 -> v1 vector */
-   float closest_pt[2], /* Closest position */
-   float* s) /* Parametric coordinate of the closest point onto [v0, v1] */
+  (const float p_f[2], /* Position */
+   const float v0_f[2], /* 1st segment vertex */
+   const float v1_f[2], /* 2nd segment vertex */
+   float closest_pt_f[2], /* Closest position */
+   float* s_f) /* Parametric coordinate of the closest point onto [v0, v1] */
 {
-  float v[2]; /* Vector from v0 to p */
-  float segment_len2; /* Square length of the [v0, v1] */
-  float dst_x_seglen; /* |p' v0| x |v0 v1|
+  double v[2]; /* Vector from v0 to p */
+  double E[2]; /* v0 -> v1 vector */
+  double p[2], v0[2], v1[2], closest_pt[2];
+  double segment_len2; /* Square length of the [v0, v1] */
+  double dst_x_seglen; /* |p' v0| x |v0 v1|
                        * p' is the orthogonal projection of p onto [v0, v1] */
-  ASSERT(p && v0 && v1);
-  ASSERT(f2_eq(f2_sub(v, v1, v0), E));
+  double s;
+  ASSERT(p_f && v0_f && v1_f);
+
+  d2_set_f2(v0, v0_f);
+  d2_set_f2(v1, v1_f);
+  d2_set_f2(p, p_f);
+  d2_sub(E, v1, v0);
 
   /* Orthogonally project the point onto the segment */
-  f2_sub(v, p, v0);
-  dst_x_seglen = f2_dot(v, E);
+  d2_sub(v, p, v0);
+  dst_x_seglen = d2_dot(v, E);
 
   /* Check if the closest point is the segment vertex 'v0' */
   if(dst_x_seglen <= 0) {
-    *s = 0;
-    return f2_set(closest_pt, v0);
+    *s_f = 0;
+    return f2_set(closest_pt_f, v0_f);
   }
   /* Check if the closest point is the segment vertex 'v1' */
-  segment_len2 = f2_dot(E,E);
+  segment_len2 = d2_dot(E, E);
   if(dst_x_seglen >= segment_len2) {
-    *s = 1;
-    return f2_set(closest_pt, v1);
+    *s_f = 1;
+    return f2_set(closest_pt_f, v1_f);
   }
 
   /* The closest point is on the segment */
-  *s = dst_x_seglen / segment_len2;
-  ASSERT(*s == CLAMP(*s, 0, 1));
-  return f2_add(closest_pt, f2_mulf(closest_pt, E, *s), v0);
+  s = dst_x_seglen / segment_len2;
+  d2_add(closest_pt, d2_muld(closest_pt, E, s), v0);
+  *s_f = (float)s;
+  ASSERT(*s_f == CLAMP(*s_f, 0, 1));
+  return f2_set_d2(closest_pt_f, closest_pt);
 }
 
 static bool
@@ -92,7 +101,6 @@ closest_point_line_segments
   struct hit_filter* filter = NULL;
   const uint32_t* ids = NULL;
   float v0[2], v1[2]; /* Segment vertices */
-  float E[2]; /* Segment vector */
   float N[2]; /* Segment normal */
   float query_pos[2]; /* Submitted position */
   float closest_point[2]; /* Computed closest point */
@@ -109,21 +117,20 @@ closest_point_line_segments
   ASSERT(geom->lines->attribs_type[S2D_POSITION] == S2D_FLOAT2);
   f2_set(v0, line_segments_get_pos(geom->lines)+ids[0]*2/*#coords per vertex */);
   f2_set(v1, line_segments_get_pos(geom->lines)+ids[1]*2/*#coords per vertex */);
-  f2_sub(E, v1, v0);
 
   query_pos[0] = args->query->x;
   query_pos[1] = args->query->y;
 
   /* Compute the closest point on the segment from the submitted query_pos */
-  closest_point_segment(query_pos, v0, v1, E, closest_point, &s);
+  closest_point_segment(query_pos, v0, v1, closest_point, &s);
 
   f2_sub(vec, closest_point, query_pos);
   dst = f2_len(vec);
   if(dst >= args->query->radius) return 0;
 
   /* Compute the segment normal (left hand convention) */
-  N[0] = E[1];
-  N[1] = -E[0];
+  N[0] = v1[1] - v0[1];
+  N[1] = v0[0] - v1[0];
 
   /* Flip the geometry normal wrt the flip line_segments flag */
   if(geom->flip_contour) f2_minus(N, N);
diff --git a/src/test_s2d_closest_point.c b/src/test_s2d_closest_point.c
@@ -405,7 +405,7 @@ test_single_segment(struct s2d_device* dev)
   CHK(s2d_scene_ref_put(scn) == RES_OK);
 
   /* Check accuracy on a configuration whose analytic distance is known */
-  FOR_EACH(a, 0, 16) {
+  FOR_EACH(a, 0, 19) {
     const float amplitude = exp2f((float)a);
     const float eps = 1e-6f * amplitude;
     FOR_EACH(i, 0, 1000) {