Finding ray direction with smallest angle from a ray to a rectangle in 3D

Question

In 3D I have been given a ray and a rectangle:

$$Ray=\mathbf{r_o}+t*\mathbf{r_d}$$ $$Rect=\mathbf{p_o}+p*\mathbf{p_x}+q*\mathbf{p_y}$$

where

• $\mathbf{r_o}$ is the ray origin
• $\mathbf{r_d}$ is the ray unit direction vector
• $t>0$
• $\mathbf{p_o}$ is the rectangle center
• $-rectWidth/2 < p < rectWidth/2$
• $-rectHeight/2 < q < rectHeight/2$
• $\mathbf{p_x}$ is the rectangle unit x-vector
• $\mathbf{p_y}$ is the rectangle unit y-vector
• $\mathbf{p_x}$ and $\mathbf{p_y}$ are orthogonal, i.e. $\mathbf{p_x} \cdot \mathbf{p_y} = 0$

I'm trying to find

$$Ray'=\mathbf{r_o} + t*\mathbf{r_d'}$$

that intersects the rectangle and minimizes the angle between $\mathbf{r_d}$ and $\mathbf{r_d'}$ (in other words maximizes $\mathbf{r_d} \cdot \mathbf{r_d'}$). If $\mathbf{r_d}$ intersects the rectangle then naturally $\mathbf{r_d'}=\mathbf{r_d}$ and otherwise $\mathbf{r_d'}$ points towards the rectangle edge.

I have tried couple of solutions which unfortunately do not give correct results:

1) Calculate $p$ and $q$ at the intersection between the ray and the plane of the rectangle, and clamp $p$ and $q$ to the rectangle domain. Then calculate $\mathbf{v}=\mathbf{p_o}+p*\mathbf{p_x} + q*\mathbf{p_y} - \mathbf{r_o}$ and $\mathbf{r'_d}=\mathbf{v}/||\mathbf{v}||$. This works "mostly right", but breaks down when the ray points away from the rectangle (i.e. there's no intersection). Even if the ray points towards the rectangle, clamping $p$ and $q$ between $-rectWidth/2 < p < rectWidth/2$, and $-rectHeight/2 < q < rectHeight/2$ doesn't give correct result in all cases.

2) Project the ray to the rectangle plane and calculate intersection of the infinite line going along the projected ray and the rectangle. If the line intersects the rectangle, test which of the two intersection point gives the smallest angle. If the line doesn't intersect the rectangle, test which of the 4 corners of the rectangle gives the smallest angle. This seems to be quite wrong result and worse than option 1)

I don't need the exact solution for my application, so if there are some approximations that doesn't have too large errors I'm happy to try those out. If it helps, the rectangle parameters can also be simplified with affine transformation of the ray so that $\mathbf{p_o}=[0, 0, 0]$, $\mathbf{p_x}=[1, 0, 0]$, $\mathbf{p_y}=[0, 1, 0]$, $rectWidth=1$ and $rectHeight=1$.

Answer 1

Here is a robust 4 step solution that gives correct result for any $Ray$ and $Rect$ according to my tests.

Step 1:

Construct 4 vectors from the origin of the ray $\mathbf{r_o}$ to each corner of the rectangle $Rect$, in counter-clockwise order as seen from $\mathbf{r_o}$. Let's call these vectors $\mathbf{v_1}$...$\mathbf{v_4}$:

$$\mathbf{v_1}=\mathbf{p_o}-\mathbf{r_o}+0.5*rectWidth*\mathbf{p_x}+0.5*rectHeight*\mathbf{p_y}$$ $$\mathbf{v_2}=\mathbf{p_o}-\mathbf{r_o}+0.5*rectWidth*\mathbf{p_x}-0.5*rectHeight*\mathbf{p_y}$$ $$\mathbf{v_3}=\mathbf{p_o}-\mathbf{r_o}-0.5*rectWidth*\mathbf{p_x}-0.5*rectHeight*\mathbf{p_y}$$ $$\mathbf{v_4}=\mathbf{p_o}-\mathbf{r_o}-0.5*rectWidth*\mathbf{p_x}+0.5*rectHeight*\mathbf{p_y}$$

To ensure that these points are in counter-clockwise order as seen from $\mathbf{r_o}$, check if $$(\mathbf{p_x}\times\mathbf{p_y})\cdot(\mathbf{p_o}-\mathbf{r_o})<0$$ If the above condition doesn't hold swap vectors $\mathbf{v_1}$ and $\mathbf{v_3}$ to swap the winding order of the points to satisfy the winding requirement.

Step 2:

Construct 4 normal vectors $\mathbf{n_1}$...$\mathbf{n_4}$ for planes that lay on each edge of $Rect$ and $\mathbf{r_o}$ using the above vectors as follows: $$\mathbf{n_1}=\frac{\mathbf{v_1} \times \mathbf{v_2}}{||\mathbf{v_1} \times \mathbf{v_2}||}$$ $$\mathbf{n_2}=\frac{\mathbf{v_2} \times \mathbf{v_3}}{||\mathbf{v_2} \times \mathbf{v_3}||}$$ $$\mathbf{n_3}=\frac{\mathbf{v_3} \times \mathbf{v_4}}{||\mathbf{v_3} \times \mathbf{v_4}||}$$ $$\mathbf{n_4}=\frac{\mathbf{v_4} \times \mathbf{v_1}}{||\mathbf{v_4} \times \mathbf{v_1}||}$$

Step 3:

Out of these 4 normals pick a normal $\mathbf{n_{min}}$ and $\mathbf{v_a}$ and $\mathbf{v_b}$ that were used for cross product for the normal, that results smallest value for $\mathbf{n_i} \cdot \mathbf{r_d}$ (i.e. the most perpendicular plane to $\mathbf{r_d}$). For example if $\mathbf{n_{min}}=\mathbf{n_2}$, then $\mathbf{v_a}=\mathbf{v_2}$ and $\mathbf{v_b}=\mathbf{v_3}$.

If $\mathbf{n_{min}} \cdot \mathbf{r_d}>=0$, then $\mathbf{r'_d}=\mathbf{r_d}$ since $Ray$ intersects $Rect$. If not, continue to the last step

Step 4:

Project the ray direction $\mathbf{r_d}$ to the plane defined by $\mathbf{n_{min}}$: $$\mathbf{r_{pd}}=\mathbf{r_d}-(\mathbf{n_{min}} \cdot \mathbf{r_d})*\mathbf{n_{min}}$$ And clamp $\mathbf{r_{pd}}$ to the edge defined by $\mathbf{v_a}$ and $\mathbf{v_b}$ as follows to get the final result $\mathbf{r'_d}$: $$ \mathbf{r'_d}= \begin{cases} \mathbf{\hat{v_a}} & \text{if }(\mathbf{r_{pd}} \times \mathbf{v_a}) \cdot \mathbf{r_d}<0\\ \mathbf{\hat{v_b}} & \text{if }(\mathbf{v_b} \times \mathbf{r_{pd}}) \cdot \mathbf{r_d}<0\\ \mathbf{\hat{r_{pd}}} & \text{otherwise}\\ \end{cases} $$ Because I need to calculate steps 1 and 2 for solid angle calculation anyway, the additional cost of this algorithm is only steps 3 and 4 in my case.