Kugelflächenfunktionen (spherical harmonics) bilden eine orthogonale Basis zur Darstellung richtungsabhängiger Größen. Als Eigenfunktionen des Laplace-Operators auf der Kugel spielen sie eine wichtige Rolle bei der Lösung partieller Differentialgleichungen. Typische Beispiele sind Gravitations- und Magnetfelder, Schallfelder oder die Verteilung von Strahlung.
Einige technische Aufgaben sind jedoch nicht auf der ganzen Kugel definiert. An einem Punkt einer Oberfläche sind häufig nur die Richtungen oberhalb der Tangentialebene relevant. Diese Richtungen bilden eine Hemisphäre. Hemisphärische Kugelflächenfunktionen (hemispherical harmonics, HSH) liefern dafür eine an den Definitionsbereich angepasste orthonormale Basis.
Motivation: Warum eine eigene Basis für die Hemisphäre?
Klassische Kugelflächenfunktionen sind auf der vollständigen Kugel mit Polarwinkel \(\theta\in[0,\pi]\) definiert. Werden dagegen nur Richtungen betrachtet, die von einer ebenen Bezugsfläche aus in eine Raumhälfte zeigen, genügt die obere Halbkugel. Ein typisches Beispiel sind einfallende oder ausgehende Strahlungsrichtungen.
\[
S_+^2=
\left\{
\mathbf{x}=(x_1,x_2,x_3)\in\mathbb{R}^3
\;\middle|\;
x_1^2+x_2^2+x_3^2=1,\quad x_3\geq0
\right\}.
\]
In Kugelkoordinaten wird diese Halbkugel durch \(\theta\in[0,\pi/2]\) und \(\phi\in[0,2\pi)\) beschrieben. Eine Fortsetzung der Daten auf die untere Halbkugel wäre künstlich und kann zusätzliche Diskontinuitäten oder unnötig hohe Frequenzanteile erzeugen. Hemisphärische Kugelflächenfunktionen berücksichtigen diesen eingeschränkten Definitionsbereich unmittelbar. Im Gegensatz zu klassischen Spherical Harmonics, die auf der gesamten Oberfläche der Einheitssphäre definiert sind, beschreiben sie Funktionen nur auf der oberen Halbkugel. Dabei handelt es sich nicht lediglich um eingeschränkte Kugelflächenfunktionen: Die Legendre-Funktionen werden so verschoben, dass die Basis bezüglich des Flächenmaßes auf der Hemisphäre orthonormal ist.
Von Spherical Harmonics zu Hemispherical Harmonics
Seien \(l\in\mathbb{N}_0\) der Grad, \(m\in\mathbb{Z}\) die Ordnung mit \(|m|\leq l\), \(\phi\in[0,2\pi)\) der Azimutwinkel und \(\theta\in[0,\pi/2]\) der Polarwinkel. Für \(x=\cos\theta\) gilt auf der oberen Halbkugel \(x\in[0,1]\). Durch die lineare Transformation
\[
x\longmapsto 2x-1
\]
wird dieses Intervall auf \([-1,1]\), den Definitionsbereich der Legendre-Polynome, abgebildet. Wir definieren daher die verschobenen assoziierten Legendre-Funktionen
\[
\widetilde{P}_l^m(x)=P_l^m(2x-1),\qquad x\in[0,1].
\]
Diese Konstruktion geht auf Gautron, Křivánek, Pattanaik und Bouatouch zurück[1]. Sie wurde insbesondere für die kompakte Darstellung hemisphärischer Daten in Rendering und globaler Beleuchtung entwickelt.
Definition der Basisfunktionen
Eine reellwertige HSH-Basis ist durch
\[
H_l^m(\theta,\phi):=
\begin{cases}
\sqrt{2}\,K_l^m\cos(m\phi)\,\widetilde{P}_l^m(\cos\theta), & m>0,\\
\sqrt{2}\,K_l^{|m|}\sin(|m|\phi)\,\widetilde{P}_l^{|m|}(\cos\theta), & m<0,\\
K_l^0\,\widetilde{P}_l^0(\cos\theta), & m=0
\end{cases}
\]
gegeben. Die Normierungsfaktoren lauten
\[
K_l^m=\sqrt{\frac{2l+1}{2\pi}\frac{(l-|m|)!}{(l+|m|)!}},
\qquad K_l^0=\sqrt{\frac{2l+1}{2\pi}}.
\]
Für einen maximalen Grad \(n\) verwenden wir die Indexmenge
\[
\mathcal{I}_{\omega}(n):=\{(l,m):0\leq l\leq n,\ -l\leq m\leq l\}.
\]
Sie enthält
\[
|\mathcal{I}_{\omega}(n)|=\sum_{l=0}^{n}(2l+1)=(n+1)^2
\]
Basisfunktionen. Die Zahl der Koeffizienten wächst somit quadratisch mit dem maximalen Grad.
Orthogonalität und Normierung
Auf der Hemisphäre \(\mathcal{H}^2\) verwenden wir das Skalarprodukt
\[
\langle f,g\rangle_{\mathcal{H}^2}
=\int_0^{\pi/2}\int_0^{2\pi}
f(\theta,\phi)g(\theta,\phi)\,d\phi\,\sin\theta\,d\theta.
\]
Die HSH erfüllen die Orthonormalitätsbedingung
\[
\langle H_l^m,H_k^q\rangle_{\mathcal{H}^2}
=\delta_{lk}\delta_{mq}.
\]
Aufgrund der Orthonormalität wird jeder Koeffizient durch eine unabhängige Projektion auf die zugehörige Basisfunktion bestimmt. Der bis zum Grad \(n\) abgeschnittene Ansatz ist zugleich die bezüglich der \(L^2\)-Norm beste Approximation in diesem Ansatzraum.
Für die vollständige Basis gilt zudem die Parseval-Beziehung:
\[
\|f\|_{L^2(\mathcal{H}^2)}^2
=\sum_{l=0}^{\infty}\sum_{m=-l}^{l}|f_l^m|^2.
\]
Die Energie der verworfenen Koeffizienten entspricht dem quadratischen \(L^2\)-Fehler des abgeschnittenen Ansatzes.
Reihenentwicklung und Approximation
Für \(f\in L^2(\mathcal{H}^2)\) lautet die HSH-Reihenentwicklung
\[
f(\theta,\phi)=\sum_{l=0}^{\infty}\sum_{m=-l}^{l}
f_l^m H_l^m(\theta,\phi),
\]
wobei die Koeffizienten durch
\[
f_l^m=\int_0^{\pi/2}\int_0^{2\pi}
f(\theta,\phi)H_l^m(\theta,\phi)
\,d\phi\,\sin\theta\,d\theta
\]
bestimmt werden. In einer Rechnung wird die unendliche Reihe bei einem maximalen Grad \(n\) abgeschnitten:
\[
f^{(n)}(\theta,\phi)=\sum_{l=0}^{n}\sum_{m=-l}^{l}
f_l^m H_l^m(\theta,\phi).
\]
Niedrige Grade beschreiben den Mittelwert und grobe Winkelstrukturen. Höhere Grade erfassen feinere Änderungen. Glatte Funktionen benötigen daher meist nur wenige Bänder, während schmale Maxima oder Unstetigkeiten größere Bandgrenzen erfordern.
Ein einfaches Beispiel ist \(f(\theta,\phi)=\cos\theta\). Wegen
\[
H_0^0=\frac{1}{\sqrt{2\pi}},\qquad
H_1^0=\sqrt{\frac{3}{2\pi}}(2\cos\theta-1)
\]
ergibt sich unmittelbar die exakte Darstellung
\[
\cos\theta=\frac{\sqrt{2\pi}}{2}H_0^0
+\frac{1}{2}\sqrt{\frac{2\pi}{3}}H_1^0.
\]
Die Kosinusverteilung liegt damit vollständig im HSH-Ansatzraum ersten Grades. Alle Koeffizienten höheren Grades verschwinden exakt. Diese Darstellung eignet sich zugleich als einfacher Testfall für eine numerische HSH-Transformation: Abgesehen von Quadratur- und Rundungsfehlern dürfen nur die Koeffizienten zu \(H_0^0\) und \(H_1^0\) von null verschieden sein.
Berechnung der Koeffizienten und Rekonstruktion
Bei der Analyse werden aus den Werten einer Funktion ihre HSH-Koeffizienten bestimmt. Dafür sind Integrationen über beide Winkel erforderlich. In azimutaler Richtung hängt die Basis nur von Sinus und Kosinus des Winkels \(\phi\) ab. Liegen die Werte dort auf einem gleichmäßigen Raster, kann dieser Teil effizient mit einer eindimensionalen FFT berechnet werden.
Für die anschließende Integration über den Polarwinkel verwendet man die Substitution \(x=2\cos\theta-1 \). Sie bildet den benötigten Bereich auf \([-1,1]\) ab. Das verbleibende Integral lässt sich dort beispielsweise mit einer Clenshaw-Curtis-Quadratur auswerten.
Die umgekehrte Aufgabe heißt Synthese: Aus den HSH-Koeffizienten werden wieder Funktionswerte für vorgegebene Richtungen berechnet. Die dabei auftretenden Summen der assoziierten Legendre-Funktionen können stabil mit der Clenshaw-Rekursion ausgewertet werden. FFT, Quadratur und Rekursion übernehmen somit jeweils einen klar abgegrenzten Teil der Berechnung.
Struktur der Basisfunktionen in den Plots
Die Abbildungen der ersten fünf Bänder zeigen, wie Grad \(l\) und Ordnung \(m\) die Winkelstruktur bestimmen. Der Grad \(l\) steuert im Wesentlichen die maximal darstellbare räumliche Frequenz: Mit wachsendem \(l\) entstehen mehr Nullstellen und feinere Strukturen in polarer Richtung. Die Ordnung \(|m|\) bestimmt die Zahl der periodischen Änderungen entlang des Azimutwinkels \(\phi\).
Eine druckfähige PDF-Fassung liegt ebenfalls vor.
Für \(m=0\) sind die Funktionen rotationssymmetrisch um die Oberflächennormale. Positive Ordnungen enthalten Kosinusanteile, negative Ordnungen Sinusanteile. Die farblich getrennten positiven und negativen Funktionswerte machen Vorzeichenwechsel und Knotenlinien sichtbar. Ein vollständiges Band zum Grad \(l\) besteht aus \(2l+1\) Funktionen.
Anwendungen in Strahlung und Rendering
Wichtige Anwendungen sind:
- Radiative Wärmeübertragung: Projektion richtungsabhängiger Strahlungsintensitäten und Reflexionseigenschaften.
- Rendering und globale Beleuchtung: Darstellung von Reflexionsfunktionen, Environment Lighting und Radiance Caches.
- Richtungsverteilungen: Kompression und Auswertung winkelabhängiger Daten aus Messungen oder Simulationen.
Die Radiosity-Gleichung setzt diffuse Oberflächen voraus und besitzt deshalb keine explizite Abhängigkeit von der ausgehenden Richtung. In der allgemeineren Rendering Equation kann dagegen \(L(\mathbf{x},\omega)\) für jedes Panel durch HSH-Koeffizienten approximiert werden. Diese Richtungsdiskretisierung lässt sich mit räumlichen Verfahren wie Panel Clustering kombinieren.
SageMath-Implementierung
Die folgenden SageMath-Funktionen implementieren Normierungsfaktor, verschobene Legendre-Funktion und reelle HSH-Basis. Grad und Ordnung können anschließend unmittelbar für Auswertung oder Plot gewählt werden.
def K_lm(l, m):
m = abs(m)
return sqrt((2*l + 1)/(2*pi) * factorial(l-m)/factorial(l+m))
def P_shift(l, m, x):
return gen_legendre_P(l, abs(m), 2*x - 1)
def hsh(l, m, theta, phi):
radial = P_shift(l, m, cos(theta))
if m > 0:
return sqrt(2)*K_lm(l, m)*radial*cos(m*phi)
if m < 0:
return sqrt(2)*K_lm(l, m)*radial*sin(abs(m)*phi)
return K_lm(l, 0)*radial
var("theta phi")
H = hsh(4, 2, theta, phi)Für einen Plot auf der Halbkugel kann der positive Anteil von \(H\) beispielsweise blau und der Betrag des negativen Anteils grün dargestellt werden. Wichtig ist stets die zulässige Indexbedingung \(-l\leq m\leq l\). Ein interaktives SageMath-Skript und das Projekt Hemispherical Harmonics auf GitHub erzeugen die Basisfunktionen für frei wählbare Werte von \(l\) und \(m\).
Literatur
Quellen
- (2004): A Novel Hemispherical Basis for Accurate and Efficient Rendering. In: Keller, Alexander; Jensen, Henrik Wann (Hrsg.): Eurographics Workshop on Rendering, The Eurographics Association, 2004, ISSN: 1727-3463.