star-2d

commit 0ba38f6c41404ef1792c8d5adbf5c7d6d5b2df58
parent f99e90503db6a45a10180c3bcb1b1c987ea657ec
Author: Vincent Forest <vincent.forest@meso-star.com>
Date:   Wed, 29 Jun 2016 14:33:09 +0200

Numerically compute PI*S/P in the advanced session test

Diffstat:
M
src/test_s2d_session_advanced.c
 | 
82
+++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++

1 file changed, 82 insertions(+), 0 deletions(-)
diff --git a/src/test_s2d_session_advanced.c b/src/test_s2d_session_advanced.c
@@ -34,6 +34,45 @@
 
 #include <string.h>
 
+#define NSAMPS 10000
+
+static INLINE float*
+ran_semi_disk_cos_local(float samp[2])
+{
+  samp[0] = rand_canonic() * 2.f - 1.f;
+  samp[1] = (float)sqrt(1 - samp[0]*samp[0]);
+  return samp;
+}
+
+static INLINE float*
+ran_semi_disk_cos(const float up[2], float samp[2])
+{
+  float tmp[2], v[2];
+  CHECK(f2_is_normalized(samp), 1);
+  ran_semi_disk_cos_local(tmp);
+  v[0] = -up[1] * tmp[0] + up[0]*tmp[1];
+  v[1] =  up[0] * tmp[0] + up[1]*tmp[1];;
+  CHECK(f2_is_normalized(v), 1);
+  return f2_set(samp, v);
+}
+
+static int
+discard_self_hit
+  (const struct s2d_hit* hit,
+   const float org[2],
+   const float dir[2],
+   void* ray_data,
+   void* filter_data)
+{
+  struct s2d_primitive* prim_from = ray_data;
+  NCHECK(hit, NULL);
+  NCHECK(org, NULL);
+  NCHECK(dir, NULL);
+  CHECK((intptr_t)filter_data, (intptr_t)0xDECAFBAD);
+  if(!ray_data) return 0;
+  return S2D_PRIMITIVE_EQ(prim_from, &hit->prim);
+}
+
 static struct s2d_shape*
 create_circle_shape
   (struct s2d_device* dev,
@@ -50,6 +89,8 @@ create_circle_shape
   unsigned i;
 
   CHECK(s2d_shape_create_line_segments(dev, &shape), RES_OK);
+  CHECK(s2d_line_segments_set_hit_filter_function
+    (shape, discard_self_hit, (void*)0xDECAFBAD), RES_OK);
 
   CHECK(nsteps > 4, 1);
   NCHECK(center, NULL);
@@ -102,6 +143,9 @@ main(int argc, char** argv)
   size_t i;
   size_t nprims;
   float tmp[2];
+  float sum, sum_sqr;
+  float E, V, SE;
+  float area, length;
   int ishape;
   (void)argc, (void)argv;
 
@@ -198,6 +242,44 @@ main(int argc, char** argv)
   FOR_EACH(i, 0, sa_size(shape_prims[2])) CHECK(shape_prims[2][i], 1);
   FOR_EACH(i, 0, sa_size(scene_prims)) CHECK(scene_prims[i], 1);
 
+  /* Check the ray tracing by numerically compute PI*S/P aka 4V/S in 2D */
+  sum = sum_sqr = 0;
+  FOR_EACH(i, 0, NSAMPS) {
+    const float range[2] = { 0.f, FLT_MAX };
+    struct s2d_attrib attr;
+    struct s2d_hit hit;
+    float P[2], N[2];
+    float s;
+
+    CHECK(s2d_scene_sample(scn, rand_canonic(), rand_canonic(), &prim, &s), RES_OK);
+
+    CHECK(s2d_primitive_get_attrib(&prim, S2D_POSITION, s, &attr), RES_OK);
+    CHECK(attr.type, S2D_FLOAT2);
+    f2_set(P, attr.value);
+
+    CHECK(s2d_primitive_get_attrib(&prim, S2D_GEOMETRY_NORMAL, s, &attr), RES_OK);
+    CHECK(attr.type, S2D_FLOAT2);
+    NCHECK(f2_normalize(N, attr.value), 0.f);
+
+    f2_normalize(tmp, f2(tmp, 1, 1));
+    ran_semi_disk_cos(N, tmp);
+
+    CHECK(s2d_scene_trace_ray(scn, P, tmp, range, &prim, &hit), RES_OK);
+    CHECK(S2D_HIT_NONE(&hit), 0);
+
+    sum += hit.distance;
+    sum_sqr += hit.distance*hit.distance;
+  }
+
+  CHECK(s2d_scene_compute_contour_length(scn, &length), RES_OK);
+  CHECK(s2d_scene_compute_area(scn, &area), RES_OK);
+
+  E = sum / (float)NSAMPS;
+  V = sum_sqr / (float)NSAMPS - E*E;
+  SE = (float)sqrt(V/(float)NSAMPS);
+  printf("PI*S / P = %g ~ %g +/- %g\n",(float)PI*area / length, E, SE);
+  CHECK(eq_epsf((float)PI*area / length, E, SE), 1);
+
   CHECK(s2d_scene_end_session(scn), RES_OK);
 
   CHECK(s2d_scene_ref_put(scn), RES_OK);