star-gs

commit 801f6391227317fa38f74aee82e1366f381c2b63
parent 6f0c04f4efc298cb78ffccd492849d61e6cd8e4c
Author: Paule Lapeyre <paule.lapeyre@yahoo.fr>
Date:   Wed,  6 Jul 2022 17:19:30 +0200

Reformulating the projections steps

Diffstat:
M
src/sgs_compute_sensitivity_translation.c
 | 
76
++++++++++++++++++++++++++++++++++++----------------------------------------

1 file changed, 36 insertions(+), 40 deletions(-)
diff --git a/src/sgs_compute_sensitivity_translation.c b/src/sgs_compute_sensitivity_translation.c
@@ -49,6 +49,12 @@ static const double EMIT_S_SZ[3] = {1, 1, 0};
 /* Vecteur de déformation (à renommer en gamma ?) */
 static const double V[3] = {0, 0, 1};
 
+struct deformation_2d {
+  double alpha;
+  double beta;
+  double u[3];
+};
+
 /*******************************************************************************
  * Helper functions
  ******************************************************************************/
@@ -104,6 +110,24 @@ compute_grad_I
   * (1 - cos(2*PI*position[X]/size[X]));
 }
 
+static void
+projection
+  (const double chi[3],
+   const double normal[3],
+   const double omega[3],
+   struct deformation_2d* out_vector)
+{
+  ASSERT(chi && normal && omega && out_vector);
+  ASSERT(d3_is_normalized(normal));
+  ASSERT(d3_is_normalized(omega));
+
+  out_vector->alpha = d3_dot(chi, normal) / d3_dot(omega, normal);
+  out_vector->u[X] = chi[X] - out_vector->alpha*omega[X];
+  out_vector->u[Y] = chi[Y] - out_vector->alpha*omega[Y];
+  out_vector->u[Z] = chi[Z] - out_vector->alpha*omega[Z];
+  out_vector->beta = d3_normalize(out_vector->u, out_vector->u);
+}
+
 static res_T
 realisation
   (void* weights,
@@ -116,6 +140,10 @@ realisation
   struct sgs_hit hit = SGS_HIT_NULL;
   struct sgs_fragment frag = SGS_FRAGMENT_NULL;
 
+  struct deformation_2d proj_u_e;
+  struct deformation_2d proj_V_e;
+  struct deformation_2d proj_V_s;
+
   double recv_min[2];
   double recv_max[2];
   double emit_e_threshold;
@@ -138,20 +166,6 @@ realisation
 
   double range[2];
 
-  double u_emit_s[3];
-  double u_spec_s[3];
-  double u_emit_e[3];
-  double u_spec_e[3];
-
-  double alpha_emit_e;
-  double alpha_spec_e;
-  double alpha_emit_s;
-  double alpha_spec_s;
-  double beta_emit_e;
-  double beta_spec_e;
-  double beta_emit_s;
-  double beta_spec_s;
-
   double rho;
   double grad_rho[3], grad_rho_2d[2];
   double I;
@@ -196,7 +210,7 @@ realisation
   surf_emit_s = frag.surface;
 
   /* Sample the cosine weighted sampling of the emissive direction */
-  ssp_ran_hemisphere_cos(rng, frag.normal, dir_emit_s, NULL);
+  ssp_ran_hemisphere_cos(rng, normal_s, dir_emit_s, NULL);
 
   /* Trace the ray from the sensitivity emissive surface */
   sgs_geometry_trace_ray
@@ -217,7 +231,7 @@ realisation
   }
 
   /* Find the reflected position */
-  reflect(dir_spec_s, dir_emit_s, frag.normal);
+  reflect(dir_spec_s, dir_emit_s, normal_s);
   sgs_geometry_trace_ray
     (sgs->geom, pos_emit_s, dir_spec_s, range, surf_emit_s, &hit);
   if(SGS_HIT_NONE(&hit)) {
@@ -237,31 +251,13 @@ realisation
     goto exit;
   }
 
-  alpha_emit_s = d3_dot(V, normal_s) / d3_dot(dir_emit_s, normal_s);
-  alpha_spec_s = d3_dot(V, normal_s) / d3_dot(dir_spec_s, normal_s);
-  u_emit_s[X] = V[X] - alpha_emit_s*dir_emit_s[X];
-  u_emit_s[Y] = V[Y] - alpha_emit_s*dir_emit_s[Y];
-  u_emit_s[Z] = V[Z] - alpha_emit_s*dir_emit_s[Z];
-  u_spec_s[X] = V[X] - alpha_spec_s*dir_spec_s[X];
-  u_spec_s[Y] = V[Y] - alpha_spec_s*dir_spec_s[Y];
-  u_spec_s[Z] = V[Z] - alpha_spec_s*dir_spec_s[Z];
-  beta_emit_s = d3_normalize(u_emit_s, u_emit_s);
-  beta_spec_s = d3_normalize(u_spec_s, u_spec_s);
-
   dir_spec_e[X] = -dir_spec_s[X];
   dir_spec_e[Y] = -dir_spec_s[Y];
   dir_spec_e[Z] = -dir_spec_s[Z];
 
-  alpha_emit_e = d3_dot(u_emit_s, normal_e) / d3_dot(dir_spec_e, normal_e);
-  alpha_spec_e = d3_dot(u_spec_s, normal_e) / d3_dot(dir_spec_e, normal_e);
-  u_emit_e[X] = u_emit_s[X] - alpha_emit_e*dir_spec_e[X];
-  u_emit_e[Y] = u_emit_s[Y] - alpha_emit_e*dir_spec_e[Y];
-  u_emit_e[Z] = u_emit_s[Z] - alpha_emit_e*dir_spec_e[Z];
-  u_spec_e[X] = u_spec_s[X] - alpha_spec_e*dir_spec_e[X];
-  u_spec_e[Y] = u_spec_s[Y] - alpha_spec_e*dir_spec_e[Y];
-  u_spec_e[Z] = u_spec_s[Z] - alpha_spec_e*dir_spec_e[Z];
-  beta_emit_e = d3_normalize(u_emit_e, u_emit_e);
-  beta_spec_e = d3_normalize(u_spec_e, u_spec_e);
+  projection(V, normal_e, dir_spec_e, &proj_V_e);
+  projection(V, normal_s, dir_emit_s, &proj_V_s);
+  projection(proj_V_s.u, normal_e, dir_spec_e, &proj_u_e);
 
   pos_emit_s_2d[X] = pos_emit_s[X];
   pos_emit_s_2d[Y] = pos_emit_s[Y];
@@ -284,9 +280,9 @@ realisation
   grad_I[Z] = grad_I_2d[Y];
 
   S =
-    - beta_emit_s * d3_dot(grad_rho, u_emit_s) * I
-    - rho * beta_emit_s * beta_emit_e * d3_dot(grad_I, u_emit_e)
-    + rho * beta_spec_s * beta_spec_e * d3_dot(grad_I, u_spec_e);
+    - proj_V_s.beta * d3_dot(grad_rho, proj_V_s.u) * I
+    - rho * proj_V_s.beta * proj_u_e.beta * d3_dot(grad_I, proj_u_e.u)
+    + rho * proj_V_e.beta * d3_dot(grad_I, proj_V_e.u);
 
   w = I*emit_s_sz[X]*emit_s_sz[Y]*PI; /* Weight */
   s = S*emit_s_sz[X]*emit_s_sz[Y]*PI; /* Sensib */