How to imlement the Phong BRDF

Question

I'm really sure that this question was asked before, and I found some code snippets by searching this via google. Unfortunately these snippets are in openGL or some other high level computer graphics language. I have an input light vector I(phi,theta). I can calculate the reflected vector quite easily. The Phong BRDF assumes that this vector is the z-axis and one can generate a random direction to this R-vector. I want to know, how I can add this vector to R(phi,theta) to generate a random view direction out of the input direction I(phi,theta), say V(phi,theta). I know this can be done by rotation matrices and their transposed ones, but I'm stuck. Can you please give some algebraic formulas for this problem?

Physics behind this problem: I have a cylinder which points to the sun. This cylinder is tilted by an angle alpha. The inside of this cylinder is painted "black" and I want to know the output of this arrangement at the end of the cylinder in dependence of the tilt angle alpha.

Thanks in advance,

Uwe

Answer 1

The following are the steps you need to follow:

1. Get the normal vector for the reflecting point as $n$.
2. Calculate reflected direction from normal vector and $I$ in the world frame, as $I_r$
3. Sample Phong BRDF. This will yield a local frame vector as $r$. The next thing you need to do is to convert this local vector to world frame vector (details are below).
4. Once you converted this to world frame as $r_w$, the view direction you want will be $r_w$. Now you just need to evaluate the BRDF, return PDF, sampled direction and throughput.

I assume the difficult part for you is the 3rd step. How do I convert a local vector to the world vector? How does the frame get defined? Let's answer the second question first.

The Phong model proposes that the view direction should be sampled around the reflected ray direction. So we know that one of the axis of the coordinate frame is the reflected ray direction. The convention is that this vector usually serves as the z-axis for this local frame. That is, if you are talking about $(0, 0, 1)^T$ in this local sampling frame, the actual vector in the world frame for $(0, 0, 1)^T$ is the reflected ray direction. Since Phong model is isotropic on $\phi$ direction, the other two axes are free to choose.

Next, how to convert from local frame to world frame? Say we have this $r$ which is sampled in a frame of which the z-axis is $(0, 0, 1)^T$, and we want to convert it as if it was sampled in a frame of which the z-axis is the reflected ray direction ($I_r$). You actually need to calculate the rotation matrix between $(0, 0, 1)^T$ and $I_r$. This can be done with the following steps:

1. Calculate the cross product of $(0, 0, 1)^T$ and $I_r$. Since we know that cross product will yield a vector that is perpendicular to both input vectors, this will be our axis of rotation, that is, we will later rotate around this axis.
2. Calculate the dot product of $(0, 0, 1)^T$ and $I_r$. This will yield the cosine value for the angle between the two vectors, and actually, if we rotate one of the vectors around the axis we get in step 1 by this angle, we will get the other vector.
3. So, we rotate the vector around the cross-product axis by the dot-product angle. This can be achieved by Rodrigues' rotation formula.

The following is the code for Rodrigues' rotation formula. np_rotation_between uses the transform lib provided by scipy and rotation_between is hand-crafted (though in Taichi Lang).

from scipy.spatial.transform import Rotation as Rot

@ti.func
def skew_symmetry(vec: vec3):
    return mat3([
        [0, -vec[2], vec[1]], 
        [vec[2], 0, -vec[0]], 
        [-vec[1], vec[0], 0]
    ])
    
def np_rotation_between(fixed: Arr, target: Arr) -> Arr:
    """
        Transform parsed from xml file is merely camera orientation (numpy CPU version)
        INPUT arrays [MUST] be normalized
        Orientation should be transformed to be camera rotation matrix
        Rotation from <fixed> vector to <target> vector, defined by cross product and angle-axis
    """
    axis = np.cross(fixed, target)
    dot = np.dot(fixed, target)
    if abs(dot) > 1. - 1e-5:            # nearly parallel
        return np.sign(dot) * np.eye(3, dtype = np.float32)  
    else:
        # Not in-line, cross product is valid
        axis /= np.linalg.norm(axis)
        axis *= np.arccos(dot)
        euler_vec = Rot.from_rotvec(axis).as_euler('zxy')
        euler_vec[0] = 0                                                # eliminate roll angle
        return Rot.from_euler('zxy', euler_vec).as_matrix()

@ti.func
def rotation_between(fixed: vec3, target: vec3) -> mat3:
    """
        Transform parsed from xml file is merely camera orientation (Taichi version)
        Rodrigues transformation is implemented here
        INPUT arrays [MUST] be normalized
        Orientation should be transformed to be camera rotation matrix
        Rotation from <fixed> vector to <target> vector, defined by cross product and angle-axis
    """
    axis = tm.cross(fixed, target)
    cos_theta = tm.dot(fixed, target)
    ret_R = ti.Matrix.zero(float, 3, 3)
    if ti.abs(cos_theta) < 1. - 1e-5:
        normed_axis = axis.normalized()
        ret_R = ti.Matrix.diag(3, cos_theta) + ((1 - cos_theta) * colv3(*normed_axis)) @ rowv3(*normed_axis) + skew_symmetry(axis)
    else:
        ret_R = ti.Matrix.diag(3, tm.sign(cos_theta))
    return ret_R

With this, you can get the rotation matrix between $(0, 0, 1)^T$ and $I_r$. Let's suppose you get the rotation matrix to rotate $(0, 0, 1)^T$ to $I_r$ as $R$, then the sampled view direction for you will be: $$ v = R(0, 0, 1)^T $$ Otherwise, if the rotation is in the opposite direction, just use $R^T$ for this matrix-vector multiplication. There you go, the core mathematical formulae you need are (1) reflection (2) Phong local sampling (3) Rodrigues' rotation formula.