Perturbing a vector within a given max angle

Question

I am supposed to implement a fuzzy mirror in my raytracer that perturbs the reflection vector randomly within a maximum given angle around the perfect reflection vector. So my idea was to sample a disk that is orthogonal to the reflection vector at distance 1 to the hitpoint. Problem is: in the provided code framework there's no way of accessing the hit location in the sample code so there's no way (?) of creating a random point on a disc, cause I can't define a disc with just a normal. Can anyone tell me how to perturb a vector randomly within a given max angle? I have access to the incoming light direction and the normal at the hitpoint and the maximum derivation angle.

Answer 1

A vector is just a length in some particular direction, right? Otherwise how would someone working in your software environment produce the "perfect" reflection vector without knowing the hitting point?

So you should not care where the hitting point is. Just assume the ray is reflected at the origin. Make a vector of length $1$ in the direction of the "perfect" reflection. There is a unique point displaced from the origin by that vector. Use that point as the center of your disk.

The random vector you get in this fashion will be exactly parallel to (and the same length as) the vector you would have gotten by taking a vector from the actual hitting point to a point displaced $1$ unit along the perfect reflection vector from that hitting point and making a disk around that displaced point. That is, it will be the same vector.

Answer 2

I put together an outline for an approach, along with some test code in Python which hopefully illustrates the idea in a practical setting. I hope this helps.

Step 1: Reflect the incident ray by negating the z coordinate.
Step 2: Compute the coordinates of the reflected ray in spherical system.
Step 3: Calculate the transformation matrix $T'$ to align the $Z'$-direction with the reflected ray and the $X'$-direction to be contained in the plane of reflection. This can be accomplished using the computed $\theta$ and $\phi$ values from Step 2.
Step 4: Choose a random elevation angle, $\theta_{random}$, from a uniform distribution in $\left[0,\theta_{max}\right]$.
Step 5: Choose random azimuth angle, $\phi_{random}$, from a uniform distribution in $\left[-\pi,\pi\right]$.
Step 6: Compute the updated reflected ray as $\left(r,\theta_{random},\phi_{random}\right)$.
Step 7: Convert this back to cartesian coordinates.
Step 8: Use the inverse (which is actually the transpose) of the transformation matrix, $T'$, from Step 3 to convert this new reflected ray back to the original coordinate system.

# Simulate the random reflection vector
import numpy
import random
random.seed(1234)
from mpl_toolkits.mplot3d import axes3d
import matplotlib.pyplot as plt

# incident ray
incident_ray = numpy.array((1.0,0.0,-1.0))  # incident ray in original x-y-z coordinates
theta_max = (numpy.pi/2)/2  # max angle to consider for random vector direction.

# **Step 1:**  Reflect the incident ray by negating the z coordinate.

reflected_ray = incident_ray*numpy.array([1.0,1.0,-1.0])

print('incident ray = ', incident_ray)
print('reflected ray = ', reflected_ray)
print('theta_max = ', theta_max)

# **Step 2:**  Compute the coordinates of the reflected ray in spherical system.

r_ref = numpy.linalg.norm(reflected_ray)
phi_ref=numpy.arctan2(reflected_ray[1],reflected_ray[0])  # range is [-pi,pi]
theta_ref = numpy.arccos(reflected_ray[2]/r_ref)

print('r_ref = ', r_ref)
print('phi_ref = ', phi_ref)
print('theta_ref = ', theta_ref)

# **Step 3:**  Calculate the transformation matrix $T'$ to align the 
# $Z'$-direction with the reflected ray and the $X'$-direction 
# to be contained in the plane of reflection.  
# This can be accomplished using the computed $\theta$ and $\phi$ values from Step 2.

# We actually only need the inverse transformation for later:

T_inv = numpy.array([[numpy.cos(theta_ref)*numpy.cos(phi_ref), -numpy.sin(phi_ref), numpy.sin(theta_ref)*numpy.cos(phi_ref)],
                     [numpy.cos(theta_ref)*numpy.sin(phi_ref), numpy.cos(phi_ref), numpy.sin(theta_ref)*numpy.sin(phi_ref)],
                     [-numpy.sin(theta_ref), 0, numpy.cos(theta_ref)]])

# **Step 4:**  Choose a random angle in the range 0 to $\theta_{max}$ 

# plot in 3D 
fig = plt.figure()
ax = fig.gca(projection='3d', 
             azim=90+phi_ref*360/(2*numpy.pi), 
             elev=0*(-theta_ref*360/(2*numpy.pi)))

num_reflections=100  # number of random reflecting vectors to compute

for iRef in range(num_reflections):
#    print('----')
    rand_theta = random.uniform(0,theta_max)

# **Step 5:**  Choose random azimuth angle $\phi$ from a uniform distribution from -pi to pi$.

    rand_phi = random.uniform(-numpy.pi,numpy.pi)  # random azimuth angle

#    print('rand_theta = ', rand_theta)
#    print('rand_phi = ', rand_phi)

# **Step 6:**  Compute the updated reflected ray as $\left(r,\theta_{random},\phi_{random}\right)$.
# we can skip this and compute the cartesion coordinates directly in step 7.

# **Step 7:**  Convert this back to cartesian coordinates.

    x_prime = r_ref*numpy.sin(rand_theta)*numpy.cos(rand_phi)
    y_prime = r_ref*numpy.sin(rand_theta)*numpy.sin(rand_phi)
    z_prime = r_ref*numpy.cos(rand_theta)

    prime_vector = numpy.array([x_prime, y_prime, z_prime])

# **Step 8:**  Use the inverse (which is actually the transpose) of the transformation matrix, $T'$, 
# from Step 3 to convert this new reflected ray back to the original coordinate system.

#    prime_vector = numpy.array([0,0,r_ref])  # check to see this inverts back to original
    new_reflected_ray = numpy.matmul(T_inv, prime_vector)
#    print('new_reflected_ray = ', new_reflected_ray)
#    print('new_reflected_ray (mag) = ', numpy.linalg.norm(new_reflected_ray))
#    print('new_reflected_ray (angle relative to original) = ', 
#          numpy.arccos(numpy.dot(new_reflected_ray,reflected_ray)/
#          (numpy.linalg.norm(new_reflected_ray)*numpy.linalg.norm(reflected_ray))))

# Make the grid
    x, y, z = 0,0,0

# Make the direction data for the arrows
    u = new_reflected_ray[0]
    v = new_reflected_ray[1]
    w = new_reflected_ray[2]

    ax.quiver(x, y, z, u, v, w, length=0.038, normalize=True)

#plot original

u = reflected_ray[0]
v = reflected_ray[1]
w = reflected_ray[2]

ax.quiver(x, y, z, u, v, w, length=0.1, normalize=True) 
plt.show()