<p>The following are the steps you need to follow:</p>
<ol>
<li>Get the normal vector for the reflecting point as <span class="math-container">$n$</span>.</li>
<li>Calculate reflected direction from normal vector and <span class="math-container">$I$</span> in the <strong>world frame</strong>, as <span class="math-container">$I_r$</span></li>
<li>Sample Phong BRDF. This will yield a <strong>local frame</strong> vector as <span class="math-container">$r$</span>. The next thing you need to do is to convert this local vector to <strong>world frame</strong> vector (details are below).</li>
<li>Once you converted this to world frame as <span class="math-container">$r_w$</span>, the view direction you want will be <span class="math-container">$r_w$</span>. Now you just need to evaluate the BRDF, return PDF, sampled direction and throughput.</li>
</ol>
<p>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.</p>
<p>The Phong model proposes that the view direction should be sampled around the <strong>reflected ray direction</strong>. 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 <span class="math-container">$(0, 0, 1)^T$</span> in this local sampling frame, the actual vector in the world frame for <span class="math-container">$(0, 0, 1)^T$</span> is the reflected ray direction. Since Phong model is isotropic on <span class="math-container">$\phi$</span> direction, the other two axes are free to choose.</p>
<p>Next, how to convert from local frame to world frame? Say we have this <span class="math-container">$r$</span> which is sampled in a frame of which the z-axis is <span class="math-container">$(0, 0, 1)^T$</span>, and we want to convert it as if it was sampled in a frame of which the z-axis is the reflected ray direction (<span class="math-container">$I_r$</span>). You actually need to calculate the rotation matrix between <span class="math-container">$(0, 0, 1)^T$</span> and <span class="math-container">$I_r$</span>. This can be done with the following steps:</p>
<ol>
<li>Calculate the cross product of <span class="math-container">$(0, 0, 1)^T$</span> and <span class="math-container">$I_r$</span>. Since we know that cross product will yield a vector that is perpendicular to both input vectors, this will be our <strong>axis of rotation</strong>, that is, we will later rotate around this axis.</li>
<li>Calculate the dot product of <span class="math-container">$(0, 0, 1)^T$</span> and <span class="math-container">$I_r$</span>. This will yield <strong>the cosine value for the angle between the two vectors</strong>, 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.</li>
<li>So, we rotate the vector around the cross-product axis by the dot-product angle. This can be achieved by <a href="https://en.wikipedia.org/wiki/Rodrigues%27_rotation_formula" rel="nofollow noreferrer">Rodrigues' rotation formula</a>.</li>
</ol>
<p>The following is the code for Rodrigues' rotation formula. <code>np_rotation_between</code> uses the transform lib provided by <code>scipy</code> and <code>rotation_between</code> is hand-crafted (though in Taichi Lang).</p>
<pre class="lang-python prettyprint-override"><code>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) -&gt; Arr:
    &quot;&quot;&quot;
        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 &lt;fixed&gt; vector to &lt;target&gt; vector, defined by cross product and angle-axis
    &quot;&quot;&quot;
    axis = np.cross(fixed, target)
    dot = np.dot(fixed, target)
    if abs(dot) &gt; 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) -&gt; mat3:
    &quot;&quot;&quot;
        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 &lt;fixed&gt; vector to &lt;target&gt; vector, defined by cross product and angle-axis
    &quot;&quot;&quot;
    axis = tm.cross(fixed, target)
    cos_theta = tm.dot(fixed, target)
    ret_R = ti.Matrix.zero(float, 3, 3)
    if ti.abs(cos_theta) &lt; 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

</code></pre>
<p>With this, you can get the rotation matrix between <span class="math-container">$(0, 0, 1)^T$</span> and <span class="math-container">$I_r$</span>. Let's suppose you get the rotation matrix to rotate <span class="math-container">$(0, 0, 1)^T$</span> to <span class="math-container">$I_r$</span> as <span class="math-container">$R$</span>, then the sampled view direction for you will be:
<span class="math-container">$$
v = R(0, 0, 1)^T
$$</span>
Otherwise, if the rotation is in the opposite direction, just use <span class="math-container">$R^T$</span> 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.</p>

CC BY-SA 4.0
2023-11-04T02:06:26.690Z
2023-11-04T02:06:26.690Z
Enigmatisms
<p>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?</p>
<p>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 &quot;black&quot; and I want to know the output of this arrangement at the end of the cylinder in dependence of the tilt angle alpha.</p>
<p>Thanks in advance,</p>
<p>Uwe</p>

<raytracing><brdf>
How to imlement the Phong BRDF
Uwe Nass
# How to imlement the Phong BRDF

(Asked by: username_0 on 2023-11-03)

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

**Tags:** raytracing, brdf


---

## 1 answers

### Answer by username_1 on 2023-11-04. Score: 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](https://en.wikipedia.org/wiki/Rodrigues%27_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).

```python
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.


---

