Check if ray intersects internals of $D$-facet

164 Views Asked by At

Given a ray $\overrightarrow{r_0} + \overrightarrow{v} \cdot t, t \in [0;+\infty)$ and a $(D - 1)$-simplex, defined by $D$-tuple of its vertices $p_i = (p_i^1, p_i^2, \dots, p_i^D), i \in \{1, 2, \dots, D\}$, in $D$-space.

I want to check whether the ray intersects an internals of the above $D$-facet or not.

What I think at the moment is to find the point common for the ray and the facet, then find its barycentric coordinates within the facet and to check whether all they are positive values of not. Firstly, it requires the normalized normal equation of hyperplane, spanned on points $p_i, i \in \{1, 2, \dots, D\}$:

$\sum \limits_{i = 1}^{D}\left(x_i \cdot \frac{n_i}{|\overrightarrow{n}|}\right) - \frac{d}{|\overrightarrow{n}|} = 0$

where $\overrightarrow{n} = (n_1, n_2, \dots, n_D)$ — a normal to the hyperplane. $\frac{d}{|\overrightarrow{n}|}$ — a distance from the origin to the hyperplane.

$ n_i = \begin{vmatrix} p_1^1 & p_1^2 & \dots & p_1^{i - 1} & 1 & p_1^{i + 1} & \dots & p_1^D \\ p_2^1 & p_2^2 & \dots & p_2^{i - 1} & 1 & p_2^{i + 1} & \dots & p_2^D \\ \vdots & \vdots & \ddots & \vdots & \vdots & \vdots & \ddots & \vdots \\ p_D^1 & p_D^2 & \dots & p_D^{i - 1} & 1 & p_D^{i + 1} & \dots & p_D^D \notag \end{vmatrix} \\ d = \begin{vmatrix} p_1^1 & p_1^2 & \dots & p_1^D \\ p_2^1 & p_2^2 & \dots & p_2^D \\ \vdots & \vdots & \ddots & \vdots \\ p_D^1 & p_D^2 & \dots & p_D^D \notag \end{vmatrix} $

Hence, an intersection point $a = (a_0, a_1, \dots, a_D)$:

$a = \overrightarrow{r_0} - \overrightarrow{v} \cdot \underbrace{\frac{(\overrightarrow{r_0} \cdot \overrightarrow{n}) - d}{(\overrightarrow{v} \cdot \overrightarrow{n})}}_\text{= -t}$

If $t \ge 0$ and $(\overrightarrow{v} \cdot \overrightarrow{n}) \ne 0$, then the ray intersects the hyperplane properly. Is it right?

Further, one need to compute a barycentric coordinates $c_i, i \in \{1, 2, \dots, D\}$ of the intersection point $a$ defined with respect to the facet:

$(c_0, c_1, \dots, c_D) \cdot \left( \begin{matrix} p_1^1 & p_1^2 & \dots & p_1^D \\ p_2^1 & p_2^2 & \dots & p_2^D \\ \vdots & \vdots & \ddots & \vdots \\ p_D^1 & p_D^2 & \dots & p_D^D \notag \end{matrix} \right) = (a_0, a_1, \dots, a_D) \Rightarrow c = a \cdot P^{-1}$

Am I right in the next assertion: if $c_i \ge 0, i \in \{1, 2, \dots, D\}$, then the $a$ lies on internals of the facet? And another (not too important) problem: is it right $\sum \limits_{i = 1}^D{c_i} = 1$ (it comes true automatically, at my mind)?

1

There are 1 best solutions below

0
On BEST ANSWER

I conclude my approach is correct. To test it I wrote simple GNU Octave script for D3 case:

user@ubuntu:~$ cat > ray_plane_intersection.m
#! /usr/bin/octave -qf
pkg load geometry
lo = -1
hi = 1
facet = unifrnd(lo, hi, 3);
line = [lo lo lo hi - lo hi - lo hi - lo];#unifrnd(lo, hi, 1, 6);#
plane = [facet(1, :) (facet(2, :) - facet(1, :)) (facet(3, :) - facet(1, :))];
point = intersectLinePlane(line, plane);
barycentric = facet'\point';
in_range = all(0 <= barycentric) && all(barycentric <= 1);
coeffs_sum = sum(barycentric);
assert(sumsq(facet' * barycentric - point') < eps);

printf("set label 1 at %f, %f, %f;\n", hi, hi, hi);
printf("set label 1 '%f %f %f';\n", barycentric(1), barycentric(2), barycentric(3));
printf("set xrange [%f:%f];\n", lo, hi);
printf("set yrange [%f:%f];\n", lo, hi);
printf("set zrange [%f:%f];\n", lo, hi);
printf("set view equal xyz;\n");
printf("set xyplane at 0;\n");
if (in_range)
    printf("set title 'hit %f';\n", coeffs_sum);
else
    printf("set title 'miss %f';\n", coeffs_sum);
endif
printf("splot '-' with lines notitle, '-' with lines notitle, '-' with points notitle;\n");
printf("%f ", facet(1, :));
printf("\n");
printf("%f ", facet(2, :));
printf("\n");
printf("%f ", facet(3, :));
printf("\n");
printf("%f ", facet(1, :));
printf("\n");
printf("e\n");
printf("%f ", line(1:3));
printf("\n");
printf("%f ", line(4:6) + line(1:3));
printf("\n");

printf("e\n");
printf("%f ", point(:));
printf("\n");
printf("e\n");
user@ubuntu:~$ chmod +x ray_plane_intersection.m
    user@ubuntu:~$ ./ray_plane_intersection.m | gnuplot -p
user@ubuntu:~$