ČESKÉ VYSOKÉ UČENÍ TECHNICKÉ V PRAZE Fakulta stavební Katedra fyziky
Vyhodnocování interferogramů metodou Fourierovy transformace
Evaluation of Interferograms Using a Fourier-Transform Method
diplomová práce
Studijní program: Studijní obor:
Geodézie a kartografie Geodézie a kartografie
Vedoucí práce:
prof. RNDr. Antonín Mikš, CSc.
Bc. Petr Pokorný
Praha 2014
Čestné prohlášení Prohlašuji, že jsem předloženou práci vypracoval samostatně a že jsem uvedl veškeré použité informační zdroje v souladu s Metodickým pokynem o etické přípravě vysokoškolských závěrečných prací.
V Tehově dne 18. 12. 2013 Bc. Petr Pokorný
Poděkování Předkládaná práce by nemohla vzniknout bez rad a konzultací se členy Skupiny aplikované optiky při katedře fyziky FSv ČVUT v Praze. Zejména bych rád poděkoval mému vedoucímu diplomové práce prof. RNDr. Antonínu Mikšovi, CSc., dále doc. Ing. Jiřímu Novákovi, Ph.D., a v neposlední řadě Ing. Pavlu Novákovi, Ph.D., za jejich odborné připomínky a návrhy. Za užitečné poznámky ke korektuře práce děkuji Anežce Kabátové.
V Tehově dne 18. 12. 2013 Bc. Petr Pokorný
Abstrakt Práce představuje teoretické poznatky k pochopení principu interferometrických měření, zejména je věnována pozornost metodě vyhodnocení interferogramů pomocí Fourierovy transformace. Nejprve je ukázána matematická formulace Fourierovy transformace pro spojité a diskrétní případy, dále je představena interferometrie jako technika velmi přesného měření s náležitostmi předzpracování, vyhodnocení a závěrečného zpracování dat. Následuje část specializovaná na způsob vyhodnocení pomocí Fourierovy transformace s popisem jednotlivých dílčích kroků. Vyhodnocení simulovaných a reálných dat je prováděno ve vlastní vytvořené aplikaci. Klíčová slova Interferometrie, interferogram, metoda Fourierovy transformace, unwrapping
Abstract Thesis presents a theoretical background for comprehension of interferometrical principles. Under major attention is interferogram evaluation using a Fourier-transform method. Firstly, a mathematical formulation of Fourier transform for continuous and discrete data points is presented. Subsequently, interferometry as a technique for a precise measurement is described, in common with preprocessing, evaluation and postprocessing of measured data. The Fouriertransform method for interferogram analysis is successively and thoroughly depicted in the following part. Simulated and real data evaluation is performed in own application. Keywords Interferometry, interferogram, Fourier-transform method, unwrapping
Obsah 1
ÚVOD...................................................................................................................................................................8
2
FOURIEROVA TRANSFORMACE ..............................................................................................................10 2.1
DEFINICE FOURIEROVY TRANSFORMACE ...................................................................................................10
2.2
DISKRÉTNÍ FOURIEROVA TRANSFORMACE .................................................................................................12
2.3
SOUVISLOST FOURIEROVY TRANSFORMACE, DISKRÉTNÍ FOURIEROVY TRANSFORMACE A FOURIEROVÝCH ŘAD ............................................................................................................................................................14
3
INTERFERENCE SVĚTLA............................................................................................................................20 3.1
INTENZITA A KONTRAST ROVINNÉ MONOCHROMATICKÉ HARMONICKÉ VLNY ............................................20
3.2
PRINCIP SUPERPOZICE ................................................................................................................................24
3.3
INTERFERENCE A KONTRAST LINEÁRNĚ POLARIZOVANÝCH MONOCHROMATICKÝCH VLN STEJNÉ FREKVENCE ................................................................................................................................................25
3.4
INTERFERENCE A KONTRAST MONOCHROMATICKÝCH ELIPTICKY POLARIZOVANÝCH VLN STEJNÉ FREKVENCE ................................................................................................................................................29
3.5
KOHERENCE VLNĚNÍ ..................................................................................................................................31
3.5.1
Prostorová koherence...........................................................................................................................31
3.5.2
Časová koherence.................................................................................................................................34
3.6 4
INTENZITA A KONTRAST POLYCHROMATICKÉHO A KVAZIMONOCHROMATICKÉHO ZÁŘENÍ ........................36
INTERFEROMETRIE.....................................................................................................................................38 4.1
DVOUSVAZKOVÁ A STŘIHOVÁ INTERFEROMETRIE .....................................................................................40
4.1.1
Dvousvazková interferometrie ..............................................................................................................40
4.1.2
Střihová interferometrie .......................................................................................................................43
4.2
ZÁZNAM A VYHODNOCENÍ INTERFEROGRAMU ...........................................................................................44
4.2.1
Záznam interferenčního pole ................................................................................................................44
4.2.2
Předzpracování zaznamenaných dat ....................................................................................................47
4.2.2.1
Model šumu v měřeném interferogramu ..................................................................................................... 48
4.2.2.2
Potlačení vysokofrekvenčního a nízkofrekvenčního náhodného šumu ....................................................... 48
4.2.2.3
Eliminace nerovnoměrné intenzity pozadí interferogramu.......................................................................... 52
4.2.2.4
Úprava kontrastu interferogramu ................................................................................................................ 53
4.2.2.5
Výběr relevantní vyhodnocované oblasti – maskování ............................................................................... 56
4.2.3
Vyhodnocení interferogramu ................................................................................................................56
4.2.3.1
Vztah fázového a dráhového rozdílu, základní podmínka metodiky vyhodnocení...................................... 57
4.2.3.2
Detekce a polynomiální aproximace interferenčních proužků s následným vyhodnocením fázového rozdílu .................................................................................................................................................................... 58
4.2.3.3
Interferometrie s posuvem fáze ................................................................................................................... 60
4.2.3.4
Interferometrie se zaváděním nosné frekvence ........................................................................................... 63
4.2.4
Konečné zpracování dat – postprocessing............................................................................................69
4.2.4.1
Analytické vyjádření diskrétních rekonstruovaných dat.............................................................................. 69
4.2.4.2
Odstranění zbytkových vlivů nejednoznačně určené prostorové frekvence ................................................ 79
6
Obsah 5
VYHODNOCENÍ INTERFEROGRAMŮ METODOU FOURIEROVY TRANSFORMACE.................81
6
UNWRAPPING ................................................................................................................................................91 6.1
ÚVOD DO PROBLEMATIKY UNWRAPPINGU ..................................................................................................91
6.2
RESIDUA, MAPY KVALITY A MASKY .........................................................................................................101
6.3
ALGORITMY SLEDUJÍCÍ PŘEDEM URČENOU CESTU (PFA) .........................................................................105
6.3.1
Goldsteinův algoritmus (GA)..............................................................................................................105
6.3.2
Quality-Guided Path Following algoritmus (QGPFA) ......................................................................106
6.3.3
Mask-Cut algoritmus (MCA) ..............................................................................................................107
6.3.4
Flynnův algoritmus (FA) ....................................................................................................................108
6.4 6.4.1
Algoritmy nejmenších čtverců bez vah (ULS).....................................................................................111
6.4.2
Vážený multi-grid (WMG) ..................................................................................................................113
6.4.3
Metoda sdružených gradientů s předpodmíněním (PCG)...................................................................114
6.4.4
Metoda minima LP-normy...................................................................................................................114
6.5 7
8
METODY MINIMALIZUJÍCÍ NORMU (MNM)...............................................................................................109
DALŠÍ (KOMBINOVANÉ) ALGORITMY........................................................................................................116
ANALÝZA ALGORITMŮ UNWRAPPINGU.............................................................................................117 7.1
ZDROJE POUŽITÝCH ALGORITMŮ ..............................................................................................................117
7.2
VLIV NÁHODNÉHO ŠUMU NA ALGORITMY UNWRAPPINGU ........................................................................118
7.3
VLIV RESIDUÍ OBRAZU FÁZOVÉHO ROZDÍLU NA ALGORITMY UNWRAPPINGU ...........................................125
VYHODNOCENÍ SIMULOVANÝCH A REÁLNÝCH DAT ....................................................................129 8.1
VÝPOČETNÍ SOFTWARE K VYHODNOCENÍ INTERFEROGRAMU POMOCÍ FTM.............................................129
8.2
VYHODNOCENÍ SIMULOVANÝCH A REÁLNÝCH DAT ..................................................................................135
9
ZÁVĚR ............................................................................................................................................................141
10
POUŽITÁ LITERATURA.............................................................................................................................142
11
SEZNAM OBRÁZKŮ ....................................................................................................................................150
12
SEZNAM TABULEK .....................................................................................................................................153
13
SEZNAM PŘÍLOH.........................................................................................................................................154
7
1 Úvod V různých technických oblastech je nutné bezkontaktně a s velkou přesností určovat tvary ploch nebo jejich odchylky od nominální hodnoty. Jednou z metod takového měření může být interferometrie s následným numerickým vyhodnocením a analýzou. V současnosti se využívá vhodných algoritmů, které zpracovávají dodané měřené interferogramy, na nichž je zaznamenána intenzita vzniklá interferencí zpravidla dvou vln [1 – 9] – jedné odrážející se od referenční plochy a druhé od plochy testované. Je známo několik postupů pro měření a vyhodnocování interferogramů, přehled lze nalézt např. v [3 – 10]. Jeden z nich je založen na použití diskrétní Fourierovy transformace. Metoda je to v současnosti díky vyspělým výpočetním softwarům numericky snadno proveditelná, navíc přináší některé výhody oproti jiným způsobům. Například samotný postup vyhodnocení částečně potlačuje variace v amplitudě měřené intenzity a pro vyhodnocení nám postačuje pouze jeden obraz. Ve frekvenčním spektru jsme také schopni provádět filtrace obrazu a částečně tak eliminovat šum vzniklý náhodnými chybami použitého detektoru. I když metoda Fourierovy transformace neposkytuje nejpřesnější výsledky oproti metodám, jež využívají více měřených obrazů, může být aplikována pro velmi přesné analýzy dynamických systémů. Samotné provedení výpočtu Fourierovy transformace je totiž v praxi prováděno zakomponováním do hardwaru měřícího zařízení. Oblastmi aplikací interferometrie jsou například testování optických sférických a asférických prvků [10, 11], vysokorychlostní dynamická 3D profilometrie mikromechanických zařízení [12 – 14], analýza hoření plamenů [15], optická tomografie nestacionárního vysokorychlostního toku plynů [16], tzv. „spekl“ interferometrie pro kvantitativní měření 3D deformačních polí [17], analýzy napětí [18], proměřování tenkých vrstev [19], charakteristiky laserového plazmatu [20], analýzy mega-gaussovských magnetických polí v laserově produkovaných plazmatech [21], měření pro určování prostorové koherence a indexů lomů prostředí [22, 23], určování prostorových polarizací [24 – 26] nebo spektroskopická polarimetrie [27]. Z biologických a medicínských aplikací jmenujme použití v tomografii [28 – 29], měření dynamiky topografie trhlin [30] nebo měření deformací biologických tkání [31]. Z oblasti měření extrémních fyzikálních jevů to pak jsou například měření extrémně malých magnetických polí s použitím elektronové holografie [32 – 34] a měření deformací krystalických mřížek se sub-Ångstromovou
8
Úvod přesností [35]. Satelitní interferometrie umožňuje určovat například pohyby tektonických desek a zemětřesení [36], oceánské slapy [37], mapování přehrad a nádrží [38] a další. Práce předkládá teoretické podklady k pochopení problematiky interferometrie jako přesné bezkontaktní metody pro měření optických ploch a vyhodnocení interferogramů metodou Fourierovy transformace. První část je věnována matematické formulaci Fourierovy transformace pro spojité a diskrétní případy, následuje kapitola popisující fyzikální původ interference světla pro monochromatické, kvazimonochromatické a polychromatické záření. V další části je představena interferometrie jako technika měření, její základní principy a postupy získávání, předzpracování, vyhodnocení a konečných úprav analyzovaných obrazů. Následuje kapitola věnovaná přímo metodě Fourierovy transformace s náležitostmi, které je nutné při rekonstrukci reálných fázových polí provádět. Jednou z částí tohoto procesu je tzv. unwrapping. Tento matematický problém je stěžejní a nejkomplikovanější při rekonstrukci reálných polí, a proto je mu věnována celá další kapitola. Algoritmy provádějící unwrapping jsou následně analyzovány v předposlední kapitole. Závěrem je představen výpočetní software vytvořený pro vyhodnocení interferogramů metodou Fourierovy transformace, s jehož využitím je provedeno zpracování a analýza simulovaných i reálných měřených dat. V přílohách je uveden seznam prací, jež autor publikoval v mezinárodních impaktovaných a v českých recenzovaných časopisech po dobu svého bakalářského a magisterského studia na Fakultě stavební ČVUT v Praze. Následně také soupis soutěžních prací autora, s nimiž se účastnil mezinárodních a celofakultních soutěží. Datovou přílohou je vytvořená aplikace „VI-FTM“ pro vyhodnocení interferogramů metodou Fourierovy transformace v prostředí MATLAB R2011a přeložená pro 32bitový operační systém Windows.
9
2 Fourierova transformace Představme nyní základní matematický nástroj této práce. Je jím tzv. Fourierova transformace a inverzní Fourierova transformace [39 – 44]. Uplatnění obecně nachází všude tam, kde je třeba vyjádřit a analyzovat periodický i neperiodický signál získaný v čase. V technické praxi se Fourierova transformace využívá zejména z toho důvodu, že slouží k převodu signálu z tzv. časové oblasti, rozumějme záznam určitého jevu (signálu) v čase, do tzv. oblasti frekvenční, tj. umožňuje zjistit, jaké frekvence zaznamenaný jev (signál) obsahuje. Frekvenční oblast nám poté nese některé důležité informace, které bychom ze záznamu v čase nemohli získat. Z potřeb techniků a fyziků vznikla tzv. diskrétní Fourierova transformace (DFT) (toto označení bylo poprvé použito roku 1967 [45, 46]). Bylo totiž nutné analyzovat měřená data diskrétní povahy. Vědecká veřejnost již dříve (od roku 1815 [42]) ovšem znala tzv. Besselovy vztahy, které jsou nyní definicí DFT. Nástup rozvoje tohoto matematického aparátu nastal roku 1965, kdy Cooley a Tukey publikovali obdivuhodně rychlý algoritmus výpočtu DFT (dnes známý jako rychlá Fourierova transformace – FFT z angl. Fast Fourier Transform). Neustále se zdokonalující výpočetní technika tento rozvoj jen umocnila. V současnosti DFT (FFT) nachází uplatnění např. ve zpracování signálu nebo analýze obrazu [47, 48], ale i v dalších disciplínách.
2.1
Definice Fourierovy transformace
V literatuře různých vědeckých oblastí se používají různé způsoby definic spojité Fourierovy transformace v N dimenzích a její inverze [41]. Tuto rozdílnost lze pokrýt zavedením tří nenulových reálných konstant A, B a k. Konstanta k musí být reálná bezpodmínečně, konstanty A a B mohou být i komplexní, ale je to bezúčelné. Konstanty jsou svázány podmínkou, vyplývající z tzv. fundamentální věty o Fourierově transformaci [41, 42], AB =
k . 2π
(1)
Nechť jsou komplexní funkce f(x) a F(X) absolutně integrovatelné a po částech hladké reálných proměnných x a X ∈ EN. Potom je Fourierova transformace funkce f(x) definována jako [41] 10
Fourierova transformace ∞
∞
−∞
−∞
∞
∞
−∞
−∞
FT { f (x)} = F ( X) = A N ∫ L ∫ f (x) exp(− i kXx ) d N x
(2)
a inverzní Fourierova transformace jako FT −1 {F ( X)} = f (x) = B N ∫ L ∫ F ( X) exp( i kXx ) d N X ,
(3)
kde FT{•} značí operátor Fourierovy transformace a FT-1{•} operátor inverzní Fourierovy transformace. Jádro integrální Fourierovy transformace K(x, X) = exp(−ikXx) je symetrické, neboť platí [39, 41, 42] K ( x , X ) = K ( X, x ) ,
(4)
a dále faktorizovatelné, tj. [39, 41, 42] N
K ( x, X ) = ∏ K j ( x j , X j ) ,
(5)
j =1
což v některých případech velmi usnadňuje výpočet Fourierova integrálu (2) a (3). Ukažme si nyní tzv. fundamentální větu o Fourierově transformaci [41] a z ní důsledek pro podmínku (1). Nechť pro funkci f(x) existuje Fourierův integrál, potom v bodech spojitosti platí [41]
FT −1{FT { f (x)}} = FT {FT −1 { f (x)}} = f (x) ,
(6)
což lze snadno ukázat dosazením z (1) a (2), ∞ ∞ ∞ ∞ FT −1 {FT { f (x)}} = B N ∫ L ∫ A N ∫ L ∫ f (ρ) exp(− i kXρ ) d N ρ exp( i kXx ) d N X −∞ −∞ −∞ −∞ ∞ ∞ ∞ ∞ N = ( AB ) ∫ L ∫ f (ρ) ∫ L ∫ exp[i k (x − ρ )X] d N X d N ρ . −∞ −∞ − ∞ − ∞
(7)
Integrál ve složené závorce je tzv. Diracova distribuce [39 – 42] N
2π δ (x − ρ ) , L exp [ i k ( x − ρ ) X ] d X = ∫−∞ −∫∞ k ∞
∞
N
(8)
kde δ(x – ρ) je Diracova funkce [39 – 42]. Dosazením do (7) dostáváme
11
Fourierova transformace 2π FT {FT { f (x)}} = AB k
N
−1
∞
∞
−∞
−∞
N ∫ L ∫ f (ρ)δ (x − ρ ) d ρ .
(9)
Integrál v (9) je díky vlastnostem Diracovy distribuce roven f(x) v bodech spojitosti funkce f(x) a v bodech konečné nespojitosti je roven střední hodnotě. Aby platilo tvrzení fundamentální věty, musí díky (9) platit
2πAB = k ,
(10)
což je podmínka (1). V dalším textu budeme uvažovat hodnoty konstant A =1, B =1, k = 2π ,
(11)
čili definujeme Fourierovu transformaci jako ∞
∞
−∞
−∞
∞
∞
−∞
−∞
FT { f (x)} = F ( X) = ∫ L ∫ f (x) exp(− i 2 πXx ) d N x
(12)
a inverzní Fourierovu transformaci jako FT −1 {F ( X)} = f (x) = ∫ L ∫ F ( X) exp( i 2 πXx ) d N X .
2.2
(13)
Diskrétní Fourierova transformace
Jelikož měřená data jsou diskrétní povahy a jejich počet je konečný, je nutné představit nástroj pro jejich Fourierovskou transformaci. Podrobný vyčerpávající popis vlastností a způsobů výpočtů může čtenář najít např. v [42]. V této části představíme pouze jednodimenzionální a dvojdimenzionální tvar, protože ty jsou využívány v technické praxi nejčastěji. Velmi rychlý způsob výpočtu publikovali Cooley a Tukey roku 1965 [49], jejich metoda je dnes používána pod názvem rychlá Fourierova transformace (FFT z angl. Fast Fourier Transform). Matematické softwary (např. MATLAB) mají její algoritmy zpravidla implementované. V programovacím jazyce C jsou volně k dispozici (free software pod GNU General Public License) knihovny projektu FFTW [50] vznikajícího na Massachusetts Institute of Technology v Bostonu (USA). Tyto knihovny nabízí možnost výpočtu v jedné nebo více dimenzích,
12
Fourierova transformace volitelného rozsahu vstupních reálných i komplexních dat a jsou považovány za nejrychlejší a nejoptimalizovanější. Definici diskrétní Fourierovy transformace1 v jedné dimenzi můžeme psát jako [42 – 44, 47] N −1
Fk = ∑ f n exp(− ink 2 π / N ) , k = 0, 1, K , N − 1 ,
(14)
n =0
kde fn, n = 0, 1, ... , N – 1, je posloupnost N konečných komplexních čísel, i = √(-1). Jednodimenzionální zpětnou diskrétní Fourierovu transformací potom rozumíme fn =
1 N −1 ∑ Fk exp(ink 2π / N ) , n = 0, 1, K , N − 1 . N k =0
(15)
Jestliže je Fk ve vztahu (15) obrazem fn ze vztahu (14) (jsou-li hodnoty Fk výsledkem diskrétní Fourierovy transformace na posloupnost čísel fn), potom (15) vyjadřuje původní posloupnost a je inverzní ke vztahu (14). Toto tvrzení plyne z dosazení (14) do (15),
N −1 f v exp(− ivk 2 π / N ) exp(ink 2 π / N ) ∑ ∑ k =0 v = 0 N −1 1 N −1 = ∑ f v ∑ exp[i(n − v )k 2 π / N ] , N k =0 v =0
fn =
1 N
N −1
(16)
kde dále platí 1 N
N −1
0, pokud v ≠ n, pokud v = n,
∑ exp[i(n − v )k 2π / N ] = 1, k =0
(17)
a tudíž můžeme (16) symbolicky rozepsat jako f n = f 0 ⋅ 0 + f1 ⋅ 0 + ... + f n ⋅1 + ... + f N −1 ⋅ 0 .
(18)
Tímto jsme vyslovené tvrzení snadno dokázali. V některých literaturách může čtenář nalézt definice odlišné zpravidla o násobnou konstantu, případně znaménko v exponentu. Odlišnost je ekvivalentní zavádění konstant (11) ve spojitém případě Fourierovy transformace. My ovšem budeme uvažovat výše uvedené definice (14) a (15), jelikož ty jsou zpravidla používány v moderních výpočetních systémech.
1
Přesnější název by byl konečná (finitní) diskrétní Fourierova transformace, protože máme k dispozici konečnou posloupnost diskrétních dat (nechť čtenář srovná s definicí Fourierovy spojité transformace). Pro jednoduchost zápisu ovšem tuto skutečnost budeme zanedbávat.
13
Fourierova transformace Geometricky lze posloupnost komplexních čísel Fk, k = 0, 1, ... , N – 1, interpretovat jako obecnou soustavu vektorů se společným počátkem [43, 44]. Definiční vztah (14) pro Fk poté můžeme považovat za rozklad každého vektoru do soustavy N vektorů fn. Jelikož je v technické praxi často vyžadována analýza dvojrozměrných obrazů, představme také dvojrozměrnou diskrétní Fourierovu transformaci. Ta je nejčastěji definována jako [42 – 44, 49] M −1 N −1
Fk ,l = ∑∑ f m ,n exp(− imk 2 π / M ) exp(− inl 2 π / N ) ,
(19)
m = 0 n =0
kde fm,n je dvojrozměrná posloupnost, m = 0, 1, ... , M – 1, n = 0, 1, ... , N – 1, periodická s periodou M v indexu m a s periodou N v indexu n. Inverzní transformace je poté definována jako f m,n =
2.3
1 M −1 N −1 ∑∑ Fk ,l exp(imk 2π / M )exp(inl 2π / N ) . MN k =0 l =0
(20)
Souvislost Fourierovy transformace, diskrétní Fourierovy transformace a Fourierových řad
Abychom demonstrovali využitelnost Fourierovy transformace v technických oblastech, zabývejme se nyní její souvislostí s Fourierovými řadami. Jejich využití nachází uplatnění např. ve spektrální analýze signálů, zpracování obrazu apod. [47, 48]. Uvažujme případ jedné proměnné a definujme Fourierovu řadu. Stejně jako tomu bylo v předchozích částech, v různých literaturách se můžou definice lišit v závislosti na volbě konstant (11). Fourierova řada k dané obecně komplexní spojité funkci f(t) periodické s periodou P, P > 0, je definována jako [39, 40, 42, 43]
ξ (t ) =
∞
∑C
k = −∞
k
exp(i 2 πkt / P ) , − ∞ < t < ∞ ,
(21)
kde i = √(-1) a pro koeficienty Ck platí
1P Ck = ∫ f (t ) exp(− i 2 πkt / P ) dt , k = 0, ± 1, ± 2, K . P0
(22)
14
Fourierova transformace Jak je patrné z (22), hodnoty koeficientů Ck jsou počítány z hodnot v intervalu periodicity. Pokud funkce nebude periodická, můžeme považovat celý její definiční obor za periodu, a tak koeficienty vypočíst. I pro funkci na konečném intervalu různou od nuly lze koeficienty určit, a to snadno jako Ck =
1 P
t0 + P
∫ f (t ) exp(− i2πkt / P ) dt ,
k = 0, ± 1, ± 2, K ,
(23)
t0
kde t0 ≤ t ≤ t0 + P. Pro reálné funkce f(t) bude platit, jak vyplývá z (22), C−k = Ck∗ , k = 1, 2, K .
(24)
Dále lze vztahy (21) a (22) pro reálné funkce s užitím Eulerova vzorce [39, 40]
r (cos ϕ + i sin ϕ ) = r exp(iϕ )
(25)
přepsat jako ∞
k =1
2 πkt 2 πkt + bk sin , P P
ξ (t ) = a0 + ∑ ak cos
(26)
kde
1P a0 = ∫ f (t ) dt = C0 , P0 2P 2 πkt ak = ∫ f (t ) cos dt = 2 Re{C k } , P0 P
(27)
2P 2 πkt bk = ∫ f (t ) sin dt = −2 Im{C k } , k = 1, 2, K , P0 P kde Re{•} a Im{•} značí reálnou a imaginární složku. Vztahy (21) a (22) lze také přepsat jako [42, 43] 2 πkt + φk , P
∞
ξ (t ) = A0 + ∑ Ak cos k =1
(28)
kde
A0 = a 0 = C 0 , Ak = a k2 + bk2 = 2 C k ,
(29)
Im{C k } , k = 1, 2, K . Re{C k }
φ k = arg C k = arctan
15
Fourierova transformace V technické praxi je úloha, kdy hledáme koeficienty (29) rozvoje (28) nazývána harmonická analýza. Opačný problém, kdy z koeficientů sestavujeme daný signál, je nazýván harmonická syntéza [42]. Každý signál jsme tedy schopni rozložit do řady „kosinusovek“, určit jejich amplitudu Ak a počáteční fázi ϕk. Veličina Ck je nazývána komplexní amplitudou. Je výhodné ji používat společně s komplexním vyjádřením, i když v praxi u reálných signálů má smysl uvažovat jen reálné složky harmonického rozvoje s kladnými frekvencemi ω, kde ω = 2πk/P. Ze základních vlastností Fourierových řad jmenujme tyto [42]: •
Posloupnost částečných součtů ξn konverguje k f v průměru, čili platí: P
lim ξ n (t ) − f (t ) dt = 0 . n →∞ ∫ 2
(30)
0
•
Mezi funkcí f a koeficienty její Fourierovy řady platí tzv. Parservalova rovnost: ∞ 1P 2 2 f ( t ) d t = Ck . ∑ ∫ P0 k = −∞
•
(31)
Koeficienty Fourierovy řady konvergují k nule: lim Ck = 0 . k →∞
(32)
Z výše uvedených vlastností lze odvodit tvrzení, že má-li být funkce aproximována trigonometrickým polynomem n-tého stupně ve smyslu nejmenších čtverců (suma kvadrátů odchylek má být minimální), potom je tímto polynomem právě částečný součet Fourierovy řady [42]. Souvislost Fourierovy řady a Fourierovy transformace tedy vystihuje následující. Jak může být ukázáno [42], pro existenci Fourierovy řady i transformace finitních funkcí jsou nutné stejné podmínky. Potom lze určit jak koeficienty Fourierovy řady (22), tak i Fourierovu transformaci této funkce, která bude mít užitím (12) tvar P
FT { f (t )} = F (k / P ) = ∫ f (t ) exp(− i 2 πkt / P ) dt .
(33)
0
Porovnáním (22) a (33) dostáváme tvrzení [42]: „Koeficienty Ck Fourierovy řady finitní funkce f jsou rovny hodnotám Fourierovy transformace F této funkce děleným délkou periody.“ Platí tedy Ck =
1 F (k / P ) , k = 0, ± 1, ± 2, K . P
(34)
16
Fourierova transformace V praxi nás ovšem nejvíce bude zajímat souvislost diskrétní Fourierovy transformace s Fourierovou řadou, jelikož měřená a analyzovaná data jsou zpravidla diskrétní povahy. Lze ukázat, že přímé vyjádření funkce F (Fourierovy transformace) a čísel Ck (komplexních amplitud) pomocí diskrétních čísel Fk (diskrétní Fourierovy transformace) není obecně možné. Platí však přibližná vyjádření [42] Ck ≈
1 Fk , N
Fk ≈
1 k F , T NT
(35)
kde N značí počet diskrétních bodů, T = P/N. Vztahy (35) jsou závislé na tom, jak dobře je splněna podmínka (32), čili jak budou koeficienty Ck klesat s rostoucím k, nebo jak bude funkce F(k/NT) na základním intervalu klesat s rostoucím k/NT. Přibližné vztahy (35) přejdou v rovnosti pouze tehdy, pokud k F =0 ∀ NT
k π > , NT T
(36)
nebo pokud
Ck = 0 ∀
k > N /2.
(37)
Na následujících obrázcích je demonstrována analýza signálu a jeho rekonstrukce pomocí DFT. Nechť funkce f(t) je dána součtem tří funkcí f1(t), f2(t) a f3(t), 2 πk i f i (t ) = cos t + φi , P
(38)
kde k1 = 1, ϕ1 = 0, k2 = 5, ϕ 2 = π/3, k3 = 9, ϕ 3 = π/5 a P = 2π. Pro DFT je dále voleno N = 629. Obr. 2 zobrazuje vypočtené amplitudy Ak a počáteční fáze ϕk s využití DFT dle vztahů (35) a (29). Pro přehlednost byl zvolen maximální stupeň zobrazovaného k = 50. Ze samotné části b) obr. 2 není patrné, které počáteční fáze jsou relevantní pro daný harmonický rozvoj. V kombinaci s částí a) ovšem zjišťujeme, že jsou to pouze ty, kde k = 1, k = 5 a k = 9, což je v souladu s naším zadáním jednotlivých funkcí fi(t), které nám svým součtem definují funkci f(t). Část c) zobrazuje polární graf, kde jsou prezentovány pouze dominantní prvky v rozvoji. Pro velikost vypočtených amplitud platí Ak ≈ 1, což odpovídá našim zadaným funkcím, a velikost počátečních fází je
ϕ1 = 0.0039 ≈ 0.2° ≈ 0, ϕ5 = 1.0621 ≈ 60.9° ≈ π/3 a ϕ9 = 0.6575 ≈ 37.7° ≈ π/5.
17
Fourierova transformace
Obr. 1
Obr. 2
Funkce f(t) = f1(t) + f2(t) + f3(t)
Zobrazení a) amplitud Ak vypočtených pomocí DFT, b) počátečních fází ϕk a c) polárního zobrazení dominantních členů rozvoje
18
Fourierova transformace
Obr. 3
Zobrazení a) rekonstruované funkce ξ(t) a b) jejího rozdílu od původní funkce f(t)
Na obr. 3 následně vidíme rekonstruovanou funkci ξ(t) a její rozdíl od původní funkce f(t). Odchylky v rekonstrukci nastaly z výše popsaných důvodů při použití DFT, jelikož transformace není počítána přes neomezený interval, ale pouze pro finitní funkci a diskrétní data.
19
3 Interference světla Představme nyní teoretické poznatky nezbytné k popisu interference, čili nutné k pochopení vzniku interferenčních obrazců, jejichž prostorové rozložení je zaznamenáváno na detektorech záření a následně (dnes zpravidla analyticky) vyhodnocováno. Celá problematika by si zasluhovala mnohem větší prostor, než zde může být poskytnut, k hlubšímu studiu tedy doporučíme čtenáři dostupnou literaturu [1 – 8].
3.1
Intenzita a kontrast rovinné monochromatické harmonické vlny
Zabývejme se nyní vztahy pro rovinné monochromatické vlny, jež jsou základem a účinným nástrojem pro popis šíření světla ve formě elektromagnetických vln. Základní kostru tvoří Maxwellovy rovnice spojené s materiálovými rovnicemi a hraničními podmínkami [1, 4, 5, 51 – 53]. Pomocí nich můžeme popsat veškeré jevy, které nastávají v námi vyšetřovaných situacích. Vyjděme z Ampérova zákona celkového proudu, jak je znám v diferenciálním tvaru Maxwellových rovnic [1, 4, 5, 51 – 53], rot H (r, t ) = j(r, t ) +
∂D(r, t ) , ∂t
(39)
kde H(r, t) je vektor intenzity magnetického pole, j(r, t) vektor proudové hustoty volných nábojů (vodivých proudů), D(r, t) je vektor elektrické indukce, r značí polohový vektor od zdroje pole a t čas. Dále pak z Faradayova indukčního zákona, který zní [1, 4, 5, 51 – 53] rot E(r, t ) = −
∂B(r, t ) , ∂t
(40)
kde E(r, t) je vektor intenzity elektrického pole a B(r, t) vektor magnetické indukce. Závislost mezi veličinami uvažujme ve formě materiálových rovnic [1, 4, 5, 51 – 53]
D(r, t ) = ε E(r, t ) , B(r, t ) = µ H (r, t ) ,
(41)
kde ε značí permitivitu prostředí a μ permeabilitu. Rovnice (39) – (41) platí za předpokladu, že všechna tělesa jsou v klidu, materiálové konstanty nejsou závislé na čase a vektorech pole a
20
Interference světla v poli nejsou permanentní magnety a feromagnetická prostředí. Pro zjednodušení a přehlednost zápisu budeme v dalším textu vynechávat označení závislosti na poloze a čase, čili použijeme např. E = E(r, t). Násobme rovnici (39) skalárně vektorem intenzity elektrického pole, rovnici (40) vektorem intenzity magnetického pole a tyto od sebe odečtěme, potom s užitím (41) dostáváme ∂H ∂E , − Ej − ε E ∂t ∂t ∂ ε E 2 B2 + . H rot E − E rot H = −Ej − ∂t 2 2 µ
H rot E − E rot H = − µ H
(42)
Úprava levé strany (42) známá z vektorové algebry [39, 40] následně vede k [1, 5, 51]
∂ ε E 2 B2 = −Ej , grad(E × H ) + + ∂t 2 2 µ
(43)
což je tzv. Poyntingův teorém. Lze ho zapsat s užitím hustot energií také jako [5] grad N +
∂ (u pole + umedium ) = 0 , ∂t
(44)
kde N = E × H je tzv. Poyntingův vektor [1, 4, 5, 51 – 53], upole je objemová hustota energie elektromagnetického pole a umedium je objemová hustota energie dodaná polem do prostředí. Ukažme si nyní význam Poyntingova vektoru. Integrujme rovnici (44) přes objem V uzavřen plochou S, tedy ∂
∫ grad N dV = − ∂t ∫ (u
V
pole
+ u medium ) dV .
(45)
V
S užitím Gauss-Ostrogradského věty [39, 40] můžeme psát ∂
∫ Nn dS = − ∂t ∫ (u S
pole
+ u medium ) dV ,
(46)
V
kde n je jednotkový vektor vnější normály. Z výše uvedeného plyne, že tok Poyntingova vektoru z objemu V ohraničeného plochou S odpovídá rychlosti ztráty energie v tomto objemu. Poyntingův vektor tak definuje tok energie prostorem. Intenzita záření I je dána jako časová střední hodnota velikosti Poyntingova vektoru [1, 4, 5].
Časovou střední hodnotu musíme uvažovat proto, že žádný senzor není schopen registrovat okamžitou hodnotu vlnění. Tedy
21
Interference světla I = N = E× H =
1T 1T N d t = E × H dt . T ∫0 T ∫0
(47)
Dále definujme tzv. kontrast K jako [1, 2, 4 – 8] K=
I max − I min , I max + I min
(48)
kde Imax značí maximální intenzitu v obraze a Imin intenzitu minimální. Uvažujme nyní rovinnou monochromatickou harmonickou vlnu, pro kterou platí [1, 4, 5, 51 – 53]
E = E 0 exp[i(kr − ωt + δ )] , H =
ε (n × E) , µ
(49)
kde E0 je vektor počáteční amplitudy, k je vlnový vektor (obecně komplexní, prozatím uvažujme pouze reálný), k = kn, k = ω/v = 2π/λ je vlnové číslo [1, 4, 5], v je rychlost šíření vlnění, λ je vlnová délka, r je polohový vektor libovolného bodu na vlnoploše (polohový vektor od zdroje vlnění k dané vlnoploše), ω je úhlová frekvence, δ značí počáteční fázi, n je jednotkový normálový vektor vlnoplochy (vektor ve směru šíření rovinné vlny), i = (−1)1/2. Pro určení směru šíření energie rovinné vlny dosaďme ze vztahu (49) do vztahu pro Poyntingův vektor, dostáváme
N = E× H =
ε ε E × (n × E ) = [(EE )n − (En ) E] = ε (EE)n = ε E 2 n . µ µ µ µ
(50)
Vidíme, že Poyntingův vektor rovinné vlny je kolineární s normálovým vektorem vlnoplochy, a energie se tak šíří ve směru, který tato normála určuje. Pro aplikace, kde neuvažujeme např. absorpci prostředí, je relevantní pouze reálná část výrazů (49). Pro vektor elektrické intenzity jí získáme jako E=
E0 {exp[i(kr − ωt + δ )] + exp[− i(kr − ωt + δ )]} , 2
(51)
čili uvážením poloviny součtu daného komplexního čísla s číslem komplexně sdruženým. Užitím Eulerova vzorce dále platí E = E0 cos(kr − ωt + δ ) .
(52)
Pro intenzitu dostáváme dosazením do (47)
22
Interference světla 1T ε 2 1 I= ∫ E dt = T 0 µ T
ε E0 µ
T
2
∫ cos (kr − ωt + δ )dt , 2
(53)
0
kde předpokládáme neměnnost prostředí s časem, čili odmocninu z poměru permitivity a permeability prostředí můžeme vytknout před integrál. Po provedené integraci dostáváme I=
1 ε 1 2 E 0 1 + [sin 2(kr + δ ) − sin 2(kr − ωT + δ )] . 2 µ 2ωT
(54)
Předpokládáme-li např. detektor s odezvou T = 10-9 s a zelené světlo s frekvencí ω ≈ 4∙1015 Hz, bude ωT = 4∙105, z čehož je patrné, že druhý člen v závorce (54) můžeme zanedbat, protože bude přibližně 106krát menší než člen první. Pro intenzitu záření můžeme tedy s ohledem na výše uvedené předpoklady psát I=
1 ε 2 E0 . 2 µ
(55)
Uvažujme nyní průchod záření absorbujícím prostředím, které je v reálných aplikacích přítomno. Zaveďme tzv. komplexní index lomu [1, 4, 5] n = nω + iκ ω ,
(56)
kde nω je index lomu prostředí a κω je tzv. index absorpce. Jak indexy naznačují, index lomu prostředí a index absorpce je závislý na frekvenci záření, které daným prostředím prostupuje. Pro vlnový vektor poté dostáváme k = kn = k 0 nn =
2π
λ0
nn = k 0 (nω + iκ ω ) n = k r + ik i ,
(57)
kde k0 a λ0 je vlnové číslo a vlnová délka ve vakuu, n je komplexní index lomu prostředí vyjádřený vztahem (56). Dosazením do rovnice rovinné vlny (49) vede k
E = E 0 exp[i(k r r + ik i r − ωt + δ )] = E 0 exp[i(k r r − ωt + δ ) − k i r ]
= E 0 exp[i(k r r − ωt + δ )]exp[− k i r ] .
(58)
Je tedy patrné, že v tomto případě dochází k tlumení vlny, a to exponenciálně vzhledem k velikosti indexu absorpce. Pro intenzitu vlnění prostupujícího tímto prostředím poté platí
23
Interference světla
I=
1 ε 2 1 ε 1 ε 2 2 E = E 0 exp[− 2k i r ] = E 0 exp[− αd ] , 2 µ 2 µ 2 µ
(59)
kde α značí tzv. koeficient absorpce [1, 4, 5], α = 2ki =2k0κω, a vzdálenost d = nr značí vzdálenost bodu vlnoplochy od jejího zdroje. Z výše uvedeného vyplývá, že intenzita exponenciálně klesá se vzdáleností. Vztah (59) se nazývá Boguer-Lambertův zákon [1, 4, 5]. Obr. 4 demonstruje vliv absorpce prostředí na intenzitu záření.
Obr. 4
3.2
Závislost intenzity záření na vzdálenosti pro koeficient absorpce α = 2, I0 = 0.5√(ε/μ)|E0|2; a) závislost v rovině y = 0, potom d = x, b) závislost v rovině (x,y), kde d = (x2 + y2)1/2
Princip superpozice
Před vyslovením samotného principu superpozice musíme ujasnit jisté termíny. Je nutné rozlišovat intenzitu elektromagnetického vlnění I (časovou střední hodnotu Poyntingova vektoru) a velikosti vektorů pole – vektoru intenzity elektrického pole E a vektoru magnetické intenzity H. Důvod, proč je toto rozlišení zdůrazněno, objasní další části textu. Uvažujme nyní lineární prostředí, tedy takové, jehož vlastnosti nezávisí na intenzitě světelné vlny šířící se tímto prostorem. Tento případ nastává, pokud intenzita vlny není příliš velká. Potom platí tzv. princip superpozice pro vektory elektromagnetického pole [1 – 5]: „Výsledný vektor intenzity pole v daném místě prostoru, v němž se nachází několik zdrojů světla, je roven vektorovému součtu dílčích vektorů intenzit, které v daném místě vytváří jednotlivé zdroje světla navzájem nezávisle.“ Tedy pro případ vektoru intenzity elektrického pole E(r, t) v bodě r, čase t a k zdrojů světla můžeme psát k
E(r, t ) = ∑ Ei (r, t ) .
(60)
i =1
24
Interference světla V případě silných elektromagnetických polí princip superpozice platit nebude. Např. u výkonných laserů je velikost vektoru intenzity elektrického pole E řádově 107 – 1011 V/m [4], což je srovnatelné s intenzitou pole uvnitř atomů. Pro účely této práce ale takováto pole neuvažujeme, vystačíme tedy s principem superpozice (60). Popisem výsledné intenzity elektromagnetického vlnění I, které vznikne složením více polí splňujících princip superpozice, se budou zabývat následující části.
3.3
Interference a kontrast lineárně polarizovaných monochromatických vln stejné frekvence
Představme nyní princip superpozice pro určení intenzity výsledného pole v případě dvou lineárně polarizovaných [1 – 5] rovinných monochromatických vln stejné úhlové frekvence, jejichž vektory E1 a E2 jsou rovnoběžné, ale šíří se v různých směrech daných vlnovými vektory k1 a k2. Tento poměrně jednoduchý případ se v praxi velmi využívá. V laboratořích je možné snadno vytvořit vlny lineárně polarizované a téměř ideálně monochromatické. S užitím (49) pro takovéto dvě vlny píšeme E1 = E 01 exp[i(k 1r − ωt + δ 1 )] ,
(61)
E 2 = E 02 exp[i(k 2 r − ωt + δ 2 )] ,
kde E01 a E02 značí počáteční amplitudy a δ1 a δ2 počáteční fáze, r je polohový vektor od počátku soustavy souřadnic. Výsledné pole je dáno principem superpozice (60) jako
E = E1 + E 2 .
(62)
Pro jeho intenzitu poté s užitím (47), (50) a (62), neuvažujeme-li absorpci prostředí, píšeme I = N = E× H =
ε ε 2 E = E1 + E 2 µ µ
2
=
ε (E1 + E 2 )(E1 + E 2 )∗ . µ
(63)
Provedením naznačených početních operací dostáváme
(E
1
+ E2 )(E1 + E2 ) = E1 + E 2 + E1E*2 + E2E1* = E1 ∗
2
2
2
+ E2
2
+ E1E*2 + E 2E1* .
(64)
Pokud platí
E1E*2 + E 2 E1* ≠ 0 ,
(65) 25
Interference světla potom dochází k tzv. interferenci světla. Prostor, ve kterém dochází k interferenci, je nazýván interferenční prostor. Levá strana (65) je tzv. interferenční člen [4]. Pokud je tento člen nulový, nedochází k interferenci. Užitím vztahů (61) pro něj dostáváme E1E*2 + E 2E1* = E01E02 exp{i[(k 1 − k 2 ) r + δ1 − δ 2 ]} + E02 E01 exp{- i[(k 1 − k 2 ) r + δ1 − δ 2 ]} ,
(66)
a dále užitím Eulerova vzorce E1E*2 + E 2 E1* = E 01E 02 {cos[(k 1 − k 2 ) r + δ 1 − δ 2 ] + i sin [(k 1 − k 2 ) r + δ 1 − δ 2 ]}
+ E 02 E 01 {cos[(k 1 − k 2 ) r + δ1 − δ 2 ] − i sin [(k 1 − k 2 ) r + δ 1 − δ 2 ]}
= 2E 01E 02 cos[(k 1 − k 2 ) r + δ 1 − δ 2 ] .
(67)
Výslednou intenzitu v místě r tedy můžeme vyjádřit jako
I (r ) =
ε E1 µ
2
ε E2 µ
+
2
+
ε 2E 01 E 02 cos[(k 1 − k 2 )r + δ 1 − δ 2 ] µ
(68)
= I 1 (r ) + I 2 (r ) + 2 I 1 (r ) I 2 (r ) cos[(k 1 − k 2 )r + δ 1 − δ 2 ] , kde I1(r) a I2(r) značí intenzity první a druhé vlny v místě r a využili jsme rovnoběžnosti vektorů počátečních amplitud. Jak je z (68) patrné, bude intenzita konstantní, pokud
(k
1
− k 2 ) r + δ 1 − δ 2 = 2πm + C ,
(69)
kde m je celé číslo a C je konstanta. Maximum intenzity bude pro C = 0 a minimum pro C = π. Porovnáním s vektorovou rovnicí roviny známou z matematiky vidíme, že (69) určuje posloupnost rovin s konstantní intenzitou a tyto roviny jsou kolmé na vektor (k1 – k2). Pro určení vzdálenosti dvou po sobě následujících rovin uvažujme rovinu m a (m + 1). Dostáváme tedy
(k
(k
1
− k 2 )rm + δ 1 − δ 2 = 2 πm + C ,
− k 2 )rm+1 + δ 1 − δ 2 = 2 π(m + 1) + C . 1
(70)
Vzdálenost Δ těchto rovin je rovna projekci vektoru (rm+1 – rm) na jednotkový vektor kolmý k rovinám (normálový vektor) (k1 – k2). Tedy
∆=
(k
− k2 ) (rm+1 − rm ) . k1 − k 2 1
(71)
Odečteme-li od sebe vzájemně rovnice (70), dostaneme
(k
1
− k 2 )(rm+1 − rm ) = 2 π .
(72)
26
Interference světla Dosazení do (71) následně vede ke vztahu
∆=
2π λ = , k1 − k 2 n1 − n 2
(73)
kde jsme dále použili definice vlnového vektoru k = kn, k = 2π/λ, n1 a n2 značí jednotkové normálové vektory rovinných vlnoploch. Pro čtverec velikosti jejich rozdílu dále platí n1 − n 2 = (n1 − n 2 )(n1 − n 2 ) = n1n1 + n 2n 2 − 2n1n 2 = 2 − 2 cos α = 4 sin 2 2
α 2
,
(74)
kde α značí úhel mezi normálovými vektory. Dosazením do (73) pro vzdálenost mezi rovinami konstantní intenzity dostáváme
λ
∆=
2 sin
α
.
(75)
2
Ze vztahu (69) můžeme vyčíst, že přenos energie probíhá v periodických vrstvách kolmých na vektor (k1 – k2), také tedy na (n1 – n2). Vložíme-li do interferenčního pole detektor záření, budeme registrovat rozložení intenzity, které se projeví tzv. interferenčními proužky. Vzdálenost dvou sousedních interferenčních proužků registrovaných detektorem bude závislá na sklonu tohoto detektoru vůči vektoru (n1 + n2) (vůči ose úhlu α). Bude-li normála detektoru svírat s tímto vektorem úhel β, potom detekovaná vzdálenost mezi proužky bude
∆β =
∆ . cos β
(76)
Na obr. 5 je zobrazen ukázkový interferenční obrazec dvou rovinných harmonických vln.
Obr. 5 Interferenční obrazec pro dvě rovinné harmonické vlny zeleného světla o úhlové frekvenci ω = 4∙1015 Hz v rovině (y,z) kolmé na osu x ve vzdálenosti d = 1 od počátku soustavy souřadnic; jednotkové normálové vektory rovin voleny n1 = [1,0,0], n2 = 1/√(5)[2,0,1], δ1 = δ2 = 0, I1 = I2 = 1
27
Interference světla Pro určení kontrastu K specifikujme přesněji hodnoty Imax a Imin tím způsobem, že Imax představuje maximální hodnotu jednoho interferenčního proužku (jeden z bílých proužků na obr. 5) a Imin minimální hodnotu intenzity proužku vedlejšího (tmavý proužek). S užitím (69) tedy můžeme psát I max = I 1 + I 2 + 2 I 1 I 2 cos[2 πm] , I min = I 1 + I 2 + 2 I 1 I 2 cos[π(2m + 1)] ,
(77)
a tedy kontrast pro interferenční obrazec vzniklý interferencí dvou lineárně polarizovaných monochromatických harmonických vln stejné frekvence, jejichž vektory amplitud jsou rovnoběžné, bude dán vztahem [1, 2, 4 – 9]
K=
2 I1 I 2 I1 + I 2
.
(78)
Jak ale bylo řečeno, nemusí k interferenci dojít vždy. Podmínka (65) říká, že k interferenci dochází, pokud interferenční člen není roven nule. Pro rovinné monochromatické vlny (61) poté s (67) pro levou stranu (65) dostáváme
E 1 E *2 + E 2 E 1* = 2E 01 E 02 cos[(k 1 − k 2 ) r + δ 1 − δ 2 ] 2T E 01 E 02 cos[(k 1 − k 2 ) r + δ 1 − δ 2 ] dt T ∫0 = 2 E 01 E 02 cos[(k 1 − k 2 ) r + δ 1 − δ 2 ] , =
(79)
kde jsme využili rovnoběžnosti vektorů počátečních amplitud. Aby k interferenci došlo, musí platit podmínka 2 E01 E02 cos[(k 1 − k 2 ) r + δ 1 − δ 2 ] ≠ 0 .
(80)
Jelikož E01 a E02 nebudou nulové, bude (80) splněno právě tehdy, když
(k
1
− k 2 )r + δ1 − δ 2 ≠
π + lπ , 2
(81)
kde l je celé číslo. Lze snadno ukázat, že pro lineárně polarizovaná záření, jejichž roviny polarizace nebudou rovnoběžné, jak tomu bylo v předchozím případě, ale budou svírat úhel ψ, bude platit pro intenzitu vztah
I (r ) = I 1 (r ) + I 2 (r ) + 2 I 1 (r ) I 2 (r ) cos[(k 1 − k 2 )r + δ 1 − δ 2 ]cosψ
(82)
28
Interference světla a pro kontrast poté
K=
3.4
2 I1 I 2 cosψ . I1 + I 2
(83)
Interference a kontrast monochromatických elipticky polarizovaných vln stejné frekvence
Představme nyní případ dvou interferujících elipticky polarizovaných monochromatických vln stejné frekvence [1 – 5]. Případ kruhově nebo lineárně polarizovaných vln je jen speciálním případem tohoto. Vlny budou popsány stejnými rovnicemi jako (61), ovšem počáteční amplitudy nebudou rovnoběžné. Vlnové vektory spolu svírají úhel α. Rozložíme-li výsledný vektor intenzity elektrického pole do směru kolmého a rovnoběžného s rovinou, jež je dána vlnovými vektory k1 a k2, potom dostáváme E = E ⊥ + E|| = (E1⊥ + E 2⊥ ) + (E1|| + E 2|| ) ,
kde
⊥
(84)
značí složky kolmé, || složky rovnoběžné. Pro velikosti jednotlivých složek užitím kosinové
věty platí E⊥2 = E ⊥ = E1⊥ + E 2⊥ + 2 E1⊥ E 2⊥ , 2
2
2
(85)
E||2 = E|| = E1|| + E 2|| + 2 E1|| E 2|| cos α . 2
2
2
Pro intenzitu dostáváme I (r ) =
[
]
ε 2 ε 2 ε E (r ) = E ⊥ (r ) + E||2 (r ) = E⊥2 (r ) + E||2 (r ) = I ⊥ (r ) + I || (r ) . µ µ µ
(86)
Výsledná intenzita je tudíž dána jako součet intenzit v jednotlivých složkách. Použijeme-li pro ně vztah (68) – uvažujeme-li každou ze složek jako intenzitu rovinné harmonické vlny, potom můžeme s použitím (82) psát
I ⊥ (r ) = I1⊥ (r ) + I 2⊥ (r ) + 2 I1⊥ (r ) I 2⊥ (r ) cos[(k 1 − k 2 ) r + δ 1⊥ − δ 2⊥ ] ,
I || (r ) = I1|| (r ) + I 2|| (r ) + 2 I1|| (r ) I 2|| (r ) cos[(k 1 − k 2 ) r + δ 1|| − δ 2|| ]cos α .
(87)
29
Interference světla Zavedením substituce ∆ = (k 1 − k 2 ) r + δ 1|| − δ 2|| a Φ = δ 1⊥ − δ 2⊥ − δ 1|| + δ 2|| lze výslednou intenzitu v bodě r vyjádřit jako I (r ) = I ⊥ (r ) + I || (r ) = I 1⊥ (r ) + I 2 ⊥ (r ) + I 1|| (r ) + I 2|| (r )
(88)
+ 2 I 1⊥ (r ) I 2 ⊥ (r ) cos(∆ + Φ ) + 2 I 1|| (r ) I 2|| (r ) cos ∆ cos α . Rozepsáním a úpravou vztahu (88) můžeme intenzitu elipticky polarizované vlny vyjádřit také jako I (r ) = I1⊥ (r ) + I 2⊥ (r ) + I1|| (r ) + I 2|| (r )[1 + K cos(∆ + δ )] ,
(89)
kde K je kontrast interferenčního pole daný výrazem
K=
2
(
)
2
I1⊥ (r ) I 2⊥ (r ) cos Φ + I1|| (r ) I 2|| (r ) cos α + I1⊥ (r ) I 2⊥ (r ) sin 2 Φ I1⊥ (r ) + I 2⊥ (r ) + I1|| (r ) + I 2|| (r )
(90)
a δ je konstantní fázový člen, pro který platí
δ = arctan
I 1⊥ (r ) I 2⊥ (r ) sin Φ I 1⊥ (r ) I 2⊥ (r ) cos Φ + I 1|| (r ) I 2|| (r ) cos α
.
(91)
Výše uvedené vztahy lze snadno převést na lineárně nebo kruhově polarizované vlny. Lineární polarizace ve stejné rovině nastává tehdy, pokud I 1|| (r ) = I 2|| (r ) = 0 . Potom pro intenzitu a kontrast dostáváme
I (r ) = I1⊥ (r ) + I 2⊥ (r ) + 2 I1⊥ (r ) I 2⊥ (r ) cos[(k 1 − k 2 ) r + δ 1⊥ − δ 2⊥ ] , K=
2 I1⊥ (r ) I 2⊥ (r ) . I1⊥ (r ) + I 2⊥ (r )
(92)
(93)
Pro kruhově polarizované vlny bude platit I 1⊥ (r ) = I 1|| (r ) = I 1 (r ) a I 2 ⊥ (r ) = I 2|| (r ) = I 2 (r ) , dosazení do výše uvedených vztahů poté vede k
I (r ) = 2[I 1 (r ) + I 2 (r )][1 + K cos(∆ + δ )] ,
(94)
kde
K=
(
)
2
I1 (r ) I 2 (r ) cos Φ + I1 (r ) I 2 (r ) cos α + I1 (r ) I 2 (r ) sin 2 Φ I1 (r ) + I 2 (r )
,
(95)
30
Interference světla
δ = arctan
3.5
I1 (r ) I 2 (r ) sin Φ I1 (r ) I 2 (r ) cos Φ + I1 (r ) I 2 (r ) cos α
.
(96)
Koherence vlnění
Výše odvozené vztahy nejsou zcela kompletní, přesněji řečeno v reálných aplikacích neplatí úplně. Podíváme-li se na ně podrobněji, musíme je doplnit předpokladem, že po celou dobu průběhu interference byla vlnění tzv. koherentní. Výše uvedené vztahy platí za ideálního případu tzv. absolutní koherence. Rozumíme tím, že během celého šíření daného vlnění nedošlo ke změně frekvence, amplitudy ani počáteční fáze. Reálná záření ovšem tyto ideální vlastnosti nemají, alespoň k minimálním změnám dochází. Vztahy pro intenzitu a kontrast je tedy nutné doplnit. Koherence vlnění vyjadřuje vzájemnou korelaci (souvislost) vlnění, které vychází ze dvou míst světelného zdroje (např. jeden světelný zdroj a stínítko se dvěma štěrbinami – obr. 6), potom mluvíme o tzv. prostorové koherenci [1 – 9], nebo korelaci vlnění vycházejícího z jednoho zdroje zkoumaného s časovým odstupem, poté hovoříme o tzv. časové koherenci [1 – 9].
3.5.1 Prostorová koherence Uvažujme situaci zobrazenou na obr. 6.
pozorovací rovina (CCD senzor)
stínítko B1 (E1)
ρ1, τ1
S
ρ2, τ2
zdroj záření B2 (E2) Obr. 6
Prostorová koherence [4]
Nechť ze zdroje záření vychází vlnění, které dopadá na stínítko, v němž jsou prakticky bodové otvory B1 a B2. Vlnění v bodě B1 charakterizuje vektor elektrické intenzity E1 a obdobně v bodě B2 vektor E2. Z otvorů se vlnění šíří dál a dopadá na pozorovací rovinu (např. CCD senzor).
31
Interference světla Zhodnoťme nyní, jaká bude vzájemná korelace (koherence) vlnění jdoucích z bodů Bi do bodu S na senzoru po drahách ρi (i = 1, 2). Pro jednoduchost uvažujme lineárně polarizované vlnění, čili můžeme bez ztráty na obecnosti přistoupit ke skalárnímu popisu. Čas, který vlnění E1 potřebuje k překonání vzdálenosti ρ1 mezi body B1 a S označme τ1, obdobně pro vlnění E2 použijeme označení τ2. Aplikace principu superpozice poté pro vlnění v bodě S v čase t vede k
E (t ) = E1 (t − τ 1 ) + E2 (t − τ 2 ) .
(97)
Pro intenzitu následně dostáváme při skalární formulaci vztahu (47) I = E (t ) E ∗ (t ) = [E1 (t − τ 1 ) + E 2 (t − τ 2 )][E1 (t − τ 1 ) + E 2 (t − τ 2 )]
∗
= E1 (t − τ 1 ) E1∗ (t − τ 1 ) + E 2 (t − τ 2 ) E 2∗ (t − τ 2 )
(98)
+ E1 (t − τ 1 ) E 2∗ (t − τ 2 ) + E1∗ (t − τ 1 ) E 2 (t − τ 2 ) . Uvědomíme-li si, že posun počátku odečítání času nebude mít vliv na integrální střední hodnotu ve smyslu
E1 (t − τ 1 ) E1∗ (t − τ 1 ) = E1 (t ) E1∗ (t ) ,
(99)
I 1 = E1 (t ) E1∗ (t ) , I 2 = E 2 (t ) E 2∗ (t )
(100)
můžeme se zavedením označení
a posunutím počátku o τ2 psát I = I 1 + I 2 + E1 (t − τ 1 + τ 2 ) E 2∗ (t ) + E1∗ (t − τ 1 + τ 2 ) E 2 (t )
= I 1 + I 2 + E1 (t + τ ) E 2∗ (t ) + E1∗ (t + τ ) E 2 (t ) ,
(101)
kde τ = τ2 – τ1. U členu I1 jsme využili vlastnost (99) v tom smyslu, že platí
E1 (t − τ 1 + τ 2 ) E1∗ (t − τ 1 + τ 2 ) = E1 (t ) E1∗ (t ) = I1 .
(102)
Úpravou posledního členu ve (101) poté pro výslednou intenzitu v bodě S dostáváme
I = I1 + I 2 + 2 Re E1 (t + τ ) E2∗ (t ) .
(103)
Korelace vyšetřovaných vln bude tedy dána posledním členem 2 Re E1 (t + τ ) E2∗ (t ) . Zavedeme-li nyní tzv. funkci vzájemné koherence [1 – 9] Γ12(τ) vztahem
Γ12 (τ ) = E1 (t + τ ) E2∗ (t )
(104)
32
Interference světla a tzv. komplexní stupeň koherence γ12(τ) výrazem [1 – 9] 1 Γ12 (τ ) = I1 I 2
γ 12 (τ ) =
1 E1 (t + τ ) E2∗ (t ) , I1 I 2
(105)
potom pro intenzitu v bodě S dostáváme
I = I1 + I 2 + 2 I1 I 2 Re[γ 12 (τ )] .
(106)
Užitím Eulerova vzorce můžeme psát pro Re[γ12(τ)]
Re[γ 12 (τ )] = γ 12 (τ ) cos ϕ12 ,
(107)
kde φ12 je argument komplexního čísla γ12(τ). Porovnání se skalární formou (68) vede pro tento argument k
ρ1
ϕ12 = 2π
λ1
−
2π ρ2 + δ1 − δ 2 = (ρ1 − ρ2 ) + ∆δ12 (τ ) , λ0 λ2
(108)
kde λ0 je střední vlnová délka vlnových délek λ1 a λ2, Δδ12(τ) je fázový rozdíl mezi vlněním v otvoru B1 a v otvoru B2, ρi (i = 1, 2) značí vzdálenost mezi bodem Bi a S. Pro intenzitu v bodě S dosazením ze (107) do (106) dostáváme I = I1 + I 2 + 2 I1 I 2 γ 12 (τ ) cos ϕ12 .
(109)
Veličina |γ12(τ)| se bude měnit v intervalu od 0 do 1. Z (109) je patrné, že pro |γ12(τ)| = 0 nedojde k interferenci. O takovém vlnění vycházejícím z bodů B1 a B2 říkáme, že je nekoherentní. Pokud |γ12(τ)| = 1, je takové vlnění tzv. absolutně koherentní. Jestliže bude platit nerovnost 0 < |γ12(τ)| < 1, potom nazýváme vlnění částečně koherentní. Užitím znalostí o oboru hodnot funkce cosφ12 můžeme vyslovit závěry o maximální hodnotě intenzity v bodě Imax nebo naopak o její minimální hodnotě Imin. Platí I max = I1 + I 2 + 2 I1 I 2 γ 12 (τ ) , I max = I1 + I 2 − 2 I1 I 2 γ 12 (τ ) .
(110)
Dosazením do vztahu (48) dostáváme tedy pro kontrast
K=
I max − I min 2 I1 I 2 = γ 12 (τ ) . I max + I min I1 + I 2
(111)
33
Interference světla Výše uvedený vztah (111) nám říká, že kontrast je přímo úměrný modulu stupně koherence |γ12(τ)|. V praxi se nám toto projeví méně či více kontrastními zaznamenanými proužky v interferenčním poli. Speciální případ nastane, pokud bod S leží na ose bodů B1 a B2. Potom je časové zpoždění τ = 0 a platí
γ 12 (0) =
1 Γ12 (0) = I1 I 2
1 E1 (t ) E2∗ (t ) . I1 I 2
(112)
Všechna monochromatická vlnění jsou absolutně koherentní. Proto se nám modul stupně koherence nevyskytoval ve vztazích odvozených v předchozích částech práce. Prakticky všechny případy vlnění používaného v praxi jsou částečně koherentní. Na reálných datech tedy pozorujeme určité rozložení kontrastu interferenčních proužků.
3.5.2 Časová koherence Časová koherence popisuje korelaci (závislost) vlnění vycházejícího z jednoho zdroje zkoumaného v jednom bodě v různé časové okamžiky. Komplexní stupeň koherence γ12(τ) bude mít v tomto případě tvar2
γ 11 (τ ) =
1 1 Γ11 (τ ) = E1 (t + τ ) E1∗ (t ) . I1 I1
(113)
Předpokládejme nyní bodový zdroj světla (např. atom), jenž vyzařuje vlnění o střední frekvenci
ν0 po dobu t0. Takto vzniklé vlnové pulzy, které vysílá atom zcela nahodile, nejsou navzájem korelovány a nemůže mezi nimi dojít k interferenci [4]. Situace je vhodným modelem, která nám dobře popisuje reálnou situaci. Uvažujme nyní následující případ. Záření z bodového světelného zdroje postupuje po dvou cestách a dopadá na senzor v bodě S. První je cesta přímá a světlo urazí od zdroje vzdálenost ρ1. Druhá nechť vede přes libovolnou optickou soustavu zrcadel a světlo urazí vzdálenost ρ2. Předpokládejme dále stejné prostředí v obou cestách. Vyšetřovaná situace je zobrazena na obr. 7. Dráhový rozdíl δρ mezi oběma vlněními bude dán vztahem
δρ = ρ1 − ρ2 .
2
(114)
Index 2 přejde v index 1, druhý zdroj je totožný s prvním.
34
Interference světla
bodový zdroj záření
ρ1 S ρ2
pozorovací rovina (CCD senzor)
optická soustava Obr. 7
Časová koherence [4]
Pro rychlost šíření vlnění c potom dostáváme pro dobu δt potřebnou k uražení tohoto dráhového rozdílu
δt = δρ c .
(115)
Bude-li tato doba menší než doba t0, po kterou bodový zdroj (atom) vyzařuje vlnový pulz (δt < t0), potom se spolu setkají v bodě S dvě části stejného vlnového pulzu, které spolu mohou interferovat. V případě opačném (δt > t0) interference není možná, protože mezi nahodile vyslanými různými pulzy nemůže nastat žádná korelace. Délka trvání jednoho vlnového pulzu je nazývána koherenční čas a je označována jako τkoh [1 – 5]. Je-li doba potřebná k uražení dráhového rozdílu mezi dvěma cestami z bodového zdroje do jednoho bodu v prostoru menší než koherenční čas τkoh, může dojít k interferenci. Vzdálenost, kterou záření urazí za koherenční čas τkoh,, je nazývána koherenční délka lkoh [1 – 5]. Lze ji vyjádřit ve tvaru [1 – 5] lkoh = cτ koh =
λ20 , π∆λ
(116)
kde c je rychlost šíření vlnění, λ0 je střední vlnová délka a Δλ je spektrální šířka vlnění. Chceme-li tedy dosáhnout dobrého kontrastu v interferenčním obrazci, je třeba použít takový zdroj světla, jehož koherenční délka bude mnohem větší, než je největší dráhový rozdíl mezi interferujícími vlnami. Musí tedy platit lkoh >> δρ .
(117) 35
Interference světla Jak je vidět z rovnice (116), koherenční délka bude tím větší (a tím lepší kontrast interferenčního obrazce), čím menší bude spektrální šířka (rozsah vlnových délek). V praxi jsou používány laserové zdroje záření, jejichž koherenční délka je řádově 106 m a více [4]. Jsme s nimi tak schopni zajistit dobré podmínky pro získávání dat k dalšímu vyhodnocení.
3.6
Intenzita a kontrast polychromatického a kvazimonochromatického záření
Výše uvedené vztahy měly za úkol vnést nadhled do studované problematiky vzniku interferenčních proužků. Nyní jsme schopni vyslovovat předpoklady a závěry i u komplikovanějších (reálných) situací. Nejčastěji v praxi používané laserové zdroje jsou generátory tzv. kvazimonochromatického záření3, některé aplikace vyžadují i použití záření polychromatického4. Intenzitu v interferenčním poli dvou stacionárních polychromatických záření (neuvažujeme-li polarizační vlastnosti a můžeme tedy použít skalární popis) v bodě Q můžeme vyjádřit pomocí tzv. obecného interferenčního zákona (106) jako [1 – 9]
I (Q) = I1 (Q) + I 2 (Q) + 2 I1 (Q) I 2 (Q) Re[γ 12 (τ )] ,
(118)
kde Re[γ12(τ)] je reálná část komplexního stupně vzájemné koherence, 0 ≤ γ 12 (τ ) ≤ 1 , τ je vzájemné časové zpoždění. Pro kontrast poté platí
K=
2 I1 I 2 Re[γ 12 (τ )] . I1 + I 2
(119)
Pro kvazimonochromatické záření dostáváme pro intenzitu v bodě Q s použitím vztahů (108) a (109)
[
]
I (Q) = I1 (Q) + I 2 (Q) + 2 I1 (Q) I 2 (Q) γ 12 (τ ) cos 2 πft + ∆δ12 (τ ) ,
(120)
kde ∆δ12 (τ ) a γ 12 (τ ) jsou funkcemi vzájemného časového zpoždění τ pomalu proměnné vůči funkci cos(2 πft ) , f je střední frekvence záření. Pro kontrast následně platí
Kvazimonochromatické záření je záření velmi úzké spektrální šířky. Jinými slovy spektrální šířka Δλ je mnohem menší, než vlnová délka daného záření λ.
3
4
Polychromatické záření je záření více vlnových délek.
36
Interference světla
K=
2 I1 I 2 γ 12 (τ ) . I1 + I 2
(121)
Kvazimonochromatické záření je tzv. částečně koherentní a kontrast interferenčního pole závisí na hodnotě γ 12 (τ ) .
37
4 Interferometrie Interferometrie je bezkontaktní optická metoda využívající interferenci světla, která umožňuje ve výsledku určovat relativní vzdálenosti s velmi vysokou přesností. Vlnová pole v optice, která jsou popsatelná teorií elektromagnetického pole [1 – 9, 51 – 53], jsou charakterizována (definována) parametry, které hovoří o jejich vlastnostech a chování. Především mezi tyto parametry patří amplituda, fáze, frekvence a polarizační stav daného vlnového pole. Odraz od plochy, kterou buď přímo vyšetřujeme, nebo se nějakým způsobem podílí v procesu měření, způsobuje změnu těchto parametrů. Budeme-li schopni určit změnu fáze vlny odražené od vyšetřované plochy oproti fázi vlny odražené od plochy referenční, získáme tak kvantitativní informace o rozdílnosti tvarů těchto ploch. V praxi ovšem nejsme schopni měřit přímo frekvenci nebo fázi (čili i její změnu). Neexistují totiž tak citlivé detektory, jejichž odezva by byla porovnatelná s periodou oscilací použitého optického záření. Tyto informace (zejména fázi nebo frekvenci) tak získáváme nepřímými metodami, kdy měříme pouze intenzitu záření (časovou střední hodnotu energie, jak bylo představeno v dřívější části práce). Potom tzv. inverzními postupy zpětně určujeme fázi nebo frekvenci těchto vlnových polí ze zaznamenaných hodnot intenzity. Jednou z široké skupiny metod vyvíjených v praxi pro získávání informací o fázi je zmíněná metoda interferometrie založená na vyhodnocení interferenčních obrazců vzniklých interferencí částečně koherentních vlnových polí [1 – 9]. Omezující podmínkou pro tyto metody je právě koherence použitých vlnových polí. Tato podmínka je nutná, protože jak bylo ukázáno v předchozí kapitole, intenzita v interferenčním poli polychromatického a zejména využívaného kvazimonochromatického záření (záření laserů) je závislá na hodnotě komplexního stupně vzájemné koherence. Nevýhodou interferometrických měření je také vysoká citlivost na termomechanické vlastnosti prostředí (změna vlhkosti, teploty, tlaku, nepříznivé vibrace, apod.) a omezený rozsah vyhodnocené fáze. Pro maximální přesnost jsou dále kladeny zvýšené požadavky na realizaci měření (laboratorní podmínky). I když při realizaci měření dochází k některým poměrně závažným komplikacím, je interferometrie v technické praxi používána právě pro svou velmi vysokou přesnost, a to v řádu nanometrů – zlomků vlnových délek použitého záření [6 – 10]. S rozvojem technologií bylo možno vyvinout algoritmy a metody, které jsou náročné na výpočetní kapacity a které dovolují potlačovat nepříznivé vlivy. V důsledku toho se interferometrie stala běžnou součástí všude tam, kde je třeba znát velmi přesné prostorové vztahy [10 – 36]. Nejen v „pozemských“ aplikacích nachází tato metoda uplatnění. Jmenujme dále například SAR nebo InSAR interferometrii 38
Interferometrie [37, 38], jež slouží k získávání prostorových informacích o topografii zemského terénu. Interference zde nastává při záznamu signálu odraženého od povrchu porovnáním s referenčním signálem na družici. V další části práce se budeme zabývat zejména interferometrií používanou v praxi aplikované optiky. Rozumějme tedy měření při výrobě a kontrole optických prvků – jak rovinných (planárních), sférických, tak i asférických. Principy interferometrie popsané v této a v dalších kapitolách jsou ale převeditelné i do jiných aplikací. Bez možnosti dostatečně přesně změřit tvar vyráběné optické plochy nejsme schopni tuto plochu, a tím i další komponenty, vyrobit, protože nebudeme schopni formulovat její vliv na výsledný záměr. Nejčastěji máme ve výrobě zadán nominální tvar optické plochy a cílem kontrolních metod (tedy i interferometrie) je poté určit odchylky od tohoto tvaru. Ty vznikají jako důsledek technologických procesů výroby. Na výkresech jednotlivých součástí jsou uváděny přesnosti, s jakou mají být vyrobeny, a cílem je potvrdit nebo vyvrátit jejich dodržení v rámci povolených výrobních tolerancí. Jeden ze způsobů, využívající interferenci světla, je porovnání s tzv. kalibrem, kontrolním prvkem, který je vyroben o jeden nebo více řádů přesněji. Kalibr je přikládán na kontrolovanou optickou plochu a při dobrém osvětlení jsou patrné interferenční proužky, z jejichž tvaru lze posuzovat míru dodržení tolerancí [6 – 10]. Touto „pohledovou“ zkouškou může zkušený vyhodnocovatel odhadnout velikost výrobní chyby s rámcovou přesností až λ/8 [10], kde λ značí vlnovou délku světla použitého při kontrole. Takovýto způsob kontroly je však časově i finančně náročný a pro každý tvar vyráběné plochy musí být vyroben speciální kalibr. Navíc může dojít k poškození vyráběné plochy nebo kalibru při jejich manipulaci. V minulosti bylo snahou takovéto a podobné kontrolní metody nahrazovat způsoby bezkontaktními a plně automatizovanými a je tomu tak i dnes. Bezkontaktní kontrolu vyráběných prvků lze realizovat použitím vhodně konstruovaných interferometrů, a to buď Fizeauova nebo Twyman-Greenova typu [1 – 10], většinou doplněných o další prvky, např. difraktivní optické elementy pro měření rotačně symetrických asférických ploch [54, 55]. S vývojem techniky, zejména CCD senzorů (z angl. Charge-Coupled Device [56]) vysokého rozlišení a piezoelektrických posuvů [57], se u těchto interferometrů přešlo k plně automatizovanému analytickému vyhodnocování tvaru kontrolovaných ploch. Využívají se vhodné algoritmy umožňující určit fázi vyšetřovaného vlnového pole na základě záznamu
39
Interferometrie intenzity pole interferenčního, které vzniká v důsledku interference polí vyšetřovaného a referenčního. Tato fáze přímo souvisí s kvalitou testované optické plochy. Přesnost těchto interferometrických analytických metod je v řádu λ/10 až λ/100 a více [6 – 10]. Vyhovuje tak naprosté většině požadavků na měření topografie ploch v optickém průmyslu. V další části se seznámíme s různými praktickými variantami získávání a vyhodnocování dat pořizovaných interferometrickými metodami.
4.1
Dvousvazková a střihová interferometrie
Tato část má za úkol představit dva základní principy vzniku interferenčních polí v praxi kontroly optických prvků. Respektive způsob realizace měření takovým způsobem, aby k vzniku interferenčních polí došlo. Těmito principy jsou: •
dvousvazková interferometrie,
•
střihová interferometrie.
Poznamenejme zde, že možných metod je celá řada, ne jen tyto dvě. V práci se o nich ale zmiňovat nebudeme, pro další studium odkážeme čtenáře na příslušnou literaturu [1 – 10].
4.1.1 Dvousvazková interferometrie Dvousvazková interferometrie využívá k vzniku interferenčního pole dva svazky záření, jak už název napovídá. Zpravidla vznikají dělením počátečního svazku na dva, které po průchodu optickou soustavou společně interferují. Tímto postupem lze nastavit takové podmínky, aby interference nastala. Jeden, zpravidla laserový, zdroj totiž umožňuje snadno generovat lineárně polarizované vstupní záření a také zaručuje dostatečnou koherentnost svazků. V kontrole optické výroby se nejčastěji používají interferometry Twyman-Greenova nebo Fizeauova typu [1 – 10], jejichž principielní schéma je zobrazeno na obr. 8 a obr. 9. V případě Twyman-Greenova interferometru (představený roku 1918 [6]) svazek generovaný ve zdroji záření prochází přes objektiv na kolimátor, kde je transformován na rovinné vlny. Dále pokračuje na dělič svazků. Na něm vždy dojde k částečnému průchodu a částečnému odrazu. Jedna část vlnění tedy projde děličem a pokračuje k testované ploše, druhá část děleného svazku se odráží a jde k referenční ploše. Od testované i referenční plochy se následně paprsky odrážejí, na děliči se jeden opět láme a druhý projde přímo. Tímto způsobem prochází tyto svazky až 40
Interferometrie k pozorovací rovině, kde interferují, a my zaznamenáváme rozložení intenzity interferenčního pole.
referenční plocha dělič svazků
kolimátor zdroj záření objektiv (laser)
testovaná plocha
pozorovací rovina (CCD senzor) Obr. 8
Základní konfigurace interferometru Twyman-Greenova typu [6]
Optické komponenty musí být co nejčistší, aby se pokud možno neprojevovaly nepříznivé difrakční jevy. Dále zejména dělič svazků musí být vysoké kvality. Pro potlačení „parazitní“ interference se nejčastěji používají antireflexní vrstvy. Také materiál děliče musí být extrémně homogenní. Od roviny, na které dochází k odrazu, se očekává, že je vyrobena s minimálně dvojnásobnou přesností, než jakou požadujeme od vyhodnocovaných vlnoploch [6]. Základní konfigurace Fizeauova interferometru je zobrazena na obr. 9. Jedná se také o dvousvazkový interferometr, jak tomu bylo v případě Twyman-Greenova, poloha referenčního elementu je ale jiná, testovaná i referenční plocha leží „v jednom rameni“. Interferometry se tedy liší dráhovým rozdílem. V případě testování ideální rovinné planparalelní skleněné desky platí pro Fizeaův interferometr vztah [6 – 9]
δρ F = 2nd ,
(122)
kde n je index lomu skla desky a d je její tloušťka. Pro Twyman-Greenův interferometr poté [6 – 9] δρTG = 2(n − 1)d .
(123)
41
Interferometrie Pro vyčerpávající popis testování optických elementů odkážeme čtenáře na dostupnou odbornou literaturu, např. [6 – 9, 11].
zdroj záření (laser)
objektiv
dělič svazků
kolimátor
pozorovací rovina (CCD senzor) Obr. 9
referenční testovaná plocha plocha
Principielní schéma interferometru Fizeauova typu [6]
Referenční vlnoplocha nám tedy vzniká odrazem od referenční plochy. Porovnáním s vlnoplochou odraženou od plochy testované poté můžeme usuzovat o odchylkách tvaru testované plochy od referenční. Vyhodnocením interferogramu dostáváme informace o fázovém rozdílu vlnoplochy testované a referenční, což nám nese informace i o rozdílu dráhovém. Pomocí Twyman-Greenova a Fizeauova interferometru můžeme testovat různé optické elementy. Nejsnazší je testování planparalelních desek nebo členů se sférickým povrchem. Měření topografie asférických ploch je komplikovanější a nebyla prozatím vyvinuta obecná metoda, která by umožňovala testovat se stejnou přesností libovolnou asférickou plochu [10]. Pro testování asférických optických ploch jsou užívány speciální adaptery a difraktivní optické elementy, které jsou navrženy pro každou testovanou asférickou plochu zvlášť (pro tvar této plochy). Tyto elementy jsou ale velmi drahé a jejich výroba je časově náročná (několik týdnů). Komerčně jsou pro účely interferometrické bezkontaktní kontroly a měření topografie asférických ploch v optice dostupné přístroje např. od firmy ZYGO [58], 4D-Technology [59], ESDI [60], Schneider [61], OptoTech [62] nebo Trioptics [63]. Difraktivní optické elementy pro měření asférických ploch nabízí např. firmy Diffraction International [55] a Silios Technology [54]. Na trhu se objevují i jiné způsoby měření asférických ploch, např. tzv. Aspheric Stitching Interferometer (ASI™) firmy QED Technologies (USA) [64].
42
Interferometrie
4.1.2 Střihová interferometrie Další z metod, jak dospět k interferenci, je založena na jiném principu než výše představené interferometry Twyman-Greenova a Fizeaova typu. Jedná se o způsob, kdy necháváme interferovat vlnu s jednou nebo více vlnami, které jsou však prostorově modifikovanými kopiemi vlny původní [5 – 9]. Prakticky je vznik referenční vlny realizován pomocí různých optických prvků, např. planparalelních desek, difrakčních prvků, polarizačních prvků apod., kdy dochází k příčnému posunu nebo pootočení vlny původní. Na obr. 10 je zobrazeno principielní schéma tzv. Murtyho interferometru, který je pro střihovou interferometrii používán zejména pro svou jednoduchost. Jako dělicí člen je zde totiž použita planparalelní deska.
registrovaný obraz
+y
+x dr = (dx, 0)
pozorovací rovina (CCD senzor)
planparalelní deska zdroj záření (laser)
objektiv kolimátor Obr. 10 Murtyho střihový interferometr [6]
Přítomnost planparalelní desky v Murtyho interferometru způsobí posun vlnoplochy o dr = (dx, dy). V překrytové oblasti vzniká interferenční obrazec, který následně vyhodnocujeme. Pro optický dráhový rozdíl bude platit [6 – 9]
δρ M = W ( x + dx, y + dy ) − W ( x, y ) .
(124)
Budeme-li předpokládat, že posun dr je dostatečně malý, abychom mohli užít Taylorova rozvoje [39, 40] se zanedbáním členů vyšších než lineárních, čili že gradient fáze ∇W ( x, y ) bude na intervalu dr téměř konstantní, potom můžeme pro dráhový rozdíl δρM psát
43
Interferometrie δρ M = W ( x + dx, y + dy ) − W ( x, y ) =
∂W ( x, y ) ∂W ( x, y ) dx + dy = ∇W ( x, y )dr . ∂x ∂y
(125)
Z výše popsaného můžeme vyslovit závěr, že střihová interferometrie nevyhodnocuje fázový rozdíl a tím deformaci vlnoplochy W přímo, ale výstupem je její gradient (sklon). Vhodnou metodou (numerickou integrací) poté můžeme gradient přepočítat na samotnou deformaci vlnoplochy.
4.2
Záznam a vyhodnocení interferogramu
Průběh laboratorního měření a následného zpracování, ať se jedná o dvousvazkovou nebo střihovou interferometrii, můžeme dělit do následujících kroků:
•
záznam interferenčního pole,
•
předzpracování zaznamenaných dat,
•
vyhodnocení dat,
•
konečné zpracování, tzv. postprocessing5.
V následující části jednotlivé body podrobněji rozebereme.
4.2.1 Záznam interferenčního pole Záznam dat, interferogramu, je prvním bodem analýzy. V současnosti se v naprosté většině používá záznamu v digitální formě, nejčastěji s využitím CCD senzorů poměrně vysokého rozlišení, které jsou umístěny v rovině pozorování (rovině detekce optického záření). Výhodou digitální analýzy interferogramů je, že může být použito různých technik zpracování obrazu [47, 48] k potlačení nepříznivých vlivů, např. šumu, nerovnoměrného rozložení intenzity pozadí a podobně. Intenzita zaznamenávaného interferenčního pole I(x, y) je spojitou funkcí souřadnic. Abychom mohli
data
zpracovávat
digitální
formou,
dochází
k plošné
diskretizaci
spojitého
zaznamenávaného signálu – tzv. digitalizaci [47, 48]. Digitální obraz je ve výsledku diskrétní funkcí dvou proměnných, čili hodnoty jsou zaznamenány pro jednotlivé diskrétní body dané polohou (xi, yj), kde i = 0, 1, ... , N – 1 a j = 0, 1, ... , M – 1, M a N představují počet obrazových
Slovo „postprocessing“, jak je na první pohled patrné, není českého původu, ovšem běžně se v technické praxi používá, a tak se mu zde nebudeme vyhýbat.
5
44
Interferometrie bodů (rozměr senzoru) v jednotlivých směrech x a y. Pro obrazové body v rovině senzoru se používá název pixely [47, 48]. Při digitalizaci dochází k průměrování dopadající intenzity I(x, y) přes aktivní oblast pixelu a nezáporná finitní funkce I(x, y), pro kterou platí 0 < I ( x, y ) < ∞ ,
(126)
I d ≤ I ( xi , y j ) ≤ I h .
(127)
je převáděna do rozmezí
Hodnoty dolní a horní hranice, Id a Ih, ohraničují interval, který je nazýván škálou stupňů šedi [47, 48]. Běžně se následně interval transformuje na rozmezí 0 až L – 1, kde L je počet stupňů šedi. Ve výpočetních systémech (např. MATLAB [65]) jsou poté takovéto diskrétní funkce ukládány ve formě matic, se kterými se dále pracuje. Při digitalizaci je nutné správně zvolit množství diskrétních bodů, kterými budeme spojitý signál nahrazovat. Musíme tedy vhodně volit tzv. vzorkovací frekvenci obrazu. Tuto frekvenci určujeme pomocí tzv. Nyquist-Shannonova kritéria6, které říká [66]: „Jestliže funkce f(t) neobsahuje frekvence vyšší než W Hz, potom je kompletně určena svými hodnotami v sérii oddělených bodů, které jsou od sebe vzdáleny 1/(2W) sekund.“ Nyquist toto kritérium publikoval roku 1928 [67], Shannon ho poté upravil do výše představené přesnější podoby roku 1949 [66]. Jinými slovy, pro správné určení (záznam) spojitého signálu musíme zvolit takové vzorkování, jehož frekvence musí být minimálně 2krát vyšší, než je nejvyšší frekvence spojité funkce, která nese danou informaci. To pro případy interferometrie znamená minimálně dva obrazové body na jeden interferenční proužek. Vzhledem k dalším nepříznivým faktorům, jako je například šum apod., je ale v praxi vhodnější zvolit vzorkovací frekvenci alespoň 4krát větší, než je maximální frekvence obrazu. Kromě samotné diskretizace obrazu je další součástí digitalizace tzv. kvantování obrazové funkce [47, 48]. Jedná se o přiřazení diskrétní celočíselné hodnoty každému obrazovému bodu v závislosti na amplitudě vstupního signálu. Počet celočíselných kvantovacích úrovní musí být dostatečně velký, aby dobře reprezentoval spojitý průběh zaznamenávané intenzity a nedocházelo
6
Toto kritérium bylo nezávisle publikováno i jinými autory, proto v některých literaturách můžeme nalézt jiné označení. My se ovšem budeme držet nejpoužívanějšího – Nyquist-Shannonovo nebo jen Nyquistovo kritérium.
45
Interferometrie ke zkreslení a vniku falešných hran. Pro potřeby interferometrie postačuje použití úrovní šedi. Signál je poté kvantován do 2N hladin. V praxi se nejčastěji používá 8-bitové kvantování obrazové informace, čili přiřazení 28 = 256 hladin. Hodnoty spojitého obrazu jsou tak rovnoměrně rozloženy do intervalu od 0 do 255. Pro speciální aplikace může být použito i kvantování vyšší, 10-bitové nebo 12-bitové.
+I 2b
255 I = a + bcosΔϕ
a
0
π/2
0
π
(3/2)π
+Δϕ
Obr. 11 Nyquist-Shannonovo kritérium pro interferometrii
Z obr. 11 můžeme snadno odvodit vliv Nyquist-Shannonova kritéria pro fázový rozdíl a gradient fáze vyšetřovaného interferogramu. Uvažujme pro jednoduchost pouze jednodimenzionální případ. Čárkovaná křivka symbolizuje spojitou funkci intenzity interferenčního pole I = a + bcosΔϕ. Jak bylo řečeno, na jednu periodu (světlý a tmavý proužek) jsou třeba minimálně 2 pixely. Na obr. 11 je naznačena jejich poloha modrými úsečkami s číslem, které udává jejich přiřazenou kvantovanou hodnotu. Je patrné, že pro změnu fázového rozdílu mezi dvěma pixely bude platit
∆φ ≤ π .
(128)
V lineárním přiblížení můžeme změnu fázového rozdílu zapsat jako ∆φ =
∂φ ∆x , ∂x
(129)
kde x značí jeden ze směrů, který uvažujeme. Potom tedy můžeme psát důsledek NyquistShannonova kritéria na fázový rozdíl a gradient vyšetřované fáze interferenčního obrazce, a to ∆φ =
∂φ ∆x ≤ π, ∂x
∂φ π ≤ . ∂x ∆x
(130)
Pro určení minimálního počtu vzorků potřebných k záznamu interferogramu postupujeme následovně. Digitalizujme-li obraz v rozsahu xmin a xmax s počtem N vzorků, potom snadno
46
Interferometrie dostáváme Δx = (xmax – xmin)/N. Dosazení do výše uvedených vztahů vede k podmínce pro minimální počet vzorků ∂φ x max − x min ≤π , ∂x N ∂φ x max − x min N> . ∂x π
(131)
4.2.2 Předzpracování zaznamenaných dat Zaznamenaná digitalizovaná data jsou zatížena vlivem spousty vnějších i vnitřních faktorů. Jako příklad jmenujme vlastnosti optické soustavy, kterou svazek prochází. Může docházet k nerovnoměrným odrazům nebo propustné materiály nepropouští záření stejnoměrně. Přítomností nečistot nebo i jiných částic může docházet k parazitním difrakčním nebo interferenčním jevům. V neposlední řadě samotný detektor nemusí být rovnoměrně citlivý, při stejné intenzitě nemusí registrovat pro jednotlivé pixely stejné hodnoty. Všechny a ještě další komplikace, jako např. mechanické vlastnosti laboratorního prostředí (proudění vzduchu měnící index lomu prostředí, mechanické vibrace způsobující nestabilitu jednotlivých komponent interferometru, apod.) vedou k vzniku tzv. šumu. Důležitou myšlenkou je [68], že fázový rozdíl interferujících vln, který se snažíme vyhodnotit, není samotný zaznamenávaný signál. Je to informace, která je tímto signálem nesena. Proto je důležité provést úpravy k potlačení šumu v samotném signálu, tedy v zaznamenaném interferogramu, ne až ve vyhodnoceném fázovém rozdílu. Cílem předzpracování dat je tedy metodami zpracování obrazu co nejvíce potlačit výše nepříznivé rušivé faktory v samotném interferogramu. V obraze můžeme rozlišit šum na vysokofrekvenční, který se projevuje „zrnitostí“ větší či menší hustoty, a nízkofrekvenční, v jehož důsledku dochází k pomalým změnám v rozložení zaznamenané intenzity s ohledem k frekvencím originálního obrazu (frekvenci interferenčních proužků). Prakticky můžeme rozdělit proces předzpracování do čtyř kroků: •
potlačení vysokofrekvenčního a nízkofrekvenčního náhodného šumu,
•
určení a eliminace nerovnoměrné intenzity pozadí interferogramu,
•
úprava kontrastu,
•
výběr relevantní vyhodnocované oblasti – maskování. 47
Interferometrie
4.2.2.1 Model šumu v měřeném interferogramu Abychom mohli nějakým způsobem modelovat danou situaci, vyjděme ze vztahu pro intenzitu interferenčního pole kvazimonochromatického záření, čili záření nejvíce používaného v praxi. Tato intenzita je dána vztahem, jak bylo odvozeno dříve, I (r ) = I1 (r ) + I 2 (r ) + 2 I1 (r ) I 2 (r ) γ 12 (τ ) cos ∆φ (r ) ,
(132)
kde r značí polohu bodu, ve kterém intenzitu analyzujeme, |γ12(τ)| je modul komplexního stupně vzájemné koherence závislý na časovém zpoždění τ interferujících vln, Δϕ(r) je fázový rozdíl. Pro kontrast poté platí
K (r ) =
2 I1 (r ) I 2 (r ) γ 12 (τ ) . I1 (r ) + I 2 (r )
(133)
Zaveďme nyní pro jednoduchost zápisu substituci a (r ) = I1 (r ) + I 2 (r ) , b(r ) = 2 I1 (r ) I 2 (r ) γ 12 (τ ) .
(134)
Vztah (132) poté můžeme přepsat jako I (r ) = a (r ) + b(r ) cos ∆φ (r )
(135)
a pro kontrast (133) dostáváme
K (r ) =
b(r ) . a (r )
(136)
Původní vztah pro intenzitu (132) lze pomocí (135) a (136) zapsat jako
I (r ) = a(r )[1 + K (r ) cos ∆φ (r )] .
(137)
Vliv šumu nyní můžeme snadno modelovat zavedením tzv. aditivního šumu Na(r) a šumu multiplikativního Nm(r). Budeme-li dále značit a(r) jako tzv. nerovnoměrné rozložení intenzity pozadí, potom můžeme pro náš model psát I (r ) = a (r )[1 + K (r ) cos ∆φ (r )]N m (r ) + N a (r ) .
(138)
4.2.2.2 Potlačení vysokofrekvenčního a nízkofrekvenčního náhodného šumu Pro zpravidla vysokofrekvenční aditivní šum můžeme uvažovat normální rozdělení N(μ,σ2) [39, 40, 44] se střední hodnotou μ = 0 a variancí σ2 = konst., jelikož vznik tohoto šumu
48
Interferometrie připisujeme součtu mnoha náhodných nezávislých procesů, zejména elektronickému šumu fotodetektoru. Protože se jedná o proměnný proces, tedy při každém záznamu dostáváme jiná data (ale se stejnou střední hodnotu μ) potom lze jeho vliv zmenšit průměrem z více zaznamenaných interferogramů. Dostaneme tedy sadu n snímků Ii(r), kde i = 1, ... , n, se směrodatnou odchylkou σi = σ, kde dále předpokládáme směrodatnou odchylku nezašuměného obrazu nulovou, tudíž směrodatná odchylka obrazu se šumem bude totožná se směrodatnou odchylkou šumu samotného. Výsledný obraz získaný průměrem vyjadřuje poté vztah n
I (r ) =
∑ I (r ) i =1
i
n
(139)
a jeho směrodatnou odchylka aplikací statistické teorie (zákona hromadění směrodatných odchylek pro nezávislá měření) [44] výraz n
σ = 2 I
∑σ i =1
n2
2 i
nσ 2 σ 2 = 2 = . n n
(140)
Jak je patrné, výsledná směrodatná odchylka obrazu získaného průměrem z n snímků klesne √(n)krát. Na obr. 12 je simulován vliv průměrování obrazu. Na části a) je původní obraz, na části b) obraz se simulovaným šumem a jeho výběrová směrodatná odchylka [44], na části c) poté obraz po průměrování z 5 simulací zavedení šumu.
Obr. 12 Zobrazení a) C. F. Gausse (1777-1855) bez šumu, b) se zavedeným aditivním šumem se střední hodnotu μ = 0 a směrodatnou odchylkou σ = 15, c) průměrovaný obraz z 5 simulací zavedení šumu
49
Interferometrie Ze simulace je patrné, že výběrová směrodatná odchylka při jednom sejmutí obrazu odpovídá směrodatné odchylce šumu. Pro 5 simulací by teoretická hodnota směrodatné odchylky měla být s využitím vztahu (140) σ = 6.71, což je hodnota totožná s vypočtenou výběrovou směrodatnou odchylkou. V některých případech nemáme k dispozici možnost opakovaného snímání interferenčního obrazu a nemůžeme ideálně potlačit aditivní šum. Poté nastupují techniky zpracování obrazu. Uvedeme zde základní výčet a myšlenky, které jsou k potlačení šumu (filtrování) používány. Dá-li se předpokládat princip superpozice, čili předpokládáme-li aditivní šum, můžeme provádět tzv. filtraci digitálního obrazu [47, 48]. Zde používáme buď filtry v tzv. prostorové oblasti, nebo filtry v oblasti frekvenční. Pro náš případ se budeme zabývat pouze filtrováním obrazu daného ve stupních šedi. Takový obraz totiž bude reprezentovat zaznamenané rozložení intenzity interferenčního pole. Filtrace v prostorové oblasti znamená, že upravujeme obraz jako takový. Filtrace v oblasti frekvenční je prováděna tím způsobem, že snímek je nejdříve do frekvenční oblasti převeden vhodnou transformací (Fourierova transformace), tam je provedena filtrace a následně je prováděna transformace inverzní. Popsané schéma zobrazuje obr. 13.
PŮVODNÍ OBRAZ
TRANSFORMACE
FREKVENČNÍ FILTRACE
PROSTOROVÁ FILTRACE
INVERZNÍ TRANSFORMACE FILTROVANÝ OBRAZ Obr. 13 Schéma frekvenční a prostorové filtrace
50
Interferometrie Zabývejme se nejprve prostorovou filtrací. Pro jednoduchost zápisu můžeme použít k popisu této filtrace tzv. operátor diskrétní konvoluce [47, 48]. Pro náš případ obrazu I(x, y) můžeme psát I ( x, y ) ∗ h ( x, y ) =
a
b
∑ ∑ I ( x − s , y − t ) h( s, t ) ,
(141)
s =− a t =− b
kde h(x, y) je filtr7 o rozměrech (2a + 1) × (2b + 1) , a a b jsou sudá čísla. Filtr nám tedy představuje matici o lichém počtu sloupců a řádků. Prvky této matice jsou celá čísla v závislosti na úloze, kterou filtr má. Zápis konvoluce říká, že algoritmus prochází vstupní data pixel po pixelu. Na pozici (x, y) nastaví střed matice filtru (pozici s = 0 a t = 0), násobí pixel po pixelu hodnoty filtru s odpovídajícími hodnotami vstupních dat a výsledek uloží na polohu (x, y). Existuje velké množství rozličných typů filtrů v závislosti na úloze. Jmenujme například lineární filtr typu „dolní propust“ hD [47, 48], který odstraňuje složky vysokých frekvencí (vysokofrekvenční šum). Maticově ho můžeme pro rozměr 3× 3 zapsat jako
1 c hD = K c c 2 1 c −1
1 c , 1
(142)
kde K je konstanta, která je volena v souladu se zachováním dynamického rozsahu dat, a c je celé kladné číslo. Nevýhodou použití tohoto filtru je, že při odstranění vysokých frekvencí dochází k rozostření obrazu. Opakem je filtr typu „horní propust“ hH, který potlačuje frekvence nízké,
− 1 − 1 − 1 hD = 8 − 1 8 − 1 . − 1 − 1 − 1 −1
(143)
Kromě nízkých frekvencí obrazu potlačuje filtr horní propust i střední intenzitu pozadí. Kromě lineárních filtrů lze použít pro filtraci digitálního obrazu v prostorové oblasti i filtrů nelineárních, které na okolí bodu (x, y) aplikují jiné operace. Příkladem může být mediánový filtr [47, 48].
7
Filtr může být také zaměňováno za „okno“, „kernel“, „konvoluční matice“ apod.
51
Interferometrie Potlačení multiplikativního šumu je velmi obtížné. Zpravidla je způsoben tzv. koherentním (spekl) šumem, který se dá modelovat negativním exponenciálním pravděpodobnostním rozdělením s vysokým kontrastem [47, 48]. V praxi ale neexistuje žádný obecný postup, který by tento šum ideálně eliminoval případně částečně potlačil. Existují však různé techniky, které se o to snaží. Pro částečné potlačení multiplikativního šumu lze použít i výše popsanou filtraci digitálního obrazu. Předpokládejme, že jsme dostatečně odfiltrovali šum aditivní, čili rovnici (138) můžeme psát jako I (r ) ≈ a (r )[1 + K (r ) cos ∆φ (r )]N m (r ) .
(144)
Logaritmováním poté dostáváme ln I (r ) ≈ ln{a (r )[1 + K (r ) cos ∆φ (r )]} + ln N m (r ) .
(145)
Na rovnici (145) následně aplikujeme výše popsané metody prostorové filtrace. Jelikož byl multiplikativní šum převeden na formu aditivní, dojde k jeho odfiltrování. Zpětným odlogaritmováním poté získáme původní obraz s potlačeným multiplikativním šumem. Další vhodnou metodou úpravy obrazu je tzv. adaptivní filtrace [47, 48]. Je založena na určení průběhu interferenčních proužků (např. pomocí gradientu obrazové funkce). Filtrace je poté prováděna pouze v jejich směrech. Díky takovému přizpůsobení k obrazu nedochází k nadměrnému rozostření interferenčních proužků, což je v některých případech důležité (např. interferogramy s vysokým počtem proužků). Filtrace ve frekvenční oblasti využívá vlastností Fourierovy transformace, která nám převádí obraz z prostorové do frekvenční oblasti, čili říká nám, jaké frekvence se v obraze nacházejí. Volbou vhodného filtru poté můžeme jednotlivé frekvence potlačovat a inverzní transformací získat obraz bez nich. Samotné použití filtrů (masek) s ostrými hranami vede díky vlastnostem diskrétní Fourierovy transformace, která je při vyhodnocování používána, ke vzniku nepříznivých struktur v rekonstruovaném obraze. Využívá se proto dvojrozměrných oken používaných v digitální analýze signálu s plynulým přechodem, založených například na Hammingově okně [47, 48, 69].
4.2.2.3 Eliminace nerovnoměrné intenzity pozadí interferogramu Pro vyrovnání nerovnoměrné intenzity pozadí lze použít výše uvedené způsoby filtrace. Máme konkrétně na mysli filtr typu „dolní propust“, který eliminuje vysoké frekvence v obraze. Jestliže 52
Interferometrie zvolíme okno konvoluční matice (kernel) dostatečně velké oproti minimální prostorové frekvenci interferenčních proužků, dochází poté k jejich vyhlazování a výsledný obraz představuje odhad střední hodnoty rozložení intenzity a (r ) . Po odečtení od rovnice (135) následně dostaneme interferogram s vyváženou hodnotou intenzity pozadí. Na obr. 14 je jednoduchý případ pro reálný snímek zobrazen.
Obr. 14 a) Reálný interferogram a b) jeho podoba po odstranění nerovnoměrného rozložení intenzity pozadí s použitím průměrujícího filtru [47, 48]
Jinou možností je využití aproximace průběhu intenzitního profilu v obraze pomocí dvojrozměrných polynomů. O dvourozměrné aproximaci bude podrobně pojednáno v části 4.2.4, zde na ní tedy odkážeme. Principem je určit průběh nerovnoměrné intenzity pozadí jako funkci dvou proměnných. Poté již snadno můžeme její vliv potlačit.
4.2.2.4 Úprava kontrastu interferogramu Jak je z části b) na obr. 14 patrné, došlo ke ztrátě informace co se týče ideálního rozložení hodnot interferogramu, není tedy ideálně využito n-bitové škály obrazu. Pro vyhodnocení interferogramů je vhodné a v některých případech nutné, aby byly jednotlivé interferenční proužky od sebe dobře rozeznatelné. Přichází tedy na řadu tzv. úpravou kontrastu. Jednou z možností úpravy kontrastu je vyrovnání histogramu [47, 48]. Histogram digitálního obrazu při n-bitovém kódování, jenž obsahuje hodnoty jasu od 0 do (2n – 1), je diskrétní funkce h(rk) = nk, kde rk je hodnota k-té intenzity, tzv. intenzitní úroveň, a nk je počet pixelů v obraze, kterým je tato intenzita přiřazena. Často se používá tzv. normalizovaný histogram p(rk). Ten vzniká snadno dělením funkce h(rk) celkovým počtem pixelů. Pro obraz o rozměrech M × N tedy dostáváme p(rk ) =
h(rk ) n = k , k = 0, 1, K , 2 n − 1 . NM NM
(146)
53
Interferometrie Normalizovaný histogram nám jinými slovy představuje pravděpodobnost výskytu intenzitní úrovně rk v obraze. Součet členů normalizovaného histogramu je roven 1. Pro pochopení principu vyrovnání histogramu uvažujme nyní spojitý případ intenzitních hodnot. Cílem je získat takové hodnoty výsledného obrazu, které budou rovnoměrně rozloženy na celém zaznamenávaném intervalu. Předpokládejme spojitou funkci r, pro kterou platí 0 ≤ r ≤ L – 1, kde L = 2n. Definujme nyní transformační funkci T(r), která přiřazuje pixelu vstupního obrazu o intenzitní hodnotě r intenzitní hodnotu s v obraze výstupním. Můžeme tedy psát s = T (r ) , r = 0, 1, K , L − 1 .
(147)
Předpokládejme, že funkce T(r) bude striktně monotónní [39, 40]. To musí být dodrženo, aby stejným hodnotám r byla přiřazena jedna hodnota s. Podíváme-li se na intenzitní úrovně jako na náhodné proměnné s hustotou pravděpodobnosti pr(r) na intervalu 0 ≤ r ≤ L – 1, potom intenzitní hodnoty s budou dány transformací [47, 48] r
s = T (r ) = ( L − 1) ∫ pr (t ) d t .
(148)
0
Jak je patrné, pravá strana rovnice (148) je kumulativní hustota pravděpodobnosti (distribuční funkce) náhodné proměnné r. Pro získání informace o hustotě pravděpodobnosti ps(s) využijeme základní znalosti teorie pravděpodobnosti, která říká, že známe-li hustotu pravděpodobnosti pr(r) a vztah pro transformaci T(r), potom platí [47, 48] ps ( s ) = pr (r )
dr . ds
(149)
Diferencováním vztahu (147) a dosazením z (148) dostáváme ds dT ( r ) = dr dr d r p r (t ) d t ∫ dr 0 = ( L − 1) p r (r ) .
= ( L − 1)
(150)
Dosazením do vztahu (149) dostáváme
ps ( s ) = pr (r )
1 1 = , s = 0, 1, K , L − 1 , ( L − 1) pr (r ) ( L − 1)
(151)
54
Interferometrie kde jsme mohli bez komplikací odstranit absolutní hodnotu, protože hodnoty hustoty pravděpodobnosti pr(r) jsou vždy kladné. Je patrné, že hustota pravděpodobnosti ps(s) bude mít s užitím transformace (148) rovnoměrné pravděpodobnostní rozdělení [39, 40, 44] bez ohledu na tvar hustoty pravděpodobnosti pr(r), což je v souladu s naším zadáním, tedy že výsledný obraz musí mít rovnoměrně rozložené intenzitní hodnoty na celém svém intervalu, jinými slovy pravděpodobnost výskytu jednotlivých intenzitních hodnot musí být stejná. Pro diskrétní případ přejde transformace (148) na tvar k
sk = T (rk ) = ( L − 1)∑ pr (rj ) = j =0
( L − 1) k ∑ n j , k = 0, 1, K , L − 1 . MN j =0
(152)
Na obr. 15 je ukázán vliv vyrovnání histogramu na interferogram z obr. 14 b). Na obr. 16 poté normalizované histogramy před a po provedení vyrovnání.
Obr. 15 a) Původní interferogram a b) jeho podoba po vyrovnání histogramu
Obr. 16 Normalizované histogramy a) původního obrazu a b) obrazu po provedení vyrovnání
55
Interferometrie
4.2.2.5 Výběr relevantní vyhodnocované oblasti – maskování Tzv. maskování dat je nutné v tom případě, jestliže při vyhodnocení apriori známe oblasti obrazu, které nejsou v místech detekovaného interferenčního pole. Tato situace může snadno nastat, pokud užíváme clonu kruhového nebo eliptického tvaru a CCD senzor je obdélníkový. Jiným případem může být nedostatečný kontrast některých oblastí, který by měl fatální dopad na určité typy vyhodnocovacích metod (např. metody automatické detekce středu interferenčních proužků, viz dále). Maskou bude matice stejného rozměru jako vyhodnocovaný obraz, jejíž prvky budou 0 pro místa, která odpovídají pixelu, kde není detekováno interferenční pole, a 1, kde pole detekováno je. Určení takové oblasti lze provést manuálně nebo například užitím hodnot kontrastu interferogramu K(r). V takovém případě je stanovena prahová hodnota T, která odděluje data přijatelná a nepřijatelná k vyhodnocení. Poté tedy bude pro prvky masky M(r) platit 0 , K (r ) < T , M (r ) = 1 , K (r ) ≥ T ,
(153)
kde r značí polohu pixelu.
4.2.3 Vyhodnocení interferogramu Předzpracované obrazy metodami popsanými v předchozí části zlepšují šanci na kvalitnější a spolehlivější rekonstrukci interferometrických dat. V praxi se pro vyhodnocení a získání údajů o fázovém rozdílu interferujících vln využívají různé principy. Dle jednotlivých způsobů se pak sestavují metodiky pro záznam a firmy dle nich konstruují jednotlivé interferometry. Z klasických metod vyhodnocení interferogramu jmenujme tyto: •
detekce interferenčních proužků a jejich polynomiální aproximace s následným výpočtem fázového rozdílu,
•
interferometrie s posuvem fáze,
•
interferometrie se zaváděním nosné frekvence (FTM).
V další části textu bude představen popis obecných principů jednotlivých metod. Pro podrobné pojednání o problematice odkážeme čtenáře na dostupnou odbornou literaturu.
56
Interferometrie
4.2.3.1 Vztah fázového a dráhového rozdílu, základní podmínka metodiky vyhodnocení Než přistoupíme k charakteristice jednotlivých metod interferometrie, ukážeme základní vztah pro převod fázového rozdílu na rozdíl dráhový. Jak bylo odvozeno v předchozí části práce, vztah pro zaznamenávanou intenzitu interferenčního pole v bodě r vzniklou dvěma rovinnými monochromatickými vlněními je I (r ) = a (r ) + b(r ) cos ∆φ (r ) ,
(154)
kde a(r) a b(r) znační střední intenzitu pozadí a amplitudu intenzity, Δϕ(r) je fázový rozdíl daný vztahem
∆φ (r ) = (k 1 − k 2 ) r + ∆δ12 ,
(155)
kde k1 a k2 jsou vlnové vektory, r je polohový vektor Δδ12 je počáteční fázový rozdíl. Dosazením za vlnové vektory k = kn = (2π/λ)n, kde k je vlnové číslo, λ vlnová délka použitého záření a n určuje směr přicházejícího paprsku (normálový vektor rovinné vlnoplochy), dostáváme
∆φ (r ) = (k1n 1 − k 2 n 2 ) r + ∆δ 12 2π 2π = n 1 − n 2 r + ∆δ 12 . λ λ 2 1
(156)
Jelikož v interferometrech používáme jeden zdroj záření, můžeme si dovolit položit počáteční fázový rozdíl Δδ12 roven nule a vlnová čísla ztotožnit. Tato úprava vede na
∆φ (r ) = k ρ1 − ρ 2 =
2π
λ
W (r ) ,
(157)
kde ρ1 a ρ2 jsou průměty polohového vektoru r do směrů daných normálovými vektory n1 a n2, W(r) značí námi hledaný dráhový rozdíl. Jestliže budou směry šíření interferujících vln stejné a vektory ρ1 a ρ2 rovnoběžné, potom dráhový rozdíl přímo udává odchylku testované od referenční plochy. Snadno jí tedy získáme pomocí vztahu W (r ) =
λ 2π
∆φ (r ) .
(158)
Podíváme-li se na vztah (154), poté ze znalosti goniometrických funkcí můžeme pro kosinus fázového rozdílu psát
cos ∆φ (r ) = cos[k∆φ (r ) + 2πl ] ,
(159)
57
Interferometrie kde k = {–1, 1} určuje znaménko fázového rozdílu a l je libovolné celé číslo. Použití jednoho interferogramu a pouhá inverze vztahu (154) by vnesla do vyhodnocení nejednoznačnost. Metodiky měření tedy musí být takové, aby určení fázového rozdílu bylo jednoznačné případně s výjimkou znaménka. To umožňuje například znalost apriorní informace o průběhu fáze, užití více fázově posunutých interferogramů nebo zavedení lineární prostorové frekvence. O těchto metodách bude v následujícím textu diskutováno.
4.2.3.2 Detekce a polynomiální aproximace interferenčních proužků s následným vyhodnocením fázového rozdílu Jednou z metod vyhodnocování fázového rozdílu interferujících vln je analytické vyjádření polohy interferenčních proužků a následné vyhodnocení [6 – 9]. Pro tuto metodu postačuje jeden interferogram, ovšem je nutná další apriorní znalost o průběhu fázového rozdílu. Dále musíme předpokládat, že lokální extrémy zaznamenané intenzity odpovídají extrémům vztahu (154), jinými slovy data předpokládáme s odfiltrovaným šumem, vyrovnanou nerovnoměrnou intenzitou pozadí a upraveným kontrastem. Přesnost této metody se pohybuje okolo λ/20 [6], kde λ je vlnová délka použitého záření. Princip metody je snadný a patrný z obr. 17. Budou-li splněny předchozí předpoklady a skutečné rozložení intenzity bude odpovídat vztahu (154), potom úsek mezi maximem (nejsvětlejším místem) a minimem (nejtmavším místem) jednoho interferenčního proužku odpovídá fázovému rozdílu Δϕ = π. Budeme-li schopni určit polohu maxim interferenčních proužků, potom mezi sousedníma dvěma bude fázový rozdíl roven 2π a s apriorní znalostí sklonu (známe přibližný tvar vyšetřované plochy, jsme tedy schopni určit hodnotu konstanty k ve vztahu (159)) můžeme daný fázový rozdíl rekonstruovat.
+I 2b
Imax I = a + bcosΔϕ
a Imin 0
+Δϕ Δϕ = π
Obr. 17 Princip metody detekce interferenčních proužků a následného vyhodnocení fázového rozdílu
58
Interferometrie Jednotlivé kroky postupu vyhodnocení tedy budou: •
určení středu interferenčních proužků a jejich polynomiální aproximace,
•
číslování proužků přiřazením interferenčního řádu,
•
rekonstrukce fázového rozdílu – jeho globální prostorová aproximace (interpolace dat mezi interferenčními proužky).
Detekce středů interferenčních proužků je možná různými způsoby. Principielně jde o určení polohy lokálních maxim obrazu, čili lze použít různé techniky založené na aplikaci funkcionální analýzy [39, 40], případně jiné vhodné metody zpracování digitálního obrazu (např. segmentace obrazu, použití Yatagaiho matic apod. [6, 47]). Číslováním interferenčních proužků rozumíme přiřazením relativní výšky fázového rozdílu jednotlivým proužkům. Schematický postup by mohl být zapsán takto: •
výběr výchozího proužku a přiřazení počáteční hodnoty,
•
volba vzestupného nebo sestupného číslování vzhledem k výchozímu proužku a danému směru,
•
číslování jednotlivých proužků.
Předpokládáme-li, že vyhodnocovaný fázový rozdíl bude „dostatečně“ spojitá a hladká plocha, což můžeme u testování optických komponent předpokládat, můžeme vyslovit následující pravidla: •
sousední interferenční proužky se mohou lišit o –1, 0 nebo 1 a ne jinak,
•
nekončí-li proužek na hranici vyhodnocované oblasti interferogramu, jde o uzavřený proužek,
•
proužky stejného čísla se nesmějí vzájemně křížit.
Z výše popsaného základního principu je patrné, že je-li interferenčnímu proužku přiřazeno číslo interferenčního řádu k, potom jeho rozdíl oproti výchozímu proužku bude k2π. Výslednou rekonstrukci fázového rozdílu (interpolaci dat mezi jednotlivými interferenčními proužky) provádíme některým z postupů globální polynomické aproximace. Podrobněji o nich bude pojednáno v další části práce.
59
Interferometrie Na obr. 18 je demonstrováno číslování jednotlivých středů interferenčních proužků a následná rekonstrukce fázového rozdílu.
Obr. 18 a) Číslování středů interferenčních proužků a b) rekonstrukce fázového rozdílu
4.2.3.3 Interferometrie s posuvem fáze Interferometrie s posuvem fáze (PSI – z angl. Phase Shifting Interferometry) je metoda, která rekonstruuje fázový rozdíl s využitím více interferogramů lišících se od sebe o fázovou modulací (fázový posuv v čase) [6 – 10, 70 – 72]. Přesnost těchto metod vyhodnocení je poměrně vysoká, až λ/100 a více. Fázový posuv se nejčastěji u koherentních vlnových polí zavádí do referenčního svazku pomocí piezoelektrického posuvu, rotující planparalelní desky, posuvných hranolů, pohybující se difrakční mřížky nebo dalších komponent [6 – 9]. Velmi výhodné je použití referenčního posuvného zrcátka na piezoelektrickém posuvu, jehož pohyb se děje ve směru rovnoběžném s dopadajícím referenčním svazkem. Potom s užitím (158) dostáváme pro fázový posuv δ v závislosti na posuvu zrcátka d
δ=
2π
λ
d,
(160)
kde λ je vlnová délka použitého záření. Princip samotného vyhodnocení vychází opět ze základní rovnice I (r ) = a (r ) + b(r ) cos ∆φ (r ) ,
(161)
60
Interferometrie kde a(r) a b(r) znační střední intenzitu pozadí a amplitudu intenzity, Δϕ(r) je fázový rozdíl referenční a testované vlnoplochy. Ve vztahu se nám vyskytují tři neznámé – a(r), b(r) a Δϕ(r), obecně musíme tedy provést minimálně tři měření intenzity. Jednou z metod je zavedení známého fázového posuvu δi, čili vztah (161) pro i-tý měřený interferogram přejde na tvar I i (r ) = a (r ) + b(r ) cos[∆φ (r ) + δ i ] , i = 1, 2, K , N ,
N ≥3,
(162)
kde předpokládáme, že veličiny a(r) a b(r) jsou pro všechny interferogramy konstantní a není přítomna chyba v zavedení fázového posuvu δi. Pro i měření je poté možno výpočet provést metodou nejmenších čtverců [39, 40, 44], kde vztah (162) ještě upravíme na tvar I i (r ) = a (r ) + b(r ) cos ∆φ (r ) cos δ i − b(r ) sin ∆φ (r ) sin δ i = c 0 (r ) + c1 (r ) cos δ i + c 2 (r ) sin δ i ,
(163)
kde, jak je patrné, c0(r) = a(r), c1(r) = b(r)cosΔϕ(r) a c2(r) = –b(r)sinΔϕ(r). Vztah (163) je lineární a snadno ho můžeme přepsat maticově pro N měření jako y = Ax ,
(164)
kde
I1 (r ) 1 cos δ 1 I (r ) 1 cos δ 2 y= 2 , A= M M M I N (r ) 1 cos δ N
sin δ 1 c0 (r ) sin δ 2 , x = c1 (r ) . M c2 (r ) sin δ N
(165)
Řešení soustavy (164) je pro nadbytečný počet měření, čili kdy N > 3, dáno známým vztahem výpočtu ve smyslu metody nejmenších čtverců x = (A T A ) A T y . −1
(166)
Vyhodnocovaný fázový rozdíl poté dostaneme jako
c (r ) ∆φ (r ) = arctan − 2 . c1 (r )
(167)
Jak je z výše uvedeného patrné, pro určení fázového posunu v bodě r obrazu o rozměrech M × N je třeba řešit M∙N soustav rovnic.
61
Interferometrie V praxi bývá však obtížné měření (zavedení) přesné hodnoty fázového posuvu δi. Častěji je proto užíváno metody, kdy je zaváděn konstantní fázový posuv, jehož hodnotu ovšem nemusíme přesně znát, naopak bývá výsledkem výpočtu. Tím se ovšem do vztahu pro detekovanou intenzitu dostává čtvrtá neznámá a musíme tedy provádět minimálně čtyři měření. Výhodou této metody však je, že nejen že určíme fázový rozdíl interferujících vln bez přesné znalosti fázového posuvu jednotlivých obrazů, ale samotný fázový posuv může sloužit ke kontrole a kalibraci zařízení, jenž tuto modulaci zaznamenávaného signálu realizuje. Základním čtyřkrokovým (vyhodnocuje čtyři měření) algoritmem je tzv. Carrého algoritmus [72]. Využívá obecného vztahu, že měřenou intenzitu v každém bodě r interferogramu si můžeme zapsat jako [70 – 72] 2i − N − 1 I i (r ) = a (r ) + b(r ) cos ∆φ (r ) − α (r ) , i = 1, K , N , 2
(168)
kde fázový posuv mezi jednotlivými snímky je 2α, N = 4. Pro jednoduchost zápisu budeme v dalším vynechávat značení závislosti na poloze, čili např. I(r) = I. Vhodnou úpravou lze ukázat, že platí cos 2α =
(I
1
− I 4 ) − (I 2 − I 3 ) 2( I 2 − I 3 )
(169)
a pro fázový rozdíl Δϕ platí
∆φ = arctan
(I
1
− I 4 + I 2 − I 3 )(− I1 + I 4 + 3I 2 − 3I 3 ) . ( I 2 + I 3 − I1 − I 4 )
(170)
Mimo tohoto je známa celá řada dalších N-krokových PSI algoritmů různých autorů [6, 7], které danou problematiku vyhodnocení fázového rozdílu dvou interferujících vln řeší. Při měření je nutné vybrat takový algoritmus, který bude co nejméně citlivý na chybu měření intenzity interferujících vln, která se projeví jak ve výsledných hodnotách fázového rozdílu, tak i v hodnotách určeného fázového posuvu. Obecně čím větší počet kroků výpočtu, tím bude algoritmus méně citlivý na chyby. Na druhou stranu se tím zvyšuje náročnost na pořízení snímků a samotného vyhodnocení. V praxi jsou nejvíce rozšířeny algoritmy tříkrokové a pětikrokové. Jak ukazuje [70] nebo [71], jsou pětikrokové algoritmy vyhovující ve většině praktických případů jak v nárocích na přesnost vyhodnocení, tak i v rychlosti výpočtu.
62
Interferometrie Metoda interferometrie s posuvem fáze skýtá tedy několik výhod oproti jiným principům: je známo velké množství algoritmů různých vlastností, složky aditivního a multiplikativního šumu jsou automaticky potlačovány výpočetním postupem, s vhodnými prvky jde o velmi přesnou metodu. Nevýhodou je ovšem nutnost použití více interferogramů, čili odpadá možnost aplikace této metody u měření a analýzy rychle se vyvíjejících dynamických procesů.
4.2.3.4 Interferometrie se zaváděním nosné frekvence Tento způsob rekonstrukce fázového rozdílu bude podrobně popsán v další části práce. Zde se seznámíme zejména s možností zavádění tzv. lineární nosné frekvence do vyhodnocovaného interferogramu a se základním principem vyhodnocení. Podíváme-li se na metody posunu fáze (PSI) popsané v předchozí části, můžeme vyslovit mimo jiné tyto závěry [6]: •
pro metody PSI jsou třeba minimálně tři měřené interferogramy, ideální počet je pět,
•
abychom co nejvíce potlačili vlivy vibrací a dalších nežádoucích jevů u metod PSI, musíme měření tří nebo více snímků provádět v co nejkratším časovém intervalu,
•
znaménko u fázového rozdílu je jednoznačně určeno,
•
náklady na měřící zařízení jsou poměrně vysoké,
•
pokud je zajištěno stálé prostředí bez vibrací, fluktuací vzduchu apod., potom metody PSI dosahují velmi vysoké přesnosti, přibližně λ/100 a vyšší.
Oproti tomu o metodách zavádění nosné frekvence platí, o čemž se přesvědčíme v následujících kapitolách: •
k vyhodnocení je třeba pouze jeden interferogram,
•
vibrace celého systému nejsou závažným problémem, jelikož je využíván pouze jeden snímek,
•
z jednoho snímku nejsme schopni určit znaménko fázového rozdílu, je třeba znát apriori minimálně jednu další informaci,
•
k vyhodnocení jsou nutné sofistikované matematické aparáty,
•
v laboratorních podmínkách je přesnost odhadována na λ/20 – λ/50.
63
Interferometrie I když metody zavádění lineární nosné frekvence dosahují menší přesnosti, než metody PSI, jsou často využívány zejména pro svojí potřebu pouze jednoho snímku. Mají tedy potenciál u vyhodnocování dynamických procesů – aplikací v reálném čase. Předpokládejme, že rozdělení zaznamenané intenzity interferogramu bez zavedené nosné frekvence bude dáno známým vztahem I (r ) = a (r ) + b(r ) cos ∆φ (r ) ,
(171)
kde a(r) a b(r) znační střední intenzitu pozadí a amplitudu intenzity, Δϕ(r) je fázový rozdíl referenční
a testované vlnoplochy,
r = (x, y) je polohový vektor a předpokládáme
monochromatické lineárně polarizované interferující vlny. Referenční i testovanou vlnoplochu uvažujme pro názornost rovinnou. Princip zavedení lineární nosné frekvence naznačuje obr. 19.
+x
+x
ϕr(r)
ϕt(r)
θx ϕr(r)
ϕt(r)
+z
+z
Δϕ(r)
ϕr(r) + F(r) Obr. 19 Zavádění lineární nosné frekvence u rovinných vlnoploch [6]
Předpokládejme, že dvě rovinné vlnoplochy kolmé na směr šíření se šíří ve směru osy z. Fázi referenční vlnoplochy označujme ϕr(r) a vlnoplochy testované ϕt(r). Pro fázový rozdíl Δϕ(r) poté můžeme psát
∆φ (r ) = φt (r ) − φr (r ) .
(172)
Zaveďme nyní funkci F(r) = F(x,y) vztahem
F (r ) = k (x sin θ x + y sin θ y ) sin θ y sin θ x = 2 π x +y λ λ = 2 π(xf 0 x + yf 0 y ) = 2 πf0r ,
(173)
64
Interferometrie kde f0 = (f0x, f0y) = λ-1(sinθx, sinθy) je tzv. lineární prostorová nosná (přenosová) frekvence, θx a θy jsou sklony ve směrech os x a y, jak je naznačeno na obr. 19 pro případ θx, k = 2π/λ je vlnové číslo daného záření. Jak je vidět, funkce F(r) je rovnicí roviny nakloněné vzhledem k souřadným osám. Praktickou realizaci zavedení této funkce lze provést náklonem referenční plochy interferometru, která však nemusí být obecně rovinná. Tvar referenční vlnoplochy ϕr(r) může být různý. Zavedení náklonu (z angl. tilt) bude mít ovšem stejný vliv, funkce F(r) bude mít stejný tvar jako v (173). Pro intenzitu interferenčního pole (171) dostáváme přičtením funkce F(r) k fázovému rozdílu Δϕ(r) I (r ) = a (r ) + b(r ) cos[2 πf0r + ∆φ (r )] ,
(174)
což je výchozí rovnice pro vyhodnocení interferogramů metodou zavádění nosné frekvence. Jiný způsob zápisu, jenž bude využíván, je dán aplikací Eulerova vzorce
cos α + i sin α = exp(iα ) ,
(175)
kde i = √(–1). S jejím využitím lze snadno vztah (174) přepsat na tvar 1 1 I (r ) = a (r ) + b(r ) exp{i[2 πf 0 r + ∆φ (r )]} + b(r ) exp{− i[2 πf 0 r + ∆φ (r )]} 2 2 ∗ = a (r ) + c(r ) exp(i 2 πf 0 r ) + c (r ) exp(− i 2 πf 0 r ) ,
(176)
1 c(r ) = b(r ) exp[i∆φ (r )] , 2 1 c ∗ (r ) = b(r ) exp[− i∆φ (r )] 2
(177)
kde
jsou vzájemně komplexně sdružené funkce. Zavedení nosné frekvence se využívá při vyhodnocení interferogramů metodou Fourierovy transformace, o čemž pojednávají další části této práce. Jak bude ukázáno, aplikací Fourierovy transformace na vztah (176) dostaneme frekvenční Fourierovské spektrum, které bude právě díky zaváděné nosné frekvenci charakteristické třemi oddělenými spektrálními komponenty a můžeme ho zapsat jako FT {I (r )} = A(f ) + C (f − f 0 ) + C ∗ (f + f0 ) ,
(178)
65
Interferometrie kde jednotlivé komponenty A(f ) , C (f − f0 ) a C ∗ (f − f 0 ) příslušejí třem sčítancům na pravé straně rovnice (176). Ve frekvenční oblasti poté probíhá filtrace jednoho z bočních spekter, které nese informaci a rekonstruovaném fázovém rozdílu. Posunem tohoto spektra takovým způsobem, aby
f0 = (0,0), se zbavíme zaváděné nosné frekvence a inverzní Fourierovou transformací získáme c(r), z něhož poté snadno získáme rekonstruovaný fázový rozdíl interferujících vln jako
Im[c(r )] ∆φ (r ) = arctan . Re[c(r )]
(179)
Odvození výše uvedených vztahů bude provedeno v následující části práce. Data, rekonstruovaná fáze, získaná vyhodnocením interferogramu pomocí Fourierovy transformace, budou zpravidla v rozmezí hodnot od –π do +π, v závislosti na definici funkce „arctan“ ve výpočetním softwaru. Ideální vyhodnocený obraz (bez šumu a dalších vlivů) bude tedy obsahovat místa nespojitosti velikosti 2π. Odstraněním těchto míst nespojitosti se zabývá tzv. unwrapping. Tomuto matematickému problému byla a je v praxi věnována značná pozornost, protože na jeho vyřešení závisí přesnost celého výsledku. Pokud data neobsahují šum a další rušivé vlivy, je odstranění nespojitostí triviálním problémem. Ovšem pokud je šum přítomný, jako ve všech praktických případech, situace se komplikuje. Shrnutí procesu vyhodnocení interferogramu metodou Fourierovy transformace je zobrazeno na obr. 20.
INTERFEROGRAM SE ZAVEDENOU LINEÁRNÍ NOSNOU FREKVENCÍ FOURIEROVA TRANSFORMACE FREKVENČNÍ FILTRACE A CENTRACE BOČNÍHO SPEKTRA INVERZNÍ FOURIEROVA TRANSFORMACE ODSTRANĚNÍ NESPOJITOSTÍ FÁZOVÉHO ROZDÍLU
REKONSTRUOVANÝ FÁZOVÝ ROZDÍL INTERFERUJÍCÍCH VLN
Obr. 20 Schéma metody Fourierovy transformace
66
Interferometrie V praxi je obecně nazýván tento postup „metodou Fourierovy transformace“ (FTM – z angl. Fourier Transform Method) nebo „Fourierovská analýza interferenčních proužků“ (FFA – z angl. Fourier Fringe Analysis). Metodu poprvé představil Takeda et. al. roku 1982 [73]. Od té doby je v popředí zájmu zejména proto, že k vyhodnocení potřebuje pouze jeden zaznamenaný interferogram, a je tedy vhodným adeptem na analýzu dynamických procesů. Ovšem má i své nevýhody, především nižší přesnost vyhodnocení oproti metodám PSI, nejsme schopni určit z jednoho snímku znaménko rekonstruovaného fázového rozdílu a analýza je náročnější na výpočetní kapacity, i když se současným rychlým rozvojem výpočetní techniky přestává být tento argument proti jejímu použití relevantní. I přes tyto všechny nedostatky nachází metoda uplatnění. Roku 2004 představil Kemao obměněnou verzi metody Fourierovy transformace, jedná se o „metodu Fourierovy transformace pomocí filtračních oken“ (WFT – z angl. Windowed Fourier Transform) [74 – 76]. Výsledný obraz vyhodnocený pomocí FTM, který je definovaný na uzavřeném intervalu, bude ovlivněn díky vlastnostem Fourierovy transformace definované na intervalu od –∞ do ∞. WFT algoritmus se snaží tuto nepříjemnou vlastnost potlačit. Principem je užití modifikované transformace, která pomocí oken transformuje spektrum do omezených oblastí. Frekvenční spektrum WFT poté dává frekvenční informaci v každém pixelu obrazu, což FTM samotná neposkytuje, a je možné provádět operace vedoucí ke kvalitnější rekonstrukci fázového rozdílu. Výpočetní nároky jsou ovšem pro WFT velmi vysoké. Jinou z metod využívající Fourierovu transformaci představil Fabian et. al. [77, 78]. Jedná se o převedení problému na tzv. Pfaffovu parciální diferenciální rovnici . Zavedeme novou funkci h(r) definovanou jako h(r ) =
c(r ) = exp[i 2∆φ (r )]. c ∗ (r )
(180)
Pro gradient funkce h(r) derivováním (180) podle r dostáváme
∇h(r ) = i 2 exp[i 2∆φ (r )]∇[∆φ (r )] = i 2h(r )∇[∆φ (r )] .
(181)
Úpravou poté vyjádříme gradient fázového rozdílu jako ∇[∆φ (r )] =
1 ∇h(r ) . i 2h(r )
(182)
Užitím (177) a (180) dále dostáváme
67
Interferometrie ∇[∆φ (r )] =
1 ∇c(r ) ∇c ∗ (r ) − ∗ . i 2 c(r ) c (r )
(183)
Gradient tedy můžeme vyjádřit pomocí funkcí c(r) a c*(r) získávaných pomocí Fourierovy transformace ze zaznamenaného interferenčního obrazce. Další řešení je postaveno na následující úvaze. Gradient fázového rozdílu skalárně násoben s vektorem elementárních přírůstků ve směrech souřadnicových os x a y, čili s vektorem dr = (dx, dy), je jiným vyjádřením totálního diferenciálu, čili můžeme psát d∆φ (r ) = ∇[∆φ (r )]dr =
∂∆φ (r ) ∂∆φ (r ) dx + dy = P (r ) , ∂x ∂y
(184)
kde P(r) je diferenciál skalární funkce, kterou pro nás bude právě fázový rozdíl. Vztah (184) můžeme tedy přepsat jako ∂∆φ (r ) ∂∆φ (r ) dx + dy − d∆φ (r ) = 0 , ∂x ∂y
(185)
což je tzv. Pfaffova diferenciální rovnice. Grafické vyjádření situace je zobrazeno na obr. 21.
+Δϕ(r) +y ∂∆φ (r ) dx ∂x
+x
∂∆φ (r ) dy ∂y
dΔϕ
Obr. 21 Pfaffova diferenciální rovnice
Řešení Pfaffovy diferenciální rovnice je známé. Jednou z možností je postupná integrace. Volíme-li počáteční bod (x0, y0), integrací podél cest rovnoběžných s osami souřadnic dospějeme do bodu (x, y) podle vztahu ∂∆φ ∂∆φ dx + ∫ dy = ∆φ ( x, y ) − ∆φ ( x0 , y0 ) , ∫ ( x , y ) ∂x ( x , y ) ∂y ( x , y0 )
0
0
( x, y)
(186)
0
kde řešení je závislé na volbě počáteční podmínky Δϕ(x0, y0). To je důsledek toho, že řešení rovnice (185) generuje nekonečné množství ploch, jež jsou rovnoběžné, levá strana této rovnice
68
Interferometrie totiž bude platit pro jakoukoli z takto určených ploch. Abychom dosáhli jednoznačného řešení, je tedy nutné zvolit počáteční podmínku.
4.2.4 Konečné zpracování dat – postprocessing Výstupem vyhodnocení jsou zpravidla reálná data, jež ovšem potřebují další závěrečné úpravy. Jmenujme zejména: •
výběr a maskování oblastí, kde nebylo interferenční pole zaznamenáváno,
•
analytické vyjádření rekonstruovaných diskrétních dat rozvojem v řady polynomů (aproximace rekonstruovaných dat),
u metod zavádění nosné lineární prostorové frekvence dále •
odstranění zbytkových vlivů nejednoznačně určené prostorové frekvence.
Princip prvního kroku je totožný s výběrem relevantní vyhodnocované oblasti, o kterém bylo řečeno v kapitole předzpracování. Cílem je místům, kde nebylo zaznamenáno interferenční pole, přiřadit nulové nebo jinak uživatelem určené hodnoty. Další část práce podrobněji pojedná o aproximaci rekonstruovaných dat rozvojem v řady polynomů a o odstranění zbytkových vlivů nejednoznačně určené prostorové frekvence.
4.2.4.1 Analytické vyjádření diskrétních rekonstruovaných dat Rekonstruovaná data z vyhodnocení jsou diskrétní povahy, čili jde o matice reálných hodnot. Takováto forma je ovšem nevhodná pro další aplikace, jako například implementace do automatických procesů výroby, analytický popis zobrazovacích vlastností jednotlivých optických členů a celých optických soustav, kde se tyto komponenty vyskytují, a podobně. Vhodným způsobem je použití koeficientů rozvoje v řadu funkcí, kterými diskrétní data vyjádříme. Poté dostáváme několik málo numerických hodnot, které budou reprezentovat celou plochu. Tento způsob rekonstrukce budeme dále nazývat aproximace dat [79, 80]. Rekonstruovanou plochu reprezentovanou diskrétními měřenými daty si můžeme představit jako neznámou funkci dvou proměnných, kterou lze nahradit rozvojem v řadu jiných, námi známých, funkcí. Určením koeficientů rozvoje získáme analytické vyjádření plochy, což je naším cílem. Obecně je dobře znám rozvoj v zobecněnou Fourierovu řadu [39, 40], který vyžívá ortogonálních funkcí na určité oblasti. Díky ortogonalitě [39, 40] je výpočet koeficientů velmi snadný, čehož bylo využíváno zejména v době, kdy náročné numerické výpočty byly komplikovaným 69
Interferometrie problémem. S nástupem a neustálým vylepšováním výpočetní technicky není nezbytné využívat pouze ortogonální funkce, ale lze plochy (funkce) vyjadřovat rozvojem v libovolné funkce. K výpočtu koeficientů těchto rozvojů poté slouží numerické metody, například známá metoda nejmenších čtverců [39, 40, 44, 80, 81] nebo optimalizační algoritmy [82, 83]. Obecný postup aproximace pomocí ortogonálních systémů funkcí můžeme popsat následovně. Uvažujme ortogonální systém s vahou ρ(x, y) funkcí dvou proměnných φn(x, y) integrovatelných s kvadrátem na oblasti S [39, 40]. Dále předpokládejme, že rekonstruovaná plocha je reprezentována funkcí f(x, y) také integrovatelnou s kvadrátem na stejné oblasti. Z matematiky je dobře známo, že za výše uvedené situace lze funkci f(x, y) vyjádřit rozvojem v tzv. zobecněnou Fourierovu řadu, platí ∞
f ( x, y ) = ∑ α k ϕ k ( x , y ) ,
(187)
k =0
kde
αk =
∫∫ ρ ( x, y ) f ( x, y )ϕ
k
( x , y ) dx dy
S
∫∫ ρ ( x, y) ϕ
2
k
( x , y ) dx dy
.
(188)
S
Jestliže je systém φn(x, y) úplný, čili neexistuje žádná od nuly různá funkce ortogonální ke všem funkcím tohoto systému, potom zobecněná Fourierova řada (187) konverguje v průměru k funkci f(x, y), platí [39, 40] n
lim f ( x, y ) − ∑ α k ϕ k ( x, y ) = 0 . n →∞
(189)
k =0
Užijeme-li prvních N + 1 členů rozvoje (187), potom aproximujeme funkci f(x, y) tzv. částečným součtem zobecněné Fourierovy řady, N
f ( x, y ) ≈ ∑ α kϕ k ( x, y ) .
(190)
k =0
Pro diskrétní data, která analyzujeme v technické praxi, potom integrace analytická přejde v integraci numerickou, máme-li k dispozici dostatečně hustou síť měřených bodů. Užití rozvojů v ortogonální systémy funkcí naskýtá další výhody. Pro směrodatnou odchylku σf rozvoje (190) na oblasti S totiž dostáváme
70
Interferometrie
σ = 2 f
∫f
2
dS
S
∫ dS
=
S
=
1 ∫ dS
N N 1 N 2 ( ) (α jα kϕ jϕ k )dS α ϕ + ∑ ∑∑ k k ∫ k =0 j ≠ k ∫ dS S k = 0 S
∫ ∑ (α ϕ ) dS = ∑ σ N
S k =0
2
k
(191)
N
k
k =0
2 k
,
S
kde jsme využili ortogonálnosti funkcí na oblasti S a platí tedy [39, 40]
∫ ϕ ϕ dS = 0 . j
(192)
k
S
Z výsledku (191) je možno vyslovit tvrzení, že směrodatná odchylka aproximace pomocí jakéhokoli systému ortogonálních polynomů je minimální a přidání nebo odebrání členů neovlivní ostatní koeficienty aproximace a není je tedy nutné přepočítávat. Jak už bylo řečeno, s nástupem výpočetní techniky a jejím neustálým zlepšováním není třeba používat výhradně ortogonálních polynomů. Můžeme volit libovolné funkce, které svou lineární kombinací nejlépe vyjádří aproximovanou funkci ve smyslu námi zvolené podmínky. Odebrání nebo přidání dalších členů rozvoje je otázkou minimální námahy a výhoda snadných výpočtů koeficientů a jejich zachování tvaru při přidání nebo odebrání členů rozvoje již není tak důležitá. Volme nyní pro diskrétní měřená data označení nezávisle proměnných xi, yj (i = 0, ..., k; j = 0, ..., l). Výraz (190) lze maticově přepsat jako
f ( x0 , y0 ) ϕ0 ( x0 , y0 ) ϕ1 ( x0 , y0 ) f ( x1 , y0 ) ϕ0 ( x1 , y0 ) ϕ1 ( x1 , y0 ) ≈ L M L f (x , y ) ϕ (x , y ) ϕ (x , y ) k l 1 k l 0 k l f ≈ Aα ,
L ϕ N ( x0 , y0 ) α 0 L ϕ N ( x1 , y0 ) α1 M , L L L ϕ N ( xk , yl ) α N
(193)
kde f(xi, yj) jsou diskrétní data získaná z měření, φk(xi, yj) jsou vyčíslené funkce voleného systému v místech (xi, yj), αk jsou námi hledané koeficienty, k = 0, 1, ... , N. Chybu aproximace v jednotlivých bodech vyjádříme jako ε = Aα − f .
(194)
Na vztah (194) následně klademe podmínky a numerickým výpočtem chceme určit koeficienty α takovým způsobem, aby tyto podmínky byly splněny. Nejčastěji bývá volena podmínka minima součtu čtverců chyb, čili musí platit [44, 80, 81] M = εT ε = min .
(195)
71
Interferometrie Nalezení řešení lze odvodit například použitím funkcionální analýzy. Provedeme-li derivaci výrazu M = (Aα – f)T(Aα – f) podle koeficientů α (neboli vypočteme-li gradient funkce M(α)), dostáváme ∂M ∂ (Aα − f )T (Aα − f ) = ∂ (αT AT Aα − 2αT AT f + f T f ) = ∂α ∂α ∂α T T = 2 A Aα − 2 A f .
∇M (α ) =
(196)
Pro stacionární body αˆ poté užitím podmínky ∇M = 0 platí αˆ = (AT A ) A T f . −1
(197)
Pro Hessián funkce M [82, 83] dostáváme ∇ 2 M (α ) =
∂ 2M ∂ (2AT Aα − 2AT f ) = ∂ (2αT AT A − 2AT f ) = 2AT A . = 2 ∂α ∂α ∂α
(198)
Jelikož je matice (ATA) symetrická a tudíž pozitivně definitní [39, 40], potom koeficienty určené podle (197) představují minimum funkce M(α). O symetrii (ATA) se můžeme přesvědčit snadno dosazením z (193), kdy dostáváme pro jednotlivé prvky Mm,n, m = 1, ... , N+1, n = 1, ... , N+1, k
l
i =0
j =0
M m ,n = ∑∑ ϕm−1 ( xi , y j )ϕn−1 ( xi , y j ) .
(199)
Díky komutativnosti násobení dvou skalárních funkcí je patrné, že platí Mm,n = Mn,m, čili matice M(α) je symetrická. V některých případech je však inverze v (197) komplikovaná, výpočet může být nestabilní, problém je poté vhodné řešit optimalizačním algoritmem, například Levenberg-Marquardt [82, 83], kdy minimalizujeme tzv. meritní (cílovou) funkci (195) N + 1 proměnných. Proměnnými jsou koeficienty rozvoje a algoritmus hledá takové jejich hodnoty, pro které bude (195) minimální. Z matematiky jsou známy různé sady polynomů či funkcí vhodných k aproximaci na různých oblastech kvůli ortogonalitě. Výpočet metodou nejmenších čtverců však nevyžaduje ortogonální funkce na dané oblasti a poskytuje tak možnost volby. Dále není nezbytné volit pro aproximaci pouze podmínku metody nejmenších čtverců. Další možnou je kromě jiných např. minimalizace L1 normy M=
( k +1)( l +1)
∑ε i =1
i
.
(200)
72
Interferometrie Čili suma absolutních hodnot chyb aproximace v jednotlivých bodech má být minimální. V praxi je ovšem nejvíce využíváno zmíněné podmínky nejmenších čtverců. Lze totiž ukázat, že tato aproximace je také tzv. nejlepším nestranným odhadem a tzv. nejpravděpodobnějším odhadem, uvažujeme-li normální pravděpodobnostní rozdělení jednotlivých měření, viz např. [44]. Vhodným ukazatelem „kvality“ aproximace je aposteriorně určovaná výběrová směrodatná odchylka rozdílu rekonstruovaných a původních dat, jež je dána vztahem [44]
(Aαˆ − f ) (Aαˆ − f ) , T
s=
(201)
(k + 1)(l + 1) − N
kde čitatel představuje sumu čtverců rozdílů měřených dat od rekonstruovaných pomocí rozvoje v řadu systémů funkcí φk(x, y), kde koeficienty rozvoje jsou dané vztahem (197) nebo výsledkem optimalizačního procesu, a jmenovatel představuje počet nadbytečných měření, čili rozdíl počtu všech měřených dat (k + 1)(l + 1) a počtu určovaných N koeficientů. Jak je patrné, výběrová směrodatná odchylka může být spočtena pouze pro nadbytečný počet dat, čili kdy počet diskrétních měřených bodů bude větší, než počet členů rozvoje. To je v našem případě vyhodnocování interferometrických dat bezproblémově splněno. Obecně veškeré lineární polynomické aproximace, kdy rozvoj je dán jako lineární kombinace koeficientů a systému funkcí, lze zastřešit rozvojem v obecnou mocninnou řadu, kterou můžeme psát jako [39, 40, 80] f ( x, y ) ≈ ∑∑ α kl ( x − x0 ) ( y − y0 ) , N
M
k
l
(202)
k =0 l =0
kde x0 a y0 jsou konstanty, N a M představují stupeň rozvoje pro jednotlivé proměnné. Koeficienty αkl jsou počítány některým z výše uvedených postupů. V praxi se pro aproximaci čtvercových oblastí používají velmi často tzv. Legendreovy polynomy stupně n. Pro jednu dimenzi jsou definované jako [39, 40, 80]
Pn ( x) =
1 dn 2 (x − 1)n n n 2 n ! dx
(203)
jsou ortogonální pro −1 ≤ x ≤ 1 s vahou ρ(x) = 1 a lze je generovat přímo ze vztahu (203), nebo splňují rekurentní relaci
Pn ( x) =
2n − 1 n −1 xPn−1 ( x) − Pn−2 ( x) , n n
(204)
73
Interferometrie kde platí P0(x) = 1 a P1(x) = x. Rozšíření na dvě proměnné lze provést snadno jako N
N
f ( x, y ) ≈ ∑ ∑ α kl Pk ( x)Pl ( y ) .
(205)
k = 0 l =0
Ze čtvercové oblasti −1 ≤ x ≤ 1 a −1 ≤ y ≤ 1, na které jsou Legendreovy polynomy definovány, je lze zavedením substituce
x′ =
2 x − (b + a ) , b−a
y′ =
2 y − (d + c ) d −c
(206)
snadno transformovat na oblasti obdélníkové, kde platí a ≤ x´≤ b a c ≤ y´ ≤ d. Na obr. 22 je zobrazen vliv stupně aproximace Legendreovými polynomy na výsledná data.
Obr. 22 Aproximace původních dat (červené tečky) Legendreovým polynomem (modrá křivka) a) stupně N = 2 a c) stupně N = 6 a b), d) odchylky aproximovaných od původních hodnot
Pro aproximace oblastí kruhových se v praxi kontroly optických prvků velmi často využívají tzv. Zernikeho funkce [84, 85], které pro stupeň n a řád m vzniknou jako
Zmn (r , θ ) cos mθ m −m = R n (r ) sin mθ Zn (r , θ )
n = 0, 1, 2, K , m = − n, − (n − 1), K , (n − 1), n,
(207)
74
Interferometrie kde 0 ≤ r ≤ 1 a 0 ≤ θ ≤2π značí polární souřadnice a Rnm(r) je tzv. Zernikeho radiální polynom definovaný jako n −m 2
R mn (r ) = ∑ k =0
(− 1) (n − k )! k
n+m n−m k! − k ! − k ! 2 2
n = 0, 1, 2, K ,
r n− 2 k ,
(208)
kde (n − m) musí být sudé. Někdy bývají funkce (207) označeny jako Zn,μ(r, θ), kde μ = (n + m)/2. Rozvoj v řadu bude dán vztahem f ( x, y ) ≈ ∑∑ [α kl Zmn (r , θ ) + β kl Zn−m (r , θ )] N
N
k =0 l =0
= ∑∑ [α kl R (r ) cos mθ + β kl R (r ) sin mθ ] . N
N
m n
k =0 l =0
(209)
m n
Zernikeho funkce jsou ortogonální na jednotkovém kruhu D2, platí tedy [84, 85]
∫∫ Z
m n
(r , θ ) Zmn′′ (r , θ ) r dr dθ =
D2
∫∫ Z
n ,µ
(r , θ ) Zn′,µ′ (r , θ ) r dr dθ =
D2
1 δ nn′δ mm′ , 2n + 2 1 + δ µ ,n 2 2n + 2
(210)
πδ nn′δ µµ ′ ,
δij je tzv. Kroneckerovo delta [39, 40] (δij = 0 pro i ≠ j, δij = 1 pro i = j). Obdobně jako Legendreovy polynomy můžeme transformovat na obdélníkové oblasti, můžeme Zernikeho funkce transformovat na oblasti eliptické s hlavními poloosami a a b pomocí substituce tan θ ′ =
b tan θ , r ′ = a 2 cosθ + b 2 sin θ . a
(211)
Obr. 23 zobrazuje ukázku Zernikeho funkcí. Dva výše zmíněné případy jsou aproximacemi ortogonálními systémy funkcí. Jak bylo řečeno, se současnou technikou lze volit funkce v podstatě libovolné. Numericky jsou pak hledány takové koeficienty, jež splní požadované podmínky kladené na odchylky rekonstruovaných a původních dat. Jako příklad takovýchto aproximací ukažme tzv. Seidelovy polynomy [4, 7], které jsou definovány jako W (r ,θ ) = ∑∑ r n (α nm cos m θ + β nm sin m θ ) , N
n
(212)
n =0 m = 0
75
Interferometrie kde r a θ jsou polární souřadnice definované jako v předcházejícím případě, αnm a βnm jsou koeficienty rozvoje.
Obr. 23 Zernikeho funkce do stupně n = 3
V některých situacích apriori víme, že aproximovaná plocha bude rotačně symetrická. Potom musí platit W ( r ,θ ) = W ( − r , θ + π ) ,
(213)
a tedy (n – m) musí být sudé číslo, n ≥ m a vymizí všechny členy se sinθ. Velkou výhodou těchto polynomů je jejich jednoduchost a snadná fyzikální interpretace koeficientů rozvoje u rotačně symetrických oblastí. Oproti Zernikeho funkcím ovšem Seidelovy polynomy nejsou ortogonální. Zobecněním polynomických aproximací je využití racionálních funkcí [86, 87], kdy v čitateli i jmenovateli vystupuje polynom. Funkci, kterou chceme popsat, potom vyjádříme například pomocí vztahu N
f ( x, y ) ≈
i
∑∑ α i =0 j =0 M
k
xi− j y j
i
∑∑ β i =0 j =0
k
x
i− j
= y
j
i PN ( x, y ) , k = j + ∑t . Q M ( x, y ) t =1
(214)
76
Interferometrie Tento rozvoj již ale nelze rozepsat maticově tak snadno jako v (193) a nelze užít přímo odhadu koeficientů rozvoje metodou nejmenších čtverců (197), ve které byl předpokládán lineární rozvoj. Jednou z možností je použít optimalizační algoritmus [82, 83], který bude vyhledávat minimum meritní funkce M, jež bude mít pro podmínku metody nejmenších čtverců tvar k
l
M = ∑∑ ε ( xi , y j ) , 2
(215)
i =0 j =0
kde
ε ( xi , y j ) =
PN ( xi , y j ) Q M ( xi , y j )
− f ( xi , y j ) .
(216)
Optimalizační algoritmus najde lokální minimum meritní funkce (215), tedy hodnoty koeficientů ve smyslu nejmenších čtverců. Jedna z podmínek správného určení minima meritní funkce je ovšem nutnost zadání přibližných (počátečních) hodnot, které budou dostatečně blízko hledanému řešení. Samozřejmě bude záviset na tvaru meritní funkce, pro jednoduché průběhy bude optimalizace úspěšná i při velmi nepřesné volbě přibližných parametrů. V praxi takto jednoduché případy ovšem nastávají zřídkakdy. Pro určený přibližných hodnot počátečních parametrů můžeme u aproximace racionální funkcí využít tvaru, jenž budeme nazývat linearizovaný. Zvolíme-li si v rozvoji (214) β0 = 1 a získáme tak
f ( x, y ) ≈
α 0 + α 1 x + α 2 y + α 3 x 2 + α 4 xy + α 5 y 2 + K , 1 + β 1 x + β 2 y + β 3 x 2 + β 4 xy + β 5 y 2 + K
(217)
můžeme vyjádřit chybu aproximace v linearizované podobě jako N
i
M
i
ε ( x, y ) = ∑∑ α k x i − j y j − f ( x, y ) − f ( x, y )∑∑ β k x i − j y j . i =0 j =0
(218)
i =1 j = 0
Vztah (218) lze už snadno zapsat ve formě maticového zápisu pro množinu měřených dat jako
ε ( x0 , y 0 ) 1 x0 ε ( x0 , y1 ) 1 x0 = M M M ε ( x , y ) 1 x k l k
y0 y1
M yl
α 0 α L − x0 f ( x0 , y0 ) − y0 f ( x0 , y0 ) L 1 f ( x0 , y0 ) α L − x0 f ( x0 , y1 ) − y1 f ( x0 , y1 ) L 2 f ( x0 , y1 ) M − , (219) M M M M M β L − xk f ( xk , yl ) − yl f ( xk , yl ) L 1 f ( xk , yl ) β1 M ε = Aα − f .
77
Interferometrie Po této úpravě určíme hodnoty počátečních koeficientů pomocí (197) a můžeme přistoupit k minimalizaci meritní funkce (215) s opravami (216). Tento postup určení přibližných hodnot koeficientů a následné optimalizace je nutný, protože koeficienty z (219) neodpovídají správným koeficientům pro rozvoj (214) ve smyslu nejmenších čtverců. Výsledné řešení by nebylo „doslovnou“ aproximací ve smyslu nejmenších čtverců. Místo minimalizace sumy čtverců oprav (216) by byly totiž použity opravy
ε ( xi , y j ) = PN ( xi , y j ) − f ( xi , y j )Q M ( xi , y j ) .
(220)
Pro správné určení hodnot αk a βk je tedy nutné vzít výsledek z (219) pouze jako počáteční stav a dále minimalizovat meritní funkci (215) s opravami (216). Další možnost výpočtu koeficientů pro aproximaci racionální funkce je založena na úvaze rozvoje modelu (215) v Taylorovu řadu [39, 40]. Předpokládejme měřená data f(xi, yj) uspořádána do sloupcové matice Z. Rozvoj je závislý na volbě parametrů α (sloupcová matice koeficientů αk a βk). Označíme-li počáteční parametry jako α0 a jejich přírůstek jako dα, můžeme psát pro první
členy Taylorova rozvoje Z(α ) = Z(α 0 + dα ) ≈ Z(α 0 ) +
∂Z(α 0 ) dα , ∂α
(221)
kde Z(α0) vyjadřuje vyčíslení modelu (215) pro počáteční hodnoty parametrů (určené výše zmíněným postupem z linearizovaného modelu), Z(α) měřené hodnoty. Chybu linearizace (221) následně vyjádříme jako ε=
∂Z(α 0 ) dα − [Z(α ) − Z(α 0 )] = Adα − f ′ . ∂α
(222)
Tím je vyřešen problém nelineárnosti modelu (214), protože vyčíslíme-li hodnoty matice A (Jakobiho matice, matice derivací modelu podle počítaných parametrů) a f' pomocí měřených dat a funkčních hodnot pro počáteční α0, můžeme pomocí (197) určit vyrovnané hodnoty přírůstků a následně vyrovnané hodnoty parametrů, αˆ = α 0 + dαˆ .
(223)
Proces je ale třeba díky linearizaci iterativně opakovat, dokud se hodnoty vyrovnaných koeficientů nepřestanou lišit v jednotlivých řešeních o tolerovanou hodnotu.
78
Interferometrie Funkcí a polynomů, jež je možné využít k aproximaci interferometricky rekonstruovaných dat, je celá řada. V článku [80]8 lze nalézt příklady a analýzu aproximace asférických ploch i dalšími funkcemi,
např.
rozvojem
v trigonometrickou
Fourierovu
řadu
[39, 40, 42, 43],
řadu
Čebyševových polynomů [39, 40, 88] nebo sféroidálních vlnových funkcí [89, 90].
4.2.4.2 Odstranění zbytkových vlivů nejednoznačně určené prostorové frekvence Při rekonstrukci fázového rozdílu interferujících polí pomocí metod využívajících zaváděnou prostorovou lineární frekvenci nedojde při vyhodnocení k jejímu dostatečně přesnému určení. Rekonstruovaná data budou touto zbytkovou frekvencí ovlivněna. Zadání ovšem požaduje určit fázový rozdíl bez jejího vlivu. Seznámení s možností vyřešení tohoto problému má za úkol následující část. Předpokládejme, že rekonstruovaný fázový rozdíl Δϕr(x, y) bude ovlivněn zbytkovou lineární prostorovou frekvencí fr = (frx, fry). Potom můžeme psát ∆φr ( x, y ) = ∆φ ( x, y ) + 2 π( f rx x + f ry y ) ,
(224)
kde Δϕ(x, y) je fázový rozdíl bez vlivu nosné frekvence, již chceme určit. Jeden ze způsobů, jak určit velikost zbytkové nosné frekvence, je použití rozvoje v řadu polynomů nebo funkcí představených v předcházející kapitole. Každý z těchto rozvojů si můžeme totiž obecně zapsat jako f ( x, y ) = a0 + a1 x + a2 y + O ( x, y ) ,
(225)
kde ai, i = 0, 1, 2, značí koeficienty rozvoje až do lineárních členů, O(x, y) značí ostatní členy uvažovaného rozvoje libovolného stupně. Sloučením vztahů (224) a (225) dostáváme a0 + a1 x + a2 y + O ( x, y ) = ∆φ ( x, y ) + 2 π( f rx x + f ry y )
(226)
a porovnáním členů stejných mocnin poté
a0 + O( x, y ) = ∆φ ( x, y ) , a1 = 2 πf rx ,
(227)
a2 = 2 πf ry .
8
Jedná se o článek publikovaný autorem této práce v časopise Jemná mechanika a optika v říjnu roku 2013.
79
Interferometrie Je tedy patrné, že členy a1 a a2 nesou informaci o zbytkové nosné frekvenci, kterou se nepodařilo eliminovat při vyhodnocení. Rekonstruovaný fázový rozdíl bez vlivu této frekvence tedy dostaneme jako
∆φ ( x, y ) = ∆φr ( x, y ) − a1 x − a2 y ,
(228)
čímž je daný problém řešen.
80
5 Vyhodnocení interferogramů metodou Fourierovy transformace Zabývejme se nyní důkladněji vyhodnocením interferogramu se zavedenou prostorovou lineární nosnou frekvencí metodou Fourierovy transformace (FTM). Poprvé tuto metodu roku 1982 představil Takeda [73]. Od té doby se na jejím vývoji a zdokonalování podílelo mnoho vědců. Uvažujme již předzpracovaný interferogram, čili obraz filtrovaný od šumu a s vyrovnaným kontrastem. Potom můžeme zaznamenané rozdělení intenzity popsat vztahem, jenž jsme odvodili v předchozích kapitolách, I (r ) = a (r ) + b(r ) cos[2 π f 0 r + ∆φ (r )] ,
(229)
kde r = (x, y) je polohový vektor, a(r) představuje intenzitu pozadí společně s nechtěným zbytkovým šumem v intenzitě vznikajícím nerovnoměrným odrazem nebo přenosem světla od objektu, který se nepodařilo dostatečně odfiltrovat, b(r) je variace v amplitudě, jež také zahrnuje část nedostatečně odfiltrovaného šumu, f0 je prostorová přenosová frekvence, Δϕ(r) je fázový rozdíl dvou interferujících vln. Dále předpokládejme, že a(r), b(r) a Δϕ(r) se mění pomalu vzhledem k frekvenci f0. Naším cílem je nalézt výraz pro fázový rozdíl Δϕ(r). Užitím Eulerova vzorce můžeme (229) přepsat do pro nás vhodného tvaru jako 1 1 I (r ) = a (r ) + b(r ) exp{i[2 πf 0 r + ∆φ (r )]} + b(r ) exp{− i[2 πf 0 r + ∆φ (r )]} 2 2 ∗ = a (r ) + c(r ) exp(i 2 πf 0 r ) + c (r ) exp(− i 2πf 0 r ) ,
(230)
1 c(r ) = b(r ) exp[i∆φ (r )] , 2 1 c ∗ (r ) = b(r ) exp[− i∆φ (r )] 2
(231)
kde
jsou vzájemně komplexně sdružené funkce, i = √(–1). Proveďme nyní dvourozměrnou Fourierovu transformaci vztahu (230) definovanou jako FT {I (r )} =
∞ ∞
∫ ∫ I (r) exp[− i2πrf ] dr ,
(232)
−∞−∞
81
Vyhodnocení interferogramů metodou Fourierovy transformace čili FT {I (r )} =
∞ ∞
∫ ∫ a(r) exp[− i2πrf ] dr
−∞ − ∞
+
1∞∞ b(r ) exp[i∆φ (r )]exp[i 2 π f 0r ]exp[− i 2 πrf ] dr 2 −∫∞−∫∞
+
1∞∞ b(r ) exp[− i∆φ (r )]exp[− i 2 π f 0r ]exp[− i 2 πrf ] dr 2 −∫∞−∫∞
= +
∞ ∞
∞ ∞
−∞ − ∞
−∞ −∞
(233)
∫ ∫ a(r) exp[− i2πrf ] dr + ∫ ∫ c(r) exp[− i2πr(f − f )] dr 0
∞ ∞
∫ ∫ c (r) exp[− i2πr(f + f )] dr . ∗
0
− ∞ −∞
Výraz (233) můžeme dále zapsat jako FT {I (r )} = A(f ) + C (f − f0 ) + C ∗ (f + f 0 ) ,
(234)
kde jsme užili A(f ) =
∞ ∞
∫ ∫ a(r ) exp[− i2πrf ] dr ,
− ∞ −∞
C (f − f0 ) = =
1∞∞ b(r ) exp[i∆φ (r )]exp[− i 2 πr (f − f 0 )] dr 2 −∫∞−∫∞ ∞ ∞
∫ ∫ c(r ) exp[− i2πr(f − f )] dr , 0
(235)
− ∞ −∞
C ∗ (f + f0 ) = =
1∞∞ b(r ) exp[− i∆φ (r )]exp[− i 2 πr (f + f0 )] dr 2 −∫∞−∫∞ ∞ ∞
∫ ∫ c (r ) exp[− i2πr(f + f )] dr . ∗
0
− ∞ −∞
Bude-li se a(r), b(r) a Δϕ(r) měnit pomalu vzhledem k frekvenci f0, potom budou Fourierova spektra A(f ) , C (f − f0 ) a C ∗ (f + f 0 ) v rovnici (234) oddělena nosnou frekvencí f0 a půjde je snadno filtrovat. Maxima C (f − f0 ) a C ∗ (f + f 0 ) se budou nacházet právě v pozici –f0 a f0, jak je patrné z (235), kde tyto pozice určují maxima exponenciely. Symbolické znázornění rozložení frekvenčního spektra (234) pro jednodimenzionální případ je zobrazeno na obr. 24. Frekvenční spektrum interferogram bez přítomnosti šumu je následně na obr. 25.
82
Vyhodnocení interferogramů metodou Fourierovy transformace
+ FT {I (r )}
–f0x
+f0x
+fx
Obr. 24 Symbolické znázornění Fourierovských spekter se zavedenou prostorovou nosnou frekvencí pro jednodimenzionální případ
Obr. 25 a) Interferogram bez přítomnosti šumu a b) jeho reálná složka frekvenčního spektra
Nyní proveďme inverzní Fourierovu transformaci spektra C ∗ (f + f 0 ) definovanou jako FT −1{C ∗ (f + f 0 )} =
∞ ∞
∫ ∫ C (f + f ) exp[i2πrf ] df , ∗
0
(236)
−∞ −∞
tedy po dosazení dostáváme
∞ ∞ ∗ ∫−∞−∫∞−∫∞−∫∞c (ρ) exp[− i2πρ(f + f0 )] dρ exp[i2πrf ] df ∞ ∞ ∞ ∞ = ∫ ∫ ∫ ∫ c∗ (ρ) exp[− i 2 πρf0 ] dρ exp[i 2 π(r − ρ )f ] df −∞ − ∞ −∞ −∞
FT −1{C ∗ (f + f0 )} =
=
∞ ∞
∞ ∞
(237)
∫ ∫ c (ρ) exp[− i2πρf ]δ (r − ρ) dρ , ∗
0
−∞ − ∞
kde δ(r – ρ) je Diracova funkce [39, 40]. Z jejích vlastností plyne, že v bodech spojitosti bude integrace (237) rovna
83
Vyhodnocení interferogramů metodou Fourierovy transformace 1 FT −1{C ∗ (f + f0 )} = b(r ) exp{− i[2 πrf 0 + ∆φ (r )]} 2 = c ∗ (r ) exp[− i 2 πf 0r ] .
(238)
Logaritmujeme-li výraz (238), dostáváme
1 ln FT −1{C ∗ (f + f0 )} = ln b(r ) + ln exp{− i[∆φ (r ) + 2πrf 0 ]} 2 1 = ln b(r ) − i[∆φ (r ) + 2πrf0 ] . 2
(239)
Hledaná funkce Δϕ(r) je tak získávána z imaginární části výrazu ln FT −1{C ∗ (f + f 0 )} , zatímco nechtěné variace v amplitudě z části reálné. Snadnou úpravou (239) dostáváme
∆φ (r ) = − Im[ln FT −1{C ∗ (f + f0 )}] − 2 πrf0 .
(240)
Jiným a zejména používaným [73, 91, 92] způsobem výpočtu funkce Δϕ(r) je užití vlastností komplexních čísel a Eulerova vzorce, kdy platí
c = K exp(− iα ) = K (cos α − i sinα ) ,
(241)
a dále tedy
Im(c) sinα = − arctan . cos α Re(c)
α = arctan
(242)
S (238) a (242) tak můžeme psát
Im[FT −1{C ∗ (f + f0 )}] 2πrf0 + ∆φ (r ) = − arctan . −1 ∗ Re[FT {C (f + f0 )}]
(243)
Hledanou funkci Δϕ(r) tudíž snadno získáme jako
Im[FT −1{C ∗ (f + f 0 )}] ∆φ (r ) = − arctan − 2πrf 0 . −1 ∗ Re[FT {C (f + f 0 )}]
(244)
Ve vztahu (244) jsme se nepříznivých variací v amplitudě b(r) zbavili podílem imaginární a reálné části. Obdobně bychom postupovali při určení funkce Δϕ(r) užitím spektra C (f − f0 ) , stejným postupem jako v předchozím získáváme
Im[FT −1{C (f − f0 )}] ∆φ (r ) = arctan − 2πrf0 . −1 Re[FT {C (f − f0 )}]
(245)
84
Vyhodnocení interferogramů metodou Fourierovy transformace Ze vztahů (244) a (245) je patrné, že pro jednoznačné určení fázového rozdílu Δϕ(r) musíme přesně určit prostorovou nosnou frekvenci f0, jelikož ji musíme odečítat, a tak nám její hodnota ovlivní výsledek. V praxi se nepoužívá výše zmíněný postup doslova. Neprovádí se inverzní Fourierova transformace přímo na C (f − f0 ) nebo C ∗ (f + f 0 ) . Funkce „arctan” je ve výpočetních softwarech definována zpravidla na intervalu od –π do π, a tak v rekonstruovaném obraze dostáváme místa nespojitosti, mezi nimiž by díky přítomnosti prostorové nosné frekvence měla rekonstruovaná plocha „velké sklony“ (velké gradienty). Při odstraňování nespojitostí by poté docházelo ke značným komplikacím, jež jsou způsobeny právě vlivem prostorové nosné frekvence. Postup je volen tedy jiný [73, 91, 92] a může být zapsán v následujících bodech: •
provedení Fourierovy transformace zaznamenaného obrazu,
•
odfiltrování střední a jedné z bočních složek frekvenčního spektra,
•
vyhledání maxima druhé boční složky frekvenčního spektra,
•
centrace této boční složky, čili f 0 → 0 , C ∗ (f + f0 ) → C0∗ (f ) , C (f − f 0 ) → C0 (f ) ,
•
inverzní Fourierova transformace na centrovanou složku,
•
odstranění nespojitostí rekonstruovaného obrazu (tzv. unwrapping).
Jelikož vektor prostorové nosné frekvence tímto způsobem (centrací) přejde ve vektor nulový, po inverzní transformaci bočních spekter dostáváme úpravou (238) výrazy
1 FT −1{C0∗ (f )} = b(r ) exp{− i∆φ (r )} = c∗ (r ) , 2 1 FT −1{C0 (f )} = b(r ) exp{i∆φ (r )} = c(r ) , 2
(246)
rekonstruovaný fázový rozdíl tedy snadno získáme jako
Im[FT −1{C0∗ (f )}] Im[FT −1{C0 (f )}] ∆φ (r ) = − arctan = arctan . −1 ∗ −1 Re[FT {C 0 (f )}] Re[FT {C 0 (f )}]
(247)
Tento postup ovšem naráží na nutnost správného určení maxima bočních spekter a tím i prostorové nosné frekvence. Jestliže data nebudou obsahovat šum, poté je to problém triviální a ve výpočetních softwarech snadno proveditelný. Pro reálná data to ovšem tak jednoduché není. Přítomnost šumu může značně ovlivňovat polohu maxima bočního spektra. Získávaná maximální hodnota z dat nebude odpovídat poloze nosné frekvence, a centrace tak nebude správná. 85
Vyhodnocení interferogramů metodou Fourierovy transformace Problém určení polohy maxima lze řešit následujícím snadným způsobem. Představme si situaci znázorněnou na obr. 26. Levá a střední složka byly odfiltrovány. V pravé složce šum ovšem způsobí přítomnost „nesprávného“ maxima. Polohu odpovídající nosné frekvenci lze poté nalézt váženým průměrem, a to sice jako
∑ f FT {I (r )} = ∑ FT {I (r )} i
f0
i
i
,
(248)
i
i
kde fi znační frekvence, pro které jim odpovídající hodnoty FT {I (r )}i přesáhly volenou konstantu T.
odfiltrované složky
+ FT {I (r )}
maximum bočního spektra za přítomnosti šumu T +f0x
+fx
Obr. 26 Určení maxima bočního spektra
Tímto způsobem lze odhadnout polohu, jež odpovídá zavedené prostorové nosné frekvenci. Výsledek není ale také zcela stoprocentní. Pokud by hodnota nesprávného maxima byla velká, jeho váha by velmi ovlivnila vážený průměr. Jestliže hodnota prostorové nosné frekvence nebude určena přesně, lze vliv nejistoty potlačit úpravou rekonstruovaného obrazu. Tento postup byl představen v části konečného zpracování dat. Nabízí se také varianta vyhledání polohy maxima robustními způsoby s detekcí odlehlých vzorků [44]. Po určení prostorové nosné frekvence přichází na řadu samotná filtrace spektra. Použití masek s ostrým přechodem (kdy dochází ke skoku mezi hodnotami masky nespojitě) vede k nechtěným jevům ve výsledném obraze. Používají se tedy filtry založené na pozvolném potlačování nežádoucích hodnot, např. filtry založené na tzv. Hammingově okně a dalších technikách běžně používaných ve zpracování obrazu [47, 48].
86
Vyhodnocení interferogramů metodou Fourierovy transformace Jak už bylo výše zmíněno, hodnoty funkce „arctan“ jsou ve výpočetních softwarech zpravidla definovány v rozmezí od –π do π. Ideální vyhodnocený obraz (bez šumu a dalších vlivů) bude tím pádem obsahovat místa nespojitosti velikosti 2π. Odstraněním těchto míst nespojitosti se zabývá tzv. unwrapping. Tomuto matematickému problému byla, a stále je, v praxi věnována značná pozornost, protože na jeho vyřešení závisí přesnost celého výsledku. Odstranění nespojitostí u ideálního obrazu je triviálním problémem a lze ho řešit snadným algoritmem. Pro data obsahující šum, čili reálná data, je situace komplikovanější. Algoritmům unwrappingu se věnuje další kapitola této práce. Obr. 27 shrnuje výše popsaný postup vyhodnocení interferogramu pomocí FTM.
Předzpracovaný interferogram
Fourierova transformace, filtrace a centrace bočního spektra
Inverzní Fourierova transformace
Odstranění nespojitostí (unwrapping) a zbytkového vlivu nosné frekvence
Rekonstruovaný dráhový rozdíl interferujících vln Obr. 27 Schéma vyhodnocení interferogramu pomocí FTM
87
Vyhodnocení interferogramů metodou Fourierovy transformace Užitím metody Fourierovy transformace jsme získali námi hledanou funkci Δϕ(r). Pro určení prostorového rozdílu této funkci odpovídajícímu užijeme již dříve odvozený vztah W (r ) =
λ 2π
∆φ (r ) .
(249)
Pro určení fázového rozdílu nám tedy při použití FTM postačuje jeden obraz interferenčního pole. Ze vztahu pro výpočet Δϕ(r) je vidět, že dojde k potlačení vlivu variací v amplitudě. Dále tím, že při samotném výpočtu užíváme frekvenčního spektra, můžeme provádět filtrace ve frekvenční oblasti, a získávat tak kvalitnější výsledky s odfiltrovaným šumem. FTM tedy sama o sobě umožňuje potlačit nepříznivé vlivy. Nastávají při ní ale některé komplikace, jako nejednoznačné určení prostorové nosné frekvence, které vyhodnocení ztěžují. Použití FTM dále naráží na to, že teoreticky je při odvozování používáno Fourierovy transformace na neomezeném intervalu, ovšem při reálném zpracování je použito její finitní diskrétní formy. Ta ovlivňuje zejména oblasti v hraničních částech obrazu, kde dochází k odchylkám od skutečných hodnot. Existují metody, jimiž lze tento nepříznivý jev částečně potlačit, ovšem není prozatím známý jednoznačný způsob, který by problém řešil úplně. Jednou z metod je maskování pruhu pixelů po kompletním vyhodnocení v hraničních oblastech. Čili výsledná data nebudou vůbec hodnoty o fázovém rozdílu v těchto částech obrazu obsahovat. Další metodou je extrapolace interferenčních proužků mimo oblast zaznamenaného interferenčního pole. Iterativní algoritmus představuje Gerchberg v [93], Roddier a Roddier v [91] a jeho modifikovanou verzi poté Massig a Heppner v [92], kde je použito speciálního frekvenčního filtru s „měkkými“ hranami. Jinou možnost extrapolace představuje Dai a Wang v [94]. Princip Gerchbergova algoritmu je založen na následující myšlence. Uvažujme zaznamenané rozložení intenzity interferenčního pole dané vztahem I (r ) = M (r )[a (r ) + c(r ) exp(i 2 πf 0 r ) + c ∗ (r ) exp(−i 2 πf 0 r )] ,
(250)
kde M(r) je funkce masky, jež pomocí hodnot 0 a 1 určuje, zda se jedná o prostor interferenčního pole nebo ne, c(r) a c*(r) jsou dány vztahy (231). Aplikací Fourierovy transformace na (250) dostáváme [91, 93] FT {I (r )} = FT {M (r )} ∗ [ A(f ) + C (f − f 0 ) + C ∗ (f + f 0 )] ,
(251)
kde * značí konvoluci. Pokud M(r) = 1 pro všechna r, potom při dostatečně velké zaváděné nosné
frekvenci
bude
frekvenční
spektrum
(251)
rozděleno
do
tří
jednoznačně
88
Vyhodnocení interferogramů metodou Fourierovy transformace identifikovatelných částí – středního a dvou středově symetrických bočních. Boční spektra budou situována v kruhových oblastech s poloměrem fc a pro jejich jednoznačné rozlišení od středního musí platit |f0| > fc. Jestliže M(r) bude obsahovat hodnoty 0 a 1, poté její nespojitý charakter ovlivní frekvenční spektrum (251). Uvažujme nyní vztah (250) ve formě I (r ) − M (r ) a (r ) = M (r ) K (r ) cos[2 πf 0r + ∆φ (r )] ,
J (r ) =
(252)
kde jsme využili vztahu pro kontrast K (r ) =
b(r ) . a (r )
(253)
Fourierovu transformaci aplikovanou na (252) lze poté vyjádřit jako [91, 93] FT {J (r )} = FT {M (r )} ∗ [C J (f − f0 ) + C J* (f + f 0 )] ,
(254)
kde CJ a CJ* jsou boční spektra upraveného interferogramu J(r). Kdyby nebyl omezený maskou M(r), byla by spektra koncentrována ve dvou kruzích o středech +f0 a –f0 s poloměrem fc. Pokud tedy položíme všechny hodnoty mimo oblast těchto kruhů rovnu nule, potom inverzní Fourierovou transformací dostaneme upravený obraz s interferenčními proužky i mimo místa původního zaznamenaného interferenčního pole. Jelikož ale neznáme poloměr fc přesně, je jeho hodnota volena experimentálně a algoritmus probíhá iterativně. Při každé iteraci jsou navíc na místa původního obrazu dosazovány původní hodnoty, aby bylo zachováno dané měření beze změny. Schéma extrapolačního algoritmu je ukázáno na obr. 28.
PŮVODNÍ UPRAVENÝ OBRAZ J(r)
FOURIEROVA TRANSFORMACE
NULOVÁNÍ SPEKTER MIMO VHODNÉ OBLASTI
N iterací EXTRAPOLOVANÝ OBRAZ
NAVRÁCENÍ ORIGINÁLNÍCH HODNOT
INVERZNÍ FOURIEROVA TRANSFORMACE
Obr. 28 Schéma extrapolace interferenčních proužků
89
Vyhodnocení interferogramů metodou Fourierovy transformace Vztah (252) předpokládá znalost průběhu členu a(r), čili funkci střední hodnoty intenzity pozadí zaznamenaného interferenčního pole. Její určení je možno provést různými způsoby. Jedním z nich je např. postup popsaný v kapitole o předzpracování interferogramu, kdy použitím dostatečně širokého průměrujícího filtru tuto střední hodnotu intenzity můžeme odhadnout.
90
6 Unwrapping 6.1
Úvod do problematiky unwrappingu
Data, rekonstruovaná fáze interferenčního pole9, získaná vyhodnocením interferogramu pomocí Fourierovy transformace budou zpravidla v rozmezí hodnot od –π do +π v závislosti na definici funkce „arctan“ ve výpočetním softwaru. Ideální vyhodnocený obraz (bez šumu a dalších vlivů) bude tedy obsahovat místa nespojitosti velikosti 2π. Odstraněním těchto míst se zabývá tzv. unwrapping10. Matematický problém odstraňování nespojitosti rekonstruovaného obrazu fáze je odborníky řešen po několik desetiletí [68]. Prvotní myšlenka je velmi snadná, ovšem je aplikovatelná pouze na ideální data, čili data vznikající z interferogramů nezatížených šumem, diskontinuitami nebo residui. Tyto jevy ovšem na reálných datech jsou. Šum vzniká přirozeně už při samotném snímání obrazu. V současnosti záznam probíhá nejčastěji na CCD (z angl. Charge-Coupled Device) senzory [56], nebo jim obdobné, kde není stoprocentně zaručena stejná citlivost všech pixelů. Snímaný obraz homogenního rovnoměrného osvětlení potom nemusí být homogenní, bude zatížen náhodným šumem. Diskontinuity a residua mohou při vyhodnocování optických ploch vznikat např. v důsledku náhodných prachových částic (nečistot) na snímači nebo vyhodnocovaném objektu nebo drobnými poškozeními povrchů. Při radarové interferometrii (SAR a InSAR Interferometry [37, 38] atp.) samotný zkoumaný povrch není spojitou plochou, uvážíme-li přítomnost lesů, stavebních objektů a dalších. Zejména radarová interferometrie a vyhodnocování fází při ní přispělo velkou mírou k rozvoji algoritmů unwrappingu, které jsou dost sofistikované na to, aby alespoň částečně eliminovaly výše popsané komplikace. Základní přehled algoritmů přináší Ghiglia a Pritt [68], kteří ve své práci detailně popisují jejich principy, aplikují je na případech SAR interferometrie a poskytují jejich zdrojové kódy. Následuje celá řada autorů, např. [95 – 98], kteří publikovali různé další algoritmy a postupy, jimiž lze tuto poměrně komplikovanou úlohu řešit.
9
Fáze interferenčního pole získávaná rekonstrukcí interferogramu odpovídá fázovému rozdílu dvou interferujících vln (případně s aditivně přidanou nosnou frekvencí).
10
Překlad slova „unwrapping“ z angličtiny by zněl asi jako „odstranění nespojitostí obrazu“. V textu se ovšem budeme držet anglického tvaru, zejména pro jednoduchost zápisu. Od toho se také odvíjí použití spojení „unwrapped fáze“ a „wrapped fáze“, čili fáze s odstraněnými nespojitostmi a fáze s nespojitostmi doposud neodstraněnými.
91
Unwrapping Pro pochopení problému nejprve uvažujme ideální obraz fáze bez šumu, diskontinuit (původní fáze je tedy spojitá) a residuí. Odstranění nespojitostí potom můžeme provést následujícím jednoduchým způsobem [73]. Vytvoříme novou funkci ϕc(x,y) (c – z angl. correction), se kterou přičtením k rekonstruované fázi ϕw(x,y) (w – z angl. wrapped) dostaneme výslednou fázi bez nespojitostí ϕu(x,y) (u – z angl. unwrapped), čili platí
φu ( x, y ) = φw ( x, y ) + φc ( x, y ) .
(255)
Funkce ϕc(x,y) bude v počátečním kroku nulová. Následuje výpočet fázové diference Δϕw(x,y). Nejprve provedeme unwrapping pouze v jednom směru – ve směru osy x. Jednotlivé pixely obrazu fáze označujme indexem i a j, kde 0 ≤ i ≤ N, 0 ≤ j ≤ M. Potom pro fázovou diferenci počátečního řádku obrazu, y = 0, píšeme ∆φw ( xi ,0) = φw ( xi ,0) − φw ( xi −1 ,0) .
(256)
V místech, kde není nespojitost, bude diference malá, což vyplývá z našeho úvodního předpokladu spojitosti původní fáze. Tam, kde nespojitost je, se ovšem bude v absolutní hodnotě blížit k 2π. Na i-té pozici potom přičteme, resp. odečteme, k funkci ϕc(x,y) hodnotu 2π. Přičtení nebo odečtení 2π bude záviset na znaménku diference. Jestliže bude záporná, budeme 2π přičítat, a naopak. Pro první řádek funkce ϕc(x,y) tedy platí +0, pokud ∆ wφ ( xi ,0) < T , φc ( xi ,0) = φc ( xi−1 ,0) + 2π , pokud ∆φw ( xi ,0) ≥ T ∧ ∆φw ( xi ,0) < 0 , − 2 π , pokud ∆φ ( x ,0) ≥ T ∧ ∆φ ( x ,0) > 0 , w i w i
(257)
kde T značí testovací hodnotu, např. 0.9∙2π, rozhodující o tom, zda na i-té pozici je diskontinuita nebo není. Uvážíme-li Nyquistovo kritérium [66, 67], potom pro jeden proužek interferogramu potřebujeme minimálně 2 pixely obrazu fáze. Z toho plyne, že fáze vlnoplochy se nezmění mezi 2 pixely o více než π. Následně můžeme volit testovací hodnotu například T = 0.9∙π. Od prvního řádku se bude odvíjet unwrapping pro sloupce stejným způsobem. Nejprve výpočtem diference, ∆φw ( xi , y j ) = φw ( xi , y j ) − φw ( xi , y j −1 ),
j = 1, K , M ,
(258)
následně rozhodnutím o znaménku pro funkci ϕu(x,y) a přičtením nebo odečtením 2π,
92
Unwrapping +0, pokud ∆φw ( xi , y j ) < T , φc ( xi , y j ) = φc ( xi , y j −1 ) + 2π , pokud ∆φw ( xi , y j ) ≥ T ∧ ∆φw ( xi , y j ) < 0 , − 2 π , pokud ∆φ ( x , y ) ≥ T ∧ ∆φ ( x , y ) > 0 . w i j w i j
(259)
Následující obr. 29 zobrazuje testovací data ϕ(x,y) (použita byla funkce „peaks“ v prostředí MATLAB) v reálné podobě a ve wrapped formě, která je pro testování získána snadno jako
φ w ( x, y ) = arctan[sin φ ( x, y ) cos φ ( x, y )] .
(260)
Jiná možnost simulace wrapped dat je dána vztahem φ ( x, y ) , 2π
φw ( x, y ) = φ ( x, y ) − 2πR
(261)
kde R{•} značí operátor zaokrouhlení k nejbližšímu celému číslu. Obr. 30 ukazuje data rekonstruovaná za použití výše uvedeného jednoduchého algoritmu a rozdíl rekonstruovaných dat od dat původních.
Obr. 29 Zobrazení a) testovacích dat (funkce „peaks“) v původní podobě a b) ve wrapped formě
93
Unwrapping
Obr. 30 Zobrazení a) rekonstruovaných dat po unwrappingu, b) rozdílu rekonstruovaných dat po unwrappingu od skutečných hodnot
Jak je z výše uvedeného obrázku patrné, byla původní spojitá fáze snadno zrekonstruována z wrapped formy s přesností v řádu 10-15, čili s přesností odpovídající zaokrouhlovacím chybám počítače. Než přejdeme k uvedení algoritmů využívaných u reálných dat, která nejsou bez šumu a residuí, zmíníme se o unwrappingu ideálního obrazu, který ovšem nebude dostatečně vzorkován. Jak již bylo vícekrát zmíněno, v důsledku Nyquist-Shannonova kritéria nesmí mezi dvěma pixely být větší změna fáze než π, čili |Δϕ| ≤ π. Diferenciální přírůstek fáze získáme snadno jako ∆φ ( x) =
∂φ ∆x , ∂x
(262)
kde Δx představuje vzdálenost mezi dvěma pixely. Jestliže zaznamenáváme obraz v rozsahu xmin a xmax s počtem N vzorků, potom snadno dostáváme Δx = (xmax – xmin)/N. Dosazení do výše uvedeného vztahu vede k podmínce pro minimální počet vzorků
∂φ xmax − xmin ≤π , ∂x N ∂φ xmax − xmin N> . ∂x π
(263)
Pro ukázku vlivu nedostatečného vzorkování (z angl. under sampling) na unwrapping uvažujme fázi danou vztahem ϕ(x) = 10∙cos 2x na intervalu 0 ≤ x ≤ 1. Užitím (263) dostáváme N>
20 sin 2 x . π
(264)
94
Unwrapping Jelikož maximální hodnota funkce sin 2x je 1, dostáváme tedy podmínku pro dostatečné vzorkování N > 20/π ≈ 6.37. Obraz bude tím pádem vhodně zaznamenán při N ≥ 7. Na obr. 31 a obr. 32 je daná situace zobrazena.
Obr. 31 a) Funkce ϕ(x) = 10∙cos 2x (modrá čárkovaná) a vzorkování při N = 10 (červeně) a N = 5 (zeleně); b) wrapped fáze a její vzorkování
Obr. 32 a) Skutečná fáze (modře čárkovaně) a unwrapping při dostatečném vzorkování (červeně); b) skutečná fáze (modře čárkovaně) a unwrapping při nedostatečném vzorkování (zeleně)
Jak je na první pohled patrné, při nedostatečném vzorkování unwrapping selhává. Důvod je jednoduchý. Mezi prvním a druhým bodem na obr. 31 b) při nedostatečném vzorkování (zeleně) algoritmus určí rozdíl hodnot větší než π, což je správě a unwrapping mezi těmito body proběhne korektně, jak je patrné na obr. 32 b). Ovšem u dalších bodů již tomu tak není. Mezi druhým a třetím je rozdíl větší než π, čili algoritmus přičetl hodnotu 2π, i když neměl. Následně je rozdíl menší než π, ale body leží za zlomem, čili korekce měla proběhnout, ale neproběhla. Výše popsaný jednoduchý algoritmus unwrappingu tedy při nedostatečném vzorkování selhal.
95
Unwrapping U reálného obrazu fáze úloha unwrappingu nebude tak snadná i při dostatečném vzorkování. Obraz bude obsahovat místa, pixely, která nebudou odpovídat skutečné hodnotě, jakou bychom dostali pro reálný obraz. Tato místa nazýváme residua. Je tedy nutné hodnotit kvalitu obrazu a unwrapping už není tak snadným jako v ideálním případe. Vytváříme poté tzv. masky nebo mapy kvality, které jsou vstupním údajem pro dané algoritmy. Pro pochopení sofistikovanějších algoritmů uveďme myšlenku, kterou představil Itoh (1982) [99]. Nechť W je tzv. wrap operátor, kterým získáme wrapped fázi v hodnotách mezi –π a +π, čili v jedné dimenzi pišme
φw ( xi ) = W {φ ( xi )} = φ ( xi ) + 2 πk ( xi ) ,
(265)
kde xi značí polohu v obraze a k(xi) je diskrétní funkce celých čísel, 0 ≤ i ≤ N. Tato funkce může být reprezentována např. pomocí vztahu (261). Definujme dále operátor diference Δ jako ∆{φ ( xi )} = φ ( xi ) − φ ( xi −1 ) , ∆{k ( xi )} = k ( xi ) − k ( xi −1 ) .
(266)
Použitím operátoru diference na wrap operátor dostáváme s (265) a (266) výraz ∆{W {φ ( xi )}} = ∆{φ ( xi )} + 2 π∆{k ( xi )}.
(267)
Opětovné použití wrap operátoru na rovnici (267) vede k W {∆{W {φ ( xi )}}} = ∆{φ ( xi )} + 2 π[∆{k1 ( xi )} + k 2 ( xi )] ,
(268)
kde k1(xi) a k2(xi) značí funkce celých čísel příslušných k prvnímu a druhému wrap operátoru. Předpokládáme-li dále, že bude platit tzv. podmínka hladkosti − π < Δ{φ ( xi )} ≤ π ,
(269)
čili že absolutní hodnota rozdílu dvou sousedních bodů fáze nebude větší než π, potom díky definici wrap operátoru, pro který také platí rozsah oboru hodnot od –π do +π, vyplývá 2 π[∆{k1 ( xi )} + k 2 ( xi )] = 0 ,
(270)
∆{φ ( xi )} = W {∆{W {φ ( xi )}}} = W {∆{φ w ( xi )}} .
(271)
a tedy
Je patrné, že fázi, jež splňuje podmínku hladkosti (269), můžeme rekonstruovat pomocí vztahu
96
Unwrapping
φ ( xi ) = φ ( x0 ) + ∑W {∆{φ w ( x j )}}. i −1
(272)
j =1
Ve dvou dimenzích lze postupovat následovně. Nejprve se například vyhodnotí první řádek (x = 0),
φ ( x0 , y ) = φ ( x0 , y0 ) + ∑ W {∆ y {φw ( x0 , y j )}}, M
(273)
j =1
kde W značí wrap operátor a Δy operátor diference ve směru řádku. Sloupce poté můžou být unwrappovány od tohoto úvodního řádku,
φ ( x, y ) = φ ( x0 , y j ) + ∑ W {∆ x {φ w ( xi , y j )}}, N
(274)
i =1
kde Δx je operátor diference ve směru sloupce. Principielní schéma ukazuje obr. 33
(x0,y0)
(x0,yj)
(xN,yM)
Obr. 33 Principielní schéma jednoduchého 2D unwrappingu Itohovým algoritmem
Na obr. 34 a obr. 35 je tento Itohův algoritmus demonstrován pro jednodimenzionální případ. Itohův algoritmus je speciálním případem N-dimenzionálního vztahu
φ (r ) = ∫ ∇φ d N r + φ (r0 ) ,
(275)
C
který říká, že unwrapping fáze po křivce C z bodu r0 do bodu r lze provést jako integraci fázového gradientu po křivce C, pokud známe hodnotu gradientu podél této křivky a hodnotu fáze v bodě r0.
97
Unwrapping
Obr. 34 a) Wrap operátor použitý na funkci ϕ(xi) = 10∙cos(2xi), b) operátor diference použitý na výsledek z části a)
Obr. 35 Zobrazení a) rekonstruované funkce pomocí Itohova algoritmu a b) rozdíl od původní hodnoty
Vztah (275) je základem pro detekování residuí ve vyhodnocovaném obrazu fáze. Pokud by byla fáze (interferogram) ideální, bude platit nejen vztah (275), ale také
∫ ∇φ d
N
r = 0,
(276)
C
který říká, že integrace fázového gradientu po libovolné uzavřené křivce C je nulová, není-li v cestě integrace diskontinuita (residuum). Pokud tomu tak není, může být residuum v daném místě přítomno. V případě, že fázové gradienty ∇φm , které můžeme získat z měření, jsou zatíženy šumem, může algoritmus řešit obecný problém, kdy dochází k minimalizaci čtverce chyby ε [68],
98
Unwrapping
ε 2 = ∫ P(∇φ m − g ) dS , 2
(277)
S
ve smyslu nejmenších čtverců. Ve vztahu (277) představuje P váhovou funkci, dS elementární plochu a g gradient ovlivněný šumem, g = ∇φ + n ,
(278)
n je vektor šumu, který ovlivňuje skutečný fázový gradient ∇φ .
Dle detekovaných poloh residuí jsou jednotlivým pixelům přiřazovány váhy v rozmezí intervalu hodnot 0 a 1. Použijeme-li pouze hodnoty vah 0 nebo 1, potom používáme pro unwrapping tzv. masky. Jestliže jsou vybírány hodnoty z celého intervalu, používáme tzv. mapy kvality (QM – z angl. Quality Maps). Masky i QM mají stejný rozměr jako analyzovaný obraz fáze. Jednotlivé algoritmy je používají zároveň nebo pouze některou z těchto variant. Před výpočtem masek nebo map kvality může výsledek výrazně zlepšit filtrování wrapped obrazů fáze. Lze tak potlačit šum a další nepříznivé jevy. Kdybychom filtrovali wrapped obraz celý, došlo by k „rozmazání“ poloh 2π diskontinuit. Ovšem filtrování čitatele a jmenovatele argumentu funkce „arctan“ při výpočtu fáze potlačíme šum a poloha 2π diskontinuit se nezmění. Symbolicky tedy můžeme s užitím vztahu odvozeného v dřívější části psát
FILTR{Im[FT −1{C ∗ (f + f0 )}]} − 2πrf0 , −1 ∗ FILTR{Re[FT {C (f + f0 )}]}
φ (r ) = − arctan
(279)
kde FILTR značí filtrační algoritmus. Ghiglia a Pritt dělí algoritmy unwrappingu obecně do dvou skupin [68]:
•
Algoritmy sledující předem určenou cestu (PFA – z angl. Path-Following Algorithms),
•
Metody minimalizující normu (MNM – z angl. Minimum-Norm Methods).
PFA definují cestu unwrappingu pomocí QM nebo masky, pozice reziduí jsou pak definovány různými způsoby v závislosti na použitém algoritmu. MNM minimalizují rozdíl mezi vyhodnocovanou parciální derivací a parciální derivací řešení ve smyslu LP-normy M − 2 N −1
M −1 N − 2
J = ε P = ∑∑ φi +1, j − φi , j − ∆xi , j + ∑∑ φi , j +1 − φi , j − ∆yi , j , i =0 j = 0
P
P
(280)
i =0 j = 0
kde ϕij je hledaná hodnota fáze v bodě o souřadnicích (xi, yj), φi+1,j – φi,j je diference (diskrétní derivace) hledané skutečné fáze ve směru osy x, φi,j+1 – φi,j je diference ve směru osy y,
99
Unwrapping ∆xi , j = W {φi +1, j } − W {φi , j } je diference wrap fáze ve směru osy x, ∆yi , j = W {φi , j +1} − W {φi , j } ve směru osy y, W je wrap operátor, M a N jsou rozměry obrazu fáze ve směru x a y. Narozdíl od PFA definují MNM globální řešení na celém obrazu fáze pro redukování negativních vlivů. Možnou obtíž při vyhodnocení představuje ale záměna skutečných lokálních deformací za rezidua a přiřazení nevhodných vah, což může ovlivnit výslednou rekonstruovanou fázi. V odborné literatuře existuje velké množství algoritmů pro unwrapping. Jako příklad volme pro skupinu PFA:
•
Goldsteinův algoritmus (GA),
•
Quality-Guided Path Following algoritmus (QGPFA),
•
Mask-Cut algoritmus (MCA) ,
•
Flynnův algoritmus (FA).
Do skupiny MNM můžeme řadit např.:
•
algoritmy minimalizující sumu čtverců bez zavádění vah (ULS – z angl. Unweighted Least Squares),
•
vážený multi-grid (WMG – z angl. Weighted Multi-Grid),
•
metodu sdružených gradientů s předpodmíněním (PCG – z angl. Preconditioned Conjugate Gradient),
•
metodu minima LP-normy.
Z algoritmů, které kombinují výše jmenované metody, ale navíc používají další nástroje, jmenujme:
•
PUMA (z angl. Phase Unwrapping MAx-flow/min-cut – unwrapping na základě optimalizace grafů),
•
SRNCP (z angl. Sorting by Reliability, following a Non-Continuous Path),
•
Iso-phase (unwrapping na základě dekompozice obrazu),
•
unwrapping s využitím IIR (z angl. Infinite Impulse Response).
100
Unwrapping
6.2
Residua, mapy kvality a masky
Než přistoupíme k samotnému popisu algoritmů, nastíníme teorii residuí pro 2D unwrapping, a dále bude uveden přehled základních map kvality a masek používaných v algoritmech. Podrobný popis a analýza překračuje rámec této práce, čtenáře tedy odkážeme na příslušnou literaturu. Představme základní myšlenky [68]: •
Křivkový integrál I = ∫ F (r )dr = ∫ ( f dx + g dy ) C
(281)
C
je nezávislý na cestě integrace C, pokud je splněna některá z následujících podmínek: 1. F(r )dr = f dx + g dy je totální diferenciál, 2. F(r ) = ∇φ (r ) , čili F je gradient funkce ϕ, 3.
∫ F(r)dr = 0 , integrace po libovolné uzavřené křivce je nulová,
4. ∇ × F (r ) ≡ 0 , rotace F je nulová, pokud jsou si jednotlivé parciální derivace rovny. •
Tzv. teorém residuí pro dvojdimenzionální unwrapping zní:
∫ ∇φ (r) dr = 2πk ,
(282)
kde integrace probíhá „kolem residua“ a k je celé číslo11. Jinými slovy uzavřený křivkový integrál kolem residua je roven k-násobku 2π. Uvedené myšlenky, jmenovitě třetí podmínka – tj. integrace po uzavřené křivce je rovna nule, jsou základem ke všem PFA. Pokud je součet diferencí wrapované fáze po libovolné uzavřené křivce v rámci měřených dat nulový, jsou data tzv. konzistentní, a lze použít libovolnou integrační cestu – integrál nezávisí na zvolené cestě. Jsou-li v obraze přítomna residua, je unwrapping závislý na cestě integrace. Detekce residuí pro 2D unwrapping může probíhat s užitím výše uvedeného teorému (282). Vyberou se nejmenší možné uzavřené cesty (křivky), tj. 2 × 2 sousední pixely, a provede se integrace diferencí wrapované fáze po této elementární cestě (kolem vybraných 4 pixelů proti směru hodinových
11
Toto platí pouze pro reálná data – reálný obraz fáze. Pokud by byla data komplexní, nemusí být k jednoznačně celé číslo. Pro podrobné vysvětlení odkážeme čtenáře např. na [68].
101
Unwrapping ručiček). Je-li výsledek nulový, jsou data konzistentní. Bude-li ale roven ±2π, je detekováno tzv. kladné, resp. záporné, residuum. Hovoříme poté o tzv. polaritě residua. Na obr. 36 a obr. 37 je zobrazen vliv residuí na jednoduchý algoritmus popsaný v úvodní částí této kapitoly. Je na první pohled patrné, že při nepřítomnosti residuí není problém unwrapping provést, kdežto jsou-li nepříznivé jevy obsaženy v obraze rekonstruované fáze, algoritmus selhává a ve směru vyhodnocení dochází k závažným chybám. Za přítomnosti residuí je tedy třeba použít některý ze sofistikovanějších algoritmů popsaných v následující kapitole.
Obr. 36 Zobrazení a) simulované wrapped funkce „peaks“ bez residuí a šumu a b) její rekonstrukce jednoduchým algoritmem popsaným v úvodní části unwrappingu
Obr. 37 Zobrazení a) wrapped funkce „peaks“ s náhodně generovanými residui a b) její rekonstrukce jednoduchým algoritmem popsaným v úvodní části kapitoly
102
Unwrapping Mapy kvality (QM – z angl. Quality Maps) jsou nezbytné pro provádění sofistikovanějších algoritmů unwrappingu, které budou popsány v následujících částech. Jedná se o datová pole stejných rozměrů jako vyhodnocovaný obraz, kde jednotlivé prvky obsahují hodnocení kvality (váhu) jednotlivých pixelů. Vytváří se pro několik PFA, ale také pro některé algoritmy založené na určení vah LP-normy. Ghiglia a Pritt ve své knize [68] publikují 4 základní druhy map, a to: •
korelační mapy,
•
pseudokorelační mapy,
•
mapy variance fázové derivace,
•
mapy maximálního fázového gradientu.
První z nich, korelační mapa kvality, je získávána z IFSAR dat (z angl. Interferometric Synthetic Aperture Radar). Pro tuto práci nejsou takováto data k dispozici, tudíž se nemá smysl mapami zabývat. Algoritmy používané pro unwrapping generují mapy kvality jedním ze způsobů popsaných dále. Pseudokorelační mapy kvality mohou být generovány z libovolných fázových dat, čili pro naší úlohu jsou použitelné. V praxi ovšem nejsou dobrým hodnocením kvality [68]. Například pokud obraz obsahuje velké změny v hodnotách fáze, které ovšem nejsou residui nebo diskontinuitami, pseudokorelace označí tato místa malou vahou, což ale vede k nejistotě ve kvalitě vyhodnocení. Proto bylo definováno nové kritérium hodnocení kvality, a to mapy variance fázové derivace. Variance fázové derivace zmn je definována jako [68]
∑ (∆ k −1
z mn =
i , j =0
x i, j
− ∆xm ,n ) + 2
k
2
∑ (∆ k −1
i , j =0
y i, j
− ∆ym,n )
2
,
(283)
kde každá variance je počítána pro okno velikosti k × k centrované na pozici (m,n); ∆xi , j , ∆yi , j jsou parciální derivace obrazu ve směrech x a y, resp. diference wrapped fáze; ∆xm ,n , ∆ym ,n jsou průměrné hodnoty těchto derivací pro jednotlivá okna. Při generování map kvality jsou tyto mapy považovány za nejspolehlivější [68]. Mapy maximálního fázového gradientu mají podobnou nevýhodu jako pseudokorelační mapy,
tedy pro velké změny v hodnotách fáze, kdy bude gradient velký, přiřadí dané pozici menší váhu,
103
Unwrapping i když se může jednat o reálný stav fáze nikoli o residuum nebo diskontinuitu. Maximální fázový gradient je definovaný jako
max ∆xi , j , max ∆yi , j ,
(284)
kde ∆xi , j , ∆yi , j jsou obdobně jako v předchozím případě diskrétní parciální derivace obrazu v jednotlivých směrech x a y – fázové diference. Maxima jsou vyhodnocována v oblasti k × k kolem daného pixelu. Masky jsou speciálním případem map kvality, a to v tom smyslu, že obsahují pouze hodnoty 0 a
1. Nulou ohodnocené části obrazu představují místa s nízkou vahou a snažíme se jejich vliv na výsledek co nejvíce potlačit, případně úplně ignorovat. Použití masek se nachází u některých algoritmů minimalizujících LP-normu, ale i v PFA (např. Flynnův algoritmus), kde na ně způsoby unwrappingu kladou speciální požadavky. Tvorbu masek můžeme rozdělit do dvou skupin:
•
explicitní tvorba uživatelem,
•
generace masky z map kvality.
První případ nastává například v situaci, kdy je při pořizování obrazu interferogramu část plochy zacloněna pupilou (např. v případě kontroly optických komponent se tak stává velmi často, pupily mívají kruhový příp. eliptický tvar). Je-li tato poloha přesně známa, resp. ví-li uživatel při vyhodnocování, že daná data v určitých místech nemá smysl uvažovat, lze velmi snadno v softwarech definovat plochy, které mají být maskovány. Těm jsou poté přiřazeny váhy 0 a ostatním váhy 1. Generaci masek z map kvality lze provést volbou hodnot tolerance T, hraničních hodnot (z angl. threshold value), kdy při překročení této meze je přiřazena ploše váha 1 a naopak. Volba ovšem musí být provedena obezřetně, aby nedošlo k odmaskování nadbytečného počtu pixelů, nebo naopak nebyla zbytečně do výpočtu zahrnována residua a diskontinuity. Ghiglia a Pritt [68] představili algoritmus založený na tzv. re-mapování map kvality takovým způsobem, kdy histogram dané mapy získá tzv. U-tvar. Automaticky je pak volena hodnota poblíž minima histogramu re-mapované mapy kvality.
104
Unwrapping
6.3
Algoritmy sledující předem určenou cestu (PFA)
V této části budou podrobněji představeny jednotlivé použité PFA. Popis do maximální hloubky by ovšem překročil rámec práce, tudíž odkážeme na dostupnou odbornou literaturu.
6.3.1 Goldsteinův algoritmus (GA) Goldsteinův algoritmus (GA) generuje tzv. větvené řezy (BC – z angl. branch cuts), které spojují jednotlivá residua a vytváří tak oblasti, ve kterých následně probíhá unwrapping. Publikovali ho autoři Goldstein, Zebker a Werner [100] a je klasickým představitelem z rodiny PFA. Je velmi rychlý a obecně dosahuje poměrně dobrých výsledků [68]. Algoritmus je dělen do tří základních kroků: •
identifikace residuí,
•
generace hranic (BC),
•
integrace okolo hranic.
Algoritmus najde první residuum a označí ho jako „aktivní“ a „vybalancované“. Umístí na toto residuum vyhledávací okno o rozměru 3× 3 pixely a hledá, zda existuje v tomto okně další residuum. Není-li přítomno, zvětší se oblast vyhledávání na 5× 5 pixelů, hledá se znovu a případně opakuje zvětšení okna, dokud algoritmus nenajde další residuum nebo dokud nenarazí na okraj oblasti. Pokud algoritmus najde další residuum (a toto residuum není součástí jiné hranice, tj. není označeno jako „vybalancované“), přičte jeho polaritu (hodnotu +1 nebo –1) k polaritě aktivního residua (pokud je „vybalancované“, přičte se nula), nově nalezené reziduum se označí za „vybalancované“ a obě rezidua se propojí pomocí hranice. Hranice se dále rozrůstá tím způsobem, že je vyhledávací okno posunuto do nové polohy (tj. vycentruje se na poloze nově nalezeného residua a toto residuum se označí jako „aktivní“) a opět se hledá další residuum. Pokud je nalezeno (a není dosud „vybalancované“), přičte se jeho polarita k celkovému „náboji“ (součtu polarit) dané cesty a obdobným způsobem se pokračuje dále. Hranice (branch cut) je dokončena, pokud je výsledný náboj cesty rovný nule, nebo pokud se algoritmus dostane k okraji oblasti. V tom případě je poslední residuum propojeno s tímto okrajem. Po dokončení se přejde k dalšímu „nevybalancovanému“ residuu a začne se tvořit další hranice, dokud algoritmus neprojde všechny dosud nevybalancovaná residua.
105
Unwrapping Ghiglia a Pritt publikují také vylepšení Goldsteinova algoritmu tvorbou tzv. dipole cuts [68] (aplikací jednoho z Huntleyho algoritmů), kdy jsou nejprve propojena residua vzájemně opačné polarity, která leží velmi blízko vedle sebe. Na obr. 38 je schematicky znázorněno vhodné i nevhodné propojení dipólů. Nemá-li některé residuum svůj protějšek, je propojeno nejkratší cestou s hranicí oblasti. Následně je aplikován klasický Goldsteinův algoritmus popsaný výše.
a)
b)
c)
Obr. 38 a) Rozložení residuí v obraze (černé – záporná polarita, bílé – kladná polarita), b) nevhodné spojení dipólů, c) vhodné propojení dipólů [68]
6.3.2 Quality-Guided Path Following algoritmus (QGPFA) Quality-guided path following algoritmus (QGPFA) se dívá na problém unwrappingu z jiného pohledu. Základní myšlenkou, na které je tato metoda postavena, je, zda-li existuje nějaká informace ve fázových datech (mimo polohy residuí), která by mohla být použita k vedení integrační cesty. Jinými slovy – cílem je získat takové informace o kvalitě jednotlivých obrazových bodů, které posléze pomohou definovat „nejkvalitnější“ cestu, po které bude probíhat integrace – unwrapping. Základem je definice mapy kvality a začátek výpočtu u pixelu, jehož váha je největší. Residua tedy nejsou konkrétně identifikována, ale využívá se toho, že residua se obvykle nacházejí v oblastech méně kvalitních dat. Mapy kvality mohou být různé v závislosti na oblasti použití algoritmu. QGPFA je velmi používaný zejména pro svojí robustnost a rychlost provedení unwrappingu. Prvním, kdo pravděpodobně použil kvalitativní hodnocení obrazu, byl Bone [101]. Fázová data vyhodnocoval výpočtem druhé derivace (resp. druhé diference, uvážíme-li diskrétní povahu digitálních dat) a volbou toleranční hodnoty, jež utvářela hranici pro definování vah vznikající masky. Tímto postupem byly unwrappovány pouze ty pixely, jejichž hodnota druhé diference
106
Unwrapping nepřekročila hodnotu tolerance. Následovala celá řada prací různých autorů, které rozšiřovaly myšlenku maskování a hodnocení pomocí kvalitativních vlastností obrazu fáze. Ghiglia a Pritt [68] prezentují algoritmus, který využívá různé druhy map kvality. Zmínka o nich byla v předchozí části, čili předpokládejme, že taková mapa buď existuje, nebo ji lze generovat ze zaznamenaného obrazu fáze, viz např. mapy variance fázové derivace. Princip algoritmu je takový, že vhodným postupem jsou nejdříve vyhodnocovány oblasti s největší vahou danou mapou kvality, až jako poslední přijdou na řadu „nejméně kvalitní“ části obrazu. Kvalita vyhodnocení tak ale závisí na vhodném výběru a generování kvalitativní mapy.
6.3.3 Mask-Cut algoritmus (MCA) Mask-Cut algoritmus (MCA) je vhodnou kombinací dvou předchozích algoritmů – mapy kvality jsou používány pro nalezení polohy hranic unwrappingu, větvených řezů, v obrazu fáze. Jako první přišli s myšlenkou definování řezů pomocí kvalitativních údajů Prati, Giani a Leuratti [102], ale jejich algoritmus ještě nedokázal nalézt polohu pro nejméně kvalitní pixely. Dle [68] byl prvním, kdo publikoval detailní popis takového algoritmu definování řezů, Flynn. Princip spočívá v nalezení residua, a poté nechání na algoritmu samém, aby v souladu s kvalitativními údaji umístil hranice pro unwrapping. Struktura algoritmu je obdobná QGPFA, ovšem unwrapping začíná z „opačné strany“, respektive v regionech s nejmenší vahou. Maska poté narůstá (rozšiřuje se), dokud není polarita residuí balancována, nulová, nebo dokud maska nedosáhne okraje obrazu. Ve zkratce můžeme kroky algoritmu zapsat jako: •
identifikace residuí,
•
generace řezů pomocí mapy kvality (MC – z angl. mask cuts),
•
tenčení řezů,
•
integrace okolo řezů.
Třetí bod (tenčení řezů) je nutné provádět, jelikož při „rozšiřování“ masky dochází k zesilování, rozšiřování, hranic. Přichází tedy na řadu morfologická operace [47, 48], která danou úlohu provede. MCA dává obdobné výsledky jako QGPFA [68], pokud máme k dispozici dostatečné kvalitativní údaje – mapy kvality. Jestliže nejsou, nebo je z nějakého důvodu nejsme schopni generovat, potom Goldsteinův algoritmus, co se týče přesnosti vyhodnocení, dosahuje výsledků lepších. 107
Unwrapping
6.3.4 Flynnův algoritmus (FA) Flynnův algoritmus (FA) vhodným způsobem využívá myšlenku, která je patrná z obr. 39. Hranice mezi nejsvětlejšími a nejtmavšími místy (místa nespojitosti) nám definují oblasti, ke kterým musíme přičíst nebo odečíst celočíselný násobek 2π, abychom provedli unwrapping. Z jiného pohledu můžeme tedy říci, že pro unwrapping musíme minimalizovat velikost nespojitosti mezi jednotlivými regiony.
Obr. 39 Zobrazení a) simulované fáze, b) fáze ve wrapped formě a c) polohy nespojitostí v obraze (detekované hrany)
Rozdělení wrap fáze do regionů a jejich unwrappování se ovšem komplikuje u reálných dat. Přítomnost šumu nedovoluje jednoznačně definovat hranice. Také residua způsobí přerušení linií nespojitostí. Oba dva tyto aspekty je nutné řešit, na první pohled by se mohlo zdát, že je to problém až neřešitelný. Flynn ovšem publikoval práci [103], kde předvedl algoritmus minimalizující váženou sumu diskontinuit. Iterativně zde hledá minimum L1-normy. Na rozdíl od výše představené myšlenky se ale neomezuje na definování hranic oblastí pouze v pozici zlomů. Volí však tyto hranice oblastí a přičítá regionům vhodný násobek 2π, pokud toto přičtení nezpůsobí vznik více diskontinuit, než jich bylo předtím. Proces je iterativně opakován a redukuje počet zlomů v každém kroku. Pokud není nalezeno žádné další možné dělení, je algoritmus ukončen a fáze považována za zbavenou nespojitostí. Takto vzniká rekonstruovaný obraz ve smyslu minimálního součtu diskontinuit, jinými slovy zbývající diskontinuity jsou minimální v souladu s příslušnou mírou (normou).
108
Unwrapping
6.4
Metody minimalizující normu (MNM)
MNM se dívají na problematiku unwrappingu z jiného úhlu pohledu. Zatímco PFA vytvářely hranice, které bránily cestám integrace, MNM zavádí matematicky formální hledisko. Problém odstraňování nespojitostí fáze je řešen ve smyslu minimalizace LP-normy M −2 N −1
M −1 N − 2
J = ε P = ∑∑ φi+1, j − φi , j − ∆xi , j + ∑∑ φi , j +1 − φi , j − ∆yi , j , i =0
P
j =0
i =0
P
(285)
j =0
kde ϕij je hledaná hodnota fáze v bodě o souřadnicích (xi, yj), φi+1,j – φi,j je diference (diskrétní derivace) hledané skutečné fáze ve směru osy x, φi,j+1 – φi,j je diference ve směru osy y, ∆xi , j = W {φi +1, j } − W {φi , j } je diference wrap fáze ve směru osy x, ∆yi , j = W {φi , j +1} − W {φi , j } ve směru osy y, W je wrap operátor, M a N jsou rozměry obrazu fáze ve směru x a y. Řešení je tak hledáno v globálním měřítku. Lze také definovat váhy, potom bychom mluvili o minimalizaci vážené LPnormy. Na minimalizační problém (285) lze nahlížet dvěma způsoby, a to: •
spojitým variačním přístupem, nebo
•
diskrétním variačním přístupem.
Pro první případ, spojitý variační přístup, můžeme rovnici (285) přepsat ve smyslu [68]
J = ∫∫ f (φ x , φ y , x, y ) dxdy ,
(286)
f = φ x −ψ x + φ y −ψ y ,
(287)
kde funkce f je definována jako P
P
kde ϕx, ϕy jsou diference (diskrétní derivace) hledané skutečné fáze ve směru x a y, ψx, ψy jsou jim odpovídající diference wrapped fáze ve směrech x a y. Fázové diference jsou následně nahrazeny spojitými parciálními derivacemi. Poté je použit variační počet pro nalezení takových podmínek, při kterých je J konstantní. Funkce ϕ, která minimalizuje řešení (286), musí splňovat EulerLagrangeovu rovnici [68]
∂ ∂f ∂ ∂f + = 0. ∂x ∂φx ∂y ∂φ y
(288)
Výpočtem parciálních derivací funkce f dostáváme
109
Unwrapping P −2 ∂f = p (φx − ψ x )φx − ψ x , ∂φx P −2 ∂f = p (φ y − ψ y )φ y − ψ y , ∂φ y
(289)
protože platí ∂/∂z(|z|P) = ∂/∂z[(z2)P/2] = pz|z|P-2. Abychom dostali řešení pro ϕ ve smysly LPnormy, musíme řešit parciální diferenciální rovnici, která vznikne po dosazení (289) do (288), ∂ [U ( x, y)(φx −ψ x )] + ∂ [V ( x, y )(φy −ψ y )] = 0 , ∂x ∂y
(290)
kde U(x,y) = |ϕx – ψx|P-2 a V(x,y) = |ϕy – ψy|P-2 mohou být uvažovány jako zobecněné váhy fázových derivací závislé na vstupních datech (nejsou to tedy váhy nezávisle volené). Pro případ P = 2 (případ metody nejmenších čtverců, minimalizuje se J = ε2), jsou váhy V a U jednotkové, potom dostáváme tzv. Poissonovu diferenciální rovnici [68]
φxx + φ yy = ψ xx + ψ yy ,
(291)
∇ 2φ = ρ ,
(292)
nebo
kde ∇ 2 je Laplaceův operátor a ρ = ψxx + ψyy, ψii je ve spojitém případě druhá derivace wrapped fáze (v diskrétním případě druhá diference wrapped fáze) ve směru i, ϕii je ve spojitém případě druhá derivace hledané skutečné fáze (v diskrétním případě druhá diference hledané skutečné fáze) ve směru i. Případ, kdy váhy jsou jednotkové, tedy případ řešení ve smyslu nejmenších čtverců, nás vede k řešení Poissonovy rovnice (292). Pokud by váhy jednotkové nebyly, musíme řešit nelineární parciální diferenciální rovnici (290). Lze ukázat [68], že v případě diskrétního variačního přístupu musíme obecně řešit problém
(φ
i +1, j
− φ i , j − ∆xi , j )U (i, j ) + (φ i , j +1 − φ i , j − ∆ yi , j )V (i, j ) −
− (φ i , j − φ i −1, j − ∆xi −1, j )U (i − 1, j ) − (φ i , j − φ i , j −1 − ∆ yi , j −1 )V (i, j − 1) = 0 ,
(293)
kde
110
Unwrapping U (i, j ) = φi +1, j − φi , j − ∆xi , j
P −2
U (i, j ) = 0
⇔ i = 0, K , M − 2, ⇔
V (i, j ) = φi , j +1 − φi , j − ∆
y i, j
P −2
V (i, j ) = 0
jinak,
⇔ i = 0, K , M − 1, ⇔
j = 0, K , N − 1, j = 0, K , N − 2,
(294)
jinak.
Pro P = 2 poté dostáváme diskrétní formu Poissonovy rovnice
(φ
i +1, j
− 2φi , j + φi−1, j ) + (φi , j +1 − 2φi , j − φi , j −1 ) = ρi , j ,
(295)
kde
ρ i , j = (∆xi , j − ∆xi−1, j ) + (∆yi , j − ∆yi , j−1 ) .
(296)
Algoritmy unwrappingu řeší právě výše uvedené diferenciální rovnice a problémy. Pro podrobnější popis konkrétních postupů odkážeme na příslušnou literaturu.
6.4.1 Algoritmy nejmenších čtverců bez vah (ULS) Využití principu metody nejmenších čtverců je výhodné zejména proto, že vede na soustavu lineárních rovnic, pro které je vyvinuto několik vhodných metod řešení. Nevýhodou však je, že může dojít k „vyhlazování“ rekonstruovaných ploch, čili můžeme získat nevhodné, nedostatečně přesné, výsledky. Ovšem v aplikacích, jako např. měření topografie optických ploch, kde chceme potlačit vliv šumu a residuí na výsledek, mohou najít tyto algoritmy uplatnění. V [68] může čtenář najít několik způsobů řešení ULS, zde budou představeny v krátkosti dva algoritmy, které jsou analyzovány a použity v další části práce. Jedná se o: •
algoritmus založený na rychlé Fourierově transformaci (FFT – z angl. Fast Fourier
Transform), •
multi-grid algoritmus (MG).
První z algoritmů využívá toho, že Poissonova rovnice (292) má až na aditivní konstantu jednoznačné řešení a nepotřebuje hraniční podmínky, pokud je funkce ϕ periodická [68]. K dosažení periodicity wrapped obrazu fáze je použito zrcadlení podle jednotlivých okrajů. Obr. 40 schematicky způsob zrcadlení naznačuje.
111
Unwrapping
Obr. 40 J. B. J. Fourier (1768-1830) jako funkce fi definovaná na 1 ≤ i ≤ 500 a 1 ≤ j ≤ 612 zrcadlená podle osy x = 500 a y = 612
Takto vzniklou periodickou wrap funkci značme ψ~i , j a jí odpovídající unwrappovanou funkci ~
φi , j . Fázové diference získáme jako ∆xi , j = W {ψ~i +1, j − ψ~i , j } , ∆yi , j = W {ψ~i , j +1 − ψ~i , j } ,
(297)
kde W je wrap operátor. Poté diskrétní Poissonova rovnice (295) přejde na tvar
(φ~
i +1, j
) (
)
~ ~ ~ ~ ~ − 2φi , j + φi −1, j + φi , j +1 − 2φi , j − φi , j −1 = ρ~i , j .
(298)
Její pravá strana bude dána vztahem (296), kde za fázové diference dosazujeme z (297). Aplikací dvojdimenzionální Fourierovy transformace na (298) dostaneme [68]
Φ m ,n =
Ρm ,n , 2 cos(πm / M ) + 2 cos(πn / N ) − 4
(299)
kde M a N značí rozměry wrapped periodického obrazu, m a n pozici v tomto obraze, Φm,n a Ρm,n ~ ~ dvojdimenzionální Fourierovu transformaci φi , j a ρ~i , j . Řešení φi , j je tedy získáváno inverzní Fourierovou transformací rovnice (299), neperiodický unwrapovaný obraz fáze φi, j následně uvážením jen těch bodů, které odpovídají pozici v neperiodickém obraze.
112
Unwrapping Druhý z algoritmů, tzv. multi-grid12 (MG), transformuje problém do mřížek (sítí) různých rozlišení, od nejhrubších k nejjemnějším, a hledá pak nejlepší řešení vyhovující podmínce nejmenších čtverců. Obecně MG metody jsou technikami pro poměrně rychlé řešení parciálních diferenciálních rovnic (PDR) na velkých sítích. Jsou založeny na myšlence aplikace Gauss-Seidelovy iterační metody [104] na hrubší (menší) mřížky. Oproti způsobům řešení PDR na základě diskrétní Fourierovy transformace, kterým se vyrovnávají co se týče rychlosti, ovšem umožňují řešit mnohem širší pole problémů, jako i nelineární parciální diferenciální rovnice.
6.4.2 Vážený multi-grid (WMG) Metoda váženého multi-gridu (WMG – z angl. Weighted Multi-Grid) je rozvinutá metoda MG, která používá mapy kvality k určení vah, kde „nejhorší“ pixely obrazu fáze jsou váženy hodnotou 0, naopak „lepší“ jsou váženy vyššími hodnotami. Jedná se tedy také o metodu nejmenších čtverců, kde ovšem zavedení vah umožňuje potlačit vliv residuí a dalších nepříznivých jevů. Hodnota, kterou algoritmus minimalizuje, je s přihlédnutím k (293) a (294) [68]
ε 2 = ∑ U (i, j )(φi+1, j − φi , j − ∆xi , j ) + ∑ V (i, j )(φi , j +1 − φi , j − ∆yi , j ) , 2
i, j
2
(300)
i, j
a váhy gradientu, U(i,j) a V(i,j), jsou definovány jako
U (i, j ) = min(wi2+1, j , wi2, j ) , V (i, j ) = min(wi2, j +1 , wi2, j ) ,
(301)
kde wi,j je hodnota mapy kvality na pozici (i, j), 0 ≤ wi,j ≤ 1. Vliv residuí tím pádem může být potlačen úplně zavedením váhy 0, pokud je ovšem jejich poloha detekována správně. Výpočet fáze ve smyslu nejmenších čtverců je následně definován rovnicí [68]
U (i, j )(φ i +1, j − φ i , j ) + U (i − 1, j )(φ i , j − φ i −1, j ) +
+ V (i, j )(φ i , j +1 − φ i , j ) − V (i, j − 1)(φ i , j − φ i , j −1 ) = c i , j ,
(302)
ci , j = U (i, j )∆xi , j − U (i − 1, j )∆xi −1, j + V (i, j )∆yi , j − V (i, j − 1)∆yi , j −1 .
(303)
kde
Metoda nejmenších čtverců s váhami je řešena iterativně, MG využívá Gauss-Seidelovu metodu.
12
z angl. grid – mřížka, mříž
113
Unwrapping
6.4.3 Metoda sdružených gradientů s předpodmíněním (PCG) Metody sdružených gradientů jsou numerické metody iterativně řešící soustavy lineárních rovnic [105], ke kterým patří také diskretizovaná Poissonova rovnice (295). Tento způsob využívá optimalizační (minimalizační) techniky a dosahuje mnohem rychlejší konvergence než např. metody největšího spádu (z angl. steepest descent) [68]. Pro určení odhadu vstupních parametrů („předpodmínění“ sdruženého gradientu) je použito nevážené metody nejmenších čtverců. Další iterace optimalizačního výpočtu (sdružených gradientů) už zavádějí váhy, které jsou určeny pomocí mapy kvality (generované podle metod pseudokorelace, variance fázových derivací nebo maximálního fázového gradientu). Důležitým krokem při použití tohoto algoritmu, ale i u algoritmů MG a všech ostatních využívajících metodu nejmenších čtverců, je po unwrappování vytvořit takové řešení, které je tzv. kongruentní – odpovídá řešení před unwrappingem. Jinak řečeno, chceme získat takové řešení, na které když použijeme wrap operátor, dostaneme stejné (odpovídající) hodnoty, jaké byly vstupem do unwrapping algoritmu. Nejjednodušší řešení může být následující. Mějme unwrappovanou plochu φ(x,y) a jí odpovídající wrap plochu ψ(x,y). Potom operaci „kongruence“ definujme jako [68]
)
ϕ ( x, y ) = ϕ ( x, y ) + W {ψ ( x, y ) − ϕ ( x, y )} ,
(304)
kde W je wrap operátor. Tento postup získání kongruentního řešení je snadný, avšak při vyhodnocování zašuměných dat nebo obrazů obsahujících residua může selhávat. Ghiglia a Pritt představují varianty, jak daný problém řešit [68].
6.4.4 Metoda minima LP-normy Minimalizace LP-normy diferencí gradientu fáze před a po odstranění nespojitostí (unwrappingu) vede na řešení nelineární parciální diferenciální rovnice, viz (290) a (293), jelikož váhy U a V jsou závislé na vstupních datech a na řešení samotném. Obecně je možné tyto rovnice řešit iterativně, jak bylo nastíněno výše, a jak je podrobně popsáno v [68]. Z rovnic (290) a (293) lze určit vlastnosti algoritmů pro volby jednotlivých P. Pro velké hodnoty, P >> 2, budou gradienty výsledku lišící se od gradientů vstupních dat hodnoceny malými váhami. Tudíž dojde k velkému vyhlazování ploch a residua značně ovlivní globální řešení, což může být v některých aplikacích nežádoucí. Výsledné gradienty tudíž nebudou shodné se vstupními na žádných místech rekonstruovaných dat. 114
Unwrapping Naopak pro P < 2 budou váhy velké tam, kde se gradienty rekonstruované fáze blíží vstupním. Pro P = 0 dostáváme s (285) řešení, kde J odpovídá sumě počtu vzorků, ve kterých gradienty řešení neodpovídají přesně vstupním. Jinými slovy použití L0-normy vede k výsledným gradientům totožným se vstupními ve všech místech, kde to je možné (vyjma residuí a diskontinuit). Porovnáme-li metody minimalizující L0-normu s Goldsteinovým algoritmem představeným výše, vidíme, že oba dva zachovávají gradienty fáze před a po odstranění nespojitostí. GA ovšem nezachová gradienty v místech přechodu řezů (branch cuts), kdežto globální pojetí optimalizačních algoritmů dá „shodné“ gradienty i v těchto místech. Zdá se tedy, že hodnota P určuje míru vyhlazení výsledku. Velké vyhlazení ale nepřináší výhodu, v praxi chceme znát skutečné gradienty, čehož dosáhneme volbou P < 2. Jak ukazuje [68], je vhodné volit P = 0. Jednak z vhodných vlastností zachování gradientu, a také protože P = 1 nepřináší žádnou další výhodu.
115
Unwrapping
6.5
Další (kombinované) algoritmy
V praxi je vyvíjeno a používáno velké množství dalších algoritmů, jejichž principy jsou založeny např. na teorii grafů, dekompozici obrazu a dalších možnostech, které nejlépe vyhovují aplikaci, pro níž jsou algoritmy sestavovány. Podrobný rozbor a popis všech jednotlivých metod je tedy v podstatě nemožný z hlediska rozsahu této práce. Čtenáře tak odkážeme na dostupnou odbornou literaturu.
Základ
kombinovaných
algoritmů
ale
zpravidla
tvoří
principy
popsané
v předcházejících částech. V práci budou používány tyto kombinované algoritmy: •
PUMA (z angl. Phase Unwrapping MAx-flow/min-cut), algoritmus založený na optimalizaci grafů pomocí max-flow/min-cut kalkulací [95]. Funkce, které mají být minimalizovány (meritní funkce) jsou Markovova náhodná pole prvního stupně. Podrobný popis těchto polí lze nalézt v [106]. PUMA je schopen díky svým vlastnostem exaktně minimalizovat LP-normu pro p ≥ 1.
•
SRNCP (z angl. Sorting by Reliability, following a Non-Continuous Path), algoritmus z rodiny QGPFA, který unwrapuje pixely v pořadí závislém na jejich kvalitě [96]. Cesta, po které probíhá integrace, je zde nespojitá.
•
Iso-phase algoritmus, algoritmus založený na dekompozici obrazu [97], provádí unwrapping po oblastech, které mají podobnou hodnotu fáze – po homogenních částech. Algoritmus tak velmi rychle a kvalitně reaguje na šum a diskontinuity ve vyhodnocovaném obraze.
•
Unwrapping s využitím IIR, algoritmus založený na využití IIR (z angl. Infinite Impulse Response) prochází vyhodnocovanou oblast ve směru fázových isolinií [98]. S použitím lokálních oken pak vylepšuje citlivost k šumu.
116
7 Analýza algoritmů unwrappingu V předchozí části bylo představeno několik algoritmů pro odstranění nespojitostí při rekonstrukci fázového rozdílu. Tato kapitola má za cíl provést srovnání jednotlivých metod, a nabídnout tak přehled vlastností pro ideální volbu v jednotlivých aplikacích. Nejprve se budeme zabývat analýzou citlivosti jednotlivých algoritmů na náhodný šum přítomný v nespojitém obraze fázového rozdílu interferujících vln získaného při vyhodnocení pomocí FTM. Následně se zaměříme na chování algoritmů unwrappingu u takových dat, jež obsahují residua. Pro některé praktické aplikace je důležité takovouto situaci prověřit.
7.1
Zdroje použitých algoritmů
Pro analýzu chování algoritmů unwrappingu, které publikovali Ghiglia a Pritt [68], byly použity zdrojové kódy napsané v jazyce C dostupné online na ftp serveru vydavatelství Wiley •
ftp://ftp.wiley.com/public/sci_tech_med
v adresáři phase_unwrapping. Dále byly získány kódy pro algoritmus PUMA dostupné na •
http://www.lx.it.pt/~bioucas/code.htm,
kódy algoritmů SRNCP a Iso-phase získané online na •
http://www.ljmu.ac.uk/GERI/90225.htm,
a pro algoritmus využívající systému IIR zdrojové kódy dostupné na •
http://www.cio.mx/~jestrada/phase_unwrapping2.html.
Z jazyka C nebo C++ byly výše uvedené kódy upraveny pro použití v prostředí MATLAB pomocí MEX funkcí [107]. Využitím MEX funkcí lze elegantním způsobem implementovat kódy napsané v jazyce C nebo C++ do prostředí MATLAB. Pro úspěšné použití tohoto postupu je předpokládána alespoň základní znalost jazyka C/C++. Pro důkladnější studium odkážeme čtenáře na dostupné zdroje, např. [107, 108].
117
Analýza algoritmů unwrappingu
7.2
Vliv náhodného šumu na algoritmy unwrappingu
Pro analýzu vlivu náhodného šumu na algoritmy unwrappingu uvažujme vztah pro wrapped fázový rozdíl Δϕw(r) dvou interferujících vln ∆φw (r ) = W {∆φ (r )} + ν (r ) ,
(305)
kde W je wrapped operátor a ν(r) je aditivní náhodný šum, pro nějž předpokládejme normální rozdělení o střední hodnotě μ a směrodatné odchylce σ a r = (x, y) určuje polohu v obraze. Pro hodnocení reakcí algoritmů na náhodný šum zaveďme dvě kritéria, a to: •
míru šumu v rekonstruovaném obraze definovanou vztahem
∑ (∆φ M
N=
i =1
− ∆φui )
2
i
M
,
(306)
kde M je počet obrazových bodů interferogramu, Δϕi je hodnota skutečného fázového rozdílu interferujících vln se zavedeným šumem ν(r) v i-tém bodě a Δϕui je hodnota fázového rozdílu po provedeném unwrappingu, a dále •
dobu t provádění unwrappingu.
Multiplikací výše uvedených kritérií získáme testovací hodnotu pro jednotlivé algoritmy T = Nt .
(307)
Pomocí představených kritérií a testovací hodnoty můžeme dále vyvozovat závěry o jednotlivých metodách odstranění nespojitostí v rekonstruovaném obraze. Na obr. 41 je zobrazena testovací funkce se zavedeným náhodným šumem, pro který uvažujeme normální rozdělení se střední hodnotou μ = 0 a směrodatnou odchylkou σ = π/4, a dále wrap podoba této funkce a rekonstrukce pomocí jednoduchého algoritmu, který můžeme popsat v následujících bodech: •
čtvrcení obrazu,
•
tvorba korekční funkce Δϕc(r) pro každou čtvrtinu způsobem popsaným v úvodní části předchozí kapitoly, výchozí bod odpovídá středu původního obrazu,
•
přičtení korekční funkce k původním hodnotám obrazu,
118
Analýza algoritmů unwrappingu •
„sešití“ obrazu s odstraněnými nespojitostmi pomocí mediánu rozdílu vzájemně sousedících řad pixelů jednotlivých čtvrtin.
Obr. 41 a) Testovaná funkce, b) její wrap podoba a c) rekonstrukce pomocí jednoduchého algoritmu
Jak je patrné, výše popsaný algoritmus si s takovýmto případem neporadil. Pro data obsahující šum je tedy nutné volit sofistikovanější metody. Algoritmy získané z výše uvedených zdrojů byly upraveny pro použití v prostředí MATLAB. Po této modifikaci je možné zadávat u jednotlivých metod různé volitelné parametry. Výsledky rekonstrukce nespojitého obrazu závisí na jednotlivých parametrech, a proto byla provedena analýza jejich vlivu na míru šumu definovanou vztahem (306). Tab. 1 představuje souhrn volitelných parametrů jednotlivých algoritmů. Podrobné vysvětlení jednotlivých parametrů může čtenář nalézt na výše uvedených adresách v části 7.1.
Tab. 1
Volitelné parametry algoritmů unwrappingu
Algoritmus Parametr(y)
Popis
Goldsteinův algoritmus (GA) C
Maximální povolená délka řezu (cut-length)
D
Eliminace „dipól“ residuí před unwrappingem (1 pro „ano“, 0 pro „ne“)
Quality-Guided Path Following algoritmus (QGPFA) Mask-Cut algoritmus (MCA) MK
Typ generované mapy kvality (1 pro „minimální gradient“, 2 pro „minimální variance“, 3 pro „maximální pseudokorelace“)
W
Velikost průměrujícího okna použitého na mapu kvality
119
Analýza algoritmů unwrappingu Flynnův algoritmus (FA) MK
Typ generované mapy kvality (1 pro „minimální gradient“, 2 pro „minimální variance“, 3 pro „maximální pseudokorelace“)
W
Velikost průměrujícího okna použitého na mapu kvality
TH
Automatické prahování – tj. automatické stanovení prahové hodnoty a binarizace mapy kvality (1 pro „ano“, 0 pro „ne“)
F
Počet pixelů, o které může být mapa kvality tenčena pomocí morfologických operací digitální analýzy obrazu
ULS algoritmus pomocí rychlé Fourierovy transformace (ULS FT) ---
---
ULS multi-grid algoritmus (ULS MG) I
Počet iterací Gauss-Seidelovy metody
C
Počet cyklů
Vážený multi-grid (WMG) MK
Typ generované mapy kvality (1 pro „minimální gradient“, 2 pro „minimální variance“, 3 pro „maximální pseudokorelace“)
W
Velikost průměrujícího okna použitého na mapu kvality
TH
Automatické prahování – tj. automatické stanovení prahové hodnoty a binarizace mapy kvality (1 pro „ano“, 0 pro „ne“)
F
Počet pixelů, o které může být mapa kvality tenčena pomocí morfologických operací digitální analýzy obrazu
I
Počet iterací Gauss-Seidelovy metody
C
Počet cyklů
Metoda sdružených gradientů s předpodmíněním (PCG) MK
Typ generované mapy kvality (1 pro „minimální gradient“, 2 pro „minimální variance“, 3 pro „maximální pseudokorelace“)
W
Velikost průměrujícího okna použitého na mapu kvality
TH
Automatické prahování – tj. automatické stanovení prahové hodnoty a binarizace mapy kvality (1 pro „ano“, 0 pro „ne“)
F
Počet pixelů, o které může být mapa kvality tenčena pomocí morfologických operací digitální analýzy obrazu
I
Počet iterací
E
Tolerance konvergence
Metoda minima LP-normy MK
Typ generované mapy kvality (1 pro „minimální gradient“, 2 pro „minimální variance“, 3 pro
120
Analýza algoritmů unwrappingu „maximální pseudokorelace“) W
Velikost průměrujícího okna použitého na mapu kvality
TH
Automatické prahování – tj. automatické stanovení prahové hodnoty a binarizace mapy kvality (1 pro „ano“, 0 pro „ne“)
F
Počet pixelů, o které může být mapa kvality tenčena pomocí morfologických operací digitální analýzy obrazu
I
Počet iterací
IPCG
Počet PCG iterací
N
Normalizační parametr
Algoritmus PUMA P
Exponent ve vztahu pro potenciál
Algoritmus SRNCP ---
---
Iso-phase algoritmus NG
Počet skupin
Algoritmus s využitím IIR τ
Hranice omezení IIR systému
σ
Vyhlazovací parametr
W
Velikost okna
Analýza, jež ukázala vliv výše uvedených parametrů na míru šumu, byla provedena v prostředí MATLAB formou simulací. Obr. 42, obr. 43 a obr. 44 zobrazují vliv volitelných parametrů na míru šumu při použití vybraných algoritmů.
Obr. 42 Vliv volitelných parametrů QGPFA algoritmu na míru šumu rekonstruovaného obrazu (modrá – MK minimálního gradientu, zelená – MK minimální variance, červená – MK maximální pseudokorelace)
121
Analýza algoritmů unwrappingu
Obr. 43 Vliv volitelných parametrů MC algoritmu na míru šumu rekonstruovaného obrazu (modrá – MK minimálního gradientu, zelená – MK minimální variance, červená – MK maximální pseudokorelace)
Obr. 44 Vliv volitelných parametrů Flynnova algoritmu na míru šumu rekonstruovaného obrazu (světlá místa znační malou hodnotu míry šumu)
Pomocí výše uvedené analýzy byly voleny takové parametry, pro které míra šumu dosahovala nejmenších hodnot. Tab. 2 porovnává míry šumu jednotlivých algoritmů s původním obrazem a relativní čas průběhu výpočtu vztažený ke Goldsteinově algoritmu. Na obr. 45, obr. 46 a obr. 47 jsou následně zobrazeny rekonstruované obrazy fázových rozdílů. Pro poslední čtyři algoritmy nebyly použity zdrojové kódy napsané ve stejném „stylu“ jako metody předchozí. Není zaručena jejich optimalizace stejného charakteru jako pro ty, které popisují Ghiglia a Pritt v [68]. Musíme tedy jejich čas výpočtu brát s dostatečnou rezervou.
122
Analýza algoritmů unwrappingu Tab. 2 Míra šumu rekonstruovaného obrazu a její relativní chyba vzhledem k původnímu obrazu a relativní čas výpočtu vztažený ke Goldsteinově algoritmu Algoritmus
N [rad]
trel [s]
0.79
---
GA
0.90 (15%)
1
QGPFA
1.33 (69%)
6.3
MCA
1.00 (28%)
5.8
FA
0.91 (16%)
23.8
ULS FT
1.69 (116%)
3.3
ULS MG
1.88 (140%)
1.82
WMG
0.89 (13%)
4.9
PCG
1.27 (62%)
22.4
Minimalizace LP-normy
0.79 (<1%)
71.2
PUMA
0.79 (<1%)
264.2
SRNCP
6.36 (712%)
1.5
Iso-phase
0.83 (6%)
6.5
IIR
1.22 (56%)
242.2
původní obraz
Obr. 45 Rekonstruovaný fázový rozdíl pomocí algoritmů sledujících předem definovanou cestu – a) GA, b) QGPFA, c) MCA, d) FA
123
Analýza algoritmů unwrappingu
Obr. 46 Rekonstruovaný fázový rozdíl pomocí algoritmů minimalizujících normu – a) ULS MG, b) PCG, c) WMG, d) Minimalizace LP-normy
Obr. 47 Rekonstruovaný fázový rozdíl pomocí kombinovaných algoritmů – a) PUMA, b) SRNCP, c) Iso-phase, d) IIR
124
Analýza algoritmů unwrappingu Z výše uvedené analýzy je patrné porovnání jednotlivých algoritmů co se týče přesnosti vyhodnocení velmi zašuměných dat a jejich relativní doby výpočtu. Pro praktické aplikace je třeba volit z představných algoritmů s ohledem na míru zašumění reálných dat. Simulovaná situace ukazovala extrémní případ, který v současné praxi testování optických ploch prakticky nenastává. Při použití Fourierovy transformace můžeme ve frekvenční oblasti navíc provádět filtrace, kterými dokážeme šum potlačit, a wrap obraz má poté příznivější charakter.
7.3
Vliv residuí obrazu fázového rozdílu na algoritmy unwrappingu
Pro analýzu vlivu residuí (diskontinuit) na algoritmy unwrappingu byl simulován fázový rozdíl Δϕw(r) dvou interferujících vln vztahem ∆φw (r ) = W {∆φ (r ) + δ KL (r )},
(308)
kde W je wrapped operátor a δKL(r) funkce náhodného rozdělení residuí maximální velikosti K × L s rovnoměrným pravděpodobnostním rozdělením.
Pro hodnocení reakcí algoritmů na přítomnost diskontinuit byla zaváděna dvě kritéria, a to: •
míra eliminace residuí vztahem Mr
R=
∑ (∆φ i =1
− ∆φuri )
2
ri
Mr
,
(309)
kde Mr je počet obrazových bodů interferogramu, v nichž nejsou přítomna residua, Δϕri a Δϕuri jsou hodnoty skutečného fázového rozdílu interferujících vln v i-tém bodě a fázového rozdílu po provedeném unwrappingu, jestliže na i-tém místě v původním obraze není přítomna diskontinuita. Tímto způsobem je hodnocena kvalita rekonstrukce míst, kde jsou zaznamenána data odpovídající skutečnému fázovému rozdílu. Dalším kritériem je poté •
doba t provádění unwrappingu.
Multiplikací výše uvedených kritérií získáme testovací hodnotu pro jednotlivé algoritmy T = Rt .
(310)
125
Analýza algoritmů unwrappingu Obr. 48 zobrazuje testovanou funkci a její rekonstrukci pomocí jednoduchého algoritmu, tab. 3 poté hodnoty jednotlivých testovacích kritérií při použití různých metod odstranění nespojitostí fázového rozdílu a obr. 49, obr. 50 a obr. 51 rekonstruované fázové rozdíly jednotlivými algoritmy.
Obr. 48 a) Testovaná funkce, b) její wrap podoba a c) rekonstrukce pomocí jednoduchého algoritmu
Tab. 3 Míra eliminace residuí v rekonstruovaném obraze a její relativní chyba vzhledem k původnímu obrazu a relativní čas výpočtu vztažený ke Goldsteinově algoritmu Algoritmus
R [rad]
trel [s]
0
---
GA
8.35∙1e-14 (0%)
1
QGPFA
4.49∙1e-15 (0%)
4.0
MCA
0.08 (8%)
1.6
FA
0.02 (2%)
14.0
ULS FT
0.11 (11%)
1.3
ULS MG
0.11 (11%)
47.4
WMG
4.04 (404%)
109.4
PCG
0.19 (19%)
26.7
Minimalizace LP-normy
1.06 (106%)
139.4
PUMA
1.60∙1e-15 (0%)
425.3
SRNCP
2.77∙1e-7 (<1%)
1.4
Iso-phase
0.13 (13%)
5.8
IIR
0.43 (43%)
221.9
původní obraz
126
Analýza algoritmů unwrappingu
Obr. 49 Rekonstruovaný fázový rozdíl pomocí algoritmů sledujících předem definovanou cestu – a) GA, b) QGPFA, c) MCA, d) FA
Obr. 50 Rekonstruovaný fázový rozdíl pomocí algoritmů minimalizujících normu – a) ULS MG, b) PCG, c) WMG, d) Minimalizace LP-normy
127
Analýza algoritmů unwrappingu
Obr. 51 Rekonstruovaný fázový rozdíl pomocí kombinovaných algoritmů – a) PUMA, b) SRNCP, c) Iso-phase, d) IIR
Analýza ukazuje, že 70% algoritmů si poradilo se simulovanými daty, a dokázalo rekonstruovat nespojitý obraz fázového rozdílu za přítomnosti residuí s chybou menší než 15%.
128
8 Vyhodnocení simulovaných a reálných dat 8.1
Výpočetní software k vyhodnocení interferogramu pomocí FTM
Pro vyhodnocení interferogramů metodou Fourierovy transformace byl vytvořen software v prostředí MATLAB R2011a principem „switched board programming“ [109], který poskytuje elegantní možnost konstrukce grafického rozhraní. Na obr. 52 je zobrazeno schéma a funkce vytvořené aplikace. Jako vstupní data aplikace jsou předpokládány předzpracované obrazy interferenčního pole. Na obr. 53 je zobrazeno hlavní okno programu, ve kterém probíhá koordinace celého vyhodnocení. Jsou zde vidět označená místa na hranici zájmového území (pupily) a otevřena volba „Oprava polohy bodů“ pro možnou editaci již zadaných bodů. Tímto způsobem lze pro kruhové pupily eliminovat oblasti, kde je zacloněno interferenční pole. Z vybraných bodů je oblast určována metodou nejmenších čtverců. Obr. 54 znázorňuje okno vytvořené k zadávání parametrů algoritmu pro extrapolaci interferenčních proužků, které slouží k potlačení jevů vznikajících u omezených oblastí v důsledku definice Fourierovy transformace na intervalu neomezeném. Gerchbergův algoritmus, představený dříve, vyžaduje znalost rozložení nerovnoměrné intenzity pozadí v obraze. Program umožňuje zadat tuto intenzitu konstantní, je-li známa z dřívější analýzy, nebo zvolit automatické určení pomocí průměrujícího filtru. Následně je zadáván počet iterací algoritmu a poloměr ořezu bočního spektra. Na obr. 55 je zobrazeno okno aplikace pro nastavení parametrů filtrace bočního spektra ve frekvenční oblasti. V příkladu je zadána volba určení maxima pomocí váženého průměru z oblasti, ve které hodnoty spektra přesahují zadaný limit. Tento limit představuje relativní část maximální hodnoty reálného spektra. Dále jsou voleny poloměry kruhových filtrů pro jednotlivé složky a poloměr útlumu pro eliminaci ostrých hran ve výsledném filtrovaném spektru.
129
Vyhodnocení simulovaných a reálných dat
PŘEDZPRACOVANÝ OBRAZ
HLAVNÍ OKNO
Výběr a editace bodů na hranici zájmové oblasti
Gerchbergův algoritmus • způsob určení intenzity pozadí • počet iterací • poloměr ořezu spektra
VÝBĚR ZÁJMOVÉ OBLASTI (PUPILY)
EXTRAPOLACE INTERFERENČNÍCH PROUŽKŮ
NASTAVENÍ FILTRACE SPEKTRA
VYHODNOCENÍ INTERFEROGRAMU
Kruhový filtr spekter s útlumem • způsob určení maxima bočního spektra • poloměr filtrace • poloměr útlumu
Algoritmus unwrappingu • volba parametrů jednotlivých algoritmů Aproximace rekonstruovaných dat • typ a stupeň aproximace
Obr. 52 Schéma programu k vyhodnocení interferogramu metodou Fourierovy transformace
130
Vyhodnocení simulovaných a reálných dat
Obr. 53 Hlavní okno aplikace
Obr. 54 Okno „Extrapolace interferenčních proužků”
131
Vyhodnocení simulovaných a reálných dat
Obr. 55 Okno pro nastavení parametrů filtrace frekvenčního spektra
Obr. 56 Dialogová okna nastavení parametrů algoritmů unwrappingu a aproximace
132
Vyhodnocení simulovaných a reálných dat Na obr. 56 jsou zobrazena dialogová okna, která slouží k volbě typu algoritmu unwrappingu a zadání jeho případných parametrů, a dále k volbě typu aproximačních polynomů (funkcí). Příklad ukazuje volbu algoritmu minimalizace LP normy pro odstranění nespojitostí fáze interferenčního obrazce s parametry popsanými v předchozí kapitole a výběr aproximace pomocí racionálních funkcí, kde jsou koeficienty počítány optimalizačním algoritmem Levenberg-Marquardt [82, 83]. Následující obr. 57, obr. 58 a obr. 59 ukazují jednotlivé kroky vyhodnocení. Jako první je představeno dialogové okno s wrapped fází a její formou po provedení unwrappingu. Následuje zobrazení rekonstruované fáze bez úprav a fáze po eliminaci vlivu zbytkové prostorové nosné frekvence. Posledním výstupem programu je aproximovaná forma fáze a její tvar při odstranění konstantního a lineárních členů rozvoje.
Obr. 57 Vyhodnocení fáze interferenčního obrazce
133
Vyhodnocení simulovaných a reálných dat
Obr. 58 Dialogové okno „Rekonstrukce vlnoplochy“
Obr. 59 Dialogové okno „Aproximace vlnoplochy“
134
Vyhodnocení simulovaných a reálných dat
8.2
Vyhodnocení simulovaných a reálných dat
Výše představená aplikace umožňuje snadno vyhodnocovat interferogramy metodou Fourierovy transformace. V této kapitole bude představeno vyhodnocení a analýza výsledků pro simulovaná a reálná data, která byla získána měřením v laboratoři Skupiny aplikované optiky FSv ČVUT v Praze [110], a dále data jež byla získána v laboratořích firmy Meopta – optika, s.r.o. [111]. Na obr. 60 a) je zobrazena simulovaná asférická plocha, která reprezentuje dráhový rozdíl dvou interferujících vln W(x, y). Její vyjádření bylo zvoleno formou rozvoje v řadu Zernikeho funkcí W (r , θ ) = −0.5 ⋅ Z 20 (r , θ ) + 0.1 ⋅ Z 3−1 (r , θ ) + 0.1 ⋅ Z 42 (r , θ ) ,
(311)
kde 0 ≤ r ≤ 1 a 0 ≤ θ ≤ 2π jsou polární souřadnice a Z nm (r , θ ) je Zernikeho funkce stupně n a řádu m. V části b) téhož obrázku je poté ukázán simulovaný interferogram daný vztahem I ( x, y ) = a ( x, y ) + b( x, y ) cos[2 πxf 0 x + kW ( x, y )] + N ( x, y ) ,
(312)
π 2 a ( x, y ) = 2 cos x + y2 2
(313)
kde
je nerovnoměrné rozložení intenzity pozadí, b(x, y) = 2 je variace amplitudy, f0x je prostorová nosná frekvence ve směru osy x, k = 2π/λ je vlnové číslo a N(x, y) je aditivní šum, u kterého předpokládáme normální rozdělení se střední hodnotou μ = 0 a směrodatnou odchylkou σ = 0.2.
Obr. 60 a) Simulovaná asférická plocha představující dráhový rozdíl dvou interferujících vln a b) simulovaný interferogram se zavedeným náhodným šumem a nerovnoměrným rozložením intenzity pozadí
135
Vyhodnocení simulovaných a reálných dat Na obr. 61 je zobrazen dialog okna pro extrapolaci interferenčních proužků. Je patrné, že nerovnoměrné rozložení intenzity pozadí bylo automaticky odhadnuto korektně. Dále extrapolované proužky v oblasti hranice pupily správně navazují, mělo by tak dojít k potlačení vlivu vyhodnocení interferogramu metodou FTM na omezené oblasti. Obr. 62 zobrazuje rekonstruovanou vlnoplochu opravenou o zbytkovou chybu z nejednoznačného určení prostorové nosné frekvence a její odchylku od plochy simulované. Obr. 63 následně ukazuje plochu aproximovanou rozvojem v řadu Zernikeho funkcí stupně n = 8 a odchylku této aproximace od původní simulované asférické plochy dané vztahem (311).
Obr. 61 Extrapolace interferenčních proužků simulovaných dat
Obr. 62 a) Rekonstruovaná vlnoplocha vyhodnocením interferogramu metodou Fourierovy transformace a b) její odchylka od simulované asférické plochy
136
Vyhodnocení simulovaných a reálných dat
Obr. 63 a) Aproximovaná rekonstruovaná vlnoplocha rozvojem v řadu Zernikeho funkcí stupně n = 8 a b) její odchylka od původní simulované asférické plochy
Z výše uvedených výsledků lze graficky porovnat kvalitu vyhodnocení simulované vlnoplochy. Pro kvalitativní hodnocení rekonstrukce fáze interferenčního obrazce použijme dvě kritéria, a to: •
PV (z angl. peak-to-valley) PV = max(W ) − min(W ) ,
•
(314)
RMS (z angl. Root Mean Square)
∑ (W M
RMS =
i =1
i
− Wri )
M
2
,
(315)
kde Wi a Wri značí simulovanou a rekonstruovanou vlnoplochu v i-tém bodě a M je počet všech bodů v zájmové oblasti. Přehled výsledků ukazuje tab. 4.
Tab. 4
Kvalitativní hodnocení rekonstrukce fáze interferenčního obrazce PV [λ]
RMS [λ]
rekonstruovaná vlnoplocha
0.104 (≈ λ/10)
0.018 (≈ λ/56)
aproximovaná vlnoplocha
0.056 (≈ λ/18)
0.018 (≈ λ/56)
Z výše představeného kvalitativního hodnocení můžeme vyslovit závěr, a to jestliže bychom předpokládali výše zmíněnou asférickou plochu a k vytvoření interferogramu by bylo použito
červeného kvazimonochromatického záření (He-Ne laser) o střední hodnotě vlnové délky
137
Vyhodnocení simulovaných a reálných dat λ = 633nm, potom 95% chyba rekonstruované fáze by měla hodnotu RMS95% = 1.96∙RMS = = 22.3nm. V praxi je snahou zajišťovat co nejlepší podmínky pro získávání interferometrických dat, a přibližovat se tak k ideálním případům, které můžeme simulovat a následně kvalitativně hodnotit. Vyhodnocení reálného interferogramu získaného měřením v laboratoři Skupiny aplikované optiky FSv ČVUT v Praze představuje obr. 64.
Obr. 64 Vyhodnocení reálného interferogramu – a) interferogram, b) a c) rekonstruovaná vlnoplocha
Pro porovnání vyhodnocení pomocí vytvořené aplikace VI-FTM v rámci této práce s komerčně nasazenými programy byly získány interferogramy měřené na interferometru ESDI Intellium Z100 HR a vyhodnocené softwarem IntelliWave firmy ESDI [60] v laboratořích firmy Meopta – optika, s.r.o. [111]. Pro jednoduchost označme analyzované interferogramy jako I1 a I2. Následující obrázky nabízí možnost porovnání jednotlivých výsledků.
Obr. 65 Měřený interferogram I1
138
Vyhodnocení simulovaných a reálných dat
Obr. 66 Vyhodnocená vlnoplocha pro interferogram I1 a) pomocí softwaru IntelliWave firmy ESDI, b) pomocí aplikace VI-FTM a c) rozdíl výsledků
Obr. 67 Měřený interferogram I2
Obr. 68 Vyhodnocená vlnoplocha pro interferogram I2 a) pomocí softwaru IntelliWave firmy ESDI, b) pomocí aplikace VI-FTM a c) rozdíl výsledků
139
Vyhodnocení simulovaných a reálných dat Tab. 5
Porovnání vyhodnocení interferogramu pomocí softwaru IntelliWave firmy ESDI a aplikace VI-FTM PV [λ]
RMS [λ]
interferogram I1
0.211 (≈ λ/5)
0.013 (≈ λ/77)
interferogram I2
0.554 (≈ λ/2)
0.019 (≈ λ/53)
Jak ukazuje tab. 5, pomocí vytvořené aplikace VI-FTM lze úspěšně vyhodnocovat reálné interferogramy. Ve výše představených dvou situacích není RMS odchylka mezi vyhodnocením dat pomocí VI-FTM a softwaru IntelliWave firmy ESDI větší než λ/53. Hodnota PV je v obou případech poměrně velká. Z obr. 66 a obr. 68 je patrné, že se největší odchylky nacházejí v hraničních oblastech. Tento jev vyplývá z použití Fourierovy transformace, která předpokládá data periodická, kdežto interferogramy jsou prostorově omezené. Je to nepříznivá typická vlastnost metody FTM a způsoby jejího potlačení jsou neustále vyvíjeny. Jednou z možností redukce těchto vlivů je maskování hraničních oblastí výsledných rekonstruovaných polí. Tab. 6 ukazuje kvalitativní hodnocení rekonstrukce výše představených interferogramů při omezení poloměru pupily výsledných dat na 85% původní velikosti. Jak je vidět, omezíme-li se na menší oblast, potom odchylky výrazným způsobem (několikanásobně) poklesnou. Tab. 6 Porovnání vyhodnocení interferogramu pomocí softwaru IntelliWave firmy ESDI a aplikace VI-FTM při omezení poloměru pupily výsledných dat na 85% původní velikosti PV [λ]
RMS [λ]
interferogram I1
0.111 (≈ λ/9)
0.010 (≈ λ/100)
interferogram I2
0.151 (≈ λ/7)
0.016 (≈ λ/63)
140
9 Závěr V technické praxi je častým cílem velmi přesně a bezkontaktně měřit topografii ploch rozličných tvarů. Je známo několik způsobů takového měření. Velmi často bývá používána metoda interferometrie, která pomocí interference vln dokáže poskytnou výsledky s přesností ve zlomcích vlnové délky použitého záření. Práce představila teoretické podklady k pochopení problematiky interferometrie, speciální pozornost byla poté věnována principu vyhodnocení interferogramů pomocí diskrétní Fourierovy transformace. Tato metoda nachází uplatnění při analýze dynamických procesů a všude tam, kde je používán pouze jeden snímek interferenčního pole. První část práce byla věnována matematické formulaci Fourierovy transformace pro spojité a diskrétní případy, následovala kapitola popisující fyzikální původ interference světla pro monochromatické, kvazimonochromatické a polychromatické záření. V další části byla představena interferometrie jako technika měření, její základní principy a postupy získávání, předzpracování, vyhodnocení a konečných úprav analyzovaných obrazů. Následovala kapitola věnovaná přímo metodě Fourierovy transformace s náležitostmi, které je nutné při rekonstrukci reálných fázových polí provádět. Jednou z částí procesu vyhodnocení je tzv. unwrapping. Tento matematický problém je stěžejní a nejkomplikovanější při rekonstrukci reálných polí, a proto mu byla věnována celá další kapitola. Algoritmy provádějící unwrapping byly následně analyzovány v předposlední kapitole a závěrem byl představen výpočetní software vytvořený k vyhodnocení interferogramů metodou Fourierovy transformace a vyhodnocována simulovaná a reálná data. Aplikace pro vyhodnocení interferogramů metodou Fourierovy transformace poskytuje možnost výběru relevantní části vyhodnocovaného obrazu (oblast kruhové pupily), extrapolaci interferenčních proužků mimo vybranou oblast, vyhodnocení pomocí filtrace a centrace bočního spektra a následnou aproximaci získaných dat analytickým vyjádřením rozvojem v řady funkcí nebo polynomů. Tento software byl napsán v prostředí MATLAB a využívá grafického rozhranní, tudíž samotné vyhodnocení je snadným úkolem. Jak ukazuje srovnání rekonstrukce interferogramů s komerčně dostupným softwarem v závěru poslední kapitoly, program může najít další uplatnění v technické praxi.
141
10 Použitá literatura [1]
BORN, Max, Emil WOLF a A. B. BHATIA. Principles of optics: electromagnetic theory of propagation, interference and diffraction of light. 7 ed. Cambridge: Cambridge University Press, 1999, 952 s. ISBN 05-216-3921-2.
[2]
BASS, Michael, ed. a Virendra N. MAHAJAN, ed. Handbook of optics Volume I , Geometrical and physical optics, polarized light, components and instruments. 3 ed. New York: McGraw-Hill, ©2010, 1195 s. ISBN 978-0-07-163598-1.
[3]
BASS, Michael, ed., Virendra N. MAHAJAN, ed. a Eric W. VAN STRYLAND, ed. Handbook of optics Volume II , Design, fabrication and testing, sources and detectors, radiometry and photometry. 3 ed. New York: McGraw-Hill, ©2010. 1178 s. ISBN 978-007-163600-1.
[4]
MIKŠ, Antonín. Aplikovaná optika. Vyd. 1. V Praze: České vysoké učení technické, 2009, 230 s. ISBN 978-80-01-04254-0.
[5]
PEATROSS a Michael WARE. Physics of Light and Optics. Brigham Young University, 2012.
[6]
MALACARA, Daniel, Manuel SERVÍN a Zacarias MALACARA. Interferogram analysis for optical testing. 2 ed. Boca Raton, FL: Taylor, 2005, 550 s. Optical engineering (Marcel Dekker, Inc.), v. 84. ISBN 15-744-4682-7.
[7]
MALACARA, Daniel. Optical shop testing. 3 ed. Hoboken, N.J.: Wiley-Interscience, ©2007, 862 s. ISBN 04-714-8404-0.
[8]
HARIHARAN, P. Optical interferometry. 2 ed. San Diego: Academic Press, 2003, 351 s. ISBN 978-0-12-311630-7.
[9]
YOSHIZAWA, Toru, ed. Handbook of optical metrology: principles and applications. Boca Raton, Fla.: CRC, ©2009. 730 s. ISBN 9780849337604.
[10] POKORNÝ, Petr, Jan OPAT a Antonín MIKŠ. Metody měření topografie ploch. Jemná mechanika a optika. 2013, č. 10, s. 301-303. ISSN 0447-6441. [11] BRAUNECKER, Bernhard, Rüdiger HENTSCHEL a Hans J TIZIANI. Advanced optics using aspherical elements. Bellingham, Wash.: SPIE Press, ©2008, 414 s. ISBN 08-1946749-9. [12] SU, Xianyu a Qican ZHANG. Dynamic 3-D shape measurement method: A review. Optics and Lasers in Engineering. 2010, vol. 48, issue 2, s. 191-204. DOI: 10.1016/j.optlaseng.2009.03.012.. [13] PETITGRAND, S., R. YAHIAOUI, K. DANAIE, A. BOSSEBOEUF a J.P. GILLES. 3D measurement of micromechanical devices vibration mode shapes with a stroboscopic interferometric microscope: A review. Optics and Lasers in Engineering. 2001, vol. 36, issue 2, s. 77-101. DOI: 10.1016/S0143-8166(01)00040-9.
142
Použitá literatura [14] GREN, P. Bending wave propagation in rotating objects measured by pulsed TV holography. Applied Optics. Washinghton: Optical Society of America, 2002, č. 34, 7237 7240. [15] LIU, J. B. a P. D. RONNEY. Modified Fourier transform method for interferogram fringe pattern analysis. Applied Optics. Washinghton: Optical Society of America, 1997, č. 25, 6231–6241. [16] SNYDER, Ray a Lambertus HESSELINK. Measurement of mixing fluid flows with optical tomography. Optics Letters. Washinghton: Optical Society of America, 1988, vol. 13, issue 2, s. 87-. DOI: 10.1364/OL.13.000087. [17] FRICKE-BEGEMANN, T a J BURKE. Speckle interferometry: three-dimensional deformation field measurement with a single interferogram. Applied Optics. Washinghton: Optical Society of America, 2001, č. 28, 5011 - 5022. [18] QUAN, C., P.J. BRYANSTON-CROSS a T.R. JUDGE. Photoelasticity stress analysis using carrier fringe and FFT techniques. Optics and Lasers in Engineering. Washinghton: Optical Society of America, 1993, vol. 18, issue 2, s. 79-108. DOI: 10.1016/01438166(93)90014-C. [19] KARAALIOGLU, Canan, P.J. BRYANSTON-CROSS a T.R. JUDGE. Fourier transform method for measurement of thin film thickness by speckle interferometry. Optical Engineering. Washinghton: Optical Society of America, 2003-06-01, vol. 42, issue 6,. DOI: 10.1117/1.1572498. [20] BORGHESI, M., A. GIULIETTI, D. GIULIETTI, L. A. GIZZI, A. MACCHI a O. WILLI. Characterization of laser plasmas for interaction studies: Progress in time-resolved density mapping. Physical Review E. 1996, č. 6, 6769 - 6773. [21] KALAL, M., B. LUTHER-DAVIES a K. A. NUGENT. Phaseamplitude imaging: the fully automated analysis of mega-gauss magnetic field measurements in laser-produced plasmas. J. Appl. Phys. 1988, s. 3845-3850. [22] CHANG, Chang, Patrick NAULLEAU, Erik ANDERSON a David ATTWOOD. Spatial coherence characterization of undulator radiation. Optics Communications. 2000, vol. 182, 1-3, s. 25-34. DOI: 10.1016/S0030-4018(00)00811-7. [23] CHANG, Chang, Erik ANDERSON, Patrick NAULLEAU, Eric GULLIKSON, Kenneth GOLDBERG a David ATTWOOD. Direct measurement of index of refraction in the extreme-ultraviolet wavelength region with a novel interferometer. Optics Letters. 2002, vol. 27, issue 12, s. 1028-. DOI: 10.1364/OL.27.001028. [24] OHTSUKA, Y., K. OKA a David ATTWOOD. Contour mapping of the spatiotemporal state of polarization of light. Applied Optics. Washinghton: Optical Society of America, 1994, č. 13, s. 2633-2636. [25] DROBCZYNSKI a KASPRZAK. Application of space periodic variation of light polarization in imaging polarimetry. Applied Optics. Washinghton: Optical Society of America, 2005, č. 16, s. 3160-3166. [26] NAIK, D. N., R. K. SINGH, H. ITOU, M. M. BRUNDAVANAM, Y. MIYAMOTO a M. TAKEDA. Single-shot full-field interferometric polarimeter with an integrated calibration scheme. Optics Letters. 2012, č. 15, s. 3282-3284. 143
Použitá literatura [27] OKA, Kazuhiko, Takayuki KATO, H. ITOU, M. M. BRUNDAVANAM, Y. MIYAMOTO a M. TAKEDA. Spectroscopic polarimetry with a channeled spectrum. Optics Letters. 1999, vol. 24, issue 21, s. 1475-. DOI: 10.1364/OL.24.001475 [28] MOMOSE, A., W. YASHIRO, H. MAIKUSA a Y. TAKEDA. Highspeed X-ray phase imaging and X-ray phase tomography with Talbot interferometer and white synchrotron radiation. Opt. Express. 2009, s. 12540-12545. [29] YASUNO, Y., V. D. MADJAROVA, S. MAKITA, M. AKIBA, A. MORISAWA, C. CHONG, T. SAKAI, K.-P. CHAN, M. ITOH a T. YATAGAI. Three-dimensional and high-speed swept-source optical coherence tomography for in vivo investigation of human anterior eye segments. Opt. Express. 2005, s. 10652-10664. [30] DUBRA, A., C. PETERSON a C. DAINTY. Study of the tear topography dynamics using a lateral shearing interferometer. Opt. Express. 2004, 6278–6288. [31] SCHEDIN, S., G. PEDRINI a H. J. TIZIANI. Pulsed digital holography for deformation measurements on biological tissues. Applied Optics. Washinghton: Optical Society of America, 2000, č. 16, 2853–2857. [32] TONOMURA, A. Electron holography. TrAC Trends in Analytical Chemistry. 1989, vol. 8, issue 6, s. 226-229. DOI: 10.1016/0165-9936(89)87007-4. [33] TONOMURA, A. Applications of electron holography. Reviews of modern physics. Mineapolis: The American Physical Society, 1987, č. 3, s. 639-669. [34] YOSHIDA, K., T. OKUWAKI, N. OSAKABE, H. TANABE, Y. HORIUCHI, T. MATSUDA, K. SHINAGAWA, A. TONOMURA a H. FUJIWARA. Observation of recorded magnetization patterns by electron holography. IEEE Transactions on Magnetics. Mineapolis: The American Physical Society, 1983, vol. 19, issue 5, s. 1600-1604. DOI: 10.1109/TMAG.1983.1062674. [35] TAKEDA, M. a J. SUZUKI. Crystallographic heterodyne phase detection technique for highly-sensitive lattice-distortion measurements. Journal of the Optical Society of America A: Optics, Image Science, and Vision. Washington: Optical Society of America, 1996, č. 7, s. 1495-1500. [36] CHANG, Wen-Yen, Chih-Tien WANG, Chih-Yuan CHU a Jyun-Ru KAO. Mapping GeoHazard by Satellite Radar Interferometry. Proceedings of the IEEE: Optics, Image Science, and Vision. Washington: Optical Society of America, 2012, vol. 100, issue 10, s. 28352850. DOI: 10.1109/JPROC.2012.2194629. [37] RIGNOT, E., L. PADMAN, D. R. MACAYEAL a M. SCHMELTZ. Observation of ocean tides below the Filchner and Ronne Ice Shelves, Antarctica, using synthetic aperture radar interferometry: Comparison with tide model predictions. Journal of geophysical research oceans. 2000, C8, s. 19615-19630. [38] VASCO, D. W., Charles WICKS, Kenzi KARASAKI a Osni MARQUES. Geodetic imaging: reservoir monitoring using satellite interferometry. Geophysical Journal International. 2002, vol. 149, issue 3, s. 555-571. DOI: 10.1046/j.1365246X.2002.01569.x.
144
Použitá literatura [39] KORN, Granino Arthur and Theresa M. KORN. Mathematical handbook for scientists and engineers: definitions, theorems, and formulas for reference and review. New York: Dover, 2000. 925 s. ISBN 0-486-41147-8. [40] REKTORYS, Karel. Přehled užité matematiky. 2. vyd. Praha: SNTL, 1968. 1140 s. [41] KOMRSKA, Jiří. Fourierova transformace. In: [online]. [cit. 2013-12-17]. Dostupné z: http://physics.fme.vutbr.cz/~komrska/Fourier/KapF02.pdf [42] ČÍŽEK, Václav. Diskrétní Fourierova transformace a její použití. Praha: SNTL, 1981, 160 s. [43] WILSON, Raymond G a Sean M MCCREARY. Fourier series and optical transform techniques in contemporary optics: an introduction. New York: Wiley, ©1995, xvii, 325 s. ISBN 04-713-0357-7. [44] MERVART, Leos and Zdenek LUKES. Adjustment Calculus. Vyd. 1. Praha: Nakladatelství ČVUT, 2007, 165 s. ISBN 978-80-01-03593. [45] COCHRAN, W.T., J.W. COOLEY, D.L. FAVIN, H.D. HELMS, R.A. KAENEL, W.W. LANG, G.C. MALING, D.E. NELSON, C.M. RADER a P.D. WELCH. What is the fast Fourier transform?. Proceedings of the IEEE. 1967, vol. 55, issue 10, s. 1664-1674. DOI: 10.1109/PROC.1967.5957. [46] COOLEY, J.W., P.A.W. LEWIS a P.D. WELCH. Historical notes on the fast Fourier transform. Proceedings of the IEEE. 1967, vol. 55, issue 10, s. 1675-1677. DOI: 10.1109/PROC.1967.5959. [47] GONZALEZ, Rafael C, Richard E WOODS a Steven L EDDINS. Digital Image processing using MATLAB. 2nd ed. Natick: Gatesmark Publishing, 2009, xviii, 826 s. ISBN 978-0-9820854-0-0. [48] GONZALEZ, Rafael C a Richard E WOODS. Digital image processing. 3rd ed. Upper Saddle River: Pearson, c2008, xxii, 954 s. ISBN 01-316-8728-X. [49] COOLEY, James W. a John W. TUKEY. An Algorithm for the Machine Calculation of Complex Fourier Series. Mathematics of Computation. 1965, č. 90, s. 297-301. [50] FFTW. [online]. [cit. 2013-12-17]. Dostupné z: http://www.fftw.org/ [51] STRATTON, Julius Adams. Electromagnetic theory. New York, N.Y: McGraw-Hill Book Company, 1941. ISBN 978-144-3730-549. [52] HAŇKA, L. Teorie elektromagnetického pole. 1. vyd. Praha: SNTL, 1982, 218 s. [53] MIKŠ, Antonín. Fyzika 2: elektromagnetické pole. Vyd. 1. Praha: Vydavatelství ČVUT, 2005, 162 s. ISBN 80-010-3164-0. [54] SILIOS Technologies diffractive optical element [online]. [cit. 2013-12-17]. Dostupné z: http://www.silios.com/ [55] Diffraction International [online]. [cit. 2013-12-17]. Dostupné z: http://www.diffraction.com/
145
Použitá literatura [56] WALDEN, R. H., N. L. SCHRYER, G. E. SMITH, R. J. STRAIN, J. MCKENNA a R. H. KRAMBECK. Buried channel charge coupled devices. The Bell System Technical Journal. New York: American Telephone and Telegraph Company, 1972, č. 7, s. 1635. [57] Piezoelectric transducer. Ultrasonics. 1971, vol. 9, issue 2, s. 124-. DOI: 10.1016/0041624X(71)90170-3. [58] Zygo Corporation [online]. [cit. 2013-12-17]. Dostupné z: http://www.zygo.com/ [59] 4D Technology :: The Leader in Innovative Optical Metrology [online]. [cit. 2013-12-17]. Dostupné z: http://www.4dtechnology.com/ [60] ESDI - Surface & Wavefront Metrology [online]. [cit. 2013-12-17]. Dostupné z: http://esdimetrology.com/ [61] Schneider GmbH & Co.KG - Fascination for Innovation [online]. [cit. 2013-12-17]. Dostupné z: http://www.schneider-om.com/ [62] OptoTech [online]. [cit. 2013-12-17]. Dostupné z: www.optotech.de [63] Trioptics - Optical Measurement, Autocollimators and visual optical instruments, MTF, EFL [online]. [cit. 2013-12-17]. Dostupné z: http://www.trioptics.com/ [64] QED Technologies [online]. [cit. 2013-12-17]. Dostupné z: http://qedmrf.com/ [65] MathWorks - MATLAB and Simulink for Technical Computing [online]. [cit. 2013-12-17]. Dostupné z: http://www.mathworks.com/ [66] SHANNON, C.E. Communication in the presence of noise. Proceedings of the IEEE. 1984, vol. 72, issue 9, s. 1192-1201. DOI: 10.1109/PROC.1984.12998. [67] NYQUIST, H. Certain Topics in Telegraph Transmission Theory. Transactions of the American Institute of Electrical Engineers. 1928, vol. 47, issue 2, s. 617-644. DOI: 10.1109/T-AIEE.1928.5055024. [68] GHIGLIA, Dennis C. a Mark. D. Pritt. Two-dimensional phase unwrapping: theory, algorithms, and software. 1 ed. New York: John Wiley, 1998, xiv, 493 s. ISBN 0-47124935-1. [69] WEBSTER, R. A generalized hamming window. IEEE Transactions on Acoustics, Speech, and Signal Processing. 1978, vol. 26, issue 2, s. 176-177. DOI: 10.1109/TASSP.1978.1163067. [70] NOVÁK, Jiří, Pavel NOVÁK a Antonín MIKŠ. Multi-step phase-shifting algorithms insensitive to linear phase shift errors. Optics Communications. Washinghton: Optical Society of America, 2008, vol. 281, issue 21, s. 5302-5309. DOI: 10.1016/j.optcom.2008.07.060. [71] NOVAK, Jiri. Five-step phase-shifting algorithms with unknown values of phase shift. Optik - International Journal for Light and Electron Optics. 2003, vol. 114, issue 2, s. 6368. DOI: 10.1078/0030-4026-00222.
146
Použitá literatura [72] CARRÉ, P. Installation et utilisation du comparateur photoélectrique et interférentiel du Bureau International des Poids et Mesures. Metrologia. 1966-01-01, vol. 2, issue 1, s. 1323. DOI: 10.1088/0026-1394/2/1/005. [73] TAKEDA, Mitsuo, Hideki INA a Seiji KOBAYASHI. Fourier-transform method of fringepattern analysis for computer-based topography and interferometry. Journal of the Optical Society of America. 1982, vol. 72, issue 1, s. 156-. DOI: 10.1364/JOSA.72.000156. [74] KEMAO, Q. Windowed Fourier transform for fringe pattern analysis. Applied Optics. 2004, č. 13, s. 2695-2702. [75] KEMAO, Qian. Two-dimensional windowed Fourier transform for fringe pattern analysis: Principles, applications and implementations. Optics and Lasers in Engineering. 2007, vol. 45, issue 2, s. 304-317. DOI: 10.1016/j.optlaseng.2005.10.012. [76] KEMAO, Q., H.X. WANG a W.J. GAO. Windowed Fourier transform for fringe pattern analysis: theoretical analyses. Applied Optics. Washinghton: Optical Society of America, 2008, č. 29, s. 5408-5419. [77] LARA-CORTEZ, Francisco, Cruz MENESES-FABIAN a Gustavo RODRIGUEZZURITA. Pfaff equation and Fourier analysis to phase extraction from an interferogram with carrier frequency. Journal of Physics: Conference Series. 2011-01-01, vol. 274, s. 012028-. DOI: 10.1088/1742-6596/274/1/012028. [78] MENESES-FABIAN, Cruz, Gustavo RODRIGUEZ-ZURITA, Alberto CORDERODAVILA, Carlos ROBLEDO-SANCHEZ, Pramod K. RASTOGI a Erwin HACK. Solving differential equations for phase retrieval in Fourier-transform methods. s. 118-123. DOI: 10.1063/1.3426096. [79] MEINARDUS, G. Aproximace funkcí: Teorie a numerické metody. SNTL, 1968. [80] POKORNÝ, Petr. Aproximace asférických ploch v optice. Jemná mechanika a optika. 2013, č. 10, s. 294-300. ISSN 0447-6441. [81] AHN, Sung Joon. Least squares orthogonal distance fitting of curves and surfaces in space. New York: Springer, c2004, xx, 125 p. ISBN 35-402-3966-9. [82] SCALES, L. Introduction to non-linear optimization. New York: Springer-Verlag, ©1985, xi, 243 p. ISBN 0333325524. [83] NOCEDAL, Jorge. a Stephen J. WRIGHT. Numerical optimization. 2nd ed. New York: Springer, ©2006. xxii, 664 s. ISBN 0-387-30303-0. [84] FRICKER, Paul. Analyzing LASIK Optical Data Using Zernike Functions. [online]. [cit. 2013-12-17]. Dostupné z: http://www.mathworks.com/company/newsletters/articles/analyzing-lasik-optical-datausing-zernike-functions.html [85] Periodic Table of the Disc Polynomials of Zernike. [online]. [cit. 2013-12-17]. Dostupné z: http://www.lumetrics.com/documents/wavefront/Zernike_table.pdf [86] NOVAK, J. a A. MIKS. Least-squares fitting of wavefront using rational function. Optics and Lasers in Engineering. Vol. 43, Issue 7, p. 776-787. ISSN 01438166. DOI: 10.1016/j.optlaseng.2004.08.002. 147
Použitá literatura [87] Numerical recipes: the art of scientific computing. 3rd ed. Cambridge: Cambridge University Press, 2007, xxi, 1235 s. ISBN 978-0-521-88068-8. [88] LIMPOUCH, Jiří. Čebyševovy polynomy [online]. 2000 [cit. 2013-12-17]. Dostupné z: http://kfe.fjfi.cvut.cz/~limpouch/numet/aprox/node16.html [89] STRATTON, Julius Adams. Spheroidal Wave Functions Including Tables of Separation Constants and Coefficients. 1. vyd. London: Chapman & Hall, 1958. 613 s. [90] MOORE, Ian C. a Michael CADA. Prolate spheroidal wave functions, an introduction to the Slepian series and its properties. Applied and Computational Harmonic Analysis. roč. 16, č. 3, s. 208-230. ISSN 10635203. DOI: 10.1016/j.acha.2004.03.004. [91] RODDIER, C. a F. RODDIER. Interferogram analysis using Fourier transform techniques. Applied Optics. Washinghton: Optical Society of America, 1987, s. 1668-1673. [92] MASSIG, J.H. a J. HEPPNER. Fringe-pattern analysis with high accuracy by use of the Fourier-transform method: theory and experimental tests. Applied Optics. Washinghton: Optical Society of America, 2001, č. 13, s. 2081-2088. [93] GERCHBERG, R.W. Super-resolution through Error Energy Reduction. Optica Acta: International Journal of Optics. 1974, vol. 21, issue 9, s. 709-720. DOI: 10.1080/713818946. [94] DAI, Mingjun a Yan WANG. Fringe extrapolation technique based on Fourier transform for interferogram analysis. Optics Letters. 2009, vol. 34, issue 7, s. 956-. DOI: 10.1364/OL.34.000956. [95] BIOUCAS-DIAS, Jos M. a Gonalo VALADAO. Phase Unwrapping via Graph Cuts. IEEE Transactions on Image Processing. 2007, vol. 16, issue 3, s. 698-709. DOI: 10.1109/TIP.2006.888351. [96] HERRAEZ, M. A., D. R. BURTON, M. J. LALOR a M. A. GDEISAT. Fast twodimensional phase-unwrapping algorithm based on sorting by reliability following a noncontinuous path. Applied Optics. 2002, č. 35, s. 7437-7444. [97] HERRAEZ, M. A., M. A. GDEISAT, D. R. BURTON a M. J. LALOR. Robust, fast, and effective two-dimensional automatic phase unwrapping algorithm based on image decomposition. Applied Optics. 2002, č. 35, s. 7445-7455. [98] ESTRADA, Julio C., Javier VARGAS, J. Mauricio FLORES-MORENO a J. Antonio QUIROGA. Windowed phase unwrapping using a first-order dynamic system following iso-phase contours. Applied Optics. 2012, vol. 51, issue 31, s. 7549-. DOI: 10.1364/AO.51.007549. [99] ITOH, Kazuyoshi. Analysis of the phase unwrapping algorithm. Applied Optics. 1982, vol. 21, issue 14, s. 2470-. DOI: 10.1364/AO.21.002470. [100] GOLDSTEIN, Richard M., Howard A. ZEBKER a Charles L. WERNER. Satellite radar interferometry: Two-dimensional phase unwrapping. Radio Science. 1988, vol. 23, issue 4, s. 713-720. DOI: 10.1029/RS023i004p00713. [101] BONE, Donald J. Fourier fringe analysis: the two-dimensional phase unwrapping problem. Applied Optics. 1991, vol. 30, issue 25, s. 3627-. DOI: 10.1364/AO.30.003627. 148
Použitá literatura [102] PRATI, C., M. GIANINI a N. LEURATTI. SAR Interferometry: A 2-D Phase Unwrapping Technique Based On Phase And Absolute Values Informations. Proceedings of the 1990 International Geoscience and Remote Sensing Symposium. IEEE, 1990, s. 2043-2046. [103] FLYNN, Thomas J. Two-dimensional phase unwrapping with minimum weighted discontinuity. Journal of the Optical Society of America A. 1997, vol. 14, issue 10, s. 2692-. DOI: 10.1364/JOSAA.14.002692. [104] LIMPOUCH, Jiří. Gauss-Seidelova metoda [online]. 2000 [cit. 2013-12-18]. Dostupné z: http://www-troja.fjfi.cvut.cz/~limpouch/numet/linalg/node24.html [105] TAKAJO, Hiroaki a Tohru TAKAHASHI. Noniterative method for obtaining the exact solution for the normal equation in least-squares phase estimation from the phase difference. Journal of the Optical Society of America A. 1988, vol. 5, issue 11, s. 1818-. DOI: 10.1364/JOSAA.5.001818. [106] KINDERMANN, Ross a J SNELL. Markov random fields and their applications. Providence, R.I.: American Mathematical Society, 1980, ix, 142 p. Contemporary mathematics (American Mathematical Society), v. 1. ISBN 08-218-5001-6. [107] Build MEX-function from C/C++ or Fortran source code - MATLAB mex [online]. [cit. 2013-12-18]. Dostupné z: http://www.mathworks.com/help/matlab/ref/mex.html [108] GETREUER, Pascal. Writing MATLAB C/MEX Code. [online]. [cit. 2013-12-18]. Dostupné z: http://classes.soe.ucsc.edu/ee264/Fall11/cmex.pdf [109] ZAPLATÍLEK, Karel a Bohuslav DOŇAR. MATLAB: tvorba uživatelských aplikací. 1. vyd. Praha: BEN, 2004, 215 s. ISBN 80-730-0133-0. [110] Applied Optics Group | CTU: Faculty of Civil Engineering, Department of Physics [online]. [cit. 2013-12-18]. Dostupné z: http://departments.fsv.cvut.cz/aog/ [111] Meopta :: Lepší pohled na svět [online]. [cit. 2013-12-18]. Dostupné z: http://www.meopta.cz/
149
11 Seznam obrázků Obr. 1
Funkce f(t) = f1(t) + f2(t) + f3(t) .......................................................................................18
Obr. 2
Zobrazení a) amplitud Ak vypočtených pomocí DFT, b) počátečních fází ϕk a c) polárního zobrazení dominantních členů rozvoje ...........................................................18
Obr. 3
Zobrazení a) rekonstruované funkce ξ(t) a b) jejího rozdílu od původní funkce f(t)......19
Obr. 4
Závislost intenzity záření na vzdálenosti pro koeficient absorpce α = 2, I0 = 0.5√(ε/μ)|E0|2; a) závislost v rovině y = 0, potom d = x, b) závislost v rovině (x,y), kde d = (x2 + y2)1/2 ...........................................................................................................24
Obr. 5
Interferenční obrazec pro dvě rovinné harmonické vlny zeleného světla o úhlové frekvenci ω = 4∙1015 Hz v rovině (y,z) kolmé na osu x ve vzdálenosti d = 1 od počátku soustavy souřadnic; jednotkové normálové vektory rovin voleny n1 = [1,0,0], n2 = 1/√(5)[2,0,1], δ1 = δ2 = 0, I1 = I2 = 1 .......................................................................27
Obr. 6
Prostorová koherence [4] ................................................................................................31
Obr. 7
Časová koherence [4]......................................................................................................35
Obr. 8
Základní konfigurace interferometru Twyman-Greenova typu [6] ................................41
Obr. 9
Principielní schéma interferometru Fizeauova typu [6]..................................................42
Obr. 10 Murtyho střihový interferometr [6].................................................................................43 Obr. 11 Nyquist-Shannonovo kritérium pro interferometrii ........................................................46 Obr. 12 Zobrazení a) C. F. Gausse (1777-1855) bez šumu, b) se zavedeným aditivním šumem se střední hodnotu μ = 0 a směrodatnou odchylkou σ = 15, c) průměrovaný obraz z 5 simulací zavedení šumu ..................................................................................................49 Obr. 13 Schéma frekvenční a prostorové filtrace.........................................................................50 Obr. 14 a) Reálný interferogram a b) jeho podoba po odstranění nerovnoměrného rozložení intenzity pozadí s použitím průměrujícího filtru [47, 48]...............................................53 Obr. 15 a) Původní interferogram a b) jeho podoba po vyrovnání histogramu ...........................55 Obr. 16 Normalizované histogramy a) původního obrazu a b) obrazu po provedení vyrovnání.55 Obr. 17 Princip metody detekce interferenčních proužků a následného vyhodnocení fázového rozdílu .............................................................................................................................58 Obr. 18 a) Číslování středů interferenčních proužků a b) rekonstrukce fázového rozdílu ..........60 Obr. 19 Zavádění lineární nosné frekvence u rovinných vlnoploch [6] ......................................64 Obr. 20 Schéma metody Fourierovy transformace ......................................................................66 Obr. 21 Pfaffova diferenciální rovnice.........................................................................................68 Obr. 22 Aproximace původních dat (červené tečky) Legendreovým polynomem (modrá křivka) a) stupně N = 2 a c) stupně N = 6 a b), d) odchylky aproximovaných od původních hodnot .............................................................................................................................74 Obr. 23 Zernikeho funkce do stupně n = 3...................................................................................76 Obr. 24 Symbolické znázornění Fourierovských spekter se zavedenou prostorovou nosnou frekvencí pro jednodimenzionální případ .......................................................................83 150
Seznam obrázků Obr. 25 a) Interferogram bez přítomnosti šumu a b) jeho reálná složka frekvenčního spektra ...83 Obr. 26 Určení maxima bočního spektra .....................................................................................86 Obr. 27 Schéma vyhodnocení interferogramu pomocí FTM .......................................................87 Obr. 28 Schéma extrapolace interferenčních proužků .................................................................89 Obr. 29 Zobrazení a) testovacích dat (funkce „peaks“) v původní podobě a b) ve wrapped formě .........................................................................................................................................93 Obr. 30 Zobrazení a) rekonstruovaných dat po unwrappingu, b) rozdílu rekonstruovaných dat po unwrappingu od skutečných hodnot ...............................................................................94 Obr. 31 a) Funkce ϕ(x) = 10∙cos 2x (modrá čárkovaná) a vzorkování při N = 10 (červeně) a N = 5 (zeleně); b) wrapped fáze a její vzorkování ..........................................................95 Obr. 32 a) Skutečná fáze (modře čárkovaně) a unwrapping při dostatečném vzorkování (červeně); b) skutečná fáze (modře čárkovaně) a unwrapping při nedostatečném vzorkování (zeleně).........................................................................................................95 Obr. 33 Principielní schéma jednoduchého 2D unwrappingu Itohovým algoritmem .................97 Obr. 34 a) Wrap operátor použitý na funkci ϕ(xi) = 10∙cos(2xi), b) operátor diference použitý na výsledek z části a) ...........................................................................................................98 Obr. 35 Zobrazení a) rekonstruované funkce pomocí Itohova algoritmu a b) rozdíl od původní hodnoty ...........................................................................................................................98 Obr. 36 Zobrazení a) simulované wrapped funkce „peaks“ bez residuí a šumu a b) její rekonstrukce jednoduchým algoritmem popsaným v úvodní části unwrappingu.........102 Obr. 37 Zobrazení a) wrapped funkce „peaks“ s náhodně generovanými residui a b) její rekonstrukce jednoduchým algoritmem popsaným v úvodní části kapitoly.................102 Obr. 38 a) Rozložení residuí v obraze (černé – záporná polarita, bílé – kladná polarita), b) nevhodné spojení dipólů, c) vhodné propojení dipólů [68] ..........................................106 Obr. 39 Zobrazení a) simulované fáze, b) fáze ve wrapped formě a c) polohy nespojitostí v obraze (detekované hrany) .........................................................................................108 Obr. 40 J. B. J. Fourier (1768-1830) jako funkce fi definovaná na 1 ≤ i ≤ 500 a 1 ≤ j ≤ 612 zrcadlená podle osy x = 500 a y = 612 ..........................................................................112 Obr. 41 a) Testovaná funkce, b) její wrap podoba a c) rekonstrukce pomocí jednoduchého algoritmu .......................................................................................................................119 Obr. 42 Vliv volitelných parametrů QGPFA algoritmu na míru šumu rekonstruovaného obrazu (modrá – MK minimálního gradientu, zelená – MK minimální variance, červená – MK maximální pseudokorelace) ..........................................................................................121 Obr. 43 Vliv volitelných parametrů MC algoritmu na míru šumu rekonstruovaného obrazu (modrá – MK minimálního gradientu, zelená – MK minimální variance, červená – MK maximální pseudokorelace) ..........................................................................................122 Obr. 44 Vliv volitelných parametrů Flynnova algoritmu na míru šumu rekonstruovaného obrazu (světlá místa znační malou hodnotu míry šumu) ..........................................................122 Obr. 45 Rekonstruovaný fázový rozdíl pomocí algoritmů sledujících předem definovanou cestu – a) GA, b) QGPFA, c) MCA, d) FA............................................................................123 Obr. 46 Rekonstruovaný fázový rozdíl pomocí algoritmů minimalizujících normu – a) ULS MG, b) PCG, c) WMG, d) Minimalizace LP-normy.....................................................124 151
Seznam obrázků Obr. 47 Rekonstruovaný fázový rozdíl pomocí kombinovaných algoritmů – a) PUMA, b) SRNCP, c) Iso-phase, d) IIR ....................................................................................124 Obr. 48 a) Testovaná funkce, b) její wrap podoba a c) rekonstrukce pomocí jednoduchého algoritmu .......................................................................................................................126 Obr. 49 Rekonstruovaný fázový rozdíl pomocí algoritmů sledujících předem definovanou cestu – a) GA, b) QGPFA, c) MCA, d) FA............................................................................127 Obr. 50 Rekonstruovaný fázový rozdíl pomocí algoritmů minimalizujících normu – a) ULS MG, b) PCG, c) WMG, d) Minimalizace LP-normy.....................................................127 Obr. 51 Rekonstruovaný fázový rozdíl pomocí kombinovaných algoritmů – a) PUMA, b) SRNCP, c) Iso-phase, d) IIR ....................................................................................128 Obr. 52 Schéma programu k vyhodnocení interferogramu metodou Fourierovy transformace 130 Obr. 53 Hlavní okno aplikace ....................................................................................................131 Obr. 54 Okno „Extrapolace interferenčních proužků”...............................................................131 Obr. 55 Okno pro nastavení parametrů filtrace frekvenčního spektra .......................................132 Obr. 56 Dialogová okna nastavení parametrů algoritmů unwrappingu a aproximace...............132 Obr. 57 Vyhodnocení fáze interferenčního obrazce ..................................................................133 Obr. 58 Dialogové okno „Rekonstrukce vlnoplochy“ ...............................................................134 Obr. 59 Dialogové okno „Aproximace vlnoplochy“..................................................................134 Obr. 60 a) Simulovaná asférická plocha představující dráhový rozdíl dvou interferujících vln a b) simulovaný interferogram se zavedeným náhodným šumem a nerovnoměrným rozložením intenzity pozadí ..........................................................................................135 Obr. 61 Extrapolace interferenčních proužků simulovaných dat...............................................136 Obr. 62 a) Rekonstruovaná vlnoplocha vyhodnocením interferogramu metodou Fourierovy transformace a b) její odchylka od simulované asférické plochy .................................136 Obr. 63 a) Aproximovaná rekonstruovaná vlnoplocha rozvojem v řadu Zernikeho funkcí stupně n = 8 a b) její odchylka od původní simulované asférické plochy................................137 Obr. 64 Vyhodnocení reálného interferogramu – a) interferogram, b) a c) rekonstruovaná vlnoplocha.....................................................................................................................138 Obr. 65 Měřený interferogram I1...............................................................................................138 Obr. 66 Vyhodnocená vlnoplocha pro interferogram I1 a) pomocí softwaru IntelliWave firmy ESDI, b) pomocí aplikace VI-FTM a c) rozdíl výsledků..............................................139 Obr. 67 Měřený interferogram I2...............................................................................................139 Obr. 68 Vyhodnocená vlnoplocha pro interferogram I2 a) pomocí softwaru IntelliWave firmy ESDI, b) pomocí aplikace VI-FTM a c) rozdíl výsledků..............................................139
152
12 Seznam tabulek Tab. 1
Volitelné parametry algoritmů unwrappingu................................................................119
Tab. 2
Míra šumu rekonstruovaného obrazu a její relativní chyba vzhledem k původnímu obrazu a relativní čas výpočtu vztažený ke Goldsteinově algoritmu............................123
Tab. 3
Míra eliminace residuí v rekonstruovaném obraze a její relativní chyba vzhledem k původnímu obrazu a relativní čas výpočtu vztažený ke Goldsteinově algoritmu .....126
Tab. 4
Kvalitativní hodnocení rekonstrukce fáze interferenčního obrazce..............................137
Tab. 5
Porovnání vyhodnocení interferogramu pomocí softwaru IntelliWave firmy ESDI a aplikace VI-FTM ..........................................................................................................140
Tab. 6
Porovnání vyhodnocení interferogramu pomocí softwaru IntelliWave firmy ESDI a aplikace VI-FTM při omezení poloměru pupily výsledných dat na 85% původní velikosti.........................................................................................................................140
153
13 Seznam příloh Příloha 1
Seznam publikací autora
Příloha 2
Seznam soutěžních prací autora
Datová příloha
Aplikace „VI-FTM“ pro 32bitový operační systém Windows
154
Příloha 1 Seznam publikací autora Impaktované mezinárodní časopisy [1]
A. Miks and P. Pokorny, "Analytical expressions for the circle of confusion induced by plane-parallel plate", Opt. Laser Eng. 50, 1517-1521 (2012) IF: 1.576
Recenzované české časopisy [1]
P. Pokorný, “Aproximace asférických ploch v optice“. Jemná mechanika a optika. 2013, 10, p. 294-300. ISSN 0447-6441.
[2]
P. Pokorný, J. Opat, A. Mikš, "Metody měření topografie ploch". Jemná mechanika a optika. 2013, 10, p. 301-303. ISSN 0447-6441.
[3]
A. Mikš and P. Pokorný, "Vliv disperse skla planparalelní desky na kvalitu zobrazení objektivu". Jemná mechanika a optika. 2013, 7-8, p. 214-218. ISSN 0447-6441.
[4]
A. Mikš, J. Novák, J. Opat, P. Pokorný, "Metoda rekonstrukce tvaru měřené optické plochy pomocí deflektometrického principu". Jemná mechanika a optika. 2013, vol. 58, no. 3, p. 89-93. ISSN 0447-6441.
[5]
A. Mikš, J. Novák, P. Novák, P. Pokorný, "Imaging properties of the lenses with paraboloidal surfaces". Jemná mechanika a optika. 2012, vol. 57, no. 11-12, p. 330-332. ISSN 0447-6441.
[6]
A. Mikš and P. Pokorný, "3D optické skenery". Jemná mechanika a optika. 2012, no. 5, p. 137-141. ISSN 0447-6441.
[7]
A. Mikš and P. Pokorný, "Vliv tloušťky planparalelní desky na zkreslení objektivu". Jemná mechanika a optika. 2011, 7-8, p. 220-223. ISSN 0447-6441.
Příloha 2 Seznam soutěžních prací autora Mezinárodní soutěže [1]
P. Pokorný, "Teoretické základy jedno-zrcadlových a dvou-zrcadlových optických skenerů". SVOČ 2012, sekce Geodézie a kartografie, 2. místo, http://www.fsv.cvut.cz/svoc/2012/vysledky.php [cit. 5.12.2013].
Celofakultní soutěže [1]
P. Pokorný: "Srovnávací studie aproximace funkcí". Vyčichlova soutěž 2013, 2. místo, http://mat.fsv.cvut.cz/vycichlo/2013/vysledky.html [cit. 5.12.2013].
[2]
P. Pokorný, "Analytical expressions for the circle of confusion induced by plane-parallel plate". Vyčichlova soutěž 2012, 3. místo, http://mat.fsv.cvut.cz/vycichlo/2012/vysledky.html [cit. 5.12.2013].
[3]
D. Bartošová, L. Nevoralová, P. Pokorný, "Výpočet potenciálu ve vnějším bodě Země z dat měřených na povrchu". Vyčichlova soutěž 2011, 1. místo, http://mat.fsv.cvut.cz/vycichlo/2011/vysledky.html [cit. 5.12.2013].
[4]
P. Pokorný, "Gibbsův jev". Vyčichlova soutěž 2011, 1. místo, http://mat.fsv.cvut.cz/vycichlo/2011/vysledky.html [cit. 5.12.2013].
