I wanted to improve the realism of ray traced images.
Real light waves are not conveniently split into RGB components, but rather a spectrum. Further, light exhibits polarisation and diffraction.
Spectral data for light sources is readily available as samples plotted as intensity against wavelength. Storing that as a simple lookup table seems to be relatively simple.
Spectral data for surfaces of materials is somewhat harder to find, but available for some materials. Usually this data is measured using a particular light source, which probably means it has to be divided by that light source sample by sample to get the actual reflectivity. The spectral data can then multiplied sample by sample with the spectrum data for the incoming ray of light.
Missing sample points in spectral data can be interpolated or extrapolated.
Light from an LED is not polarised but is monochromatic. Its spectrum would be zero everywhere except for the one wavelength it has. Treating a spectrum as a histogram is possibly more useful than treating it as a line.
class Spectrum {
val sampleWavelengths: Array<Double>
val sampleValues: Array<Double>
}It is common to model the polarisation of light using Stokes vectors and Mueller calculus.
This project defines its light in C as Stokes values for R,G,B along with a reference X vector.
This points to a structure something like this:
class Light(
val sampleWavelengths: List<Double>,
val sampleStokes: List<Vector4>,
// Equal to the tangent vector at the last hit point
val referenceX: Vector3
)Fresnel equations describe the behaviour of light polarised in two directions.
- s-polarised light (perpendicular to the plane of incidence)
- p-polarised light (parallel to the plane of incidence)
For complex index of refraction ñ = n + κi and angle of incidence θ,
a² + b² - 2acosθ + cos²θ
F_s = ────────────────────────
a² + b² + 2acosθ + cos²θ
a² + b² - 2asinθtanθ + sin²θtan²θ
F_p = ───────────────────────────────── F_s
a² + b² + 2asinθtanθ + sin²θtan²θ
⎧ 2bcosθ ⎫
δ_s = arctan⎪───────────────⎪
⎩cos²θ - a² - b²⎭
⎧2cosθ((n² - κ²)b - (2nκ)a⎫
δ_p = arctan⎪─────────────────────────⎪
⎩(n² + κ²)²cos²θ - a² - b²⎭
Where
__________________________
2a² = ⎷(n² - κ² - sin²θ)² + 4n²κ² + n² - κ² - sin²θ
__________________________
2b² = ⎷(n² - κ² - sin²θ)² + 4n²κ² - n² + κ² + sin²θ
As a Mueller matrix,
⎡ A B 0 0 ⎤
M = ⎢ B A 0 0 ⎥
⎢ 0 0 C S ⎥
⎣ 0 0 -S C ⎦
Where
F_s + F_p
A = ─────────
2
F_s - F_p
B = ─────────
2
_______
C = cos(δ_s - δ_p)⎷F_s F_p
_______
S = sin(δ_s - δ_p)⎷F_s F_p
Determining the new X at the hit:
val incomingRefX = computeX(normal, ray.direction)
Vector3 computeX(y: Vector3, z: Vector3) {
return (y cross z).normalize()
}How to rotate Stokes data once you have the new reference frame:
// second parameter is `referenceX` computed earlier
// third parameter is `reflectionRay.direction`
fun rotateReferenceFrame(light: Light, newReferenceX: Vector3, dir: Vector3): Light {
val dotX = light.referenceX dot newReferenceX
val detX = dir dot (light.referenceX cross newReferenceX)
val phi = atan2(detX, dotX)
val c2p = cos(2.0 * phi)
val s2p = sin(2.0 * phi)
val newSampleStokes = light.sampleStokes
.map { stokes -> rotateStokes(stokes, c2p, s2p) }
return Light(light.sampleWavelengths, newSampleStokes, newReferenceX)
}
fun rotateStokes(stokes: Vector4, c2p: Double, s2p: Double): Vector4 {
return Vector4(
stokes.x,
c2p * stokes.y + s2p * stokes.z,
-s2p * stokes.y + c2p * stokes.z,
stokes.w)
}Materials stuff: the index of refraction is a complex number, n + i k
- n is the classical refractive index we're most familiar with
- k is the extinction coefficient, indicating the amount the light is attenuated passing through the material
How to compute a Mueller Matrix for a hit:
fun fresnelMuellerMatrix(refractiveIndex: Complex, sinTheta: Double, cosTheta: Double, tanTheta: Double): Matrix4x4 {
val n = refractiveIndex.real
val k = refractiveIndex.imaginary
val n2 = n * n
val k2 = k * k
val st2 = sinTheta * sinTheta
val left = sqrt((n2 - k2 - st2) * (n2 - k2 - st2) + 4.0 * n2 * k2)
val right = n2 - k2 - st2
val a2 = 0.5 * (left + right)
val b2 = 0.5 * (left - right)
val a = sqrt(a2)
val b = sqrt(max(b2, 0.0))
val ct2 = cosTheta * cosTheta
// orthogonal
val ortA = a2 + b2 + ct2;
val ortB = 2.0 * a * cosTheta;
// parallel
val parA = a2 + b2 + st2 * tanTheta * tanTheta
val parB = 2.0 * a * sinTheta * tanTheta
// Fresnel parameters
val F_ort = (ortA - ortB) / (ortA + ortB)
val F_par = ((parA - parB) / (parA + parB)) * F_ort
val D_ort = atan2(2.0 * b * cosTheta, ct2 - a2 - b2)
val D_par = atan2(2.0 * cosTheta * ((n2 - k2) * b - 2 * n * k * a), (n2 + k2) * (n2 + k2) * ct2 - a2 - b2)
val phaseDiff = D_ort - D_par
// Matrix components
val A = 0.5 * (F_ort + F_par)
val B = 0.5 * (F_ort - F_par)
val C = cos(phaseDiff) * sqrt(F_ort * F_par)
val S = sin(phaseDiff) * sqrt(F_ort * F_par)
return Matrix4x4( A, B, 0.0, 0.0,
B, A, 0.0, 0.0,
0.0, 0.0, C, S,
0.0, 0.0, -S, C)
}So when the ray hits a new surface:
- Rotate the Stokes data to the hit position
- Create the Mueller matrix for the hit
- Multiply M x S to get the new Stokes vector for the outgoing reflected ray.
One source had the transmission angle for complex index of refraction:
⎧ sinθi ⎫
θt = arctan⎪──────────────────⎪
⎩ Re √(ñ²-sin²θi) ⎭
This source has:
⎧ Re[ñ₂ sinθt] ⎫
θt_real = arctan⎪──────────────⎪
⎩ Re[ñ₂ cosθt] ⎭
(where complex angle sinθt = ñ₁/ñ₂ sinθi)
In another source, we have:
θt = arcsin(sinθi/n′)
n′²= ½a[1 ± √(1 - 4b/a²)] (choose positive result)
where a = n² - k² + sin²θi
b = (n² - k²)sin²θi - n²k²
k′ = nk/n′cosθt
Still not sure:
- What about the refracted ray? Intuitively, you would usually subtract the amount of light reflected from the amount of incoming light, but is this simple for Stokes vectors?