Wir gehen von einer Menge
Γ ={x1, . . . , xN}
gestreuter Punkte xj ∈ Rd mit d ∈ N≥2 aus. Dabei sei Γ in einer Menge Ω ⊂ Rd enthalten, welche sich als Vereinigung endlich vieler beschr¨ankter, zusammenh¨angender Teilmengen des Rd schreiben l¨asst. Weiterhin sind uns f¨ur j = 1, . . . , N die gest¨orten Abtastwerte
f˜(xj) = f(xj) +zj
gegeben. Dabei seien f : Ω → R eine st¨uckweise stetige Funktion und zj unabh¨angige normalverteilte Zufallsvariablen mit Erwartungswert 0 und unbekannter Varianz σ2. Die
abgetasteten Funktionswerte sind demnach mit additivem Gaußschen Rauschen belegt.
Ziel ist eine Entst¨orung der Daten, d.h. eine m¨oglichst gute Rekonstruktion der originalen Abtastwerte f(xj) f¨urj = 1, . . . , N.
Wir nehmen an, dass Γ quasi-uniform gem¨aß der folgenden Definition (siehe [10]) ist.
Definition 2.61 Erf¨ullen die maximale Dichte (engl. maximal density) δ(Γ) := max
x∈Ω min
xk∈Γkx−xkk2 und der minimale Abstand (engl. minimal spacing)
µ(Γ) := min
mit einer a priori fixierten KonstantenC passender Gr¨oße (siehe auch Bemerkung 2.62), so heißt Γ quasi-uniform in Ω.
Bemerkung 2.62 1. Die Quasi-Uniformit¨at stellt sicher, dass die Menge Ω ¨uberall
”gen¨ugend dicht“ mit Abtastpunkten xj ∈ Γ abgedeckt ist. Die maximale Dichte δ(Γ) gibt den maximalen Fehler an, den man bei der Bestapproximation eines Punktes x∈Ω durch einen der Abtastpunkte xj ∈Γ begeht.
2. Betrachten wir im Fall √d
N ∈N etwa den d-dimensionalen W¨urfel Ω = [1,√d
3. Die Autoren in [10] verwenden C = 2 unabh¨angig von der Raumdimension d in der Definition der Quasi-Uniformit¨at. In diesem Fall ist ein uniformes Gitter f¨ur d ≥ 16 per Definition nicht mehr quasi-uniform (siehe vorherigen Punkt dieser Bemerkung). Hier sei daher C > √
d/2. Es kann sinnvoll sein, auch f¨ur d < 16 Werte f¨ur C zuzulassen, die etwas gr¨oßer als 2 sind.
Der im Folgenden diskutierte und von uns in [66] neu vorgeschlagene Entst¨ orungsalgo-rithmus bedient sich einer klassischen Shrinkage-Prozedur, wobei die Wavelet-Transformation entlang zu konstruierender Pfade erfolgt, und stellt somit eine Verall-gemeinerung der EPWT [101] dar. In Anlehnung an das Cycle-Spinning [31] verwenden
wir eine Durchschnittsbildung ¨uber mehrere mittels unterschiedlicher Pfade entst¨orte Datens¨atze, um das Resultat der Entst¨orung zu verbessern.
Zusammenfassen l¨asst sich unser neuer Entst¨orungsalgorithmus wie folgt (vgl. [66]).
Algorithmus 2.63 Entst¨orung mittels Wavelets entlang von Pfaden
Gegeben: Γ ={x1, . . . , xN}=:{x01, . . . , x0N}=: Γ0 ⊂Rd, f˜(xj) =:c0j f¨ur j = 1, . . . , N, wobei N = 2t f¨ur ein t ∈ N und eine biorthogonale Wavelet-Filterbank mit Zerlegungs-filtern h, g und Rekonstruktionsfiltern h,˜ g˜ mit P
k∈Zhk =√
2 sowie ein Tiefpass-Filter γ mit P
k∈Zγk= 1.
Zerlegung: F¨uhre die folgenden vier Schritte f¨ur l = 0,1, . . . , L−1 mit L < t durch:
1. Konstruiere einen geeigneten Pfadvektor pl ∈ RN/2
l, d.h. eine Permutation von {1, . . . , N/2l}, welche eine Reihenfolge der Punktexlpl(j) und der zugeordneten Wer-te clpl(j) darstellt.
2. Wende den (periodischen) Tiefpass-Filter (¯h−k)k∈Z bzw. den (periodischen) Hoch-pass-Filter (¯g−k)k∈Z auf den Vektor (clpl(j))N/2j=1l an und erhalte nach Downsampling die Tiefpass-Koeffizienten (cl+1j )N/2j=1l+1 bzw. die Hochpass-Koeffizienten(dl+1j )N/2j=1l+1 mit
cl+1j =X
k∈Z
clpl(k)¯hk−2j+1
dl+1j =X
k∈Z
clpl(k)¯gk−2j+1.
3. Wende den Tiefpass-Filter (¯γ−k)k∈Z auf den Vektor (xlpl(j))N/2j=1l (separat f¨ur jede der d Komponenten) an und erhalte nach Downsampling neue gestreute Punkte (xl+1j )N/2j=1l+1 mit
xl+1j =X
k∈Z
¯
γk−2j+1xpl(k). Setze Γl+1 :={xl+11 , . . . , xl+1N/2l+1}.
4. F¨uhre ein Hard-Shrinkage mit einem Parameter θ >0 an den Hochpass-Koeffizi-enten (dl+1j )N/2j=1l+1 resultierend in
d˜l+1j =
(dl+1j falls |dl+1j | ≥θ 0 falls |dl+1j |< θ durch.
Rekonstruktion:Initialisiere(˜cLj)N/2j=1L := (cLj)N/2j=1L und f¨uhre die folgenden drei Schritte f¨ur l =L, L−1, . . . ,1 durch:
5. Wende Upsampling und anschließend den Tiefpass-Filter ˜h auf (˜clj)N/2j=1l an.
6. Wende Upsampling und anschließend den Hochpass-Filter g˜auf ( ˜dlj)N/2j=1l an.
7. Addiere die Ergebnisse von Schritt 5 und 6, um
(˜cl−1pl−1(j))N/2j=1l−1 = X
zu erhalten. Invertiere die Permutation mit dem Resultat (˜cl−1j )N/2j=1l−1.
Ausgabe:Den gestreuten Punkten xj ∈Γzugeordnete gegl¨attete Funktionswerte(˜c0j)Nj=1. Iteration: Wiederhole die Schritte 1 bis 7 (unter Verwendung neuer Pfadvektoren) T -mal (z.B. T = 64) und bilde den Durchschnitt der Ausgaben (˜c0j)Nj=1.
Bemerkung 2.64 1. Die Konstruktion geeigneter Pfade ist wesentlich f¨ur den Er-folg des Algorithmus. Wir werden im n¨achsten Abschnitt 2.4 erl¨autern, wie wir die Pfade festlegen. Man beachte, dass wir in jedem Zerlegungslevel einen neuen Pfad bestimmen. Dadurch ist die im obigen Entst¨orungsalgorithmus verwendete Transformation nicht echt eindimensional.
2. Schritt 3 des Algorithmus aktualisiert in jedem Zerlegungslevel die zugrundeliegen-de Punktmenge Γl durch Γl+1. Durch das Downsampling findet eine Ausd¨unnung in Form einer Halbierung der Kardinalit¨at der Punkmenge statt, d.h.
|Γl+1|= 1
2|Γl|= N 2l+1. Die Wahlγ = √1
2h ist naheliegend, aber andere Tiefpass-Filter sind m¨oglich. Man beachte, dass, wenn Γ = Γ0 ⊂ Ω ⊂ Rd ist, f¨ur l ≥ 1 trotz der Normierung P
k∈Zγk = 1 nicht notwendigerweise Γl ⊂ Ω folgt. F¨ur die numerischen Resul-tate in Abschnitt 2.7 hat es sich als ausreichend erwiesen, f¨ur die Filterung der Punktmenge Γl :={xl1, . . . , xlN/2l} ein reines Downsampling entlang des Pfades zu verwenden. Den zugeh¨origen trivialen
”Filter“ bezeichnet man zuweilen auch als
”Lazy-Filter“.
3. Man beachte, dass wir im Gegensatz zur urspr¨unglichen EPWT in jedem Zerle-gungslevel mit gestreuten Punkten xlj ∈ Rd arbeiten. Die
”Punkte“, entlang derer
Pfade konstruiert werden, sind bei der EPWT in [101] ab dem zweiten Zerlegungs-level sogenannte Indexmengen, welche sich als Vereinigung zweier im Pfad des vorherigen Level aufeinanderfolgender Punkte bzw. Indexmengen ergeben. In [130]
wird von Tenorth auch eine sogenannte
”Mittelpunkt-EPWT“ diskutiert, bei der eine Indexmenge durch ihren Mittelpunkt identifiziert wird. Dieses Vorgehen ent-spricht Algorithmus 2.63 unter Verwendung eines entsprechend normierten Haar-Filters f¨ur die Aktualisierung der Punktmenge Γl.
4. F¨ur das Thresholding in Schritt 4 sind andere Shrinkage-Funktionen wie Soft-, Firm- oder Garrote-Shrinkage (siehe etwa [96]) denkbar. Die Wahl des Shrinkage-Parameters θ hat entscheidenden Einfluss auf das Resultat der Entst¨orung und muss daher sorgf¨altig erfolgen.
5. Die einzelnen Iterationen des Algorithmus sind vollkommen unabh¨angig voneinan-der. Der Gesamtaufwand des Algorithmus steigt folglich linear mit der Anzahl der verwendeten Iterationen. Andererseits erm¨oglicht die Unabh¨angigkeit der Iteratio-nen eine Parallelisierung. Damit kann der zeitliche Aufwand bei der Durchf¨uhrung des Algorithmus erheblich reduziert werden.
6. Der Algorithmus 2.63 ist grunds¨atzlich f¨ur beliebige Dimension d ∈ N≥2 anwend-bar. Außer der Voraussetzung, dass die Punktmenge Γ quasi-uniform (siehe De-finition 2.61) in einer Menge Ω, welche sich als Vereinigung endlich vieler be-schr¨ankter, zusammenh¨angender Teilmengen des Rd schreiben l¨asst, ist, bestehen weiterhin f¨ur die Anwendbarkeit von Algorithmus 2.63 keinerlei Einschr¨ankungen an die Form und Verteilung der Punktmenge. Insbesondere sind wir nicht auf zweidimensionale, regul¨are und rechteckige Gitter angewiesen. Der Pfad (im Sin-ne eiSin-nes Polygonzuges von xlpl(1) nach xlpl(N/2l)) muss, wenn Ω etwa nicht zusam-menh¨angend oder nicht konvex ist, im Allgemeinen nicht vollst¨andig in Ωenthalten sein, was f¨ur die Durchf¨uhrung des Algorithmus 2.63 grunds¨atzlich jedoch nicht von Bedeutung ist.
7. Algorithmus 2.63 l¨asst sich wie in Abbildung 2.4 als modifizierte Wavelet-Filterbank interpretieren. Dabei wird vor der Anwendung der Zerlegungsfilter ein Permu-tationsoperator Pl eingef¨ugt, welcher die Tiefpass-Daten (clj)N/2j=1l zu einer Folge (clpl(j))N/2j=1l entlang des konstruierten Pfades pl umordnet. Analog wird nach der Rekonstruktion der Tiefpass-Werte der zugeh¨orige inverse Permuationsoperator (Pl)−1 erg¨anzt. Man beachte, dass bei der Darstellung in Abbildung 2.4 jedoch die Adaptivit¨at der Permutationsoperatoren nicht deutlich wird. Weiterhin sind die im Hintergrund ablaufenden Aktualisierungen der Punktmengen Γl nicht dargestellt.
c0 P0
Abbildung 2.4: Algorithmus 2.63 als modifizierte Filterbank (vgl. Abbildungen in [108]).