Disclaimer: This article is mainly about my personal notes on this black dog hair paper in SIG23. There should be no threshold for reading, because I myself have not entered the graphics field. There are formulas, but not many. The formula tags are the same as the original paper. All the content is just my tinkering. If there are any misunderstandings in the formulas, please correct me. Thank you very much!
Mengqi Xia, Bruce Walter, Christophe Hery, Olivier Maury, Eric Michielssen, and Steve Marschner, “A Practical Wave Optics Reflection Model for Hair and Fur,” ACM Transactions on Graphics (TOG), vol. 42, no. 4, article 39, pp. 1-15, Jul. 2023.
1. Related Work
Ray-based fiber model 在人类毛发的研究中，Marschner[2003]的模型广泛应用于业界，通过分析介电圆柱和圆锥中的光线路径，将散射分为R、TT和TRT。Zinke[2009]加入了漫反射成分。Sadeghi[2010]提出了艺术家控制参数化方法。d’Eon[2014]和Huang[2022]提出了非可分离的表征方法，也就是说方位角和纵向角度之间有耦合效应，不可简单分离。Chiang[2016] 进一步优化了该模型，使其适用于生产级渲染。 除了人类毛发，Khungurn和Marschner[2017]探讨了椭圆形毛发的建模。Yan[2015, 2017]研究了带有内部髓质的动物毛发。Aliaga[2017]将模型推广至有更加复杂截面的纺织纤维。
Fiber Model Based on Wave Optics Linder [2014] resolved cylindrical fibers with perfectly circular cross-sections and investigated the scattering behavior. Xia [2020] demonstrated several important differences compared to the geometric optics model by performing two-dimensional wave simulations on cylinders with arbitrary cross-sections. However, this model is a regular circle and does not have microscopic geometric structures such as hair scales. Benamira and Pattanaik [2021] proposed a faster hybrid model that uses wave optics only for forward scattering and relies on geometric optics for the rest.
Plane model simulation based on physical optics Physical optics planar models use physical optics approximations to simulate the scattering behavior of light on nearly planar, rough surfaces that can be represented as height fields, including Gaussian random, periodic, precomputed, and scratched surfaces. They combine Kirchhoff scalar diffraction theory with path tracing methods to handle scattering and reflection, and calculate the diffraction of light on different rough surfaces through various models such as Beckmann-Kirchhoff and Harvey-Shack. Although these models are effective on planar surfaces, the closed nature of fiber geometry and complex light interactions require more complex treatments.
Computational Electromagnetics,Spickle Effect and Stylized noise This has been mentioned in the previous article, so I will skip it here.
Hair modeling is based on scanning electron microscope (SEM) images of hair, which can accurately restore the microstructure of hair, such as hair scales and roughness.
Next, use "WAVE SIMULATION WITH 3D FIBER MICROGEOMETRY" to calculate the reflection and diffraction of the rough fiber surface.
On this basis, speckle theory is introduced to analyze the statistical characteristics of scattering patterns, and noise is used to describe these speckles, which greatly optimizes the model.
Then, by comparing with the actual measured data, reasonable fiber parameters (such as size, skin angle and surface roughness) are derived and finally integrated into the rendering system.
3. OVERVIEW
3.1 Fiber scattering models
I translate this into fiber scattering model. This describes the interaction between single fibers. The key here isBCSDF (Bidirectional Curvilinear Scattering Distribution Function)Unlike the bidirectional scattering distribution function (BSDF) commonly used in surface reflection and refraction, the BCSDF is designed specifically for curved fibers. The following formula states that when a given wavelength of light is irradiated on a fiber, it will be reflected or transmitted from another direction after entering from one direction. $$ L_r(\omega_r, \lambda) = \int L_i(\omega_i, \lambda) S(\omega_i, \omega_r, \lambda) \cos \theta_i d\omega_i \tag{1} $$ The left side of the formula represents the given wavelength $\lambda$,EmissionThe radiant brightness in the direction $\omega_r$. $ L_i(\omega_i, \lambda)$ isIncidentThe radiance in the direction $\omega_i$. $S(\omega_i, \omega_r, \lambda)$ is the bidirectional curvilinear scattering distribution function, which describes how the light is "scattered" by the fiber. $\cos \theta_i$ is to take into account the effect of the angle of incidence. If the light hits the fiber at a very flat angle, its effect will be smaller than when it hits it vertically.
Writing the above formula in spherical coordinates, and treating each different interaction of light with the hair fiber as a different mode, and then summing up the different scattering terms: $$ S(\theta_i, \theta_r, \phi_i, \phi_r, \lambda) = \sum_{p=0}^{\infty} S_p(\theta_i, \theta_r, \phi_i, \phi_r, \lambda) \tag{2} $$ The first scattering term $S_0$ describes the surface reflection, which is often referred to as the direct reflection term $R$ in the past. This term generally represents the statistical average of the reflection from the smooth fiber or the rough fiber surface. This term can be calculated more accurately here in the paper.
Recall that previous rendering methods (e.g. Marschner [2003]) typically decompose each scattering pattern $S_p$ into two separate functions: the longitudinal function $M_p$ and the anazimuthal function $N_p$ . This approach has been criticized as being inaccurate and should be avoided using the following separable approximation. $$ S_p(\theta_i, \theta_r, \phi_i, \phi_r, \lambda) = M_p(\theta_i, \theta_r)N_p(\theta_i, \phi_i, \phi_r, \lambda) \tag{3} $$ Therefore, XIA[2023] samples the scattering parameters of multiple rough fibers and takes the average, denoted as $S_0,avg$ . $f(\theta_h, \phi_h, \lambda)$ is the noise component, represented by two half-range vector angles and wavelength, which is used to correct the deviation of the current specific fiber instance from the mean. $$ S_{0,\text{sim}}(\theta_i, \theta_r, \phi_i, \phi_r, \lambda) \approx S_{0,\text{avg}}(\theta_i, \theta_r, \phi_i, \phi_r , \lambda) f(\theta_h, \phi_h, \lambda) \tag{4} $$ Currently, the complete fiber scattering model is as follows. The first term represents the scattering mode of reflection from the fiber surface, combined with 3D wave simulation. The subsequent summation terms are the sum of other higher-order scattering modes. $$ S_{\text{prac}}(\theta_i, \theta_r, \phi_i, \phi_r, \lambda) = S_{0,\text{prac}}(\theta_i, \theta_r, \phi_i, \phi_r, \lambda ) + \sum_{p=1}^{\infty} S_p(\theta_i, \theta_r, \phi_i, \phi_r, \lambda) \tag{5} $$ In practice, the final scattering formula is much simpler to accommodate surfaces with more complex geometry.
3.2 Speckle theory
Speckle theory describes the random light intensity distribution phenomenon produced when light interacts with a rough surface. The $A$ in Goodman's [2007] formula represents the superposition of all phase vectors, i.e. the resulting phase vector (Phasor): $$ \mathbf{A} = \frac{1}{\sqrt{N}} \sum_{n=1}^{N} a_n = \frac{1}{\sqrt{N}} \sum_{n=1} ^{N} a_n e^{i\phi_n} $$ The above formula is a comprehensive expression of the speckle intensity, in other words, it is used to express the intensity of scattered light. The reason for this phenomenon is that light reflects, refracts and interferes with each other between many tiny surfaces. The phase difference in the propagation of light leads to uneven light and dark patterns. Through this formula, it is possible to statistically calculate how light interferes with each other on the fiber surface to obtain the speckle pattern.
4. WAVE SIMULATION WITH 3D FIBER MICROGEOMETRY
In wave optics simulations, light is considered an electromagnetic wave. It consists of magnetic and electric fields that are perpendicular to each other. The interaction of light with hair (such as scattering) can be transformed into an analysis of how the electromagnetic field is affected by the object.
4.1 Wave optics
First of all, it is important to understand that when the electromagnetic field changes with a time period, it is called a time-harmonic field. In a time-harmonic field, the electric and magnetic fields can be represented by complex numbers as phase vectors (phasors). The clever thing is that the electric and magnetic fields are perpendicular to each other. $$ E_{\text{inst}} = \Re(E e^{j\omega t}), \quad H_{\text{inst}} = \Re(H e^{j\omega t}) \tag{6} $$ They are the electric field and magnetic field, respectively, but the real part is taken separately. The complex part contains the amplitude and phase information of the field.
The following two sets of equations areMaxwell's equationsIn the time-harmonic field form: $$ \nabla \times \mathbf{E} = -\mathbf{M} – j\omega \mu \mathbf{H} \ \nabla \times \mathbf{H} = \mathbf{J} + j\omega \varepsilon \mathbf{E} \tag{7} $$ Among them, $\varepsilon$ is the dielectric constant (affects the electric field), and $\mu$ is the magnetic permeability (affects the magnetic field). These two terms describe how the material affects the propagation of electromagnetic fields.
When an object (such as a fiber) is illuminated by an incident wave, the incident electric and magnetic fields are represented by $\mathbf{E}_i$ and $\mathbf{H}_i$, respectively. But the presence of the object changes these fields, so that what we observe isTotal Field. $$ \mathbf{E}_1 = \mathbf{E}_i + \mathbf{E}_s, \quad \mathbf{H}_1 = \mathbf{H}_i + \mathbf{H}_s \tag{8} $$ In short, total field = incident field + scattered field.
The light energy propagates outward along with the scattered fields $\mathbf{E}_s$ and $\mathbf{H}_s$, so calculating the scattering function of the fiber is the key. How to calculate it?
Full-wave simulation is the most accurate. Full-wave simulation requires discretizing the object into a grid and requires high resolution (generally, at least 10 grid cells per wavelength), which requires processingMillions of grid cells. It was also mentioned in the previous article.
That is to say, even simulating a short fiber segment that is only tens of microns long requires processing millions of grid cells.
Therefore, using Physical Optics Approximation (PO)In PO, the electric current and magnetic current on the surface of the object (respectively denoted as $J$ and $M$) can be regarded as the secondary source of the scattered field. The electromagnetic current generates secondary radiation, forming scattered waves. PO assumes that only a single reflection occurs on the surface of the object, ignoring multiple reflections and complex diffraction effects. After obtaining the electric current and magnetic current on the surface, the scattered waves generated by them are calculated. The properties of these far-field waves are derived from the surface electric current and magnetic current.
Just one PO is not enough, an octree algorithm must be incorporated to accelerate far-field calculations.
4.2 Physical Optics Approximation
Specifically, PO makes two simplifying assumptions: single scattering and local plane assumption. This method is also general enough.
As shown in the figure above, multiple points are sampled on the surface of an object, which is approximated as a plane. The current and magnetic field on the tangent plane are calculated to generate a scattering field. Through this scattering field, the far-field scattered waves generated by these currents and magnetic currents are calculated. At the same time, the Octree structure is used to divide it into smaller voxels. The calculation results of each leaf node are aggregated to the parent node to obtain the total radiation contribution.
Surface current calculation.
The calculation of surface current and magnetic current is a critical step in PO simulation. For each sampling point, its normal vector $n(r{\prime})$ and area information are stored. Next, the interaction between the incident wave and the surface current is calculated for each small plane.
According to the direction of the incident wave $e_i$ and the surface normal $n(r{\prime})$, the incident electric field is decomposed into parallel polarization and perpendicular polarization components. The reason for the decomposition is that parallel polarization and perpendicular polarization are completely different expressions in the Fresnel equation.
Parallel polarization component: The electric field of light is parallel to the reflection plane of the incident light. Vertically polarized component: The electric field of light is perpendicular to the reflection plane of the incident light.
Therefore, the incident field $E_i$ is decomposed into a parallel polarization component $E_i^p$ and a perpendicular polarization component $E_i^s$, which looks like this: $E_i = E_i^p + E_i^s$.
The reflected field $E_r$ is expressed as the sum of parallel and perpendicular polarization components, and the coefficients next to the incident field represent the reflection from the Fresnel equations: $$ E_r = E_r^p + E_r^s = F^p E_i^p + F^s E_i^s \tag{9-1} $$ Then, the total electric field $E_1$ is the sum of the incident and reflected fields: $$ E_1 = E_i + E_r \tag{9-2} $$ According to high school physics, we know that electricity can generate magnetism, and magnetism can also generate electricity. Therefore, we have the following expression: $$ M = -n \times E_1, J = n \times H_1 \tag{10} $$ Therefore, the induced current and magnetic flux on the surface can be calculated based on the incident field and the reflected field. After obtaining the current and magnetic flux, the scattered wave in the far field is calculated.
In theory, the incident field can be anything. But here the author uses Gaussian-windowed plane waves. The amplitude of this wave follows a normal distribution and is easy to calculate.
To summarize, here we decompose the incident light into two components, and use the Fresnel equation to calculate the reflected field. Then we add the reflected field to the incident field to get the total field. Through the relationship between the electric/magnetic field and the electric/magnetic current, we calculate the induced electric/magnetic current on the surface. In this way, we can simulate the scattering of the fiber.
That is to say,Calculation of electric and magnetic currents on the fiber surfaceThe scattering behavior of light can indeed be obtained.
Far-field radiation in 3D
In the previous section, we used the surface currents to obtain the surface currents and magnetic fluxes. This section uses this information to calculate the far-field scattering.
based onHuygens's PrincipleConvert the original scattering problem into a radiation problem.
Huygens' principle states that the electricity at each wavefront can be considered a new secondary wave source. It feels like one thing leads to another.
The electric current $J$ and magnetic current $M$ on the surface of the hair fiber are regarded as secondary sources of light, which are then reradiated to obtain the scattered field.
This is a bit difficult, as it involves the Method of Moments (MoM) in electromagnetism. Readers can study Gibson's [2021] book in depth. If you learn it, you must teach me well. But it doesn't matter. I found a picture from Professor Yan's 2021 paper, and I guarantee you can understand it.
The red ball in the figure representsSecondary radiation source, each radiates outward, similar to the secondary radiation waves generated by surface currents and magnetic currents. Secondary radiation sourceEmit waves of equal intensity in all directionsThis is consistent with the idea that each surface point in the scattering formula in this paper contributes to scattering in all directions. The figure shows a small area ( $\delta \mathbf{r}$ ) far from the light source (the distance is $r$ ).Far field areaIn , the behavior of the analytical wave at different locations $\mathbf{r}_1$ and $\mathbf{r}_2$ is observed along the two directions $\hat{r}_1$ and $\hat{r}_2$ respectively. This is very similar to the behavior of the scattered electric field in the far field in the scattering formula. The light wave on the right side of the figure willFar field superpositionComplex interference fringes are formed, which is similar to the integral term in the scattering formula, superimposing the electromagnetic waves of the secondary radiation source at a distance.
The formula is: $$ E_s(\mathbf{r}) = j \omega \mu_0 \frac{e^{-jk_0 R}}{4\pi R} \hat{r} \times \int_\Gamma \left[ \hat{r} \times \mathbf{J}(r{\prime}) + \frac{1}{Z_0} \mathbf{M}(r{\prime}) \right] e^{jk_0 r{\prime} \cdot \hat{r}} d\mathbf{r{\prime}} \tag{11} $$ Combined with the figure, the formula describes the performance of the scattered electric field $E_s(\mathbf{r})$ at a certain point $\mathbf{r}$ far away (far field), and decays with the distance in a relationship of $1/R$. This formula is a two-dimensional complex plane, that is, the solution of Maxwell's equations is regarded asTime-harmonic electromagnetic wavesin the form of.
看这个上式结构，形式上跟 $\mathbf{E}(\mathbf{r}, t) = \mathbf{E}_0 e^{j(\omega t – k \cdot \mathbf{r})}$ 很像，叫电场的常见时谐解。
$ j \omega \mu_0$ This imaginary term describesEffect of Magnetic Field on Electric Field, which is related to the frequency of the electromagnetic field$\omega$ and the magnetic permeability$\mu_0$ in the vacuum. The higher the frequency of the light wave, the stronger the electric field.$e^{-jk_0 R}$ is the phase factor, which represents the phase change of the light wave during propagation. The wave number$k_0 = \frac{2\pi}{\lambda}$, $\lambda$ is the wavelength of the electromagnetic wave. It means that when the wave propagates to a distance R, the phase of the electric field will change.$\hat{r}$ is the unit vector pointing from the scattering object to the observation point, indicating the direction of wave propagation. The cross product$\times$ operator indicates that the direction of the electric field is calculated. Make sure that the calculated electric field is consistent with the direction of wave propagation.
The integral term $ \int_\Gamma \left[ \hat{r} \times \mathbf{J}(r{\prime}) + \frac{1}{Z_0} \mathbf{M}(r{\prime}) \right] e^{jk_0 r{\prime} \cdot \hat{r}} d\mathbf{r{\prime}} $ is the core of the formula. By integrating $\Gamma$ on the surface of the object, we can get the contribution of each point on the hair surface to the scattered electric field. The surface current and surface magnetic current at each point are summed, and then multiplied by the phase factor. Furthermore, $\mathbf{J}(r{\prime})$ is the surface current density, $\hat{r} \times$ ensures that the scattered electric field generated by the current is orthogonal to the wave propagation direction$\hat{r}$. $Z_0$ is the free space impedance, and the specific value$Z_0 = \sqrt{\frac{\mu_0}{\varepsilon_0}} \approx 377 \, \Omega$. I personally understand that this constant is the proportional relationship between the electric field and the magnetic field in a vacuum. The electric field and the magnetic field have different impedances, so they need to be scaled to be linearly added, that is, the magnetic current is normalized to a form similar to the current. The exponential term is also a phase factor, so I won't explain it in detail.
Careful readers may find out why the surface current density $\mathbf{J}(r{\prime})$ has a cross product $\hat{r} \times \mathbf{J}(r{\prime})$ , while the surface magnetic current density $\mathbf{M}(r{\prime})$ does not.
According to Maxwell's equations, the electric field and magnetic field are orthogonal to each other. The current $\mathbf{J}$ will generate a magnetic field, and the changing magnetic field will in turn generate an electric field.Electric and magnetic fields are coupled. However, let us revisit the following two parts of Maxwell's equations. $$ \nabla \times \mathbf{E} = -\frac{\partial \mathbf{B}}{\partial t} \ \nabla \times \mathbf{B} = \mu_0 \mathbf{J} + \mu_0 \varepsilon_0 \frac{\partial \mathbf{E}}{\partial t} $$ The generation of magnetic field is achieved directly by displacing electric field, and the direction is orthogonal to the electric field. The generation of electric field is the rate of change of magnetic field, and the directionality needs to be adjusted.
To understand it from another perspective, the left side of the integral term describesThe strength of the electromagnetic field decreases with distance, phase variation and directivity are also taken into account. The integral term describesElectromagnetic current and phase contribution at each surface point.
In summary, formula (11) is an accurate expression that describes how the electromagnetic field on the surface of an object is scattered at a distance of $r$. With the formula for the electric field, the magnetic field is easy to calculate.
When the scattering distance $R$ is large enough, the simplified expression of the far field is as follows, and the near-field effect can be directly ignored. In other words, in the far field, the propagation characteristics have become stable, so the integral term can be simplified to a term related to the scattering direction $\hat{r}$, that is, $E_s^{\text{far}}(\hat{r})$. $$ E_s(r) = \frac{e^{-jk_0 R}}{R} E_s^{\text{far}}(\hat{r}), \quad H_s(r) = \frac{e^{- jk_0 R}}{R} H_s^{\text{far}}(\hat{r}) \tag{12} $$ In addition, if the contribution of each point is calculated directly, the time complexity is $O(MN)$, that is, the number of discrete points is $M$ $*$ and the number of scattering directions is $N$. Here, the octree is introduced to reduce the number of discrete points on the calculation surface through space partitioning, and the time complexity is reduced to $O(M+log(M)N)$. The specific implementation method is in 4.3.
4.3 Multilevel fast Physical Optics
Accelerate far-field scattering calculations using Multilevel Fast Physical Optics (MFPO).
The multi-layer fast multipole algorithm (MLFMA) proposed by Chew [2001] is used to accelerate the solution of electromagnetic field scattering problems. This algorithm first constructs an octree structure for the surface of hair, for example, where each leaf node represents a sampling point. After the Octree is constructed, the surface current and magnetic current of each leaf node are calculated. Then, starting from the leaf node, the data are accumulated layer by layer. The original $O(MN)$ is reduced to $O(M + \log(M)N)$.
The author then introduces the three key parts of the octree acceleration algorithm. The first is the far-field scattering kernel. $$ e^{jk_0 r{\prime} \cdot \hat{r}} = e^{jk_0 (r{\prime} – c_L) \cdot \hat{r}} \prod_{i=1}^{L} e^{jk_0 (c_i – c_{i-1}) \cdot \hat{r}} \tag{13} $$ The ultimate goal is to calculate the scattering contribution of the electromagnetic wave at each point on the surface of the object to a far-field observation area (distance $R$, direction $\hat{r}$). Specifically, the scattered electric field $E_s(r)$ and magnetic field $H_s(r)$ are calculated, that is, the contribution of the electromagnetic wave emitted from each point $r{\prime}$ on the hair surface in the far field. The left side of the formula is the phase change from a surface point $r{\prime}$ to the far-field observation direction $\hat{r}$. The final time complexity is $O(MN)$. The author divides the surface into different regions, each of which specifies a reference center point. The $c_0, c_1, …, c_L$ in the formula are the node centers of different levels in the octree.
This alone cannot reduce the amount of calculation. Therefore, it is necessary to merge the contributions of surface points that are closer to each other through the high-level nodes of the octree because their phase changes are very small.
Let's use Professor Yan's picture to explain. For a certain area $\delta \mathbf{r} $ in the far field, the contribution of all sampling points to the reference point $\mathbf{r}$ will be accurately calculated first. However, other points in the area are approximated using the parent node of the octree. In other words, the nearby $\mathbf{r}_1$ and $\mathbf{r}_2$ no longer need to consider so many sampling points, but directly use the total contribution obtained by the octree.
The author defines a direction set, and each parent node of the octree stores the cumulative contribution data about different direction sets. Therefore, the parent node contains not only spatial information, but also cumulative information in multiple scattering directions. Finally, a complete 360-degree scattering field distribution is obtained at the root node.
Next, starting from the second level of the tree, merge upwards in sequence. The specific method of upsampling is to perform forward FFT on the scattering contribution of the child node in the direction, then expand the frequency domain data with zero-padding, and finally convert it to the spatial domain through inverse FFT. Ultimately, the parent node can obtain more accurate scattering information in different directions.
Performance
Octrees have a good acceleration effect. Triple trees have the best effect. The more complex the fiber, the better the optimization effect.
Fiber microgeometry and scattering patterns
The cross section of a hair fiber is generally not a perfect circle, but an ellipse. The geometric parameters of the fiber are defined by the major radius $r_1$ and the minor radius $r_2$ of the ellipse. In order to simulate the microscopic roughness of the fiber surface, the author superimposed aGaussian random height field, simulating the real fiber surface. Furthermore,Cuticle tilt, the simulated flakes are arranged obliquely on the fiber surface.
By comparing with traditional light-based hair models, wave optics simulations found that in addition to the significant forward-scattering phenomenon predicted by XIA[2020], complex wavelength-dependent granular patterns were observed, which generated rich color effects when converted to RGB colors.
The study pointed out some regularities observed from the simulations:
Regardless of their actual position, fibers with the same geometric parameters (such as radius, roughness, cuticle tilt, etc.) will generateSimilar grain patterns.
If the fibers have different geometrical parameters, they will generateParticle patterns with different statistical properties, that is, the scattering patterns are obviously different.
The position of the scattered spots depends on the incident angle of the light, and the direction of the deviation followshalf vectordirection.
The size of the speckle pattern increases with the wavelength of the light, which is consistent with the results of Goodman [2007].
5. A PRACTICAL FIBER SCATTERING MODEL
A PRACTICAL FIBER SCATTERING MODELDesigned to account for microscopic geometric variations and complex scattering behavior of fibers.
Previous studies generally used Lut to store the scattering distribution function. This method consumes a lot of space and requires simultaneous recording of the longitudinal and azimuthal scattering distributions.
Now, the authors propose aWavelet-based noise representationThe compact fiber scattering model can represent more geometric complexity and achieve better scattering effect.
In short, the authors want to compact the statistical speckle phenomenon so that it can beMean,variance,Autocorrelation function(ACF) and other statistics. Therefore, the author usedSpeckle TheoryTo describe the patterns produced by the random interference of light.
5.1 Speckle statistics
Here the author mentionedFully Developed SpeckleThis concept. The following is my personal understanding. When light hits a rough surface (such as a fiber/hair surface), each small area on the surface will scatter the light. Due to the tiny irregularities of the surface, the scattered light will interfere with each other, producing a complex light intensity distribution. This distribution is manifested as a series ofBright spots and dark spots, which we callSpeckleThe tiny features of the surface (such as roughness) become sufficiently irregular across the illuminated area (relative to the wavelength of the light) that the phase and intensity of the scattered light at each point is random, which is calledFully developed speckle.
At this time, the fully developed speckle can be described by the complex Gaussian distribution of Goodman [2007]. That is, the real part $\mathcal{R}$ and the imaginary part $\mathcal{I}$ of the electromagnetic field obey the complex Gaussian distribution in space. $$ p_{\mathcal{R},\mathcal{I}}(\mathcal{R}, \mathcal{I}) = \frac{1}{2 \pi \sigma^2} \exp\left( – \frac{\mathcal{R}^2 + \mathcal{I}^2}{2\sigma^2} \right) \tag{14} $$ The real and imaginary parts of the field are independent and normally distributed, with zero mean and the same variance.
The electromagnetic field intensity $I$ and the probability density function of the intensity distribution obeying the exponential distribution: $$ I = \mathcal{R}^2 + \mathcal{I}^2 \ p_I(I) = \frac{1}{2\sigma^2} \exp\left( -\frac{I}{2\sigma ^2} \right) \tag{15} $$ There is nothing much to say about these formulas. In short, the speckle field is very random!
The author makes the light intensity follow an exponential distribution to ensure the statistical characteristics in a single direction. Secondly, by studying the statistical relationship of the light intensity between two points, the ensemble average of the two points is measured. Here, the autocorrelation function (ACF) is used. $$ C(I_{p_1}, I_{p_2}) = \frac{\overline{(I_{p_1} – \overline{I_{p_1}})(I_{p_2} – \overline{I_{p_2}})}}{\sigma(I_{p_1}) \sigma(I_{p_2})} \tag{16} $$ The value of the autocorrelation function is between -1 and 1.When the value is close to 1, indicating that the scattering behaviors of the two light intensities are very similar. This is the key to efficiently reproducing the particle structure in the fiber scattering field.
5.2 Wavelet noise representation of the speckles
The author introduces wavelet noise to represent the noise component of speckle $f(\theta_h, \phi_h, \lambda)$. The specific formula is as follows: $$ f(\mathbf{x}) = \sum_{b=0}^{n-1} w_b(\mathbf{x}) I\left(2^b g_{\lambda}(\mathbf{x})\ right) \tag{17} $$ The core idea of the formula is to decompose the light intensity of the speckle into noises of different frequency levels and perform weighted combination.
By adjusting the weights of different frequency bands, the generated autocorrelation function $ C_f(\mathbf{x}_1, \mathbf{x}_2)$ is close to the final noise with the target autocorrelation function.
According to the Wiener-Khinchin theorem, the autocorrelation function can be expressed byFourier TransformTo calculate. The specific formula is as follows:
$$ C_f(\mathbf{x}_1, \mathbf{x}2) = \mathcal{F} \left( \mathcal{F}^{-1} \left( \sum{b=0}^{n-1 } w_b I_b \right)^2 \right) \tag{18} $$ This means that we can obtain the autocorrelation function by calculating the Fourier transform of the wavelet noise. The autocorrelation function and the power spectral density function are a pair of Fourier transforms, which is amazing.
Here the original paper gives a proof of approximating the ACF by weighted summation of frequency bands.
Traffic saving: By weighting each frequency band of the noise, the autocorrelation function of the entire noise can be approximated without dealing with the interactions between different frequency bands.
The non-negative weights $v_b$ are found by the least squares method so that the weighted sum of the autocorrelation functions $C_b(\mathbf{x}1, \mathbf{x}_2)$ of each frequency band can approach the target autocorrelation function $C_t(\mathbf{x}_1, \mathbf{x}_2)$ .
$$ C_t(\mathbf{x}_1, \mathbf{x}2) \approx \sum{b=0}^{n-1} v_b C_b(\mathbf{x}_1, \mathbf{x}_2) \tag{19} $$ After calculating the weights, we need to ensure energy conservation. That is, we adjust the expected value of the noise function $f(\mathbf{x})$ to $\mathbb{E}[f(\mathbf{x})] = 1$ . $$ \begin{aligned} & \mathrm{E}\left[S{\text {avg }}\left(\theta_i, \theta_r, \phi_r, \phi_r, \lambda\right) f\left(\theta_h, \phi_h, \lambda\right)\right] \\\ & =S_{\text {avg }}\left(\theta_i, \theta_r, \phi_r, \phi_r, \lambda\right) \mathrm{E}[f(\mathbf{x})] \ & \approx S_{\text {avg }}\left(\theta_i, \theta_r, \phi_r, \phi_r, \lambda\right) \end{aligned} \tag{20} $$
Although the effect is good, it has limitations. For example, at grazing incidence angles (when the light is almost parallel to the surface), the fitting accuracy decreases, resulting in a decrease in scattering accuracy. Degeneration in the Forward Direction problem, when the light is in the forward direction (the light direction is consistent with the surface normal), the half-vector direction will degenerate.
6. Validation
6.1 Wave simulation validation
The scattered intensity in the xy plane was calculated and calculated through 3600 azimuth angles $\phi_r$, and finally these angles were averaged into 360 directions.
First, let's compare it to Mie scattering. It is not as good as Mie scattering at grazing angles. The PO approximation is more accurate when the radius of the object is large relative to the wavelength. It is inaccurate when the curvature of the object is small.
Then compare BEM. To simulate the wave scattering of a small ellipsoid, BEM took three hours, while PO took two seconds.
Finally, let’s compare 2D BEM. 1D Gaussian height fields are wrapped around circular and elliptical cross sections. PO wins hands down.
6.2 Measurement
A HeNe laser with a wavelength of 633 nm was used, and the beam spot size was 0.7mm (along the length of the hair) × 3mm (perpendicular to the hair direction)The laser is shone through a tiny hole onto a small area of the human hair sample.
In short, the effect is good!
6.3 Noise representation validation
In one word, good!
7. Rendering
Integrated into PBRT-v3. The original ray-based model is denoted as $S_{\text{ray}}$ , while the diffraction model is denoted as $S_{\text{diffract}}$ .
For diffraction, it is approximated here as single-slit diffraction, where the width of the slit is equal to the diameter of the cylinder (fiber). $$ f_{\text{diffract}}(\theta_i, \phi_d, a) = a \cos \theta_i \cdot \text{sinc}^2(a \cos \theta_i \sin \phi_d) \tag{22} $$ This diffraction model is used to combine with the longitudinal function to obtain a completeBidirectional scattering distribution function (BCSDF)Precompute the diffraction factors in the form of a table of size $50 \times 50 \times 200$ and use a table of the same size for importance sampling.
The extinction cross section can be simply understood as the "effective area" over which light interacts with an object. The extinction cross section of a fiber will be larger than its actual geometric cross section. The extinction cross section will be close to twice the geometric cross section. Light is both reflected and diffracted in a fiber, so we need to divide the total energy between these two phenomena. As a rule of thumb, half of the energy can be used for diffraction and the other half for reflection. $$ S_{\text{diffract}}(\theta_i, \phi_i, \theta_r, \phi_r, \lambda) = \frac{1}{2} \left[S_{\text{ray}}(\theta_i, \phi_i , \theta_r, \phi_r, \lambda) + f_{\text{diffract}}(\theta_i, \phi_d, D/\lambda)\right] \tag{23} $$ Through importance sampling, reflection and diffraction of light are fairly considered.
Next, the final rendering formula describing the light scattering phenomenon is in a nice form: $$ S_{0,\text{prac}}(\theta_i, \phi_i, \theta_r, \phi_r, \lambda) = S_{0,\text{avg}}(\theta_i, \phi_i, \theta_r, \phi_r, \lambda) f(\theta_h, \phi_h, \lambda)\ \\ = \frac{1}{2} \left[ S_{\text{ray}}(\theta_i, \phi_i, \theta_r, \phi_r, \lambda) + f_{\text{diffract}}(\theta_i, \ phi_d, D/\lambda) \right] f(\theta_h, \phi_h, \lambda)\ \tag{24} $$ Reflection + Diffraction + NoiseThe noise function $f(\theta_h, \phi_h, \lambda)$ is at the end, so that the scattered light is no longer concentrated in a few directions.
Focus on this part: $$ S_{0,\text{prac}}(\theta_i, \phi_i, \theta_r, \phi_r, \lambda) = S_{0,\text{avg}}(\theta_i, \phi_i, \theta_r, \phi_r, \lambda) f(\theta_h, \phi_h, \lambda) $$ $S_{0,\text{avg}}$ represents the average reflection/refraction and diffraction behavior of the hair surface, which is calculated by preprocessing (BSDF table).
Next, let’s talk about how to extract BCSDF from PO.
The scattered electric and magnetic fields ($E_s^{\text{far}}$ and $H_s^{\text{far}}$) are obtained through Maxwell's equations.
Then the Poynting vector is used to calculate the energy flow, which is equivalent to calculating the intensity of light.
In order to use the simulation results for rendering, the scattered intensity is related to the incident power. The scattered intensity is calculated by the formula, where $R^2$ is the far field spherical area:
The scattered power $P_s$ and the absorbed power $P_a$ are calculated by integrating the intensities of the scattered light and the absorbed light, respectively.
The 1D dimension describes the wavelength. 2D describes the direction of incident light. 2D describes the direction of the outgoing light.
The storage requirements can be reduced through the memory effect of scattering.
When rendering, for each incident light direction, query the adjacent table data and apply the corresponding angle offset (angle offset in memory effect).
The final table size used is $25 \times 32 \times 72 \times 180 \times 360$ , which is very large, requiring about 15GB of memory.
8. Result
Compared with XIA[2020], the new model is more vivid. The colorful glints are better. It reflects the glints produced by multiple wavelengths on the hair.
9. DISCUSSION AND CONCLUSION
The first practically applicable 3D wave optics fiber scattering model.
Capable of generating speckle patterns commonly seen in optical scattering.
Although the 3D simulation is highly accurate, its memory usage and computational cost are very large. A noise utility model is proposed to reduce the computational cost by capturing the statistical properties of the fiber scattering spot (such as the autocorrelation function).
Currently, this model is mainly used to simulate the first-order reflection mode. In the future, it can be extended to higher-order scattering modes.
Further development of new techniques may be needed to predict or measure the statistical properties of higher-order modes.
EN 
CN
Leave a Reply