The notebook is divided into three parts: elementary, intermediate and advanced.
The beginner's chapter mainly introduces one-dimensional standard and non-standard FDTD theories.
The intermediate chapter will introduce the two-dimensional FDTD theory. The two-dimensional NS-FDTD theory improves the accuracy of solving Maxwell's equations by combining different finite difference models, which is one of the core of this method.
The advanced version will introduce the three-dimensional FDTD theory and apply it to conductive media (such as metals, plasmas, etc.) to study the propagation characteristics of electromagnetic waves in these materials.
Advantages of the NS-FDTD algorithm
Finite-difference time-domain (FDTD) is one of the most famous numerical algorithms for calculating electromagnetic wave propagation. The FDTD method can simulate structures of arbitrary shapes, nonlinear media, and can calculate the propagation of broadband electromagnetic waves. NS-FDTD reduces computing resources by optimizing the calculation of single-frequency waves, allowing higher-precision structures to be calculated with the same resource consumption. Using the NS-FDTD method, researchers have successfully and accurately simulated a special type of electromagnetic mode - Whispering Gallery Modes (WGM).
The whispering gallery mode refers to the phenomenon that electromagnetic waves or sound waves propagate near the inner wall of a circular, spherical or ring-shaped structure and circle around it many times without being easily scattered to the outside.
An intuitive example is that under the circular dome of St. Paul's Cathedral in London, if one person whispers softly on one side, another person on the other side, tens of meters away, can still hear it. This is because the sound waves propagate along the circular wall and stay on a specific path, thus reducing energy loss.
The traditional FDTD method often has large errors when calculating WGM, but the NS-FDTD method has higher accuracy, so the calculation results are closer to the theoretical values, as shown in the figure above.
(a) Figure: Analytical solution of Mie theory, used as a reference standard.
(b) Figure: Simulation results using the traditional FDTD method on a coarse grid deviate greatly from the theoretical values.
(c) Figure: Simulation using the NS-FDTD method on the same coarse grid. The results are highly consistent with the Mie theory and are significantly better than the traditional FDTD calculation results.
Beginner's Edition
First, we introduce the one-dimensional standard/non-standard FDTD theory. In computer simulation, additional considerations are required for numerical stability and boundary conditions.
1. Finite Difference Model
Many partial differential equations (PDEs) do not have analytical solutions, so they can only rely on numerical simulation. Finite difference method (FDM) is the most common numerical calculation method. The basic idea is to use differential expressions to approximate the solution of differential equations.
1.1 Forward Finite Difference (FFD)
To find the derivative of a one-dimensional function, first use Taylor Series Expansion to expand the function $f(x)$ at $x+\Delta x$: $$ f(x + \Delta x) = f(x) + \Delta x \frac{df(x)}{dx} + \frac{\Delta x^2}{2!} \frac{d^2 f(x)}{dx^2} + \cdots , \tag{1} $$ When $\Delta x$ is small enough, we can ignore the higher-order terms (i.e. $ \frac{\Delta x^2}{2!} \frac{d^2 f(x)}{dx^2} + \cdots $ ), and only keep the first-order derivative terms. Finally, we can get: $$ \frac{df(x)}{dx} \approx \frac{f(x+\Delta x) – f(x)}{\Delta x} . \tag{2} $$ Formula (2) is called the forward finite difference (FFD) approximation.
The term $\frac{\Delta x^2}{2!} \frac{d^2 f(x)}{dx^2} + \cdots$ is ignored in the Taylor expansion, so the truncation error is a first-order error, that is, $O(\Delta x) $ .
1.2 Backward Finite Difference (BFD)
Similarly, the expansion of the function $f(x)$ at $x – \Delta x$ is: $$ f(x – \Delta x) = f(x) – \Delta x \frac{df(x)}{dx} + \frac{\Delta x^2}{2!} \frac{d^2 f(x)}{dx^2} + \cdots . \tag{3} $$ Easy to get: $$ \frac{df(x)}{dx} \approx \frac{f(x) – f(x- \Delta x)}{\Delta x} . \tag{4} $$ The error is the same as above.
1.3 Central Finite Difference (CFD)
The previous two algorithms only use the current point and one adjacent point for calculation, so the accuracy is low. The center difference uses the two adjacent points before and after to improve the accuracy.
For the function $f(x)$, subtract formula (1) from formula (3), eliminate the second-order derivative term, and obtain: $$ f(x + \Delta x) – f(x – \Delta x) = 2\Delta x \frac{df}{dx} + O(\Delta x^3) . \tag{5} $$ Further sorting gives: $$ \frac{df}{dx} \approx \frac{f(x + \Delta x) – f(x – \Delta x)}{2\Delta x} . \tag{6} $$ Some textbooks use $ \Delta x/2 $ to replace $ \Delta x $ to obtain formula (7). The two are actually equivalent. $$ \frac{df}{dx} \approx \frac{f(x + \Delta x/2) – f(x – \Delta x/2)}{\Delta x} . \tag{7} $$ Both of the above formulas are called second-orderCentral finite differenceformula.
In addition, for the function $f(x)$, adding formula (1) to formula (3) and eliminating the first-order derivative term, we can obtain: $$ f(x+ \Delta) + f(x- \Delta) = 2f(x) + \Delta x^2 \frac{d^2 f}{dx^2} + O(\Delta x^4). \tag{8} $$ After sorting, we get: $$ \frac{d^2 f}{dx^2} \approx \frac{f(x + \Delta x) – 2f(x) + f(x – \Delta x)}{\Delta x^2} \tag{9} $$
1.4 Higher-Order Finite Difference
In the finite difference method, the calculation accuracy is improved by increasing the number of sampling points to obtain a higher-order difference formula.
In order to improve the accuracy, two additional sampling points $x + 2\Delta x, x – 2\Delta x$ are added on the basis of the second-order central finite difference method, and Taylor expansion is performed here:
For the right-hand side points, the Taylor expansion is: $$ f(x+2\Delta x) = f(x) + 2\Delta x\frac{df(x)}{dx} + 2\Delta x^2 \frac{d^2f(x)}{dx^2}+O(\Delta x^3) , \tag{10} $$ Similarly, $$ f(x – 2\Delta x) = f(x) – 2\Delta x \frac{df(x)}{dx} + 2\Delta x^2 \frac{d^2 f(x)}{dx^2} + O(\Delta x^3) . \tag{11} $$ By linear combination (for example, inverting formula (11) and halving it, and adding the two formulas together), we can eliminate the high-order error terms and obtain: $$ \frac{d^2 f(x)}{dx^2} \approx \frac{1}{\Delta x^2} \left[ \frac{4}{3} (f(x + \Delta x) + f(x – \Delta x)) – \frac{1}{12} (f(x + 2\Delta x) + f(x – 2\Delta x)) – \frac{5}{2} f(x) + \cdots \right] . \tag{12} $$ This model uses four points and reduces the error by weighted averaging, but it brings potential problems such as numerical instability.
When we replace differential equations with higher-order finite differences, we introduce additional numerical solutions, called spurious solutions, which do not actually exist but are caused by the discretization properties of higher-order difference equations.
For example, a first-order differential equation such as $ \frac{df}{dx} = f(x) $ usually has one independent solution, a second-order differential equation such as the wave equation $\frac{d^2f}{dx^2}=k^2f(x)$ usually has two independent solutions (two free parameters), a fourth-order differential equation will have four independent solutions, and so on.
In numerical calculations, if a fourth-order finite difference method is used to approximate a second-order differential equation, the end of the difference equation will be forcibly raised, resulting in additional spurious solutions. These additional solutions do not correspond to the original physical system.
In particular, the second-order differential wave equation describing electromagnetic wave propagation should be solved using the second-order central finite difference method rather than higher-order methods.
2. One-dimensional wave equation
To simulate the propagation of waves on a computer, numerical methods are needed for approximate calculations, among which the finite difference time domain method plays a powerful role.
The mathematical expression of the one-dimensional wave equation is: $$ \frac{\partial^2 \psi}{\partial t^2} = v^2 \frac{\partial^2 \psi}{\partial x^2} , \tag{1} $$ Where $\psi(x,t)$ represents the amplitude of the wave, $v$ is the speed of wave propagation (light in a vacuum is $3 \times 10^8$ m/s).
The physical meaning of this equation is that the change of waves is related to the change of time and space.
For an infinitely long wave medium, the solution is usually a traveling wave of a very simple form: $$ \psi(x,t) = A \cos(kx – \omega t) + B \sin(kx – \omega t) . \tag{2} $$ This situation applies to electromagnetic waves propagating in free space. For example, if the boundary is fixed, the traveling wave can be expanded by a sine function: $$ \psi(x,t) = \sum_{n} A_n \sin(k_n x) \cos(\omega_n t). \tag{3} $$ However, if the boundary is irregular, such as electromagnetic waves reflecting on the ground, water waves hitting the coastline, sound waves propagating in multiple rooms, etc., it cannot be simply expanded into a sine or exponential form, and the analytical solution will be extremely complicated. For example, sound waves will have different reflection paths when they contact walls of different media, making it impossible to directly solve the propagation path of the sound waves.
Analytical solutions usually assume that the wave speed $v$ is constant, that is, the wave propagates in a homogeneous medium. But in reality, waves travel through inhomogeneous media, such as sound waves that have different propagation speeds in rooms of different temperatures. In these cases, the wave equation becomes a variable coefficient partial differential equation, such as: $$ \frac{\partial^2 \psi}{\partial t^2} = v(x)^2 \frac{\partial^2 \psi}{\partial x^2} , \tag{4} $$ The speed $v$ will change with the position $x$, which makes it impossible to directly use Fourier transform to obtain the analytical value.
Or there is an excitation source on the right side of the wave equation, that is, multiple waveforms may have different sources interfering with each other. In this case, it is also difficult to find an analytical solution.
And in the real world, the propagation of energy will cause losses, so it is necessary to introduce loss terms into the wave equation. At this time, the wave equation will become a nonlinear equation or a high-order differential equation, making it impossible to solve directly.
Therefore, we introduce the finite difference time domain method to solve these problems.
2.1 Standard FDTD algorithm
Since computers cannot directly calculate continuous mathematical equations, they needDiscretization, which divides space and time into grids, and then lets the computer gradually calculate the fluctuations of each grid.
In continuous form, the one-dimensional wave equation $$ \frac{\partial^2 \psi}{\partial t^2} = v^2 \frac{\partial^2 \psi}{\partial x^2} , \tag{1*} $$ The equivalent form is $$ \left( \frac{\partial^2}{\partial t^2} – v^2 \frac{\partial^2}{\partial x^2} \right) \psi(x,t) = 0. \tag{5} $$ Therefore, using the interval $\Delta x$ to discretize space, each point is represented by the index $i$$x = i\Delta x$; using the interval $\Delta t$ to discretize time, each moment is represented by the miniature $n$$t = n\Delta t$, where $n$ are all integers. So the wave function is written as: $$ \psi_i^n = \psi(i\Delta x, n\Delta t) . \tag{6} $$ The derivative of the wave function is then calculated using the central finite difference method.
The central finite difference approximation of the second-order partial derivative functions in time and space directions is $$ \frac{\partial^2 \psi}{\partial t^2} \approx \frac{\psi_i^{n+1} – 2\psi_i^n + \psi_i^{n-1}}{\Delta t^2} , \tag{7} $$
where $\psi_{i+1}^{n}$ is the fluctuation value of the adjacent grid point on the right.
Substituting formula (7) and formula (8) into the wave equation (formula (1)), we obtain formula (9). $$ \frac{\psi_i^{n+1} – 2\psi_i^n + \psi_i^{n-1}}{\Delta t^2} = v^2 \frac{\psi_{i+1}^{n} – 2\psi_i^n + \psi_{i-1}^{n}}{\Delta x^2} , \tag{9} $$ After sorting out, we get the standard 1D FDTD calculation formula.1D standard finite difference time domain (FDTD) algorithm: $$ \psi_i^{n+1} = 2\psi_i^n – \psi_i^{n-1} + \left( \frac{v^2 \Delta t^2}{\Delta x^2} \right) (\psi_{i+1}^{n} – 2\psi_i^n + \psi_{i-1}^{n}) . \tag{10} $$ This formula depends on the calculation of the two step points before and after in time and space. In other words, the state of the current point is affected by the adjacent points.
2.2 Non-standard FDTD algorithm
In the standard FDTD method, the wave equation is discretized using central differences. However, this method hasNumerical dispersionProblem: There will be large errors when the grid is coarse.
The non-standard FDTD algorithm solves the above problems. NS-FDTD makes special treatment for monochromatic waves, ensuring accuracy even if the values are discrete.
First, suppose that we are studying a certain monochromatic wave $$ \psi(x,t) \;=\; e^{\,i(kx – \omega t)}, \tag{11} $$ where $k$ is the wave number ($k = 2\pi / \lambda$) and $\omega$ is the angular frequency ($\omega = 2\pi f$).
The spatial differential of formula (11) in the continuous case is $$ \frac{\partial \psi}{\partial x} \;=\; i\,k\, \psi(x,t). \tag{12} $$ For example, if we use standard differences directly to calculate $\Delta_x \psi(x,t)$, we get $$ \Delta_x \psi(x,t) = \psi(x+\Delta x,t) – \psi(x,t) = e^{\,i(kx – \omega t)}\bigl(e^{\,i\,k\,\Delta x} – 1\bigr) , \tag{13} $$ It is not consistent with the real differential. It can only be approximated when $\Delta x$ is small enough. In other words, if the grid is coarse, the so-called error will occur.
NS-FDTD introduces the correction factor $s(\Delta x)$ to minimize the error of the discrete operator $$ \Delta_x \psi(x,t) \approx s(\Delta x)\,\frac{\partial \psi(x,t)}{\partial x}. \tag{14} $$ This $s(\Delta x)$ is determined by the phase change at $\Delta x$. For example, our goal is $$ \Delta_x \psi(x,t) = \psi(x+\Delta x,t) – \psi(x,t) . \tag{15} $$ Then the error factor is derived from formulas (12), (13), (14), and (15): $$ s(\Delta x)
=\frac{\Delta_x \psi(x,t)}{\partial_x \psi(x,t)}
\frac{e^{\,i\,k\,\Delta x} – 1}{i\,k\,\Delta x}
\times \Delta x
\frac{e^{ik\Delta x} – 1}{ik}. \tag{16} $$ By using Euler's formula, the error factor is simplified to $$ s(\Delta x) = \frac{2}{k}\,\sin!\Bigl(\frac{k\,\Delta x}{2}\Bigr)\,e^{\,i\,\frac{k\,\Delta x}{2}} . \tag{17} $$ Note that when $\Delta x \to 0$, the nonstandard difference degenerates back to the case of "standard central difference", which is almost consistent with the true derivative approximation. In fact, formula (17) contains two parts: phase factor and amplitude, the former corresponds to the exponential part.
Similarly, the correction function in the time direction is derived $$ s(\Delta t) = \frac{2}{\omega}\,\sin\Bigl(\frac{\omega\,\Delta t}{2}\Bigr)e^{-\,i\omega\,\Delta t/2}. \tag{18} $$ To summarize, we start with the wave equation, discretize space and time, approximate it with the central difference, and then introduce the correction function (the exponential term of the correction function is actually handled implicitly) $$ \frac{\psi_i^{n+1} – 2\psi_i^n + \psi_i^{n-1}}{\left(\frac{2}{\omega} \sin(\omega \Delta t / 2)\right)^2} = v^2 \frac{\psi_{i+1}^{n} – 2\psi_i^n + \psi_{i-1}^{n}}{\left(\frac{2}{k} \sin(k \Delta x / 2)\right)^2}. \tag{19} $$ After sorting, we get $$ \psi_i^{n+1} = -\psi_i^{n-1} + \left(2 + u_{\text{NS}}^2 d_x^2 \right) \psi_i^n , \tag{20} $$ in $$ u_{\text{NS}} = \frac{\sin(\omega \Delta t / 2)}{\sin(k \Delta x / 2)} . \tag{21} $$
3. One-dimensional FDTD stability
When doing numerical simulation, we hope that the time step $\Delta t$ is as large as possible to reduce the total number of calculation steps, but it cannot exceed a certain limit, otherwise the numerical solution will diverge. How is this limit derived?
Hours (Jikan): When it comes to hours, add "かん" after "hour" to indicate the time period, such as "いちじかん".
Minutes (ふん/ぷん）: The pronunciation of "minute" is "ふん" or "ぷん", which varies depending on the number. For example, "いっぷん" (1 minute), "にふん" (2 minutes).
The verbs in Japanese areRespectandSimplifiedThere are two forms. In verb sentences, the honorific form is mainly used in formal or polite situations, while the simplified form is used in casual or daily conversation.
Comparison table of four quadrants of the formal and simplified forms of verbs
In Japanese, "give" needs to be distinguished according to different subjects. This series of verbs is used to express the relationship of giving and receiving objects, and is called "giving and receiving verbs."
Note: When using "あげる", you need to consider the speaker's psychological tendency. If there is a significant difference in psychological distance, you may not be able to use "あげる".
In addition to "あげる", "やる" also means "to give", and its usage is basically the same as "あげる".
Get it from others.1From N2Item N is obtained there.
N1(Receiver/わたし) + は + N2(giver) + に/から+ N (item) + を+ もらう, the subject is the recipient, and 「から」 is used to emphasize the source. 「もらう」 means to get.
When the verb isOther verbsWhen expressing a wish, the object is usually indicated by the particle "を". However, when using the "たい" form, the particle "が" can be used.
Generally speaking, "Vますshaped＋たいです" is used for the first person. It is more polite to use "V ますshaped + たがっています" to indirectly ask or describe the wishes of another person (usually in the third person).
Word endings are "う, つ, る" → become "って" Word endings are "む, ぶ, ぬ" → become "んで" Words ending with 「く」 → become「いて」 Words ending with 「ぐ」 → become「いで」 Words ending with 「す」 → become「して」
Word endings are "う, つ, る" → become "った" Word endings are "む, ぶ, ぬ" → become "んだ" Words ending with 「く」 → become「いた」 Words ending with 「ぐ」 → become「いだ」 Words ending with 「す」 → become「した」
V1て、V2て、V3Masu.Indicates a series of actions that occur in sequence or continue. Usually the last action uses the "ます" form to indicate the completion of the entire action sequence. Only the last place can reflect the tense.
V (て形) + ください, which means to ask the other person to do something. It can be used in daily conversations, service places, and communication with superiors and elders.
Verb continuation: V る+ predetermined だ, used to express planning to do something. Noun continuation: N + の+ predetermined だ, indicating the scheduled plan for an event or matter.
It is more formal and objective than 「〜つもりだ」 and is often used for arrangements or plans that have already been finalized.
"middle（ちゅう）" is often used after a Chinese word with an action meaning to indicate that the action is in progress. For example, "会会议中（かいぎちゅう）", " "Reluctantly(べんきょうちゅう)"、" in the official career（しごとちゅう）”、”Business(えいぎょうちゅう)"、"Moving（It's so nice）"、"Preparing（じゅんびちゅう）" wait.
Indicates the continuation of a state (static verb), used to describe the continuation of a certain state, that is, the action has been completed, but the state remains unchanged.
Indicates that an action has not yet occurred or is not yet completed.
The honorific form is often used to answer questions such as "もうV ましたか." (Has it already been completed?) when asking whether a certain action has been completed.
Vru (verb basic form/dictionary form) + 前 (に), before doing something. Usually indicates that another action needs to be completed before the main action occurs.
(A third party) wants to do something. / (A third party) wants something.
Third person + V (verb form) + たがる, indicating that the third person "wants to do something".
「-たがる」 is a derived verb form of 「-たい」, which is used exclusively for the third person. It should be noted that in the 「V たがる」 structure, the object of the verb needs to be indicated by 「を」, and cannot be changed to 「が」.
Moreover, the て form of 「たがる」 is commonly used by Japanese people, which is 「たがる→たがっています」.
N (noun) + が見える, meaning "to be able to see..." or "to be able to see...". Describes visually visible situations and is often used to express natural scenery, buildings, specific places, etc.VisibilityUnlike "N を見る" which means to take the initiative to see something, "N が見える" is more of a description.Nature in sightsituation.
V (verb form) + てはいけない, which means prohibiting or not allowing something to be done. It is suitable for conveying discipline, rules and regulations, or orders in the relationship between superiors and subordinates.
V (verb form) + なければならない, emphasizing the necessity or obligation of a certain behavior.
If not + prohibited = have to, obligation, must
It's like this + it's like this Nakute wa + ikenai Nakerezawa + Nara Nai Nakerezha + Ikenai なくちゃ (spoken) + ならない なくちゃ（spoken）＋いけない なきゃ (spoken) + ならない なきゃ（spoken）＋いけない
N＋だけ＋で、に、は、を、が, etc. When "だけ" is followed by "を" or "が", these two particles can usually be omitted. It emphasizes that it is limited to a certain object or range.
da ke de: Indicates that a certain condition alone is sufficient. “Only…”
V (verb form) + ことがある, indicating that one has had a certain experience, equivalent to "once...been...been..." V(Verb form) + ことがない, negative form, indicating that a certain experience has not been had.
Used to describe a person's past experiences, usually referring to having done something before or having had this experience before a certain time.
N (noun) + がする, expresses a sensory experience that is the result of natural perception rather than active perception.
sound（おと）がする (the sound of a certain environment or object), sound（こえ）がする (someone’s voice), 匂い（におい）がする (to give off a certain smell), taste（あじ）がする(having a certain smell), smelly（くさい）がする (smell something strange). . .
Indicates speculation or possibility. The tone carries a certain degree of uncertainty, but is based on certain information or judgment. This sentence pattern can also be used to ask questions, with the meaning of confirmation or inquiry. The simplified form is "だろう".
Verbs and adjectives: Use simplified sentences with 「かな」 Noun: Directly connected to 「かな」
A soliloquy expressing a slight exclamation or uncertainty. It is often used in soliloquy or in situations with speculation, indicating the speaker's uncertainty or hesitation about something.
expressEven in some cases, the result does not change, is often used to express a firm result or state, that is, no matter what the conditions are, the result remains the same.
『Quoted sentences + と+ 言う/hear く/answer える, used to quote the other party's original words, with a strong sense of objective description. Simplified sentence + と+ 言う/hear く/answer える, converting the original words into simplified sentences with a more natural tone, suitable for use in daily conversations.
Used to quote someone else's words or express indirect quotation, indicating information heard or obtained.
Simplified Chinese Sentences + と思う, used to express one's own thoughts, opinions or speculations. The particle "と" is used to indicate the content of the thought and to indicate the object of "思う". It is mainly used in the first person, that is, to express one's own subjective judgment or opinion. Generally, direct quotations are not used, so the content before "と" should be expressed in simplified form. Regardless of the tone of the sentence, the content before "と" must be expressed in simplified form.
Verb simplified form + かもしれない いAdjective simplified form + かもしれない na-adjective stemda+ Kamoshirena Noun N (simplified form removes "だ") + かもしれない
This tone indicates that the speaker is not completely sure about a situation, but thinks it is possible. It is more uncertain than 「〜でしょう」 and carries a sense of speculation.
N + しか + V (negative form) means "only..." or "merely...", emphasizing the small quantity or limited scope.Although the verb is in negative form, the meaning of the sentence is positive., emphasizing small quantity or limited scope.
More than...1···，···N2(also)... Not only..., but also...
N1+ Dakedenaku + N2+ も, indicating that it is not limited to N1, and N2It can be used to describe multiple aspects of objects, people, characteristics, etc. It is also suitable for use in sentences that need to emphasize the expansion of the scope.
N1 (I/giver)はN2 (receiver)にV てあげる, it expresses help to others with a sense of favor, so it should be avoided in formal occasions or to people you are not close to.
N1 (giver)はN2 (me/receiver)にV てくれる, to express that someone else does something for you (or someone on your side), the recipient (usually "I") uses "に" as a prompt, but "私に" can be omitted.
Verb ta-form + ばかりだ, emphasizing that the event has just been completed and that a short period of time has passed, with an implication of freshness or “just ended”.
Time is highly subjective and depends on the speaker's feelings.
V-shaped +: Hope that the other party will do something V ない形+でほしい: Hope that the other party will not do something (N に)V てほしい/ V ないでほしい：Indicate the object of hope
V Simplified + みたいだ. AⅠSimplified Chinese + みたいだ. AⅠⅠSimplified (～da)+ みたいだ. N Simplified (～da)+ みたいだ.
It expresses the speaker's subjective speculation about the current situation. He is not sure whether it is true, but judges a possibility from the appearance or condition.
V (simplified form) + と语っていた AⅠ(いadjective) + と语っていた AⅡ（なadjective stem +da）+ I want to say N(noun +da）+ I want to say
To paraphrase a third person's speech or opinion, it uses the form of "言っていた", which has a retrospective meaning, that is, to paraphrase information that happened in the past.
SOV→（Simplified verb + S）はOをV (BUY ったN) はパソコンです. (DRINK N) It’s good. N→こと, fill in a formal aspect here, used for abstract and objective statements, customary rules and behaviors.
In addition to being used as a general noun to mean "things", "こと" can also be used as a "formal word" to connect the conjoined form of a verb or adjective.Nounization, so that it can serve as the subject, object, etc. of the sentence. In addition to "こと", there are also "の" and "もの" in the formal aspect.
Let me briefly summarize my current understanding.
Computer graphics based on wave optics is undoubtedly a huge challenge for graphics researchers. It not only requires computational processing of electromagnetic waves, but also involves some theories of quantum mechanics. It is not easy to bridge the gap between physics and graphics. In graphics based on wave optics, the propagation of light in a medium is no longer just a straight line, but will exhibit various characteristics such as spin, deflection and diffraction due to different wavelengths. From the optical rotation caused by soap bubbles, polarizers, oil film gloss or sugar media to how radio radiation propagates between urban buildings, all are inseparable from the theory of wave optics.
There are also many experts on Zhihu who have shared very detailed basic teaching on wave optics rendering theory.
However, the reading threshold is high and it is too far from practical application. Wave optics rendering involves too much mathematical and physical knowledge, and even if you read the paper of Professor Yan and his doctoral students line by line, it is difficult to make any progress.
The holy grail of graphics is ray tracing, and global wave optics rendering is the holy grail of graphics holy grails.
At Sig21, Shlomi proposed combining path tracing with physical optics to achieve a realistic simulation of electromagnetic radiation propagation. There are many methods for calculating electromagnetic wave transmission, from high-precision but computationally intensive wave solvers to fast but inaccurate geometric optics methods. The following figure shows the current calculation methods for electromagnetic wave transmission. The left one is the most accurate, and the right one is the fastest.
Wave Solvers focus on finding exact solutions to Maxwell's equations, but are not practical for large scenes, which are usually done with FDTD, BEM or FEM.
PO is based on high-frequency approximation of electromagnetic wave calculations, but it is barely sufficient for the frequency of visible light. PS: [Xia 2023]'s black dog hair belongs to the Physical Optics method. The full-wave reference in this article should still belong to the Wave Solvers method, but it uses some ideas of PO (equivalent current, etc.) on the basis of BEM.
In addition to PO, there is another method between Wave Solvers and Geometrical Optics, called Hybrid GO-PO. I personally think it should be called a hybrid method of geometric optics and physical optics. The Uniform Theory of Diffraction (UTD) incorporates the diffraction effect into geometric optics to calculate the transmission of electromagnetic waves under high-frequency conditions. In my opinion, UTD compensates for the deficiencies of the boundary conditions of geometric optical rays by calculating the diffraction coefficient, which means that the rays of geometric optics can also turn. This operation is very practical in the field of radar detection antenna design. In addition to UTD, Hybrid GO-PO also involves a technology called Shooting and Bouncing Rays (SBR). This technology simulates multiple reflections of rays on the surface of an object, which is also based on geometric optics.
I personally think that understanding light as a plane sine wave is somewhat limited. For example, the 3b1b optics series video considers the electric field of light as a plane sine wave.
Although it can explain most phenomena and is very suitable for introductory learning, for the study of wave optics rendering, simply describing the electric field of light as a plane sine wave cannot further explain, for example, the Gaussian beam mentioned in this paper.
But I still strongly recommend readers who are not familiar with wave optics to watch this series of videos first.
The video shows the special wave optical phenomenon that light "rotates" when passing through a right-handed chiral medium, and popularizes a series of interesting phenomena in wave optics and the principles behind them.
Very thought-provoking.
Finally, it even talks about how to use matter waves to construct holographic images.
In classical electromagnetic field theory, we usually usePlane WaveExpand the electromagnetic field, and the photons of each mode are spatiallyNonlocalWith infinite spatial extension,A photon is in a mode with a well-defined frequency and wave vector., so it "fills" the entire area of the plane wave in space. This is very common for electromagnetic field patterns in free space, but this plane wave is not suitable for describing localization.
Another way to expand isLocal wave packetThis way of understanding allows us to arrange the electromagnetic field in space, that is, to form a series of wave packets with a certain width.
Simply understanding light as a sine wave actually violates the uncertainty principle. If light is regarded as a sine wave, it means that the light is completely monochromatic, but monochromatic light is impossible in reality. The effect of optical signals is basically in the frequency domain, and Fourier transform is required to convert the spatial domain to the frequency domain. The bandwidth-time uncertainty principle states that the narrower the bandwidth, the longer the signal will be in time. Conversely, the wider the bandwidth, the shorter the time length of the signal.
In classical electromagnetic theory, a wave packet can correspond to a region where energy is concentrated.
But please note that photons cannot be simply understood as a wave packet, but as a probability wave.Probability wave functionDescribes the probability of a photon's existence. The more concentrated the wave packet, the more obvious the particle nature. This is the wave-particle duality explained by quantum electrodynamics. A photon is not necessarily in only one wave packet, it can also be described as a superposition of multiple wave packets. Because the photon state is essentially the excitation of the quantum field,Allows superposition of wave packets at different locations.
That is to say,A photon can "span" multiple wave packets, that is, its wave function can exist in the form of multiple wave packets in space, rather than being confined to a specific location. When a wave packet mode has a photon, this wave packet can be regarded asLowest excited state; If there are multiple photons on the wave packet (higher-order excited states), higher energy will be reflected.
An electromagnetic beam mentioned in this articleGaussian beamIt is a specific electromagnetic wave solution and can be regarded as a wave packet.Gaussian beamIt is a structure with stable amplitude and phase.Single wave packet, which shows a Gaussian distribution on the cross section. This is different from the plane wave solution we usually use in classical electromagnetic theory, because the wavefront of the Gaussian beam has a changing curvature and is not an infinitely extended plane.
Ensemble of wavesRefers to a collection of multiple waves that may have different frequencies, phases, or propagation directions. If we consider a system in which multiple independent wave packets (such as multiple pulsed laser beams) are superimposed on each other, these wave packets can be understood as aEnsemble of wavesIn other words, if multiple uncorrelated wave packets, especially those with random phases, are superimposed on each other, they can statistically constitute a wave ensemble.
When we expand the electromagnetic field into a series of wave packets, we can regard each wave packet as a random event, and their arrival time and phase are random variables. For a collection of multiple wave packets, we can observe the characteristics of these wave packets at different times and locations. In a statistical sense, usingEnsemble meanTo analyze the energy distribution and wave behavior of light.
Since I am not a physicist, I do not approve of writing the differential before the integral formula (
As mentioned above, light can be viewed as a collection of multiple wave packets. There may be phase and time differences between different wave packets, so it is necessary to introduceCross-correlation functionandInterference Functionto describe the coherence and relative phase between wave packets at two different locations, providing aSimilarity in time and spaceIf we sayEnsemble meanUsed to describe the average behavior of a signal in a statistical sense, thenCross-correlation functionIt describes the time averageThe correlation between different locations and different timesandVolatility Correlation.
Specifically, the cross-correlation function describes the coherence of the waves at positions $\vec{r}_1$ and $\vec{r}_2$. If the two waves are coherent, then the Gamma value will be large:
The cross spectral density (CSD) matrix $W(\vec{r}_1, \vec{r}_2, \omega)$ is the cross-correlation functionFourier Transform, expressed in the frequency domainThe correlation between two positions.
Remember the autocorrelation function we discussed in Xia2023's article Black Dog Hair? The cross-correlation function requires two signals to describe, while the autocorrelation function only has one signal. From another perspective, the autocorrelation function (ACF) is a special case of the cross-correlation function. The autocorrelation function describesThe correlation between the signal and itself at different times or locations, and the cross-correlation function isCorrelation between two different signals or the same signal at different locationsThe cross-correlation function and the cross-spectral density are a pair of Fourier transform pairs, and the autocorrelation function and the power spectral density are a pair of Fourier transform pairs.
The cross spectral density (CSD) and the mutual coherence function describe the coherence of the fluctuations between different locations.Radiance Cross-Spectral Density (RCSD)It is a generalization of CSD, describing the correlation between light radiation (i.e., energy density) at different positions and directions. It can be understood that RCSD provides a coherence description similar to CSD based on radiation measurement. In the formula, $L(\vec{r}_1, \vec{r}_2, \omega)$ represents the radiation intensity coherence between positions $\vec{r}_1$ and $\vec{r}_2$ at frequency $\omega$.
Radiative cross-spectral density transfer equation (SDTE)With classicLight Transport Equation (LTE)Similar, but more suitable for wave optics. LTE describes the transmission of optical radiance from point to point, while SDTE uses the RCSD function to describe the propagation of optical radiation, which is equivalent to treating the transmission as coherent transmission between regions.
The RCSD in SDTE expresses the propagation of radiation in the form of an integral between regions, which can be understood as replacing the traditional reflection and scattering with the RCSD matrix and diffraction operator. And note that SDTE is based onRCSD function, rather than the specific light intensity value, which is quite different from LTE.
In further discussion of the useBoundary Element Method (BEM)andAdaptive Integral Method (AIM)Before we discuss the accelerated full-wave reference simulator, let's first briefly review the previous article. The previous article introduced the PO method and SDTE/RCSD theory. These methods are used for different scattering calculation needs, but their basic theories and scope of application are different. This article will discuss a method that provides high-precision surface scattering simulation by combining BEM and AIM.
To sum up in the original author's words, Wave optics is a very new branch in physically based rendering. Although wave optics phenomena can be seen everywhere in life, their impact on the picture is not very big. There is actually a lot of room for improvement in this direction.
Yunchen Yu, Mengqi Xia, Bruce Walter, Eric Michielssen, and Steve Marschner. 2023. A Full-Wave Reference Simulator for Computing Surface Reflectance. ACM Trans. Graph. 42, 4, Article 109 (August 2023), 17 pages. https: //doi.org/10.1145/3592414
Speech Report
Compared to the more commonly used ray tracing techniques, wavelight simulations are more time-consuming but also more accurate.
For example, the microscopic scratches on the metal surface and the hair fiber structure cannot be rendered by traditional optical models, and the colorful color effects we observe in real life cannot be rendered.
Wave optics-based rendering is a difficult problem because solving Maxwell's equations requires a lot of complex calculations. Existing wave-based appearance models usually adopt some approximation methods.
One approximation is the scalar field approximation.
In the wave scattering problem, the electric field and the magnetic field are different vector field quantities, and light with different polarizations consists of field quantities pointing in different directions.
Some approximate models replace these two vector field quantities with a single scalar function, thus being able to calculate the intensity of the light energy but giving up modeling polarization.
Another approximation is the first-order approximation. This assumes that the light reflects only once from each part of the model structure, ignoring multiple reflections. However, there are many cases where these approximations are not valid.
For example, Yu et al., in collaboration with Dr. Lawrence’s group at Penn State University, created surfaces with cylindrical cross-sections that cause multiple light reflections and produce structural colors, phenomena that cannot be well understood or predicted using approximate models.
The authors wanted to characterize surface scattering as accurately as possible by calculating the bidirectional reflectance distribution function (BRDF).
Existing models all use various approximations, such as ray-based, scalar, or first-order approximation models. Without a reference quality BRDF, it is difficult to see what each reflection model is missing or what scenarios it is suitable for.
The authors' solution is to build a 3D 4-way simulation to compute the BRDF of surfaces with well-defined microgeometry.
Claims to compute a reference-quality BRDF for surface samples without using any approximations.
In terms of speed, dddd.
The authors next describe how their simulation works and how it relates to the BRDF.
First, in the left picture, we input a surface sample (defined as a height field) and an incident direction (expressed on the projected hemisphere). We define an incident field propagating toward the surface and calculate a scattered field from the surface.
In the middle picture, the beam is the incident field, and the scattered field in the background is also shown in this picture.
The output is the BRDF model for a given incident direction. Each point in the hemisphere represents an outgoing direction, and the color represents the color of the reflected light in that direction. The BRDF is expressed in RGB colors, which are converted from spectral data.
For many surfaces with low roughness, the reflection pattern is symmetric around the specular direction and changes as the incident light direction moves.
The following is a discussion on how to use the boundary element method to solve the problem of signal scattering only on the surface, thereby reducing the dimensionality of the problem.
The surface signal is the surface current solved from Maxwell's equations. After discretization, the problem is constructed as a linear system, and the surface current and scattered field are solved.
To make the computation feasible, the authors symmetrize the linear system and use a minimum residual solver suitable for symmetric matrices.
In addition, the matrix-vector multiplication is accelerated using the adaptive integration method, which is an acceleration method based on the fast Fourier transform and was originally used in radar calculations.
Most of the code uses the Cuda C++ package for acceleration.
Next, some results are shown, illustrating how the BRDF it computes compares to the BRDF derived from previous methods.
[Yan 2018] use scalar field approximation BRDF models that only consider one refraction.
[Xia 2023] This paper uses vector field quantities but only considers one refraction.
The most accurate method is ours, which not only uses vector field quantities, but also takes into account reflections of all orders.
In the above figure, each incident direction corresponds to five BRDFs, representing different calculation methods.
A relatively smooth material covered with a bunch of isotropic bombs.
The first row shows the normal instance, and the reflection pattern is pretty much centered.
In the second row, the pattern shifts to the left because the incident light comes from a certain oblique direction.
Since the surface is not too rough, the results of the five methods are very similar.
Another material has some corners (corner cubes). The three faces of each corner cube will reflect the light multiple times, making the reflected light return in the direction of incidence. This is called retroreflection.
Our simulator can also simulate this phenomenon. The four methods on the left all failed.
The reason is that if only one reflection is considered, when the light hits one face of the corner cubes, it will be predicted to go down into the lower hemisphere.
The final example is a surface covered with spherical pits.
Multiple reflections occur due to the high slopes of the surface at the edges of the pits.
Different methods show obvious differences.
You can see the extra lobe predicting. (The part slightly to the right in the middle)
In addition, it is brighter overall.
These differences arise from interference between reflections of different orders.
Finally, a technique for efficiently computing the BRDF of a very large number of densely sampled directions is briefly introduced.
If the surface to be simulated is large, this will be slow and require a lot of GPU memory.
But the calculation is linear and can decompose a large surface into multiple smaller sub-areas.
The incident field is projected onto each sub-area, the smaller sub-area simulation is performed first, and then the scattered fields are integrated to obtain the BRDF of the entire large-area surface.
By applying different complex value scale factors to the scattered fields of different sub-areas, the BRDF of a large surface corresponding to different incident directions can be synthesized.
This is because applying an appropriate phase shift to the local incident field in each sub-region produces a total incident field with a different net direction on the surface. In this figure, the incident field propagates vertically. If the same field with five foci at different spatial locations is superimposed, a spatially wider field is obtained, still propagating in the vertical direction. If these fields are linearly combined and an appropriate complex-valued scaling factor is applied to the field in each sub-region, the overall field can be made to propagate in a slightly tilted direction.
Here we explain why different complex value scale factors can produce different incident directions.
This factor can adjust the amplitude and phase of the wave. For example, if two stones are thrown into the water at the same time, the waves on the water will be stronger. If one stone is thrown a little later, the ripples may cancel each other out (destructive interference). This factor controls the time of throwing the stones. By adjusting the phase to control the superposition of waves, the waves can be "guided" to propagate in different directions. Search for "Beamforming" in detail. It is widely used in radar, wireless communication, sonar and other fields.
These three pictures represent the wavefront of the light wave, i.e. the wave crest, which can be understood as the waveform cross section of the light during propagation.
In the upper figure, the incident field distribution is concentrated in the center when the light is incident perpendicular to the surface. The light field is concentrated in the center and propagates in the vertical direction (i.e. the direction of the yellow line in the middle).
At the bottom left, the effect of superposition of multiple identical incident fields, but with the same phase (i.e. no phase difference) when superimposed. The addition of multiple incident fields makes the entire field wider in space, but the propagation direction remains vertical.
In the lower right corner, the superposition effect of multiple incident fields is shown, but a phase difference is added during the superposition, which is equivalent to "offsetting" the direction of the incident field, showing an inclined direction.
In the demo, two videos show the movement of BRDF patterns. (I'm too lazy to make GIFs)
Finally, the big guys took a photo with their hands making a heart shape.
paper
mineFirst articleHaving briefly introduced the content and results of this work, we will now move directly to the theoretical derivation (Section 3-5).
3 FULL-WAVE SIMULATION
The holistic approach starts with a surface model, described by a height field and its material properties (such as the complex refractive index), and specifies a target point.
To compute the BRDF for a given incident direction, an incident field is defined that propagates from that specific direction to the target point.
This incident field is used as input and processed through a surface scattering simulation to solve for the corresponding scattered electromagnetic field.
In this section (FULL-WAVE SIMULATION), we will focus on the principles of BEM in application scenarios. The next section will explain how to efficiently implement the BEM algorithm, and the last section will explain how to combine multiple simulation results to synthesize a BRDF that is densely sampled in the incident and scattering directions.
Let me start with a symbol table to scare you.
3.1 Boundary Element Method: The Basics
The boundary element method (BEM) mainly solves the scattering problem of a single frequency, that is, how electromagnetic waves (including electric and magnetic fields) of a specific frequency are reflected and scattered at the boundaries of different media. The boundary here divides the space into two uniform regions, and the material properties of the two regions (medium parameters at the incident site) are represented by ($\epsilon_1, \mu_1$) and ($\epsilon_2, \mu_2$). Among them, $\epsilon$ represents the dielectric constant (dielectric constant), and $\mu$ represents the magnetic permeability (magnetic permeability coefficient).
In this approach, we deal with complex-valued fields that contain both amplitude and phase information (i.e., the propagation state of the wave). To simplify the formula, we assume that all waves are "time-harmonic" - that is, the waves change with a specific period over time. Throughout the text, the $e^{j\omega t}$ term is omitted to simplify the presentation.
3.1.1 Maxwell's Equations and Surface Currents
First, Maxwell's equations describe how electric fields (E) and magnetic fields (H) interact with each other, determining how light waves propagate and scatter between different materials. For simplicity, this is expressed in a "time-harmonic" form:
The left side of the equal sign describes the degree of "rotation" of the electric and magnetic fields in space. $M$ and $J$ areSurface current density(imaginary current), respectively representing the density of magnetic flow and current (electric and magnetic current densities). This formula can be understood as that when the electric field "rotates" near the boundary, changes in magnetic flow and magnetic field occur; the rotation of the magnetic field also causes changes in the electric field and current.
The core idea of the boundary element method is:Introducing surface currents on boundaries, using these currents to indirectly describe the field distribution without having to calculate all points in each region. The three-dimensional problem is reduced to a two-dimensional problem on the boundary.
3.1.2 Source-Field Relationships
How does an imaginary electric current on the surface (the "source" of the electromagnetic waves) produce the scattered electromagnetic field (the "field")?
As shown in Fig. 2., in the region $R_1$, the total field (incident field and scattered field) are represented as $E_1$ and $H_1$ respectively;
In the region $R_2$ , the total field is denoted by $E_2$ and $H_2$ . The upper scattered field is generated by the upper electromagnetic current; the lower scattered field is generated by the lower electromagnetic current.
In a homogeneous medium, Maxwell's equations can be written in integral form to describe how electric and magnetic fields are generated.
The left side of the equal sign represents the electromagnetic field strength at $r$, that is, it describes the "effect" of the field. $\mathcal{L}$ and $\mathcal{K}$ are integral operators that represent how the field is generated from the surface current and magnetization. These two operators are defined as:
This function converts the source field of the scattering surface into the distribution of the electromagnetic field in the scattering area.
The formula (11) in this paper is actually the same as that in [Xia 2023], but this paper implicitly incorporates the Green's function in the operator and the integration domain is wider. In essence, both describe how to generate the scattered electric field $E(r)$ from the current density $\mathbf{J}$ and the magnetic flux density $\mathbf{M}$.
When solving Maxwell's equations, the Green's function is used to integrate the influence of various sources (such as current and charge) on the electromagnetic field in space. Assuming that the electromagnetic field changes with time in the form of $e^{j\omega t}$ (single frequency), a form similar to the Helmholtz equation can be obtained: $(\nabla^2 + k^2) \mathbf{E} = -j \omega \mu \mathbf{J}$ , which is actually a "frequency domain" form of Maxwell's equations. The Green's function is introduced to establish a source-field relationship, that is, the current $\mathbf{J}$ and the charge $\rho$ are used as "sources" to calculate the distribution of the electromagnetic field $\mathbf{E}$ and $\mathbf{H}$. Then the Green’s function satisfies the following equation: $(\nabla^2 + k^2) G(\mathbf{r}, \mathbf{r}{\prime}) = -\delta(\mathbf{r} – \mathbf{r}{\prime})$ , $\delta$ is the Dirac delta function, which represents the standard waveform generated by a “point source” in space, and describes the wave field excited by a “point source” in space. Through the Green’s function, the influence of any current distribution $\mathbf{J}$ on the field point $\mathbf{r}$ in space can be expressed! Then, the “current” or “charge” at each source point is diffused to the entire space through the Green’s function, producing a cumulative effect on each field point. $\mathbf{E}(\mathbf{r}) = \int_V G(\mathbf{r}, \mathbf{r}{\prime}) \mathbf{J}(\mathbf{r}{\prime}) d \mathbf{r}{\prime}$ This formula is the electric field expressed as an integral superposition of source currents. Finally, combine the integral form of Maxwell's equations with the idea of Green's function. For example, the integral form of the electric field can be expressed as: $\mathbf{E}(\mathbf{r}) = -j \omega \mu \int_V G(\mathbf{r}, \mathbf{r}{\prime}) \mathbf{J}(\mathbf{r}{\prime}) d \mathbf{r}{\prime} – \nabla \times \int_V G(\mathbf{r}, \mathbf{r}{\prime}) \mathbf{M}(\mathbf{r}{\prime}) d \mathbf{r}{\prime}$ . Through convolution, we can obtain the cumulative influence of the tiny current distribution at each source point on the field point, which is the result of the propagation and superposition of the integral of the Green’s function in space. The Green function is used as a convolution kernel to propagate the current distribution at each point in space to the target field point through integration, realizing the cumulative impact of the current at each source point in space on the entire field.
The electric field in the dielectric regions $R_1$ and $R_2$ is expressed as formula (6) and (7) respectively, which include two operators.
This gives the manifestations of electric and magnetic fields in different dielectric regions.
In summary, this section transforms the Maxwell equations into integral expressions by generating scattered electromagnetic fields through hypothetical surface currents. Green's function transforms current and charge distribution into integral superposition of electromagnetic fields, showing the specific implementation of the source-field relationship. Finally, the expressions of the electric and magnetic fields in the regions $R_1$ and $R_2$ are given, showing the influence of different medium parameters on the field.
3.1.3 Boundary Conditions
When electromagnetic waves propagate at the boundary of two different media, reflection and refraction will occur. At this time, the energy of the wave cannot disappear out of thin air, but transitions smoothly at the interface. If the electric field or magnetic field is discontinuous at the boundary, an unrealistic energy jump will occur (that is, the energy suddenly disappears or increases), which violates the law of conservation of energy.
Specifically, you can search: "Interface conditions for electromagnetic fields".
Therefore, certain boundary conditions must be met at the boundary of the medium to ensure the continuity of the electromagnetic field and the conservation of energy. Specifically, when an electromagnetic wave propagates at the interface of two different media, the tangential components of the electric and magnetic fields need to maintain continuity at the boundary:
And the net electromagnetic current density on the boundary is zero, that is, the electromagnetic current densities on both sides of the boundary are opposite in direction and equal in magnitude.
These two conditions must be met at the same time to avoid violating the laws of physical conservation.
3.1.4 Integral Equations
Combining the above formulas (6), (7), (8) and (9), we get the integral equations for the electric and magnetic fields, which are called the PMCHWT (Poggio-Miller-Chang-Harrington-Wu-Tsai) equations:
EFIE（Electric Field Integral Equation）, electric field integral equation. $j \omega \mu_1 (\mathcal{L}_1 \mathbf{J})(\mathbf{r})$ and $j \omega \mu_2 (\mathcal{L}_2 \mathbf{J})(\mathbf{r})$ represent the contribution of current density $\mathbf{J}$ to the electric field in medium 1 and medium 2. $\mathcal{K}_1 \mathbf{M}$ and $\mathcal{K}_2 \mathbf{M}$ represent the contribution of magnetic current density $\mathbf{M}$ to the electric field in medium 1 and medium 2.
MFIE(Magnetic Field Integral Equation), magnetic field integral equation.
In general, this is a boundary integral equation specifically used to solve electromagnetic scattering problems caused by dielectric objects. With this PMCHWT equation, the distribution of electromagnetic fields can be accurately calculated.
3.1.5 Solving for Current Densities
In this section, we need to use the PMCHWT equation above to calculate the distribution of "electric current" and "magnetic current" on the surface of the object. These distributions determine how the electromagnetic wave will be "reflected" or "refracted" when it hits the object.
The surface current density $\mathbf{J}$ and magnetic flux density $\mathbf{M}$ are solved by discretizing the boundary element. A basis function $f_m(\mathbf{r})$ is defined for the discrete element, and the current density and magnetic flux density distribution are expressed by the basis function expansion method.
N is the total number of basis functions; $ I_{J_m}$ and $I_{M_n}$ are the unknown coefficients of the corresponding basis functions, representing the current and magnetic flux intensities on each unit.
Through this basis function expansion, the continuous surface current density and magnetic current density are decomposed into a linear combination of a series of basis functions.
In order to solve the unknown coefficients of the corresponding basis functions, the electric field integral equation (EFIE) and the magnetic field integral equation (MFIE) are transformed into a linear equation system. This is done using the Galerkin Method. The basic idea of this method is to apply the integral equation to each basis function and perform weighted averaging so that the integral equation holds true in the projection direction of each basis function. Simply put, this method involves discretization, finding the basis, and calculating the coefficients. A high-dimensional linear equation system can be simplified using linear algebra methods.
In this way, the EFIE part of the original continuous form of the PMCHWT equation can be converted into a finite number of linear equations, and the problem can be transformed into solving the following matrix equation.
$$ V_{\mathrm{E}}^m =\int_S \mathbf{f}_m(\mathbf{r}) \cdot \mathbf{E}_i(\mathbf{r}) d \mathbf{r}\tag{ 17} $$
$$ V_{\mathrm{H}}^m =-\int_S \mathbf{f}_m(\mathbf{r}) \cdot \mathbf{H}_i(\mathbf{r}) d \mathbf{r}\tag {18} $$
In formula (12), $I_J$ and $I_M$ need to be calculated.
In a loose sense, formulas (13)-(16) respectively represent the contribution of each small block's current density to the electric field, the contribution of each small block's magnetic flux density to the electric field, the contribution of each small block's current density to the magnetic field, and the contribution of each small block's magnetic flux density to the magnetic field. Formulas (17)(18) respectively represent the "driving force" of the external incident electric field on the small block's current and the "driving force" of the external incident magnetic field on the small block's magnetic flux. Emphasize the elements in the matrix, such as $A_{EJ}^{mn}$, which is actually a double integral. Since the integration of the source point $\mathbf{r}{\prime}$ has been completed in $\mathcal{L}_1$ and $\mathcal{L}_2$, the original paper looks like a single integral.
Although it is not written in the original paper, it is recommended that smart readers derive it by themselves. Try to expand formula (13) according to formula (4) mentioned above. I have tried to derive it here. Please correct me if there are any mistakes. First, substitute the two operators and pay attention to the position of the gradient operator here:
At the same time, we note the symmetry of the Green's function $\nabla G(\mathbf{r}, \mathbf{r}{\prime}) = -\nabla{\prime} G(\mathbf{r}, \mathbf{r}{\prime})$, so:
The same operation is used to obtain the remaining matrix elements. Here we directly copy the content of the additional materials of the paper:
$$ \begin{aligned} A_{\mathrm{EM}}^{mn}= & A_{\mathrm{HJ}}^{mn}=\int_{\mathbf{f}_m} \int_{\mathbf{f} _n} \mathbf{f}_m(\mathbf{r}) \cdot\left[\nabla G_1\left(\mathbf{r}, \mathbf{r}^{\prime}\right) \times \mathbf{f}_n\left(\mathbf{r}^{\prime}\right)\right] d \mathbf{r}^{\prime} d \mathbf{r} \\ & +\int_{\mathbf{f}_m} \int_{\mathbf{f}_n} \mathbf{f}_m(\mathbf{r}) \cdot\left[\nabla G_2\left(\mathbf{r}, \mathbf{r}^{\prime}\right) \times \mathbf{f} _n\left(\mathbf{r}^{\prime}\right)\right] d \mathbf{r}^{\prime} d \mathbf{r} \\ A_{\mathrm{HM}}^{mn} & =-j \omega \varepsilon_1 \int_{\mathbf{f}_m} \int_{\mathbf{f}_n} \mathbf{ f}_m(\mathbf{r}) \cdot \mathbf{f}_n\left(\mathbf{r}^{\prime}\right) G_1\left(\mathbf{r}, \mathbf{r}^{\prime}\right) d \mathbf{r}^{\prime} d \mathbf{r} \\ & +\frac{j}{ \omega \mu_1} \int_{\mathbf{f}_m} \int_{\mathbf{f}_n} \nabla \cdot \mathbf{f}_m(\mathbf{r}) G_1\left(\mathbf{r}, \mathbf{r}^{\prime}\right) \nabla^{\prime} \cdot \mathbf{f} _n\left(\mathbf{r}^{\prime}\right) d \mathbf{r}^{\prime} d \mathbf{r} \\ & -j \omega \varepsilon_2 \int_{\mathbf{f}_m} \int_{\mathbf{f}_n} \mathbf{f}_m(\mathbf{r}) \cdot \mathbf{f}_n\left(\mathbf{r }^{\prime}\right) G_2\left(\mathbf{r}, \mathbf{r}^{\prime}\right) d \mathbf{r}^{\prime} d \mathbf{r} \\ & +\frac{j}{\omega \mu_2} \int_{\mathbf {f}_m} \int_{\mathbf{f}_n} \nabla \cdot \mathbf{f}_m(\mathbf{r}) G_2\left(\mathbf{r}, \mathbf{r}^{\prime}\right) \nabla^{\prime} \cdot \mathbf{f}_n\left(\mathbf{r}^{\prime}\right) d \mathbf{r} ^{\prime} d \mathbf{r}\end{aligned} $$
After getting each element of the matrix, we introduce shift-invariant functions to help us get the Green's function and its gradient in different coordinate systems.
Then the matrix elements are expanded into a combination of different components of the basis functions, which act on the translation invariant functions. All elements of the final matrix have the following form. This form can speed up the construction and solution of the boundary element matrix.
$$ \int_{\mathbf{f}m} \int{\mathbf{f}_n} \psi_m(\mathbf{r}) g\left(\mathbf{r}-\mathbf{r}^{\prime}\ right) \xi_n\left(\mathbf{r}^{\prime}\right) d \mathbf{r}^{\prime} d \mathbf{r} $$
Finally, we get:
Here, $x, y, z, x^{\prime}, y^{\prime}, z^{\prime}$ are the Cartesian components of $\mathbf{r}, \mathbf{r}^{\prime }$. Now we have for $i=1,2$:
where $\mathbf{f}{mx}, \mathbf{f}{my}, \mathbf{f}_{mz}$ are the $x, y, z$ components of the vector basis function $\mathbf{f}_m$ . Similarly, we have:
$$ \begin{aligned} & A_{\mathrm{EM}, i}^{mn}=A_{\mathrm{HJ}, i}^{mn}=\int_{\mathbf{f}m} \int{ \mathbf{f}n} \mathbf{f}{mz}(\mathbf{r}) g_{2, i}\left(\mathbf{r}-\mathbf{r}^{\prime}\right) \mathbf{f}{ny}\left(\mathbf{r}^{\prime}\right) d \ mathbf{r}^{\prime} d \mathbf{r} \\ & -\int{\mathbf{f}m} \int{\mathbf{f}n} \mathbf{f}{my}(\mathbf{r}) g_{2, i}\left(\mathbf{r}-\mathbf{r}^{\prime}\right) \mathbf{f}{nz }\left(\mathbf{r}^{\prime}\right) d \mathbf{r}^{\prime} d \mathbf{r} \\ & +\int{\mathbf{f}m} \int{\mathbf{f}n} \mathbf{f}{mx}(\mathbf{r}) g_{3, i}\left(\mathbf{r} -\mathbf{r}^{\prime}\right) \mathbf{f}{nz}\left(\mathbf{r}^{\prime}\right) d \mathbf{r}^{\prime} d \mathbf{r} \\ & -\int{\mathbf{f}m} \int{\mathbf{f}n} \mathbf{f}{mz}(\ mathbf{r}) g_{3, i}\left(\mathbf{r}-\mathbf{r}^{\prime}\right) \mathbf{f}{nx}\left(\mathbf{r}^{\prime}\right) d \mathbf{r}^{\prime} d \mathbf{r} \\ & +\int{\mathbf {f}m} \int{\mathbf{f}n} \mathbf{f}{my}(\mathbf{r}) g_{4, i}\left(\mathbf{r}-\mathbf{r}^{\prime}\right) \mathbf{f}{nx}\left(\mathbf{r}^{\prime}\right) d \ mathbf{r}^{\prime} d \mathbf{r} \\ & -\int{\mathbf{f}m} \int{\mathbf{f}n} \mathbf{f}{mx}(\mathbf{r}) g_{4, i}\left(\mathbf{r}-\mathbf{r}^{\prime} \right) \mathbf{f}_{ny}\left(\mathbf{r}^{\prime}\right) d \mathbf{r}^{\prime} d \mathbf{r} \end{aligned} \tag{S.19} $$
Lastly, we also have:
$$ \begin{aligned} A_{\mathrm{HM}, i}^{mn} & =-j \omega \varepsilon_i \int_{\mathbf{f}m} \int{\mathbf{f}n} \mathbf {f}{mx}(\mathbf{r}) g_{1, i}\left(\mathbf{r}-\mathbf{r}^{\prime}\right) \mathbf{f}{nx}\left(\mathbf{r}^{\prime}\right) d \ mathbf{r}^{\prime} d \mathbf{r} \\ & -j \omega \varepsilon_i \int{\mathbf{f}m} \int{\mathbf{f}n} \mathbf{f}{my}(\mathbf{r}) g_{1, i}\left(\mathbf{r}-\mathbf{r}^{\prime}\right) \mathbf{f}{ny }\left(\mathbf{r}^{\prime}\right) d \mathbf{r}^{\prime} d \mathbf{r} \\ & -j \omega \varepsilon_i \int{\mathbf{f}m} \int{\mathbf{f}n} \mathbf{f}{mz}(\mathbf{r}) g_{1, i}\left( \mathbf{r}-\mathbf{r}^{\prime}\right) \mathbf{f}{nz}\left(\mathbf{r}^{\prime}\right) d \mathbf{r}^{\prime} d \mathbf{r} \\ & +\frac{j}{\omega \mu_i} \int{\mathbf{f}m} \int{\mathbf{f}n } \nabla \cdot \mathbf{f}_m(\mathbf{r}) g{1, i}\left(\mathbf{r}-\mathbf{r}^{\prime}\right) \nabla^{\prime} \cdot \mathbf{f}_n\left(\mathbf{r}^{\ prime}\right) d \mathbf{r}^{\prime} d \mathbf{r} \end{aligned} \tag{S.20} $$
This part of the code is inMVProd ClassIt doesn’t matter if you don’t understand it, because I just copied and scribbled it randomly, which is beyond the scope of graphics research.
In addition, the author also discussed the symmetry of EFIE and MFIE. It is this symmetry that makes the computational efficiency and space utilization lower. Note:
Since the matrix is symmetric, we do not need to calculate all the matrix elements, nor do we need to store all the matrix elements.
After solving for the surface current density, equation (6) can be used to calculate the scattered field propagating outward from the scattering surface.
3.2 Rough Surface Scattering: The Specifics
When simulating the interaction between electromagnetic waves and rough surfaces, the irregular geometric structure of the surface will have a significant impact on the scattering characteristics. The authors discretized the rough surface and divided the continuous surface into multiple units, each of which can be numerically calculated using the PMCHWT equation.
3.2.1 Rough Surface Samples
Represent the rough surface as a two-dimensionalHeight field, and then the height field is discretized and divided into multiple rectangular units.
Only a surface sample of size $$L_x \times L_y$$ is considered for each simulation, and a step size $d$ is chosen to define the discretization grid.
$$ x_s = s \cdot d, \quad s = 0, 1, \ldots, N_x \ y_t = t \cdot d, \quad t = 0, 1, \ldots, N_y \tag{21} $$
Among them, $ N_x = L_x / d$ and $N_y = L_y / d$, respectively, represent the number of cells divided in the $x$ and $y$ directions. At each discrete point $(x_s, y_t)$ there is a height field $h(x_s, y_t)$ function.
The authors concluded that the height variations of rough surfaces are very small, typically only a few micrometers, which is comparable to the wavelength of visible electromagnetic waves.
3.2.2 Basis Elements and Functions
Each primitive has four corners, each corner has a different height, and each corner contributes differently to the current and magnetic flux. Therefore, four basis functions are defined on each small square to approximate the situation on each small square.
In most simulations, the step size $d$ is chosen to be around $\lambda / 16$ of the wavelength to ensure accuracy.
Each primitive is parameterized by two parameters u and v, both in the range [-1, 1].
The shape of the primitive is represented by a bilinear function $\mathbf{r}(u, v)$, where $(s, t)$ represents the index of the current primitive, and the primitive is determined by the coordinates of the four vertices:
The introduction of Jacobian is used to transform the coordinate system and ensure that the basis functions have appropriate proportional relationships in different $u, v$ directions.
3.2.3 Gaussian Beam Incidence
Since I don't know much about optics, the following is my personal understanding. A Gaussian beam is a phenomenon in which the middle part of the traveling wave field appears to be concave inward when the laser is emitted. In other words, a Gaussian beam describes the energy distribution of light on the cross section. The focus of plane waves and spherical waves is to describe the direction of energy propagation. During the propagation process, the wavefront shape of a Gaussian beam is approximately a spherical wave.
Gouy phase is a phase delay effect in Gaussian beam propagation. The phase of the beam will increase after passing through the focus. On the other hand, Gaussian beams satisfy a solution of Maxwell's equations under paraxial conditions and can be approximated as non-uniform spherical waves. I feel that there is no need to study it too deeply at present.
Back to the original paper, the advantage of Gaussian beam is that the size of the incident field can be controlled, and thus the surface induced current density of an area slightly larger than the irradiated area can be controlled to be a non-zero value.
A Gaussian beam is an electromagnetic wave whose amplitude is distributed in a two-dimensional Gaussian pattern in a plane perpendicular to the propagation direction [Paschotta 2008]. Its energy is mainly concentrated near the center of the beam. Looking at the figure (a) above, a Gaussian beam can be described by the focal plane $P$ , the center point $o$ , and the beam waist $w$ . The field intensity decays with position as $e^{-r^2 / w^2}$ . When the distance from the center exceeds $2.5w$ , the field intensity decays to a very small value and can be considered to be almost zero.
However, Gaussian beams also have a certain degree of divergence. The divergence angle $\theta$ is approximately proportional to the wavelength $\lambda$ and inversely proportional to the beam waist diameter $w$. Formula:
When the beam is incident on the surface obliquely, the Gaussian beam has an elliptical cross section in the focal plane, as shown in the figure below. There are different beam waist sizes in two perpendicular directions: one parallel to the plane of the incident direction and the surface normal, and the other perpendicular to it.
In order to ensure that the irradiation area of the Gaussian beam in different directions on the surface is the same, two beam waist widths in the lateral direction are introduced here:
According to the original idea, $N$ basis functions are used to represent the current magnetic flux density, and the size of the matrix is $2N \times 2N$ . If the matrix is solved directly (LU decomposition, Cholesky decomposition, etc.), the total complexity may be $\mathcal{O}(N^3)$ . Even if the conjugate gradient method is used, the total complexity is still $\mathcal{O}(N^2)$ . In some small-scale simulations, the size of the basis function is about $960*960$ , and the storage requirement is about 29.4GB. Using the Adaptive Integral Method, AIM, the total storage requirement is about 76.8 MB based on 8 bytes. What is so magical about this method? Let me show you together!
4.1.1 Approximating Matrix Elements
AIM was first proposed by Bleszynski et al. [1996]. The core idea of AIM is to approximate the effect of each basis function as the effect of a set of point sources, avoiding the direct calculation of the exact interaction between each pair of basis functions, and at the same time spread the influence of each basis function through FFT to improve the calculation efficiency.
The calculation method of matrix elements in AIM is to approximate the calculation by linear combination of some items of matrix elements. This is why the readers are asked to further deduce after the above formulas (13)-(16) to derive the final result so that it conforms to the form of AIM.
AIM first creates a global 3D Cartesian grid in the space containing the electromagnetic fields and field sources, as shown in Figure (6) below.
In order to further simplify formula (27), the AIM algorithm approximates the original basis function as a point source on a set of grid points in this three-dimensional Cartesian coordinate system. In other words, it is a continuous-to-discrete transformation, which is convenient for subsequent FFT.
Kai Yang and Ali E Yilmaz. 2011. Comparison of precorrected FFT/adaptive integral method matching schemes. Microwave and Optical Technology Letters 53, 6 (2011), 1368–1372.
4.1.2 Base and Correction Matrices
Based on formula (29), a set of base approximation matrices $B_{EJ}, B_{EM}, B_{HJ}, B_{HM}$ are defined as approximations to specifically handle basis function pairs with a long distance. These matrices simplify the calculation by introducing $\Lambda$ matrices and convolution operations. At the same time, for basis function pairs with a short distance ($d_{near}$), correction matrices are introduced to reduce the error. $C_{EJ}, C_{EM}, C_{HJ}, C_{HM}$ is a sparse matrix defined as follows:
$A_X^{mn}$ is the exact value of the original matrix, while $B_X^{mn}$ is the approximation of the fundamental matrix. By subtracting the approximation of the fundamental matrix, a more accurate correction term is obtained to compensate for the error of the close basis function pairs.
In summary, the final approximate form of each matrix in the AIM method can be written as follows:
In other words, the original matrix can be approximated by a combination of the fundamental matrix and the correction matrix.
4.1.3 Fast Matrix-Vector Multiplication
Fast Matrix-Vector Multiplication is the core of AIM.
Since the modified matrix $C$ obtained above is a sparse matrix, $C$ has non-zero values only on close-range basis functions, so the multiplication operation of the matrix $C$ is very fast.
Using the convolution property of the basic approximation matrix $B$, the product of the matrix $B$ and the vector is calculated. The calculation process is divided into three steps:
The first step is to project the vector onto a sparse matrix grid.
The second step is also the core step, which is to propagate the data of the grid points to the entire network, that is, to calculate the influence of each point on other points. The closer the two points are, the larger the propagation function in the matrix $G$ will be. FFT is used to speed up the process.
The third step maps the result back to the original basis function space.
4.2 GPU-Accelerated Iterative Solving
The computational focus of the AIM method is shifted to fast Fourier transform (FFT) and sparse matrix operations on the GPU. In summary, the large matrix is currently divided into the basic matrix $B$ and the correction matrix $C$ to handle the long-distance and short-distance basis function pairs respectively.
cuFFT: Converts the underlying matrix multiplication operation to a convolution calculation in the frequency domain
cuSPARSE: Accelerate the correction matrix C for sparse matrix calculations
And optimize the calculation strategy:
For small-scale simulation tasks, only one GPU is needed.
For large-scale simulation tasks, distribute the tasks to 4 GPUs
4.2.1 Small-Scale Simulations
For small-scale simulation tasks (e.g. $12 \mu m \times 12 \mu m$), the Fourier transform of the propagation function (i.e. the Fourier transform of the matrix $G$) must be calculated and stored in advance. In small-scale tasks, the sparse correction matrix C occupies less than 5GB of video memory and can be handled by a single GPU.
4.2.2 Large-Scale Simulations
For large-scale simulation tasks (e.g. $24 \mu m \times 24 \mu m$), since the video memory of a single GPU is not enough to store all the data, the author distributes the computational tasks to 4 GPUs. In this scale simulation, the number of basis functions will reach 960 × 960, and storing all non-zero elements of the correction matrix (including row and column indices and complex floating-point values) requires about 20GB of video memory. The strategy is still the same as the small scale, and each GPU is allocated about 5GB of memory to store the correction matrix $C$.
The MINRES solver is executed on the host CPU, while the matrix-vector product $y = Ax$ is computed on the GPU. But don’t worry about the transfer time, the vector is only about 30MB.
4.3 FFT-Accelerated Scattered Field Evaluation
The FFT is used to accelerate the calculation of the scattered field. The field scattered from the surface is evaluated in the far field region to finally find the surface BRDF.
After solving the BEM, the current density $\mathbf{J}$ and the magnetic flux density $\mathbf{M}$ on the surface are obtained. These density distributions define the electromagnetic sources on the surface and can be used to calculate the scattered field in the far field region. The formula is very simple, and it will weaken with distance and have a certain phase change:
The right side of the formula $\mathbf{E}(\hat{r})$ and $\mathbf{H}(\hat{r})$ are the amplitudes of the specific directions $\hat{r}$ in the far field. In different directions, the intensity of the scattered field may be different.
To avoid solving these integrals directly, the author uses the point source approximation ($\Lambda$ matrix) of $\mathbf{J}$ and $\mathbf{M}$ previously (4.1.2) to discretize each integral term $F_i(\hat{r})$ and rewrite it in Fourier transform form, as shown in formula (38):
Convert the continuous field strength calculation into a discrete sum so that it can be quickly calculated using FFT. And from the above formula, it can be observed that $F_i(\hat{r})$ is actually the Fourier component of $h_i(p)$ at the spatial frequency $-k\hat{r}$.
The required spatial frequencies are not on the FFT grid but can be interpolated; we add zero padding prior to the FFT step, to ensure enough resolution in the frequency domain for the trilinear interpolation to be sufficiently accurate.
5 HIGH RESOLUTION BRDF GENERATION
After a lot of work, I finally got back to the familiar BRDF calculation. The key here is to use the linear superposition of small-scale simulations to reconstruct the far-field scattering of the large-scale incident field, rather than trying to do it all at once. A grid of $N^2$ small-scale Gaussian beams is linearly combined to approximate the large-scale incident field.
This is done using a technique called beam steering, which significantly reduces computational costs by not having to simulate each direction.
5.1 Basic and Derived Incident Directions
First, $N^2$ Gaussian beams propagating in a certain direction $\mathbf{u}$ form a grid of $N \times N$ points in the receiving plane. These beams are combined to generate a large total field.
Then, a complex scaling factor is introduced into each Gaussian beam to adjust the phase of each beam and thus adjust the propagation direction of the combined field. These directions are called desired directions.
When the angle between the target incident direction $-\omega_i$ and the basic direction $\mathbf{u}$ is close to the divergence of the small beams, aliasing artifacts begin to appear. An example is shown in Fig. 8 (d).
In our framework, we decide on a primary waist w and choose a collection of basic incident directions. In general, a smaller waist width means a larger divergence angle, so that more incident directions can be derived from each basic direction, thereby reducing the number of basic directions required. A larger waist width will reduce the divergence angle of each Gaussian beam, making the total field divergence after combination smaller, thus producing a more accurate incident direction.
The center of each hexagon corresponds to a basic direction of incidence. All the directions of incidence in the entire hemisphere are divided into several territories, and each territory belongs to a basic direction of incidence. In the hemispherical projection, this proportional relationship and the cosine factor cancel each other out, so you can useHexagons of equal sizeTo express.
For each cardinal direction, several derived directions that deviate from the cardinal direction can be generated through "beam steering".
When the incident angle of light is very small (for example, close to the surface normal direction), the range of derived directions near the basic incident direction is very small, because the light at a small angle is more concentrated and will not diffuse greatly. Conversely, the range of derived directions will be larger. In summary, the larger the incident angle (the closer the angle is to the horizontal), the larger the coverage of the derived directions. The author of the formula also mentioned: $1/\cos \theta_i$.
5.2 Individual Simulations and Synthesized Results
In order to calculate the BRDF, we need to know the scattering of this large area of incident light. However, directly simulating such a large area of light will require a high computational cost. Therefore, we:
The large-scale incident field is simulated using a superposition of small-scale simulations.
Think of it as using many small flashlights (Gaussian beams) to cover an area, rather than one giant searchlight.
First, determine the size of the flashlights (the size of the beam, i.e. the waist width $w$) and evenly distribute them on the surface. Written as a formula, this arrangement of flashlights is the grid point ${x_s}, {y_t}$, representing the center position of each Gaussian beam. The grid spacing is generally the same as the waist width to ensure uniform coverage of the light and keep the divergence angle low.
Let each Gaussian beam produce the same electromagnetic field in its central region, just repeating the effect at different locations.
Next, to get the total scattered field of a large beam, we needPhase FactorTo perform “adjustment” and “superposition”, please refer to formula (39) for details.
Finally, Combining with Eq. 39, the scattered fields in the far field region corresponding to the pair of directions $(w_i,w_o)$ are given by:
where $\mathbf{n}$ is the surface normal at the macro scale ( $+\mathbf{z}$ ). Note that Eq. 42 and Eq. 43 can also be applied in single simulations, where $\Phi_i$ is computed from a single Gaussian beam.
Disclaimer: This is a pure garbage article, which is a supplement to the previous two articles. I just sorted it out for my own review. All titles can be redirected to the paper homepage. Special terms are given in Chinese and English as much as possible. If there are any mistakes, please point them out. Thank you very much.
Wave optics is used to render hair. The surface electromagnetic field is calculated to obtain the scattered field, and then noise is added to simulate the Glints effect.
I found that the authors of this series are all very good-looking. (crossed out)
1. Background
Hair rendering has been mainly based on ray tracing technology, which cannot handle wave optics effects, such as strong forward scattering and subtle color changes on the hair surface. Previous research [Xia et al. 2020] demonstrated that diffraction effects play a key role in the color and scattering direction of fibers. However, this study did not consider surface roughness and the microstructure of the fiber epidermis (such as tilted keratin scales).
2. Motivation
In order to make up for the lack of treatment of diffraction and forward scattering (such as Glints phenomenon) in the existing light optics model.
Although full-wave simulations can produce very detailed scattering data, the computational effort is still too high and must be accelerated or simplified in some way to achieve hair or fur rendering in large-scale scenes.
We wanted to develop a model that could efficiently handle various fiber geometry variations.
3. Methods
Hair modeling is based on scanning electron microscope (SEM) images of hair.
Use "WAVE SIMULATION WITH 3D FIBER MICROGEOMETRY" to calculate the reflection and diffraction of rough fiber surfaces. That is, PO.
Speckle theory is introduced to analyze the statistical characteristics of the scattering pattern, and noise is used to accelerate it.
Wave optics, iridescence effect, Quetelet scattering model of water droplets on water surface
Combining Mie scattering, Quetelet scattering (light interference) and dynamic changes of water droplets, the rainbow-like color effect of water droplets on the water surface and in the steam is realistically rendered, surpassing the traditional single Mie scattering model.
1. Background
Iridescence is common in nature, especially in water droplets, fog and steam. It can generally be explained by Mie scattering. Mie scattering describes the scattering effect that occurs when light encounters spherical particles of the same wavelength. It is one of the important theories currently used to simulate natural phenomena such as water droplets, clouds, rain and fog.
However, while Mie scattering can explain the optical properties of isolated water droplets, it cannot fully explain phenomena such as the iridescence of water droplets on the surface of water and the complex rainbow patterns in vapor. Phenomena depend not only on how individual particles scatter light, but also on surface reflections, interference effects, and dynamic changes in particle size.
2. Motivation
Mie scattering can only deal with isolated light scattering phenomena and cannot explain more complex optical interference effects.
Accurately simulating these natural phenomena can greatly improve the realism and look and feel of image rendering.
Existing computer optical models and rendering methods are mostly limited to Mie scattering and cannot explain the interaction of light in a multi-particle environment, such as light interference and reflection between water droplets or between water droplets and surfaces.
3. Methods
The "Quetelet scattering model on water" is used to explain the rainbow effect produced by water droplets floating on the water surface. By building an empirical model, thermal imaging technology is used to relate temperature to the size and height of water droplets. Quetelet scattering phase function and BRDF (bidirectional reflectance distribution function) are used to render particle groups and water surfaces.
A water droplet growth and evaporation model was developed to simulate the dynamic changes of water droplets in steam. Combined with Mie scattering, water droplets of non-uniform size were used to simulate the rainbow color changes in steam. In order to improve rendering efficiency, an acceleration algorithm based on motion blur was used, which increased the calculation speed by 10 times compared with traditional methods.
The method of learning hair higher-order scattered radiance online combined with a small multilayer perceptron (MLP) significantly accelerates hair rendering in a path tracing framework, reducing computation time and introducing only a small amount of bias.
1. Background
The multiple scattering of hair is very complex, especially in the path tracing process, because it is necessary to simulate the multiple scattering of light between hairs, which makes it difficult to converge.
2. Motivation
Develop a method to improve computational efficiency while maintaining high-quality simulation of multiple scattering effects.
In the existing technology, some methods make simplifying assumptions about the scene or lighting. This paper hopes to propose a general method that does not make any assumptions about the scene.
3. Methods
A small multilayer perceptron (MLP) is used to learn higher-order scattered radiance online. This MLP network learns the scattering properties of hair in real time during the rendering process, without relying on pre-computed tables or simulations.
The MLP is integrated into the path tracing framework to infer and compute higher-order diffuse radiation contributions.
The renderer's bias and speedup can be adjusted in real time to find the optimal balance between computational efficiency and rendering quality.
Accelerated area light source Glints, microsurface models, real-time rendering
Rendering glints under area lights is done in real time by combining Linearly Transformed Cosines (LTC) with a microsurface count model based on the binomial distribution.
1. Background
Many real-world materials (such as metals, gemstones, etc.) have a glittering appearance, which is caused by the reflection of micro-surfaces. However, glitter is a discrete phenomenon, and the computational complexity of wave optics simulation is too large.
Previous studies have mostly focused on using infinitesimal point light sources to render flash effects, which is a reasonable simplification for distant light sources like the sun, but in reality most light sources are essentially area light sources. Existing technologies have not been able to effectively handle flash rendering under area light sources.
2. Motivation
Glint rendering under area lights. Area lights (such as the light shining into a room through a window) are a common type of light, and how to efficiently render glint effects under such lights is an unsolved problem. We hope to develop a method that can accurately render glint effects under area lights while meeting the needs of real-time rendering.
It is hoped that it can be easily integrated into existing real-time rendering frameworks without introducing significant additional overhead to existing area light shading methods.
3. Methods
Glint reflection probability estimation computes the probability that a microfacet is correctly oriented to reflect light from a light source to an observer, using Linearly Transformed Cosines (LTC) for large sources and a locally constant approximation for small sources.
The number of reflective microsurfaces is counted using a binomial distribution-based counting model.
An innovative real-time strand-based hair rendering framework is proposed, which ensures the consistent appearance of hair at different view distances and achieves significant rendering acceleration through seamless level-of-detail (LoD) transition.
1. Background
Strand-based hair rendering is becoming increasingly popular in film, television and game production for its realistic appearance, but it is computationally very expensive, especially at long viewing distances.
The current LoD method is prone to noticeable discontinuities in the transition from hair strands to cards, resulting in inconsistent appearance.
2. Motivation
Solve discontinuity in dynamics and appearance. Existing solutions for converting hair strands to hair cards have significant differences in appearance and animation performance. The goal of this paper is to achieve seamless LoD transition from far to near, eliminating appearance changes during transition while maintaining computational efficiency.
3. Methods
Encapsulates multiple hair strands within an elliptical volume using an elliptical thick hair model. The shape and overall structure of the hair cluster is maintained at different LoDs, providing a consistent look as the view distance changes.
The elliptical bidirectional curve scattering distribution function (BCSDF) simulates single and multiple scattering phenomena within hair clusters and is suitable for hair distribution scenarios ranging from sparse to dense and from static to dynamic.
Dynamic LoD adjustment and hair width calculation.
A model based on tiny example microstructures that renders glinty effects in real time, significantly reducing memory usage and computational overhead while maintaining the realism of high-frequency reflection details.
1. Background
The shimmering details produced by complex microstructures can significantly improve the realism of renderings, especially on materials such as metals and gemstones. These details usually require high-resolution normal maps to define each micro-geometry, but such methods have high memory requirements and are not suitable for real-time rendering applications.
2. Motivation
Reduce memory and computational overhead.
Leveraging material self-similarity: Many materials have independent structural features and self-similarity, and small samples are used to implicitly represent complex microstructures, thereby reducing memory requirements.
3. Methods
A tiny example-based procedural model based on the microstructure of a small sample can generate complex sparkle details by reusing a small number of samples based on the self-similarity of the material.
Precomputed Normal Distribution Functions (NDFs) Precompute and store small samples of normal distribution functions (NDFs) using 4D Gaussians. Stored in multi-scale NDF maps and called by simple texture sampling at rendering time.
A tiny example-based NDF evaluation method combines texture sampling with a small example NDF evaluation method to quickly generate the shiny appearance of complex surfaces.
A practical hair appearance aggregation model that significantly accelerates hair rendering while maintaining realistic visual effects by reducing the number of geometric hairs and combining multiple scattering of light, using neural networks to achieve real-time dynamic simplification.
1. Background
If there are too many hairs, the light scattering and reflection of each hair will greatly increase the calculation amount, especially when simulating multiple light scattering.
Most existing simplification methods improve rendering efficiency by reducing the number of hairs, but this method has great limitations. This method can cause the hair to look too rough or dry, and the reflection and scattering effects of light are not realistic.
2. Motivation
Reducing geometric complexity.
Improving rendering efficiency.
3. Methods
An aggregated fur appearance model is proposed, which uses a thick cylinder to represent the optical behavior of a group of hair clusters. By analyzing the optical properties of individual hairs (such as the incident angle of light), the model can accurately reflect the aggregated appearance of hair clusters.
A lightweight neural network is used to map the optical properties of individual hairs to parameters in the aggregate model.
A dynamic level-of-detail scheme based on view distance and number of light bounces is proposed to dynamically simplify the geometric structure of hair.
Multi-modal data is important for predictive rendering, especially in accurately modeling material reflection behavior. By combining spectral, spatial information and micro-geometric details, the realism and computational efficiency of reflection models can be improved.
1. Background
The need for predictive rendering aims to accurately simulate the appearance of materials.
Most current databases on material reflection behavior are limited to a single dimension, usually covering only the spectral domain or the spatial domain, and lack descriptions of microgeometry details.
2. Motivation
In order to address data limitations, multimodal data can not only better simulate the reflection of materials under different lighting conditions, but also reveal the influence of the microscopic geometry of the material surface on light scattering.
Multimodal reflectance data can help develop more realistic and efficient reflectance models.
3. Methods
Building a multi-modal reflection database, including spectral data, spatial distribution data and microgeometry details of the material.
Simulating microgeometry of the microgeometry of a material surface.
A bidirectional scattering distribution function (BSDF) based on free-space diffraction can efficiently simulate the diffraction phenomenon of light around the edges of objects in complex scenes through ray tracing without the need for geometric preprocessing, and is particularly suitable for path tracing technology.
1. Background
Free-space diffraction is an optical phenomenon in which light is diffracted when it encounters an edge or corner of an object, bending some of its energy into the shadowed area. This phenomenon is important for modeling the propagation of light waves, especially at long wavelengths, such as radar, WiFi, and cellular signals.
The limitations of traditional methods such as the Geometric Theory of Diffraction (GTD) and the Unified Theory of Diffraction (UTD) are the extremely high computational complexity caused by the need to deal with light rays that interfere with each other, especially in complex geometric scenes. Existing methods rely on scene simplification and specific geometric structures and cannot effectively handle complex three-dimensional scenes.
2. Motivation
Addressing diffraction rendering in complex scenes. Existing diffraction simulation methods are difficult to scale and make compatible with path tracing techniques.
Existing diffraction methods often rely on complex nonlinear interference calculations, while path tracing uses linear rendering equations. This paper hopes to design a free-space diffraction BSDF that works efficiently within the path tracing framework without requiring major modifications to the path tracer.
3. Methods
The Fraunhofer diffraction edge model is based on Fraunhofer diffraction. Near the intersection of light and geometric objects, the relevant edges are identified and the diffraction effects are calculated. When the light hits the object, the BSDF of free space diffraction is constructed through geometric analysis to quantify how the light propagates around the geometric object and how much energy is diffracted.
The importance sampling strategy evaluates the geometric edges around the points where the ray interacts with the object and samples and traces the diffracted rays.
Multi-layer oil film rendering, iridescence effect, wave optics
A multi-layer interference rendering method. It can express the iridescence effect of periodic multi-layer structures. By introducing the Huxley method from biology, it can achieve efficient calculation independent of the number of layers.
1. Background
Thin-film interference is an optical phenomenon caused by the wave properties of light waves, which produces iridescence when the viewing angle or illumination angle changes. It usually appears in single-layer or multi-layer structures in nature, such as butterfly wings, beetle shells and dielectric mirrors.
The limitations of existing methods such as recursive calculation method and transfer matrix method (TMM) are that the computational complexity increases significantly with the number of layers. Simplified methods ignore multiple reflections in thin films.
2. Motivation
Improving efficiency for multilayer structures.
Applied to physical rendering of complex materials.
3. Methods
A multilayer interference model based on Huxley's approach is proposed. It can efficiently calculate the reflection and transmission coefficients in periodic multilayer structures and supports multiple materials and absorption effects.
Based on BRDF implementation. Implemented as a BRDF (Bidirectional Reflectance Distribution Function), it can be integrated into traditional rendering systems such as PBRT-v3.
Full-wave simulator based on the boundary element method (BEM) that can calculate light scattering on rough surfaces with high accuracy. It is used to evaluate and improve approximate reflection models in computer graphics, especially when multiple scattering, interference and diffraction effects are significant.
1. Background
Surface reflection models are usually based on geometric optics, which assumes that light propagates in the form of rays. For scenes where surface features are comparable to the wavelength of light, geometric optics models cannot accurately capture wave effects such as diffraction and interference.
Based on wave optics approximations, such as Beckmann-Kirchhoff theory and Harvey-Shack model, they still produce errors under multiple scattering and complex geometric structures.
2. Motivation
Since existing reflection models have different accuracy in different situations, there is a lack of reliable benchmarks to verify their accuracy. The goal of this paper is to develop a simulator based on full-wave theory to minimize approximations and achieve high-precision surface reflection calculations through numerical discretization, thereby providing a reference tool that can be used to evaluate the accuracy of various reflection models.
Addressing multiple scattering and wave effects.
3. Methods
Boundary Element Method (BEM), accelerated by Adaptive Integral Method (AIM).
The simulator's full-wave simulation completely solves Maxwell's equations and can accurately simulate wave phenomena such as light propagation, interference, and scattering.
And it can efficiently calculate BRDF (efficient BRDF computation).
We propose an efficient Physical Light Transport (PLT) framework that exploits the principles of partially coherent light and wave optics to achieve accurate rendering of interference, diffraction, and polarization effects in complex scenes through an improved rendering algorithm, bringing its performance close to that of classic “physically based” rendering methods.
1. Background
Most existing rendering methods ignore the wave characteristics of light, especially in complex scenes, which makes it impossible to render physical phenomena such as interference and diffraction of light, which are particularly important on certain materials (such as iridescent coatings, optical discs, etc.). To solve this problem, a rendering framework based on Maxwell's electromagnetic theory is proposed.
Although PLT provides a theoretical full-wave model that can simulate the coherence, interference and diffraction of light, existing methods are very computationally difficult.
2. Motivation
Simplifying the physical light transport model.
Introducing new coherence-aware materials and developing material models that can perceive light coherence will improve the usability of PLT in practical scenarios.
3. Methods
Restricting the coherence shape of light, through thermodynamic derivation, proves that this approximation is reasonable under most natural light sources.
An extended Stokes-Mueller calculus is used to combine the radiation, polarization and coherence properties of light as new rendering primitives. The generalized Stokes parameters can fully quantify all properties of light and accurately simulate complex optical phenomena caused by these properties, such as interference and diffraction.
Wave BSDF and importance sampling.
New coherence-aware material models take full advantage of the coherence properties of light to expand the scope of application of PLT.
The first hair scattering model based on microsurface theory is proposed to accurately describe the scattering behavior of hair, including non-separable scattering lobe structure, elliptical cross section, efficient importance sampling and forward scattering spot (glint-like) effect.
1. Background
Complexity of hair rendering. Most existing hair scattering models simplify the mathematical calculations through separable scattering lobes, which are fast but not ground truth.
Most hair scattering models are based on geometric simplification, treating hair as smooth cylinders, which leads to deviations in scattering behavior.
2. Motivation
Introducing a physically-plausible microfacet model more accurately describes the scattering behavior of hair: the surface microscopic roughness, the tilted scale structure, and the non-separable scattering lobe shapes.
Improving sampling efficiency and physical accuracy.
3. Methods
The hair modeling is combined with microfacet theory, and GGX or Beckmann normal distribution is applied to describe the microscopic roughness of the surface. And it is non-separable lobes.
The bidirectional curve scattering distribution function (BCSDF) describes the complex interaction of light on the hair surface.
Support for elliptical cross-sections and efficient sampling. Support for elliptical cross-sections for hair.
The first global light transport framework based on Maxwell's electromagnetic theory that can handle partially coherent light is proposed, which accurately simulates the interference and diffraction effects of light and extends the traditional radiometric-based light transport theory to the field of wave optics.
1. Background
Existing light transport models are usually based on geometric optics and radiometry, which ignore the wave characteristics of light and cannot simulate phenomena such as interference and diffraction. They cannot accurately reproduce wave optical effects such as rainbow effect, grating, thin film interference, light polarization, etc., which is the limitation of classical radiometric light transport models.
Current models can only handle local treatment of wave effects, but cannot account for the transmission and coherence of light in global scenes.
2. Motivation
Achieving global wave-optics consistency in light transport, that is, combining Maxwell's electromagnetic theory.
Combining wave optics with classical geometric optics (integrating wave optics with classical geometric optics) can deal with the wave effect of light and be consistent with classical geometric optics in the short wavelength limit.
3. Methods
Modeling partially-coherent light is divided into two parts: two-point coherence description and light source model. Different from traditional radiance, this paper introduces a "cross-spectral density function" based on the partial coherence of light, which can capture the interference characteristics of light. The physical model of natural light sources is based on the principle of spontaneous radiation in quantum mechanics.
Generalizing the light transport equation. The spectral-density transport equation is used to calculate the interference and diffraction effects of light during propagation. This paper also proves that the framework can be simplified to classical geometric optics in the short wavelength limit, so it can be seamlessly integrated with existing light transport methods.
A generalized ray formal model is proposed for wave optics rendering. By solving the sampling problem, weak locality, linearity and completeness are simultaneously established in wave optics. Bidirectional wave optics path tracing and efficient rendering are achieved in complex scenes.
1. Background
The classical model of light transport is based on ray optics, which assumes that light propagates as a point query in a linear manner. However, ray optics cannot capture the wave nature of light and ignores interference and diffraction phenomena, such as the iridescence effect, thin film interference, and diffraction of long-wave radiation.
Although wave optics can accurately describe the interference and diffraction effects of light, traditional sampling and path tracing techniques are difficult to apply due to its nonlinear behavior.
2. Motivation
Solving the sampling problem in wave optics. In order to apply wave optics in bidirectional optical transmission, it is necessary to solve the sampling problem under weak locality.
Develop a novel formalism of wave optics that enables efficient applications in inverse path tracing and bidirectional light transport while maintaining linearity and completeness.
Improving wave-optics rendering efficiency, making the convergence speed of wave-optics rendering close to that of classical ray optics rendering systems.
3. Methods
Introduction of the generalized ray. Perform weak local linear queries. Generalized rays are no longer limited to point queries at a single location, but occupy a small spatial region. They can capture the interference and diffraction effects of light.
Weak locality and linearization. In wave optics, perfect locality and linearization cannot be achieved simultaneously. Therefore, perfect locality is abandoned. Weak locality is adopted to ensure that generalized rays can be linearly superposed.
A new framework for light-matter interaction is proposed, which unifies the formulas for scattering and diffraction by decomposing partially coherent light into the Hermite-Gauss space and modeling matter as a locally stationary random process, and enables efficient calculation and description of complex optical phenomena.
1. Background
The light observed in daily life is usually composed of many independent electromagnetic waves. Due to the complexity of partially-coherent light, the coherent properties of partially coherent light, such as reflection on microscopic geometric surfaces, the appearance of coating materials, grating effects, etc., cannot be explained by classical radiosity theory.
The limitations of existing tools only allow for rendering of specific materials and are difficult to generalize.
2. Motivation
Building a general-purpose light-matter interaction framework to efficiently process partially coherent light and simplify the complexity of existing computational tools.
Decomposing light coherence properties, the Hermite-Gauss space is introduced in the hope of decomposing and representing the coherence of light in a computationally feasible way, which is widely applicable to various optical phenomena.
Simulate the optical properties of pearlescent materials, and provide a theoretical basis for the design and reverse rendering of pearlescent materials.
1. Background
Wide Applications of Pearlescent Materials. These materials have unique gloss and color-changing effects and are widely used in packaging, ceramics, printing, cosmetics and other fields.
The complex optical processes of pearlescence are derived from multiple scattering and wave optical interference between pigment flakes. Existing models are difficult to fully describe these complex optical behaviors.
2. Motivation
Building a More Comprehensive Model for Pearlescent Materials. Existing pearlescent material models do not adequately account for the complex structure of pigments and the effects of the manufacturing process. The goal is to expand the range of pearlescent appearances that can be represented by introducing new optical simulation models.
A generic pearlescent material model can also be used in reverse rendering.
3. Methods
An optical model based on interference pigments is proposed, which takes into account the multilayer structure of pigment flakes, the directional correlation of particles, thickness variation and other characteristics.
Systematic Study of Parameter Space, exploring the effects of orientation, thickness, and arrangement of pigment flakes on the material’s appearance.
Inverse Rendering helps interpret light scattering phenomena in the real world.
A wave optics model based on non-paraxial scalar diffraction theory is used to simulate the iridescence effect on microscopic scratched surfaces, from local spots to smooth reflections at long distances.
1. Background
Optical Effects of Scratches: Under directional lighting (such as sunlight or halogen lamps), these scratched surfaces will show complex iridescent patterns, which are caused by the diffraction of incident light by the scratch structure. This cannot be reproduced in the geometric optics model.
Although existing analytical models are able to reproduce the iridescence effect of some microstructures (such as optical discs), simulation of the optical behavior of locally resolved scratches remains a challenge.
2. Motivation
Provide a Wave-Optical Scratch Rendering Framework, which can accurately simulate the optical effects caused by scratches, including light spots, iridescence and other visual phenomena.
3. Methods
Wave-Optical Model Based on Non-Paraxial Scalar Diffraction Theory: The method in this paper can accurately simulate the diffraction behavior of light on micro-scale surface features at large angles of incidence and reflection.
Vector Graphics Representation of Scratch Surfaces.
Multi-Scale BRDF Model.
Integration and Optimization in Physically-Based Rendering Systems.
Describes the material iridescence effect of 1D photonic crystals (i.e. Bragg mirrors). Simplifies to a single bounce BRDF for fast computation under certain conditions.
1. Background
Iridescence in nature is manifested in organisms, plants or gemstones. It is caused by specific microscopic geometric structures whose size is comparable to the wavelength of visible light. The most prominent example is photonic crystals, which produce structural colors by repeating in one-, two- or three-dimensional structures.
1D photonic crystals, or optical properties of Bragg mirrors. Most existing works use the classic transfer matrix method to calculate the optical effects of multilayer films, but as the number of films increases, the computational complexity increases significantly.
2. Motivation
Simplifying the computation of Bragg mirror reflectance, introducing a more concise, closed-form reflection formula and exploring fast approximation methods in RGB spectral rendering.
Investigating the effects of rough Bragg layers to explore the influence of surface roughness on optical performance.
3. Methods
Introduce the closed-form reflectance formula. Based on Yeh's formula (Yeh88 Formula), do RGB spectral approximation (RGB Spectral Approximation).
Analyze the effect of roughness on optical transmission.
The appearance of a rough Bragg layer is efficiently rendered using the Single-reflection BRDF Model.
Simulating speckle statistics under near-field imaging conditions in scattering media accelerates speckle rendering in biological tissue imaging applications and provides support for speckle-based imaging techniques.
1. Background
When performing deep imaging in biological tissues, imaging becomes very difficult due to multiple scattering of light inside the tissue. When irradiated with coherent light (such as laser), high-frequency speckle patterns are generated inside the tissue. The statistical properties of speckle patterns, especially the memory effect, provide the basis for tissue imaging techniques (such as fluorescence imaging and adaptive optical focusing).
The limitations of existing models are that they mainly focus on far-field imaging, while near-field conditions are ignored.
2. Motivation
Developing a Physically Accurate and Efficient Model for Near-Field Speckle Rendering.
Improving Computational Efficiency of Speckle Simulations. The wave equation solver is too computationally intensive.
3. Methods
Monte Carlo Path Integral Rendering Framework.
Aperture and Phase Function Approximations.
Importance Sampling.
[Kajiya and Kay 1989] Kajiya-Kay Model
The originator of hair, no need to say more
The hair is simplified as a thin and long cylinder, and the light reflection behavior of the hair surface is simulated by extending the Phong model.
1. Background
Based on the concept of Phong lighting model, it is extended into an empirical model suitable for hair rendering.
2. Motivation
Hair has very unique optical properties, such as specular reflection, subsurface scattering, etc., and the emergence of these phenomena is closely related to the geometric shape and surface structure of hair.
3. Methods
The Kajiya-Kay model is based on the idea of the Phong model and is an extension of the Phong model.
Cylindrical Hair Representation.
[Marschner 2003] Light Scattering from Human Hair Fibers
The originator of hair, no need to say more +1
It is able to capture key visual effects that existing Kajiya-Kay models cannot describe, such as multiple highlights and scattering variations associated with fiber axis rotation.
1. Background
Limitations of the Kajiya-Kay Model assumes that hair is only an opaque cylinder, ignoring key phenomena such as internal reflection and transmission.
2. Motivation
Hair is a dielectric material, and especially light-colored hair (such as blonde, brown, and red) has significant translucency. Therefore, there is a need for a More Accurate Hair Scattering Model.
3. Methods
The 3D full hemispherical light scattering of a single hair was measured.
The Transparent Elliptical Cylinder Model is proposed.
A new combined scattering and diffraction model that simulates light scattering and diffraction phenomena for hair with an elliptical cross-section without pre-calculation.
1. Background
Still with wave optics as the background, when light interacts with objects whose size is close to the wavelength of light, interference and diffraction effects become significant.
Rendering hair requires considering its geometric properties as well as the wave effects of light. While ray tracing can simulate most scattering phenomena, it falls short when it comes to diffraction in hair.
2. Motivation
Addressing Diffraction and Elliptical Cross-sections. A model combining the wave and ray properties of light is proposed to handle the light diffraction phenomenon of hair without pre-calculation. Supports hair fibers with arbitrary elliptical cross-sections.
3. Methods
The ray part (Ray Interaction with Elliptical Fibers) introduces a complete light transport model, continues the traditional ray model, and handles most scattering effects.
The Wave Diffraction by Elliptical Fibers section introduces a new diffraction scattering lobe function that captures the strong forward scattering effect that occurs when light interacts with hair.
The "dual scattering" model is widely used in real-time rendering, and there is no need to explain this classic model.
1. Background
In light-colored dense hair, multiple scattering is a key factor in determining the overall hair color.
Existing methods based on path tracing or photon mapping are too slow to render and often ignore the circular cross-section of hair fibers.
2. Motivation
Need for a Physically Accurate and Efficient Multiple Scattering Model.
3. Methods
Dual Scattering Model, global multiscattering and local multiscattering. The global multiscattering part aims to calculate the light that passes through the hair volume and reaches the neighborhood of the target point, while the local multiscattering considers the scattering events within this neighborhood.
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 In the study of human hair, the model of Marschner [2003] is widely used in the industry. It analyzes the light paths in dielectric cylinders and cones and separates the scattering into R, TT and TRT. Zinke [2009] added a diffuse reflection component. Sadeghi [2010] proposed an artist-controlled parameterization method. d'Eon [2014] and Huang [2022] proposed a non-separable characterization method, that is, there is a coupling effect between the azimuthal and longitudinal angles, which cannot be simply separated. Chiang [2016] further optimized the model to make it suitable for production-level rendering. In addition to human hair, Khungurn and Marschner [2017] explored the modeling of elliptical hair. Yan [2015, 2017] studied animal hair with an internal medulla. Aliaga [2017] generalized the model to textile fibers with more complex cross-sections.
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.
Looking at the structure of the above formula, it is very similar in form to $\mathbf{E}(\mathbf{r}, t) = \mathbf{E}_0 e^{j(\omega t – k \cdot \mathbf{r})}$, which is called the common time-harmonic solution of the electric field.
$ 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.
The historical evolution of hair rendering research
In 1989, Kajiya and Kay extended the Phong model to hair drawing and proposed an empirical model for hair drawing, the Kajiya-Kay model. This model simplifies hair into a series of elongated cylinders and assumes that light is simply reflected on the surface of hair. Specular reflection: highlights. Diffuse reflection: simulates the scattering of light inside the hair and the overall brightness of the hair.
In 2003, Marschner et al. published the paper "Light Scattering from Human Hair Fibers", proposing a physics-based hair reflection model, known as the Marschner model.
Traditional hair rendering models, such as the Marschner model and the Kajiya-Kay model, usually simplify hair into a single-layer cylinder, ignoring the influence of the medulla. In their 2015 paper "Physically-Accurate Fur Reflectance: Modeling, Measurement and Rendering", Yan et al. proposed a more accurate and efficient hair rendering model to address this problem, modeling each hair as two concentric cylinders. The modeling combines the complete hair bidirectional scattering distribution function (BSDF) to accurately describe the multipath propagation and scattering behavior of light in hair. In addition, to ensure physical authenticity, a large amount of physical measurement work was carried out, including nine different hair samples, and finally some reflection parameters of the database were opened for artists to adjust. In order to improve rendering efficiency, Yan [2015] combined the consideration of near-field scattering (R, TT, TRT) of [Zinke and Weber 2007], pre-calculated common scattering paths, stored in Lut, and realized efficient light scattering calculation of single hair fibers. For rendering pipeline optimization, Yan [2015] mainly focused on the dual-cylinder model.
From Marschner[2003] to the extended models of d'Eon[2011] and Chiang[2016], although the continuous increase of hair parameters (such as azimuthal roughness numerical integration of d'Eon[2011] and near-field azimuthal scattering of Chiang[2016]) has increased rendering accuracy, its complex scattering path and large amount of pre-calculation limit its practicality and real-time performance. The double-cylinder hair reflection model proposed by Yan[2015] also has the problems of high computational cost and low practicality. Therefore, Yan[2017] proposed a simplified version, which achieves fast integration through analytical methods and greatly reduces the number of Lobes.
The figure above clearly shows that Marschner [2003] used a longitudinal-azimuthal decomposition representation to simplify the complex three-dimensional light scattering process into two relatively independent dimensions. The longitudinal scattering function describes the propagation and scattering of light along the axis of the hair fiber. The azimuthal scattering function describes the scattering of light in the cross section of the hair fiber (the plane perpendicular to the fiber axis). This model considers T, TT and TRT. The energy conservation problem was corrected in d'Eon [2011]. Yan [2015]'s double cylinder model (hair cuticle and hair medulla) complicated the light interaction and considered R, TrT, TtT and TrRrT. Yan [2017] introduced a unified refractive index (IORs) to simplify the light path propagation and no longer distinguish the refractive indices of different materials, namely R, TT, TRT, TT^s and TRT^s (^s represents the simplified path). Yan[2017] pointed out that unifying IORs does not actually significantly affect the rendering results, and it can still maintain a high degree of realism becauseThe refractive index of the hair cortex and medulla is very close..
In order to solve the problem of high computational complexity in previous hair rendering, Yan[2017] proposedDivision of near field and far field, and introducedLayered Rendering Strategy. The near field is mainly the fine scattering and reflection of light on a single hair. This area requires a high-precision physical model to render hair, such as wave optical phenomena such as interference and diffraction. The far field describes the overall scattering effect of light under the collective action of a large number of hair fibers. In this area, the microscopic structure of a single hair can be averaged, which is more suitable for calculation using statistical/approximate methods to optimize and improve rendering efficiency.Classification criteria: Based on the distance between the light and the hair fiber and the size of the hair fiberThat is, set aThresholdWhen the distance between the light and the hair is less than the threshold, it is classified as near field; otherwise, it is classified as far field.Layered rendering processFirst, ray tracing is used to determine whether it is near field or far field. In the near field, the Mie scattering theory and Fresnel equation are used to calculate reflection and transmission in combination with the pre-calculated scattering table. In the far field, the statistical scattering function + pre-calculated scattering table + MC integral are used to reduce the complexity. Finally, the two are superimposed and normalized. Here is a detailed explanation based on the three points of the paper:
Simple Reflectance ModelAlthough Yan's model [2015] introduced the medulla and considered more physical details, it still has high computational complexity, especially in the conversion between near field and far field. Yan [2017] proposedSimplified hierarchical reflection model, retaining the key physical phenomena of reflection and reducing unnecessary complex light scattering paths. They describe reflection as three main light scattering paths (R: specular reflection of light on the surface of the hair cuticle, TT: light passes through the hair cuticle and is transmitted from the other side after internal scattering, TRT: light enters the hair, reflects once inside, and finally transmits). Finally, combined withSimplified Bidirectional Scattering Distribution Function (BSDF)The reflection of the capture path reduces the number of lobes required in the calculation (usually used to describe the distribution curve of different scattering directions). Compared with the previous model, the number of lobes in the calculation process is reduced from 9 to 5.
Improved Accuracy and PracticalityHigh-precision models (such as the Marschner model, Yan [2015], etc.) require complex numerical calculations and a large amount of pre-calculated data, and are therefore difficult to implement in real-time applications. Low-precision empirical models lack sufficient physical reality. Therefore, many improvements were made in Yan [2017]. Although the model simplifies the light scattering path,Combining physical phenomenaThe model is more accurate than the traditional empirical model by reasonably simplifying the longitudinal and azimuthal scattering of light.Transition processing between near field and far fieldTraditional models often fail to smoothly handle the optical transition between the near field and the far field. Yan et al. introducedA near-field-far-field analysisThe solution accurately simulates reflections when the light is close to the hair fiber, while quickly approximating the overall reflection behavior of the light in the far field. This makes the rendering efficient enough for real-time rendering.
Analytic Near/Far Field SolutionThere is a huge difference in the treatment of the near field (the short-range interaction between light and a single hair fiber, i.e., scattering behavior) and the far field (the long-range collective effect between light and a large number of hair fibers). In order to achieve a seamless transition between the near field and the far field, the authors used an analytical integration method instead of cumbersome numerical integration. The analytical integration can directly calculate the reflection function without the need for complex numerical solutions or pre-calculations, which greatly reduces the calculation time.
Significant Speed Up
Reduce the number of scatter paths used to describeNumber of lobes;
Combination of analytical integration and pre-calculation;
A simplified BSDF and analytical reflection calculation formula are used to combine ray tracing and reflection calculation inParallelization on GPU, the rendering speed of the model is increased by 6-8 times compared with previous methods.
To summarize briefly, the reflection model proposed by Yan [2017] has good effects and performance. By unifying the IOR of the cortex and medulla, the model only needs 5 lobes to represent the complex scattering of fur, and the tensor approximation is used to minimize the storage overhead. Based on this model, the analytical integration of the far and near fields is proposed to extend the model to multi-scale rendering. It is very simple to implement the BCSDF model in real-time rendering. There are already many implementation methods, and it has been applied to the film and television industry. [The Lion King (HD). 2017 movie] (2019 Oscar Nominee for Best Visual Effects)
XIA[2023] proposed a hair reflection model based on wave optics. Traditional hair rendering models are mostly based on geometric optics approximation. These models work well when processing larger hair fibers, but have poor performance on subtle optical phenomena (such as colored spots on hair, i.e.glints). These scattering effects, including reflection, transmission, and multiple scattering, are difficult to accurately describe with simple geometric optics models. As the diameter of the hair fiber approaches or becomes smaller than the wavelength of light (visible light), wave optics effects become increasingly important, and geometric optics models are unable to capture these effects.
Wave optical effects of hair, such asInterference and diffraction of lightThe computational complexity is very high. Wave optics simulation requires calculating the propagation of electromagnetic fields, not just the path of light. Hair and fur have highly irregular microstructures that further affect the scattering of light. Methods based on geometric optics cannot handle these wave phenomena, and full-wave simulations require high computing resources.
As early as XIA[2020], it was proposed to use wave optics to accurately describe the interaction between light and fibers, and to use the boundary element method (BEM) to simulate the fiber scattering of light at any cross section. In addition, XIA[2020] pointed out that due to the diffraction effect, the fiber exhibits an extremely strong forward scattering effect. Therefore, the wave optics effect should focus the light in the direction of forward scattering. It was also pointed out that the small fiber scattering effect depends significantly on the wavelength of light, resulting in strong wavelength scattering. In addition, the singular softening phenomenon brought by the wave field is also the key to determining the real caustic effect. In order to control the amount of calculation of the BEM simulation, the shape of the fiber is ideally a regular cross-sectional shape. However, Marschner[2003] pointed out that the irregularity of the hair surface has an important influence on the appearance of the hair. Whether such an effect is significant in wave optics is still a problem that needs to be explored and solved.
Traditional geometric optics methods are based on ray tracing, which predicts the propagation path of light by simulating the reflection and refraction of light on hair fibers. However, this method is insufficient when dealing with light waves with wavelengths comparable to the fiber size, and cannot capture the effects of diffraction.Complex optical effects are produced. In actual measurements, fiber scattering shows some sharp optical features, which are caused by the diffraction effect of light. Including the slight color shift in black dog hair, which is also caused by the interference and diffraction of light.
In order to deal with these phenomena that cannot be explained by geometric optics, XIA[2023] developed a 3D wave optics simulator based on physical optics approximation (PO) and used GPU to accelerate computational efficiency. The space is processed through an octree structure. The simulator has a certain degree of versatility and can handle arbitrary 3D geometric shapes, that is, it can handle the microstructure of the fiber surface.
However, XIA[2023] points out that it is unrealistic to directly apply this simulator to the current mainstream rendering framework due to the high computational complexity. Therefore, it is necessary to first migrate the model to the existing hair scattering model and then add aDiffraction lobe of elementary diffraction theoryFinally, aRandom ProcessThe modulation method is used to capture the optical speckle effect. Although it is procedural noise, it is still consistent with the physical simulation result, and the visual effect is close to reality.
XIA[2023] divides the current hair/fiber rendering into two types: traditionalRay-based Fiber Models, the other isWave-based Fiber Models.
Linder[2014] proposed an analytical solution to deal with the scattering behavior of cylindrical fibers, but it is only applicable to perfect circular cross-sections and cannot handle complex hair surface structures. XIA[2020] studied the scattering behavior of fibers with arbitrary cross-sectional shapes through two-dimensional wave optics rendering, showing the manifestation of diffraction effects, but the paper assumes a perfect extrusion structure, that is, the fiber surface is regular. Bennamira&Pattanaik[2021] proposed a hybrid model that uses wave optics to solve the problem of only forward diffraction and traditional geometric optics in other scattering modes. However, XIA[2023] further considered theDependence of longitudinal angle of incidence.
At the end of the paper, the procedural noise is fitted to the speckle pattern in wave optics, and a very realistic effect is produced through statistical property fitting.
XIA[2023] also mentionedComputational Electromagnetics ToolsIt plays an important role in dealing with complex interactions of rays and fibers, especially when using numerical methods such as BEM.Computational ElectromagneticsIt is a computational method used to study electromagnetic phenomena. Since light is an electromagnetic wave, many phenomena in optics can be analyzed using electromagnetic tools. CEM is often used in optics to calculate the interaction between light and the surface of an object (such as hair fibers). Common CEM algorithms include:
Finite-Difference Time-Domain (FDTD): A numerical method for solving Maxwell's equations by discretizing them in space and time, first proposed by Kane Yee in 1966. It is a direct time-domain method Kane Yee[1966], Taflove[2005].
Finite Element Method (FEM)：It is used to solve the electromagnetic field distribution in complex geometric structures by dividing the solution area into a finite number of elements Jin[2015].
Boundary Element Method (BEM): Also known as the Method of Moments (MoM), this is a numerical method that reduces the amount of computation by only dealing with the electromagnetic fields on the surface of an object Gibson[2021], Huddleston[1986], Wu[1977].
Although CEM has many acceleration algorithms, such as Song[1997]'s Multilevel Fast Multipole Algorithm (MLFMA), the improvement is still minimal in hair and fur simulation.
Since full-wave simulations are computationally expensive, XIA [2023] proposed the physical optics approximation (PO) to simplify the reflection and diffraction processes on the surface of, for example, hair fibers.
Physical Optics Plane ModelyesAn application of physical optics (PO) that is specifically designed to simulate the behavior of light on flat or nearly flat surfaces.The scattering and diffraction effects on rough surfaces are effectively calculated by Beckmann-Kirchhoff[1987] and Harvey-Shack[1979]. Gaussian random surfaces by He[1991], Kajiya[1985], periodic static surfaces by Stam[1999], and scratched surfaces by Werner[2017] all use physical optics approximations to deal with surface reflection and diffraction.
For more complex diffraction, Krywonos[2006], Krywonos[2011] proposed improved methods for processing diffraction on rough surfaces. Holzschuch, Pacanowski[2017] proposed a dual-scale microsurface model that combines reflection and diffraction to simulate rough surfaces. Recently, Falster[2020] combined Kirchhoff scalar diffraction theory and path tracing to handle secondary reflection and scattering. Yan[2018] used physical optics to render the mirror microgeometry of rough surfaces.
Unlike a flat surface, the fiber surface is a closed curved surface, and the geometric shape of the hair fiber makes the interaction with light more complex. In addition to reflection and scattering, forward diffraction scattering and large-scale shadow effects need to be dealt with.
XIA[2023] also discussedSpeckle effectandProcedural noiseApplication in hair rendering.
Speckle is a grainy image or diffraction pattern produced when light interacts with a rough surface. Its statistical properties have been extensively studied. When coherent light (such as a laser) is irradiated onto a rough surface or passes through a scattering medium, a random pattern of light and dark spots is produced. In layman's terms, it's like when you shine a laser pointer on a rough wall, you see a granular, flickering pattern instead of a smooth spot of light. Because light is scattered at tiny surface irregularities, light waves from different paths interfere with each other, some are strengthened (forming bright spots) and some cancel each other out (forming dark spots), resulting in this speckled pattern.
Previous studies have explored howMonte Carlo methodto simulate the speckle effect in volume scattering Bar[2019, 2020]. However, these models are mainly applicable toHomogeneous media, which is not applicable to heterogeneous structures such as hair fibers. Steinberg, Yan [2022] studied speckle rendering of planar rough surfaces. However, the authors pointed out thatThe speckle effect on fiber surfaces is different from that on flat surfaces, showing different statistical characteristics.
Therefore, XIA[2023] proposedAccurately capture the statistical characteristics of fiber speckle patternsBy studying the special geometric structure and speckle distribution of the fiber surface, the scattering effect of hair fibers is simulated to provide better speckle effects.
It should be noted that although both thin-film interference and speckle effect are caused by the interference of light, they have significant differences in physical mechanism, visual performance and rendering methods in computer graphics. Monte Carlo methods of thin-film interference, such as random film thickness sampling, can be used to generate random spots of speckle effect to improve the realism of rendering. Approximate algorithms such as hierarchical thickness sampling and pre-calculated interference patterns can also learn from each other. Thin-film interference often involves interference of light waves at different scales, and speckle effect also involves multi-scale scattering of microscopic surface structures.
Between the two, thin film interferometry has a relatively low rendering complexity, and pre-computation can be fully utilized to avoid the burden of real-time calculation. However, the speckle effect has highly random and statistical characteristics, and a large number of random interference paths need to be processed, especially for simulating heterogeneous structures such as hair. Current research such as XIA[2023] is working to improve its efficiency, but there is still a large gap compared to thin film interferometry.
XIA[2023] uses the Wavelet band-limited noise of Cook, DeRose[2005] to control the microscopic geometric changes of hair fibers. This noise is different from conventional procedural noise, such as Perlin[1985], Olano[2002], Perlin, Neyret[2001], etc. A significant advantage of Wavelet noise is that itStatistical distributions can be calculated and controlled.
The advantage of the practical wave optics fiber scattering model of XIA[2023] is its realisticColored highlights (glints)Previous geometric optics models usually assume that the fiber surface is a smooth dielectric cylinder, without considering the complex interaction of light waves on the surface irregular structure. In actual tests, the XIA[2023] model performs well in rendering time, can be used in production environments, and generates more delicate and realistic optical effects than traditional models.
XIA[2023] is an important breakthrough that buildsFirst 3D wave optics fiber scattering simulatorPrevious fiber models (including early wave optics models such as Xia et al. 2020) mostly assumed thatLongitudinal and azimuthal directionsThe scattering onSeparable, which greatly simplifies the calculation. However, the authors' simulation results show that the highlights areInseparable, which is a phenomenon that previous models could not accurately handle. The simulator also predictedSpeckle patternsThis is a phenomenon that has not been captured by all previous fiber scattering models based on geometric optics and wave optics.5-dimensional scattering distributionThe method is to use tabulation, which is very memory intensive. Therefore, procedural noise is used to directly replace a five-dimensional table.
XIA[2023] has only been simulated once so farSpeckle Effect in Reflection Mode, higher order reflection modes are still being studied. And light-colored hair may require higher computational requirements to simulate perfectly. The wave optics fiber scattering model used in this study canEasily combined with previous fiber models.
References
Zotero one-click generated, needs to be corrected.
[1] JT Kajiya and TL Kay, “RENDERING FUR WITIt THREE DIMENSIONAL TEXTURES,” 1989.
[2] SR Marschner, HW Jensen, and M. Cammarano, “Light Scattering from Human Hair Fibers,” 2003.
[3] A. Zinke and A. Weber, “Light Scattering from Filaments,” IEEE Trans. Visual. Comput. Graphics, vol. 13, no. 2, pp. 342–356, Mar. 2007.
[4] L.-Q. Yan, C.-W. Tseng, HW Jensen, and R. Ramamoorthi, “Physically-accurate fur reflectance: modeling, measurement and rendering,” ACM Trans. Graph., vol. 34, no. 6, pp. 1–13, Nov. 2015.
[5] L.-Q. Yan, HW Jensen, and R. Ramamoorthi, “An efficient and practical near and far field fur reflectance model,” ACM Trans. Graph., vol. 36, no. 4, pp. 1–13, Aug. 2017.
[6] M. Xia, B. Walter, C. Hery, O. Maury, E. Michielssen, and S. Marschner, “A Practical Wave Optics Reflection Model for Hair and Fur,” ACM Trans. Graph., vol. 42, no. 4, pp. 1–15, Aug. 2023.
Glints Effect Study
Traditional rendering methods based on geometric optics, such as Yan [2014, 2016], useBidirectional Reflectance Distribution Function (BRDF)To simulate the mirror reflection surface, there are certain limitations.
Yan [2014, 2016] pointed out that traditional BRDF models usually use a smooth normal distribution function (NDF), assuming that the microfacets are infinitely small. But in reality, real surfaces often have obvious geometric features, such as micron-level bumps and flakes in metallic paint, which can cause significant glints under strong directional light sources (such as sunlight). Yan et al. simulated these small-scale surface geometric features more accurately through high-resolution normal maps, and proposed a new method to effectively render these complex specular highlights.
Traditional uniform pixel sampling techniques have too large variance when capturing highlights in these small ranges, resulting in low rendering efficiency and inability to effectively handle the uneven distribution of highlights caused by the complexity of the light path. Therefore, Yan [2014, 2016] introduced a search based on normal distribution and targeted sampling.
In hair renderings, you can observe that the hair and fur will show a shimmering effect of changing color when illuminated by strong directional light sources.
XIA[2023] uses optical speckle theory to simulate highlight noise, and adds the diffraction lobe of basic diffraction theory to process the diffraction effect of light on the surface of fiber structures such as hair, thereby rendering colored highlight effects.
XIA[2023], Chapter 8, states that glints can be easily observed in sunlight. Although subtle when viewed from a distance, these color effects can significantly enhance the appearance of the hair when viewed up close, sometimes causing a slight change in the hue of the fiber.
In Figure 9, the model of XIA[2023] also produces colored shimmer effects on light-colored hair. The shimmer is more subtle on light-colored hair than on dark-colored fibers because multiple scattering averages out the colors, resulting in reduced color contrast. Compared to XIA[2020], XIA[2023] not only handles wavelength-dependent reflections better, but also improves its ability to handle the angle of the hair cuticle, capturing the shift in highlights caused by the tilt of the hair cuticle.
References
[1] L.-Q. Yan, M. Hašan, W. Jakob, J. Lawrence, S. Marschner, and R. Ramamoorthi, “Rendering glints on high-resolution normal-mapped specular surfaces,” ACM Trans. Graph., vol. 33, no. 4, pp. 1–9, Jul. 2014.
[2] L.-Q. Yan, M. Hašan, S. Marschner, and R. Ramamoorthi, “Position-normal distributions for efficient rendering of specular microstructure,” ACM Trans. Graph., vol. 35, no. 4, pp. 1–9, Jul. 2016.
This paper discusses the theoretical basis of the physical wave simulation three-dimensional wave optical fiber scattering simulator used to generate high-precision light scattering simulation data in the rendering black dog hair paper "A Practical Wave Optics Reflection Model for Hair and Fur".
Calculating light reflection from rough surfaces is an important topic. Small-scale geometric structures, such as the tiny features of hair fibers, have a significant impact on the reflection behavior of light. The BRDF describes how a surface reflects light given an incident and outgoing direction. The limitations of geometric optics have been repeated many times. This model that treats light as a straight line propagation fails to capture the wave nature of light when the microstructure is close to the wavelength of light.
Theoretical models that use wave optics to approximate diffraction include the Beckmann-Kirchhoff theory of [Beckmann and Spizzichino 1987] and the Harvey-Shack model of [Krywonos 2006]. The former describes the light reflection behavior of rough surfaces, while the latter is a series of models based on wave optics that more accurately describe the scattering behavior of light on complex surfaces.
Existing models are all aimed at the average reflection behavior of large-area surfaces, ignoring local detail changes. Yan [2016, 2018] is able to capture the changes in light reflection from microstructures in different regions of space. Even models based on electromagnetic wave propagation still require certain approximate processing due to computational complexity. These methods are not actually ground truth.
In order to accurately capture the interference effects, XIA[2023] aims to develop a reference simulation tool that simulates the propagation of light faithfully according to Maxwell's equations. The only approximation is the numerical discretization, which ultimately generates the traditional bidirectional reflectance distribution function (BRDF) as output.This simulator truly achieves ground truth.
That is to say, this simulator can accurately simulate the wave characteristics of light, including interference, diffraction, multiple scattering, etc. The approximations used in the simulator are only meshing and numerical integration errors.
Through high-precision full-wave simulation, it is possible to generateHigh angular and spatial resolution BRDF data.
At the same time, the simulator is able to handle large surface areas (such as 60 × 60 × 10 wavelengths). For example, using visible light with a wavelength of about 500 nanometers, 60 wavelengths is equivalent to 30 microns. In other words, the simulator's calculations are based on the scale of light wavelengths.Real physical sizeIn this case, light of different wavelengths will correspond to different numbers of discretized units.The larger the wavelength(For example, the wavelength of red light is longer than that of blue light). For the same physical size, the required discretization units (such as mesh division) will beRelatively less, so the amount of computation required will beRelatively smaller, processing speed may alsoFaster.
Specifically, the surface is represented as a height field, each grid point corresponds to a height value, and quadrilaterals are used as primitives.
For the scattered field, the boundary integral formulation is used to transform the scattering problem of electromagnetic waves into an integral equation that is solved only on the surface boundary. The key implementation method is the boundary element method (BEM). The adaptive integral method (AIM) based on the three-dimensional fast Fourier transform (3D FFT) is then used to accelerate the calculation process of the boundary integral.
And use GPU to accelerate the parallel processing of large-scale surface scattering problems.
And the paper uses a combined small-scale simulation results to characterize the surface bidirectional scattering behavior.
Related Work
Reflection model based on wave optics
The old-fashioned geometric optics vs. wave optics. This article mainly compares surface scattering models. Classical models of geometric optics include: Cook-Torrance model [Cook and Torrance 1982], Oren-Nayar model [Michael 1994]. In wave optics, physical optics approximations are mainly used to simplify the full wave equation. That is, the first-order approximation (single scattering) in the black dog is used to estimate surface reflection. Classical models includeBeckmann-Kirchoff theoryandHarvey-Shack Model, which use approximate equations in scalar form to model wave optics effects. They are widely used to estimate reflectance on various surface types, such asGaussian random surface,Periodic surfaceEtc. However, the calculation results of these methods are often spatial average results, and it is impossible to perform high-resolution detail reflection.
Gaussian random surface models of He et al. (1991) and Lanari et al. (2017).
Periodic surface models by Dhillon et al. (2014), Stam (1999), and Toisoul and Ghosh (2017).
Multilayer planar surface model of Levin et al. (2013).
Surface data table model of Dong et al. (2016).
Study of scratched surfaces by Werner et al. (2017).
In addition, physical optics approximation is also used to estimate theSpace changing appearance,For example:
Surface data table from Yan et al. (2018)
Random surface models of Steinberg and Yan (2022).
Some hybrid surface models apply physical optics models to some surface components (such as roughness at small scales), while using geometric optics models for larger scales. Applications of these hybrid models include:
Surface roughness models by Falster et al. (2020) and Holzschuch and Pacanowski (2017).
Thin-film interference model of Belcour and Barla (2017).
Suspended particle model by Guillén et al. (2020).
In addition, physical optics models are used to handle inter-surface effects at longer distances. For example:
Studies by Cuypers et al. (2012) and Steinberg et al. (2022) explored these long-range effects.
Scattering methods based on wave optics, such asLorenz-Mie theoryandT-Matrix Method, which is also usedVolumetric ScatteringCalculations, for example:
Theories of Bohren and Huffman (2008) and Mishchenko et al. (2002).
Application of volume scattering by Frisvad et al. (2007) and Guo et al. (2021).
In addition, complex-valued ray tracing techniques proposed by Sadeghi et al. (2012) and Shimada and Kawaguchi (2005) have been applied to rendering natural phenomena and structural color effects.Long-range effects between surfacesandVolumetric ScatteringSuch issues are currently beyond the scope of this study.
Numerical Methods in Computational Electromagnetism (CEM)
There are many methods for numerical calculations:
Oskooi et al. (2010) proposed a numerical method based on difference solution of Maxwell's equationsFinite Difference Time Domain (FDTD)FDTD has been used to predict the appearance of wavelength-scale structures (e.g. Auzinger et al. (2018), Musbach et al. (2013)). However, the overhead is considerable as the simulation area increases!
Finite Element Method (FEM)It is a widely used numerical method for solving partial differential equations, and can also be used in electromagnetic problems. It solves the problem by discretizing the simulation domain in three dimensions. Similar to FDTD, the amount of calculation is too large, not to mention for real-time rendering.
Gibson (2021) provides a detailedBoundary Element Method (BEM)The main advantage of BEM is that it reduces the dimensionality of discretization by converting the scattering problem into an integral equation on the surface of the object. In FDTD and FEM, the entire three-dimensional space needs to be discretized, while BEM only needs to discretize the surface of the object, which significantly reduces the dimensionality and complexity of the calculation.
The paper chose BEMThe main reason is its scalability, which is conducive to the processing of complex surface structures.
There are many ways to speed up BEM:
Liu and Nishimura (2006), White and Head-Gordon (1994)Fast Multipole Method (FMM).
Bleszynski et al. (1996) proposed a three-dimensional fast Fourier transform (3D FFT)Adaptive Integration Method (AIM).
Liao et al. (2016), Pak et al. (1997)Sparse Matrix Canonical Grid Method (SMCG).
Thesis selected AIMThe reason is that AIM is suitable for processing an area with a relatively small axial size.
result
BRDF value isHemisphereThe standardSpectral data to XYZ to RGB conversion, a colored BRDF map is generated.
That is, as the height field resolution increases, the BRDF output gradually stabilizes.8 samples per wavelengthThe resolution is sufficient to produce accurate results.
By comparing with existing wave optics models,The simulator in this articleIt has the highest accuracy and can handle complex optical phenomena and geometric structures. It is suitable for scenes with high precision requirements and is suitable for multiple reflections, interference, and complex surfaces. However, the computational cost is high, which is a compromise between efficiency and accuracy.
OHS and GHS ModelsThe calculation is simple and suitable for smooth surfaces and medium roughness surfaces, but the error is large at large incident angles and on complex surfaces. GHS has improved accuracy at large angles compared to OHS.
Kirchhoff modelThe accuracy is relatively high, but it can only be maintained within the first order range.
Cutting Plane MethodComputationally efficient, suitable for relatively simple surface geometries. Not so for complex ones.
The accuracy of this article is the best and is suitable for high-precision scenarios.
Comparison of coherent regions. As the illumination coherence increases, the resolution and detail of the BRDF becomes richer. When the coherence is high, the BRDF contains more high-resolution details.
In addition, the paper also shows the method used to accelerate BRDF calculation.Beam steeringAs shown in Figure 14, the surface isSpecular Reflection, while in the other directionRetroreflective effectThe paper calculates the BRDF values for a series of gradually changing incident angles, as shown in the figure below.
The BRDF image in each incident direction is reduced to a thin line segment. The comparison in Figure 15 shows that Tangent Plane cannot accurately model the surfaceSecond-order reflection(That is, light that is emitted after multiple reflections.) However, if the surface is smooth, using Tangent Plane is still very fast and accurate.
Furthermore, the paper compares the simulator results with BRDF measurements of real surfaces, especially inMultiple reflection effectAs shown in Figure 17, the top is the actual measurement, the middle is the theoretical model of the paper, and the bottom is the Tangent Plane.
You may ask why there is such a big difference? The paper uses an idealized geometric model, while the surface in the experiment may have some slight geometric deviations, which may affect the accuracy of the reflection. In a nutshell, the simulator in the paper can describe higher-order reflections!
Future Work
When dealing with more complexStructured surfaceWhen the light is scattered on the surface, the simulator can more accurately simulate the propagation and scattering behavior of light on the surface.
The future direction of work is of course to reduce computational overhead while ensuring accuracy.
Then the processing area of this BRDF approximation model is expanded.
At the same time, the simulator in this paper can be used as a benchmark reference.
References
A Full-Wave Reference Simulator for Computing Surface Reflectance
Petr Beckmann and Andre Spizzichino. 1987. The scattering of electromagnetic waves from rough surfaces. Artech House.
Andrey Krywonos. 2006. Predicting Surface Scatter using a Linear Systems Formulation of Non-Paraxial Scalar Diffraction. Ph. D. Dissertation. University of Central Florida.
Ling-Qi Yan, Miloš Hašan, Bruce Walter, Steve Marschner, and Ravi Ramamoorthi. 2018. Rendering Specular Microgeometry with Wave Optics. ACM Trans. Graph. 37, 4 (2018).
Ling-Qi Yan, Miloš Hašan, Steve Marschner, and Ravi Ramamoorthi. 2016. Positionnormal distributions for efficient rendering of specular microstructure. ACM Transactions on Graphics (TOG) 35, 4 (2016), 1–9.
EN 
CN