Cristian Grava
Vasile Buzuloiu
ELEMENTE DE PRELUCRAREA ŞI ANALIZA IMAGINILOR
2007
EDITURA UNIVERSITĂŢII DIN ORADEA
Descrierea CIP a Bibliotecii Naţionale a României GRAVA, CRISTIAN Elemente de prelucrarea şi analiza imaginilor / Cristian Grava, Vasile Buzuloiu. Oradea : Editura Universităţii din Oradea, 2007 ISBN 978-973-759-377-1 I. Buzuloiu, Vasile 621.397.3 (075.8)
EDITURA UNIVERSITĂŢII DIN ORADEA ESTE ACREDITATĂ DE CNCSIS, COD 149.
2
Cuprins: Prefaţă ………………………………………………………… 7 1. Introducere ………………………………………………………. 9 2. Reprezentarea imaginilor ………………………………………... 12 2.1. Digitizarea imaginilor ………………………………………. 12 2.2. Eşantionarea imaginilor …………………………………….. 15 2.3. Reprezentarea spaţială a imaginilor ………………………… 20 2.4. Proprietăţi ale imaginilor digitale ……………………….….. 25 2.4.1. Proprietăţi metrice ale imaginilor digitale ……….…… 26 2.4.2. Proprietăţi topologice ale imaginilor digitale ………… 27 2.4.3. Relaţii de vecinătate între pixeli ……………………… 28 2.4.4. Paradoxuri de conexitate ..……………………………. 31 2.4.5. Alte proprietăţi topologice şi geometrice ……………... 33 2.5. Reprezentarea spectrală a imaginilor ………………………... 34 2.5.1. Transformata Fourier (TF) bidimensională …………... 36 2.5.2. Proprietăţile transformatei Fourier bidimensionale …... 36 2.5.3. Proprietăţi specifice TF bidimensionale ……………… 43 3. Îmbunătăţirea imaginilor ………………………………………… 52 3.1. Calitatea unei imagini ………………………………………. . 52 3.2. Tehnici de îmbunătăţire a imaginilor ………………………… 55 3.3. Operatori punctuali de îmbunătăţire a imaginilor …………… 57 3.3.1. Operatori punctuali de modificare a contrastului ……… 58 3.3.2. Decuparea intervalelor de niveluri de gri ……………… 62 3.3.3. Modificarea histogramei ………………………………. 64 3.4. Operatori liniari de vecinătate pentru îmbunătăţirea imaginilor. Filtrarea liniară a imaginilor …………………………………. 68 3.5. Efectul în frecvenţă al operatorilor liniari de vecinătate …….. 73 3.6. Filtrarea neliniară a imaginilor ………………………………. 76 3.6.1. Filtre neliniare de ordine ………………………………. 77 3.6.2. Filtre de ordine multi-etaj ……………………………… 79 3.6.3. Proprietăţi ale filtrelor de ordine ………………………. 81 3
3.6.4. Filtre de ordine de domeniu …………………………… 83 3.6.5. L-filtre ………………………………………………… 85 4. Transformări integrale ale imaginilor ……………………………. 87 4.1. Transformări integrale unitare ………………………………. 87 4.2. Matrici unitare ………………………………………………. 93 4.3. Transformări unitare ale unor semnale unidimensionale …… 98 4.4. Transformări unitare ale unor semnale bidimensionale …….. 100 4.5. Transformata Fourier discretă unidimensională (DFT-1D) …. 104 4.6. Proprietăţi ale transformatei DFT-1D ……………………….. 105 4.7. Transformata Fourier discretă bidimensională (DFT-2D) ….... 108 4.8. Proprietăţi ale transformatei DFT-2D ……………………….. 111 4.9. Transformata Cosinus discretă unidimensională ……………. 113 4.10. Transformata Cosinus discretă bidimensională …………… 115 4.11. Transformata Sinus discretă unidimensională …………….. 120 4.12. Transformata Sinus discretă bidimensională ………………. 121 5. Restaurarea imaginilor …………………………………………… 122 5.1. Filtrarea inversă ……………………………………………… 123 5.2. Filtrul invers cu constrângeri ………………………………... 125 6. Morfologie matematică ………………………………..……. …... 130 6.1. Transformarea Hit or Miss ………………………………….. 130 6.2. Erodarea ……………………………………………………. 132 6.3. Dilatarea …………………………………………………….. 134 6.4. Proprietăţile operaţiilor morfologice ………………………... 136 6.5. Transformări morfologice derivate …………………………. 138 6.5.1. Operatori de extragere a conturului ……….…………. 138 6.5.2. Deschiderea şi închiderea …………………………….. 140 6.6. Trierea dimensională a obiectelor …………………………… 142 6.7. Caracterizarea morfologică a formelor ……………………… 143 6.7.1. Reconstrucţia după marker …………………………… 144 6.7.2. Distanţa Haussdorf …………………………………… 145 6.7.3. Extragerea skeletonului morfologic ………………….. 145 6.7.4. Skeletonul generalizat ………………………………… 149 6.8. Extinderea morfologiei matematice la imagini cu niveluri de gri.. 151 4
6.8.1. Trecerea de la mulţime la funcţie ……………………… 152 6.8.2. Trecerea de la funcţie la mulţime ……………………… 153 6.8.3. Operaţii cu funcţii …………………………………….. 154 7. Segmentarea imaginilor ….………………………………………. 157 7.1. Segmentarea orientată pe regiuni …………………………….. 157 7.1.1. Etichetarea componentelor ……………………………. 158 7.1.2. Metoda arborelui cuaternar (quad-tree) ………………. 159 7.2. Segmentarea imaginilor cu niveluri de gri …………………… 161 7.2.1. Segmentarea bazată pe histogramă ……………………. 161 7.2.2. Segmentarea bazată pe creşterea şi fuziunea regiunilor .. 165 7.3. Segmentarea orientată pe contururi …………………………. 168 7.3.1. Operatori de tip gradient ……………………………… 169 7.3.2. Operatori de tip compas ……………………………… 173 8. Compresia imaginilor ……………………………………………. 176 8.1. Compresia imaginilor binare ………………………………… 177 8.1.1. Codarea de nivel înalt …………………………………. 177 8.1.1.1. Aproximări poligonale …………………………. 177 8.1.1.2. Codul Freeman …………………………………. 178 8.1.1.3. Descriptori Fourier ……………………………… 181 8.1.2. Codarea la nivel de bloc ………………………………. 183 8.1.2.1. Metoda arborelui cuaternar (Quad-tree) ………… 183 8.1.2.2. Metoda WBS (White Block Skipping) …………. 184 8.1.3. Codarea la nivel de bit ………………………………… 185 8.1.3.1. Codarea RLE (Run Length Encoding) ………….. 185 8.1.3.2. Codarea entropică (Huffman) …………………… 186 8.2. Compresia imaginilor cu niveluri de gri ……………………… 188 8.2.1. Codarea pe plane ………………………………………. 188 8.2.2. Metode predictive de compresie ………………………. 190 8.2.3. Compresia cu transformate ……………………………. 194 Bibliografie ……………………………………………………….. 204
5
6
Prefaţă Prelucrarea şi analiza imaginilor ajunge pe zi ce trece, tot mai mult, în categoria „bunurilor de larg consum” ca urmare a aceleiaşi mutaţii suferită de calculatoarele electronice care reprezintă suportul „hard şi soft” al acestui domeniu relativ nou dar devenit indispensabil, al existenţei sociale contemporane. Numai că, fără îndoială, nu se rezumă la „suport”: analiza şi prelucrarea imaginilor – şi prin extensie, a secvenţelor video, a semnalelor multimedia şi, mai general, a semnalelor multidimensionale – înseamnă în primul rând algoritmi; algoritmi de prelucrări şi algoritmi de analiză. Şi, ajungând aici, ne dăm seama că, în fapt, avem de-a face cu „modelare matematică”. Cândva, pe vremea bunicilor bunicilor noştri, adică în zorii epocii moderne din istoria omenirii, cu arta aceasta se ocupau doar genii de prim rang – un Galilei, Newton, Leibnitz – dar numărul celor ocupaţi cu modelarea matematică a crescut vertiginos în secolul XIX şi apoi al XX-lea astfel încât astăzi, cum învăţământul universitar a ajuns „învăţământ de masă”, modelarea matematică a ajuns şi ea în învăţământul de masă; adică s-ar vrea un „bun de larg consum”. Numai că în matematică nu există „cale regală” adică pe care să înaintezi fără efort (asta o ştim încă din antichitate: i-a spus-o mentorul său împăratului Alexandru cel Mare!). Cele de mai sus se vor o justificare a faptului că această carte, care este şi un suport de curs pentru studenţii ingineri ai Universităţii din Oradea dar poate fi şi o lectură utilă pentru toţi cei interesaţi de subiect, conţine pagini întregi de matematică, constituind subiecte alese din nevoile disciplinei. În mod ideal, o carte de prelucrarea şi analiza imaginilor de sine stătătoare ar trebui să aibă un material imagistic ilustrativ, mult mai bogat. Rabatul pe care îl facem de la acest deziderat este justificat prin faptul că materia unui astfel de curs nu se rezumă la teorie ci are drept componentă principală şi partea de aplicaţii, în faţa ecranului: în felul acesta studentul „testează” pe viu ce influenţă au asupra rezultatelor modificările în algoritmi şi vede cum se modifică imaginea. De asemenea, ne-am decis să folosim în titlu un partitiv „Elemente de …” fiindcă de fapt cartea este doar o uşoară introducere într-un domeniu astăzi deja vast, care conţine printre multe altele, tehnicile de restaurare a 7
imaginilor (din care fac parte cele de reconstrucţie a imaginilor din proiecţii care stau la baza imagisticii medicale, a defectoscopiei nedistructive şi a teledetecţiei satelitare), ustensilele software pentru fotografia digitală (ne mărginim la un singur exemplu: corecţia în „timp real” a ochilor roşii) şi tot felul de sisteme de supraveghere a căror componentă principală trebuie să fie una capabilă de recunoaştere a formelor (maşini, feţe etc). În pofida faptului că „prelucrarea şi analiza imaginilor” este, cum se zice azi, un „domeniu de vârf”, noi am fost activi în el de mai bine de 30 de ani datorită cercetărilor de televiziune digitală şi prelucrare în timp real a semnalului digital în colectivul de cercetare al Catedrei de Electronică Aplicată de la Politehnica din Bucureşti; primul sistem digital de analiză a imaginilor a fost terminat în 1981 iar în următorii ani a fost reprodus în câteva exemplare la Fabrica de Calculatoare Electronice. Perioada anilor ’80 n-a fost, din păcate, propice dezvoltărilor tehnologice la noi şi astfel întâietatea şi avantajul pe care apucaserăm să îl avem „în lagărul socialist”, n-a dat roade. După 1990 ne-am trezit într-o lume care între timp progresase mult în acest domeniu. Totuşi, experienţa pe care o avem ne-a permis să reluăm ideile iar astăzi, Laboratorul de Analiza şi Prelucrarea Imaginilor LAPI al Politehnicii din Bucureşti este destul de bine cunoscut şi peste hotare; mai mult, transplantul cunoştinţelor l-am făcut şi spre alte universităţi – Braşov, Oradea – şi aşa se face că astăzi, în Catedra de Electronică a Universităţii din Oradea există un colectiv important de cadre tinere implicate în acest domeniu, cu stagii în străinătate şi colaborări internaţionale. Avem toate condiţiile preliminare pentru un învăţământ responsabil. Cartea de faţă, scrisă aproape în totalitate de Cristian Grava, adună elementele esenţiale ale unui curs introductiv dar presupune că studentul – cititorul – a fost „expus” înainte unor cursuri pregătitoare printre care cel de teoria statistică a semnalelor, cel de sisteme liniare şi cursurile fundamentale de matematică pentru ingineri. Ar mai fi de remarcat că, pentru compactizarea materialului s-a preferat o înşiruire neliniară a subiectelor, care de exemplu în capitolul 2 trimite la capitolul 4 dar, fiind vorba de subiecte reluate de la cursuri anterioare, suntem convinşi că cititorul va putea parcurge textul fără dificultăţi, cel mult eventual cu reluări. Prof. Vasile Buzuloiu 8
1. Introducere Dezvoltarea spectaculoasă din ultimii ani a tehnologiei informaţiei şi a componentelor electronice a condus la impunerea de soluţii inimaginabile până nu demult pentru numeroase probleme tehnice, ca de exemplu în industrie la conducerea proceselor de producţie sau la controlul de calitate. Tehnicile de exploatare a informaţiei vizuale ocupă o poziţie importantă şi de extremă actualitate. Domeniul prelucrării şi analizei imaginilor grupează tehnicile de achiziţie, transformare şi utilizare a informaţiei vizuale din imaginile reprezentate, transmise şi exploatate în formă digitală, în sisteme de calcul de uz general sau calculatoare specializate. Printre aplicaţiile importante ale prelucrării şi analizei imaginilor se pot aminti aplicaţiile în medicină (investigarea de organe ale corpului uman), aplicaţii în industrie şi în tehnică (cartografierea solului, prospecţiuni geologice, controlul tehnic automat al diverselor produse, robotică etc.) precum şi aplicaţii dedicate unor domenii ca arta, aplicaţiile militare etc. Principalele probleme ale prelucrării şi analizei de imagini sunt: 1. Reprezentarea şi modelarea imaginilor: • Eşantionarea şi cuantizarea imaginilor • Reprezentarea spaţială a imaginilor • Reprezentarea spectrală a imaginilor. Transformata Fourier
Imagine discretizată
Transformare
Imagine transformată
Figura 1.1. Transformarea unei imagini. 9
Prelucrarea se face asupra acestei imagini
• Modelele imaginilor pot fi: stohastice; deterministe. 2. Îmbunătăţirea imaginilor. Se face prin: • prelucrări punctuale. Exemple: modificarea contrastului şi luminanţei; modificarea histogramei. • prelucrări pe vecinătăţi. Exemple: accentuarea contururilor; reducerea zgomotului imaginilor; pseudo-colorarea. • prelucrări geometrice. • prelucrări integrale. 3. Restaurarea imaginilor. Presupunem că f(x,y) este imaginea originală care datorită procesului de captare a suferit o transformare (degradare), liniară sau neliniară, obţinându-se imaginea degradată fd(x,y). f(x,y)
transformare
fd(x,y)
transformare-1
Figura 1.2. Transformarea inversă a unei imagini. Imaginea originală se poate obţine din imaginea degradată prin aplicarea unei transformări inverse celei suferite în procesul de captare. Apar probleme când în procesul de captare intervine şi zgomot, caz în care se impune şi o etapă de restaurare a imaginii, pentru a se obţine o aproximare a imaginii originale fˆ(x,y) : 10
f(x,y)
Transformare
+
fd(x,y)+zg
Restaurare
ˆf (x,y)
zg=zgomot
Figura 1.3. Restaurarea imaginilor. • Reconstrucţia imaginilor din proiecţii.
4. Analiza imaginilor, care implică: • măsurători automate pe imagini; • segmentarea imaginilor (extragerea obiectelor).
5. Compresia imaginilor, care implică reducerea cantităţii de informaţie. Exemplu: Pentru o imagine de 512×512 pixeli (29×29), în care fiecare pixel este reprezentat pe 8 biţi, adică cu 28=256 niveluri de gri, cantitatea de informaţie este: (29×29)pixeli×23biţi/pixel=221 ≈ 2 Mb Prin compresia imaginilor se urmăreşte reducerea acestei cantităţi de informaţie. Schema bloc a unui lanţ de analiza şi prelucrarea imaginilor este prezentată în figura 1.4. Scenă
senzor
CAD
Îmbunătăţire + Restaurare
Măsurători
date
Segmentare Display
Display Compresie
CAD = convertor analog-digital
Listă de obiecte
Măsurători pe obiecte
Clasificare
Descrierea scenei
Memorie
Figura 1.4. Schema bloc a unui lanţ de analiza şi prelucrarea imaginilor. 11
2. Reprezentarea imaginilor Sistemele vizuale ale organismelor vii percep mediul înconjurător 3-dimensional, prin intermediul unor latici de senzori de lumină bidimensionale (de exemplu, retina din ochii mamiferelor) şi refac spaţiul 3D prin integrare temporală şi/sau vedere binoculară. În pofida faptului că senzaţia este de câmp continuu al imaginilor percepute, laticea senzorilor este discretă. În sistemele tehnice imaginate de om, aceste proprietăţi se păstrează: informaţia imagistică din mediul înconjurător se proiectează pe latici bidimensionale de senzori de lumină şi astfel se discretizează spaţial, iar semnalele de la fiecare senzor se discretizează în timp şi în valoare astfel încât „imaginile” ajung în sistemele de calcul sub formă digitală, pentru prelucrare şi analiză.
2.1. Digitizarea imaginilor Imaginile pot fi descrise de distribuţia spaţială a intensităţii luminoase într-un plan. Din punct de vedere matematic, distribuţia spaţială a intensităţii luminoase (I) poate fi descrisă printr-o funcţie continuă de două variabile spaţiale continue (x,y)=p: I(x,y)=I(p) Calculatoarele existente nu pot trata imaginile ca funcţii definite pe un domeniu continuu ci doar ca matrici discrete de numere. Din acest motiv este necesară transformarea şi reprezentarea imaginilor continue ca matrici bi-dimensionale de puncte, prin discretizare. Un punct al unei astfel de matrice se numeşte pixel (din engleză = picture element). Un pixel reprezintă intensitatea luminoasă sau culoarea corespunzătoare unui anumit punct din matrice. Prin urmare, un pixel este caracterizat prin poziţia şi prin valoarea sa. 12
Pentru tratări teoretice, imaginea bidimensională poate fi reprezentată ca o funcţie continuă (analogică) bidimensională I(x,y)=f(x,y), unde x şi y sunt coordonatele spaţiale. Valoarea funcţiei într-un punct oarecare (x,y) va reprezenta: • luminanţa din punctul respectiv, în cazul în care funcţia f(x,y) este o funcţie reală. În acest caz avem o imagine cu niveluri de gri, numită impropriu şi imagine alb-negru; • culoarea din punctul respectiv, în cazul în care funcţia f este o funcţie vectorială, (f1(x,y),f2(x,y),f3(x,y)) = (R,G,B). În acest caz avem o imagine color, cu componentele fundamentale (R, G, B). Trecerea de la o imagine color la o imagine cu niveluri de gri se face prin adunarea componentelor fundamentale ponderate cu anumiţi coeficienţi, adică printr-o combinaţie liniară a acestor componente. Pentru prelucrarea digitală a imaginilor analogice (de exemplu cu ajutorul unui calculator) este nevoie de discretizarea imaginilor, proces în urma căruia imaginea este transformată într-o matrice care conţine elementele de imagine (pixel). În practică, camera video de tip CCD (Charge Coupled Device) realizează discretizarea imaginilor chiar în procesul de captare. Pentru afişare, imaginile se pot converti din nou în formă analogică. Discretizarea imaginilor analogice se realizează în doi paşi: • discretizarea spaţială (eşantionarea), cu ajutorul unei reţele discrete
f(l·∆x,k·∆y), l,k∈Z. În urma acestei operaţii rezultă o imagine (matrice) cu L linii şi K coloane: l≤L, k≤K. Prin urmare, se obţine L×K pixeli, iar imaginea obţinută se va scrie printr-o expresie de forma: A={a(l,k), 1≤l≤L, 1≤k≤K}, l,k,L,K∈Z; a11
a12
…
a1,K
a21 … aL,1
a22 … aL,2
… … …
a2,K … aL,K
13
• discretizarea în valoare (cuantizarea): fq(l·∆x,k·∆y) l,k∈Z, l≤L, k≤K, cu:
fq(l·∆x,k·∆y)∈{f1,…,fn}, unde n este numărul nivelurilor de cuantizare a imaginii (numărul nivelurilor de gri). De exemplu, pentru n=2 avem o imagine binară. În acest caz, fiecare eşantion obţinut în pasul anterior (la eşantionare) este cuantizat folosind un număr finit de biţi. Astfel, fiecare pixel va avea un anumit nivel de gri (pentru imagini alb-negru) sau o anumită culoare (pentru imagini color), codificată printr-un număr constant de biţi. În reprezentare binară, un pixel oarecare al,k este codat: (al,k)binar=bn-1bn-2…b1b0 căruia îi corespunde o valoare zecimală: (al,k)zecimal=bn-1·2n-1+bn-2·2n-2+…+b1·21+b0·20=q
(2.1)
care reprezintă nivelul q din scara de 2n nivele de gri considerate. În mod uzual, negrul este considerat ca având nivelul logic 0 (în binar 00 …00), iar albul ca având nivelul logic 1 (în binar 11…11). De exemplu, în cazul unei imagini reprezentate pe 8 biţi avem un număr de 28 =256 niveluri de gri, în care negrul este codat cu nivelul q=0, iar albul este codat cu nivelul q=255. Pixelul a34 având valoarea codată binar cu octetul 00101001 este pixelul al 4-lea de pe rândul 3 şi are nivelul q=41 în scara de niveluri de gri amintită. În continuare, când se va vorbi despre imagini digitale sau simplu despre “imagini”, se va face referire la imagini eşantionate şi cuantizate, iar când se va vorbi despre discretizarea imaginilor se va face referire la discretizarea spaţială (eşantionarea) şi în valoare (cuantizarea) imaginilor.
14
2.2. Eşantionarea imaginilor Se consideră o matrice de eşantionare, cu pasul (∆x,∆y) care transformă imaginea dintr-o funcţie (continuă) într-un şir: f(x,y) →( f (l ⋅ ∆x,k ⋅ ∆y ))l,k∈Z esantionare
(2.2)
y
∆y
x ∆x
Figura 2.1. Matricea de eşantionare a unei imagini. Transformarea inversă (din şir în imagine) este posibilă în condiţiile teoremei eşantionării. Modelul matematic al semnalului eşantionat este: not . f ( x, y ), x = l ⋅ ∆x, y = k ⋅ ∆y f e ( x, y ) = 0, în rest
(2.3)
f e ( x, y ) = ∑ ∑ f ( x, y ) ⋅ δ ( x − m ⋅ ∆x, y − n ⋅ ∆y ) = f ( x, y ) ⋅ δ p ( x, y ) , m n
unde:
δ p ( x, y ) = ∑ ∑ δ ( x − l ⋅ ∆x, y − k ⋅ ∆y ) = δ periodic
(2.4)
l k
este impulsul Dirac δ periodic: 1, pt. x = 0 δ( x ) = 0, în rest
15
(2.5)
Transformata Fourier a semnalului eşantionat este: Fe (u , v) = ℑ{ f e ( x, y )} = ℑ f ( x, y ) ⋅ δ p ( x, y )
{
}
(2.6)
Se ştie că dacă funcţia f(t) este periodică (cu perioada T), seria sa Fourier este: f (t ) = ∑ C k ⋅ e
j
2π k ⋅t T
, unde coeficienţii:
(2.7)
k
2π
−j 1 C k = ∫ f (t ) ⋅ e T TT
k ⋅t
dt
(2.8)
Din acest motiv, deoarece funcţia δp(x,y) este periodică cu perioada ∆x pe x, respectiv ∆y pe y:
δ p ( x, y ) = ∑ ∑ C kl
2π 2π j k ⋅ x + l ⋅ y ∆x ∆y ⋅e
, unde:
(2.9)
k l
C kl =
1 ∫ ∫ δ p ( x, y ) ⋅ exp − ∆x ⋅ ∆y ( ∆x ) ( ∆y )
2π 2π j k⋅x+ l ⋅ y dxdy (2.10) ∆y ∆x
∆x ∆x , , respectiv Deoarece se integrează pe un interval ∆x = − 2 2 ∆y ∆y , şi dacă se presupune că într-un dreptunghi cu laturile ∆x, ∆y = − 2 2 ∆y cade un singur impuls δp şi numai unul (cel din origine):
16
z
k·∆y
l·∆x y
x
Figura 2.2. Eşantionarea imaginilor.
⇒ C kl =
∞ ∞ 1 ∫ ∫ δ p ( x, y ) ⋅ exp − ∆x ⋅ ∆y −∞ −∞
2π 2π j k⋅x+ l ⋅ y dxdy (2.11) ∆y ∆x
∞
Deoarece: ∫ δ (t ) ⋅ f (t )dt = f (0) , iar în cazul de faţă: f(0)=1 −∞
⇒ C k ,l =
1 ∆x ⋅ ∆y
(2.12) 2π
2π
j k ⋅x + l⋅ y 1 ∆x ∆y ⇒ δ p ( x, y ) = ∑ ∑ ⋅e k l ∆x ⋅ ∆y 2π 2π j k ⋅ x + l ⋅ y 1 ∆y ⇒ Fe (u , v) = ℑ ⋅ ∑ ∑ f ( x, y ) ⋅ e ∆x ∆x ⋅ ∆y k l
(2.13)
(2.14)
Pe baza proprietăţii de liniaritate a transformatei Fourier: 2π 2π j k ⋅ x + l ⋅ y 1 ∆y ⇒ Fe (u , v) = ⋅ ∑ ∑ ℑ f ( x, y ) ⋅ e ∆x ∆x ⋅ ∆y k l
17
(2.15)
⇒ Fe (u , v) =
1 2π 2π 1 k, v − l = ⋅ ∑ ∑ ℑ(u − k ⋅ ∆u , v − l ⋅ ∆v ) ⋅ ∑ ∑ ℑ u − ∆x ⋅ ∆y k l ∆x ∆y ∆x ⋅ ∆y k l not . 2π
unde: ∆u =
∆x
not . 2π
, ∆v =
(2.16)
∆y
Spectrul semnalului eşantionat F(u,v) se obţine prin periodizarea (repetarea) spectrului său Fe(u,v): FTJ
F(u,v)
v v ∆u u
u
∆v
Figura 2.3. Spectrul semnalului eşantionat. Recuperarea semnalului original din semnalul eşantionat se poate face cu un filtru trece-jos (FTJ) cu parametri adecvaţi (figura 2.4). ∆v-vmax ∆v vmax umax
∆u-umax
∆u
u Figura 2.4. Parametrii FTJ.
18
v
Parametrii umax, respectiv vmax reprezintă frecvenţele spaţiale maxime din spectrul funcţiei f, în direcţia u, respectiv v (ce corespund coordonatelor spaţiale x, respectiv y). Pentru ca semnalul original să fie corect recuperat cu ajutorul unui FTJ, trebuie ca parametrii acestuia să satisfacă condiţiile: ∆u − u max ≥ u max ∆v − v max ≥ v max
∆u ≥ 2 ⋅ u max ⇒ ∆v ≥ 2 ⋅ v max
(2.17)
Prin urmare, frecvenţele de tăiere (ξ, η) a FTJ de recuperare a semnalului original în cele două direcţii (u,v), trebuie să îndeplinească condiţiile: u max ≤ ξ ≤ ∆u − u max v max ≤ η ≤ ∆v − v max
(2.18)
Astfel, în acest caz, teorema eşantionării se poate enunţa astfel: 2π 2π , ∆v = adică frecvenţele de eşantionare pe x, respectiv dacă ∆u = ∆x ∆y pe y, sunt mai mari sau cel puţin egale cu dublul frecvenţelor maxime din spectrul lui f pe direcţia u, respectiv v, atunci recuperarea semnalului original de imagine f(x,y) se poate face exact, din eşantioane, cu un filtru ideal trece-jos cu funcţia de transfer: fe(x,y)
H(u,v)
f(x,y)
Figura 2.5. FTJ necesar pentru extragerea semnalului original. 1, pentru u ≤ ξ si v ≤ η H (u, v) = 0, în rest
(2.19)
unde frecvenţele de tăiere (ξ, η) sunt alese în mod corespunzător, adică astfel încât: 19
u max ≤ ξ ≤ ∆u − u max v max ≤ η ≤ ∆v − v max
(2.20)
Observaţie: Condiţia de recuperare enunţată de teorema eşantionării este suficientă dar nu şi necesară. Acest lucru este ilustrat de exemplul următor: v umax ∆u vmax ∆v-vmax
u F(u,v)
∆v
∆u-umax
H(u,v)
Figura 2.6. Caz particular de extragere a semnalului original. Se observă că deşi condiţia de recuperare nu este îndeplinită, recuperarea se poate face cu un FTJ ideal corespunzător, care să extragă doar zona corespunzătoare spectrului funcţiei.
2.3. Reprezentarea spaţială a imaginilor Odată digitizate (eşantionate şi cuantizate), imaginile pot fi prelucrate şi analizate cu sisteme de calcul uzuale sau dedicate. Reprezentarea imaginilor se poate face sub diverse forme (spaţială, spectrală etc.), adecvate diverselor aplicaţii. În cazul cel mai simplu, pixelii sunt localizaţi pe o reţea rectangulară. Poziţia unui pixel este dată în mod analog notaţiei utilizate 20
pentru elementele unei matrice. Primul indice (l) exprimă poziţia pe linie, iar cel de-al doilea indice (k) exprimă poziţia pe coloană (figura 2.1).
0
1
coloane
k
0
K-1 x
1
Cl,k
linii
l
L-1 y
Figura 2.7. Reprezentarea imaginilor digitale ca matrici de pixeli dispuşi într-o reţea rectangulară bi-dimensională. Dacă imaginea conţine L×K pixeli, aceasta poate fi reprezentată printr-o matrice de dimensiune L×K, unde indicele l=0…L-1, iar k=0…K-1. L reprezintă numărul de linii, iar K numărul de coloane. Ca şi în cazul matricilor, sensul pozitiv al axei verticale (y) este de sus în jos şi nu de jos în sus, cum este cazul reprezentărilor grafice bidimensionale uzuale. Sensul pozitiv al axei orizontale (x) este cel uzual, de la stânga la dreapta (figura 2.1). Rezoluţia spaţială de reprezentare a unei imagini poate fi definită ca reprezentând numărul total de pixeli (de exemplu L×K) sau poate fi definită ca fiind egală cu numărul de pixeli pe unitatea de suprafaţă (în pixeli/cm2 sau în pixeli/inch2). Rezoluţia spaţială a unui sistem de achiziţie de imagini se poate defini ca fiind egală cu numărul de pixeli pe unitatea de lungime (pixeli/mm sau pixeli/cm). Pe baza acestor noţiuni, se mai poate defini şi sensibilitatea unui sistem de vedere sau a unui sistem de 21
achiziţie de imagini, ca fiind unitatea minimă de lungime care poate fi observată într-o imagine achiziţionată. Fiecare pixel reprezintă nu numai un punct al unei imagini, ci o regiune rectangulară a acesteia, care defineşte o celulă elementară a imaginii. Valoarea asociată unui pixel reprezintă în mod adecvat media intensităţii luminoase din celula corespunzătoare. În figura 2.8 este ilustrată una şi aceeaşi imagine, reprezentată printr-un număr diferit de pixeli.
(a)
(b)
(c) (d) Figura 2.8. Imagine digitală cu diferite rezoluţii: (a) - 16×16 pixeli; (b) - 32×32 pixeli; (c) - 64×64 pixeli; (d) - 256×256 pixeli. În cazul unor pixeli de dimensiune mare (figura 2.8.a şi 2.8.b), nu numai că rezoluţia spaţială este mică, dar apar nişte artefacte (zgomote) 22
deranjante datorate discontinuităţilor de niveluri de gri de la marginile pixelilor, care distrag atenţia privitorului de la conţinutul propriu-zis al imaginii. Atunci când dimensiunea pixelilor devine mai mică, adică atunci când creşte rezoluţia spaţială (figura 2.8.c şi 2.8.d), efectele descrise mai sus devin mai puţin pronunţate, putându-se ajunge până la impresia de continuitate spaţială a imaginii. Acest lucru de întâmplă când rezoluţia spaţială a imaginii devine mai mare decât rezoluţia sistemului uman de vedere, adică atunci când dimensiunea unui pixel al imaginii devine mai mică decât dimensiunea minimă pe care o poate percepe ochiul uman. Nu există un răspuns general valabil legat de numărul optim de pixeli necesar pentru a crea senzaţia de continuitate spaţială a unei imagini. În cazul observării vizuale a unei imagini, trebuie ca dimensiunea unui pixel să fie mai mică decât dimensiunea corespunzătoare rezoluţiei spaţiale a sistemului vizual, la o distanţă nominală a observatorului. În cazul unei aplicaţii concrete, dimensiunea unui pixel trebuie să fie obligatoriu mai mică decât dimensiunea celui mai mic obiect pe care dorim să îl vizualizăm. În general, într-o aplicaţie dată, cel care impune o limită a numărului de pixeli este sistemul de achiziţie a imaginilor. De exemplu, chiar dacă se utilizează sisteme de achiziţie cu o rezoluţie ridicată, de 1000×1000 = 1 milion de elemente, rezoluţia spaţială relativă este de 10-3. Aceasta poate fi considerată o rezoluţie slabă, deoarece în cazul măsurării unei lungimi, a unei tensiuni electrice sau a unei frecvenţe, o rezoluţie sau o precizie satisfăcătoare începe de la 10-6. Însă, în cazul măsurării unor astfel de mărimi uni-dimensionale, se efectuează măsurători relativ la un singur punct, în timp ce o imagine de 1000×1000 conţine un milion de puncte. Prin urmare, o imagine poate oferi informaţii referitoare la variaţia spaţială a unui semnal. În plus, dacă se achiziţionează secvenţe temporale de imagini, se pot obţine informaţii care nu sunt accesibile dintr-o imagine statică. Astfel se pot obţine informaţii legate de variaţiile temporale ale unui semnal şi prin urmare se poate studia cinematica şi dinamica temporală a acestuia.
23
O reţea rectangulară reprezintă cea mai simplă, dar şi cea mai răspândită geometrie a unei imagini digitale. Pe lângă aceasta, mai există şi alte aranjamente geometrice ale pixelilor sau alte forme ale celulelor elementare. Aceste forme şi dispuneri geometrice sunt similare configuraţiilor cristaline posibile în cazul corpurilor solide 3D în fizică, chimie sau mineralogie. Dacă se iau în considerare doar poligoane regulate, există doar trei forme de reţele regulate posibile: triunghiulare, pătrate sau hexagonale (figura 2.9).
Figura 2.9. Forme de reţele regulate posibile în 2D: (a) - reţea triunghiulară; (b) - reţea pătrată; (c) – reţea hexagonală. În cazul imaginilor 3D, pixelul se transformă în voxel (din engleză = volume element). Într-o reţea rectangulară, fiecare pixel reprezintă valoarea medie a nivelului de gri (sau de culoare) dintr-un cub elementar. Poziţia unui voxel este indicată prin trei indici: un indice de linie (l), un indice de coloană (k) şi un indice (m) pentru ”adâncime” (figura 2.10).
24
z m x
l
y
k
Figura 2.10. Reprezentarea imaginilor digitale ca matrici de voxeli dispuşi într-o reţea rectangulară tri-dimensională.
2.4. Proprietăţi ale imaginilor digitale Imaginile digitale au unele proprietăţi, metrice sau topologice, diferite de proprietăţile funcţiilor bidimensionale continue. Pe baza celor prezentate până în acest punct, se pot trage următoarele concluzii: • o imagine digitală este formată din elemente de imagine (pixeli) de dimensiune finită; • în mod uzual, pixelii sunt aranjaţi sub forma unei reţele rectangulare; • o imagine digitală reprezintă o matrice bidimensională a cărui elemente sunt numere întregi care corespund nivelurilor de cuantizare a gamei de niveluri de gri. • unele proprietăţi ale imaginilor continue nu au o analogie directă în domeniul imaginilor digitale.
25
2.4.1. Proprietăţi metrice ale imaginilor digitale Distanţa dintre doi pixeli dintr-o imagine digitală reprezintă o mărime cantitativă. Distanţa dintre punctele de coordonate (i,j) şi (k,l) poate fi definită în diferite moduri: • distanţa euclidiană: d E [(i, j ), (k , l )] =
(i − k )2 + ( j − l )2
(2.21)
Avantajul distanţei euclidiene este faptul că este intuitivă, dar are dezavantajul unui cost mare de calcul datorită radicalului din formulă şi datorită valorii neîntregi care rezultă şi deci a interpolării necesare. Distanţa dintre două puncte poate fi exprimată şi prin numărul minim de paşi elementari de pe reţeaua discretă, dintre punctul de start şi punctul final. • Dacă sunt permise doar deplasări orizontale şi verticale, se poate defini distanţa d4 sau distanţa interbloc: d 4 [(i, j ), (k , l )] = i − k + j − l
(2.22)
Această distanţă este similară distanţei dintre două locaţii dintr-un oraş cu o reţea rectangulară de străzi şi blocuri închise de clădiri. • Dacă sunt permise şi deplasări diagonale, se poate defini distanţa d8 sau distanţa de tip şah: d 8 [(i, j ), (k , l )] = max{i − k , j − l }
26
(2.23)
2.4.2. Proprietăţi topologice ale imaginilor digitale Adiacenţa pixelilor este un concept important în prelucrarea imaginilor digitale. Oricare doi pixeli sunt vecini în sensul distanţei d4 dacă există o distanţă d4 =1 între cei doi pixeli. În mod analog, doi pixeli sunt vecini în sensul distanţei d8 dacă există o distanţă d8=1 între cei doi pixeli. Cele două tipuri de vecinătăţi sunt ilustrate în figura de mai jos:
V4
V8
Figura 2.11. Vecinătatea V4 şi V8. Pe baza adiacenţei pixelilor se pot defini regiunile, ca mulţimi conexe de pixeli adiacenţi. O cale dintre un pixel P şi Q este o secvenţă de puncte A1, A2, …, An, unde A1=P şi An=Q, iar Ai+1 este vecin cu Ai, i=1, …,n. O regiune reprezintă o mulţime de pixeli în care există o cale între oricare pereche de pixeli ai săi, iar pixelii acelei căi sunt incluşi şi ei în mulţimea respectivă. Dacă există o cale între doi pixeli ai unei imagini, aceşti pixeli sunt conecşi. Relaţia de conexitate este reflexivă, simetrică şi tranzitivă şi defineşte o descompunere a mulţimii (în cazul de faţă imaginea) în clase echivalente (regiuni). Să presupunem că Ri sunt regiuni disjuncte din imagine şi că aceste regiuni nu ating marginile imaginii (pentru a evita cazurile speciale). Fie R reuniunea tuturor regiunilor Ri. Fie RC complementara mulţimii R în raport cu imaginea. 27
Submulţimea lui RC care este conexă în raport cu marginile imaginii se numeşte fundal, iar restul mulţimii RC se numesc găuri. Dacă nu avem găuri într-o regiune, aceasta se numeşte regiune simplu conexă. O regiune cu găuri se numeşte regiune multi-conexă. Trebuie observat faptul că noţiunea de regiune implică doar proprietatea de conexitate. Regiunilor li se pot atribui proprietăţi secundare care îşi au originea în interpretarea imaginilor. Astfel, unele regiuni din imagine se numesc obiecte. Procesul prin care se determină care regiuni dintr-o imagine corespund fiecărui obiect se numeşte segmentarea imaginilor. De exemplu, nivelul de gri al unui pixel reprezintă o proprietate simplă care poate fi utilizată pentru a defini obiectele dintr-o imagine. Dacă un pixel are un nivel de gri mai mare decât anumite praguri predefinite, el aparţine unui anumit obiect. Toate punctele care satisfac această proprietate şi care sunt conexe, constituie un obiect. O gaură constă din punctele care nu aparţin unui obiect şi sunt înconjurate de obiecte. Toate celelalte obiecte constituie fundalul. Un exemplu îl constituie un text negru pe o pagină albă, în care literele reprezintă obiectele. Regiunile albe înconjurate de litere reprezintă găuri (de exemplu în interiorul literei O). Toate celelalte regiuni ale hârtiei reprezintă fundal.
2.4.3. Relaţii de vecinătate între pixeli Una din proprietăţile importante ale imaginilor discrete este reprezentată de relaţiile de vecinătate dintre pixeli, deoarece pe baza acestora se pot defini regiunile conexe şi obiectele. Într-o reţea rectangulară bidimensională se pot defini două tipuri de vecinătăţi ale pixelilor (figura 2.12.a şi 2.12.b ).
28
l-1, k-1 l-1, k l-1, k+1
l-1, k l, k-1
l, k
l, k-1
l, k+1
l, k
l, k+1
l+1, k-1 l+1, k l+1, k+1
l+1, k
Figura 2.12. Vecinătăţi definite pe o reţea rectangulară: (a) - vecinătatea V4; (b) - vecinătatea V8. Definirea acestor vecinătăţi se poate face şi pe baza unor relaţii matematice, dar pentru moment vor fi definite într-un mod simplu, pentru o mai bună înţelegere. Astfel, se poate spune despre doi pixeli că sunt vecini dacă au cel puţin o latură comună. În acest caz un pixel va avea 4 vecini, obţinându-se vecinătatea V4 (figura 2.12.a). Se poate defini şi o altă vecinătate, în cadrul căreia doi pixeli sunt vecini dacă au cel puţin un colţ comun. În acest caz, un pixel va avea 8 vecini, obţinându-se vecinătatea V8 (figura 2.12.b). Ambele tipuri de vecinătate sunt necesare pentru a defini obiectele şi regiunile conexe. Se spune despre o regiune (sau un obiect) că este conexă atunci când se poate ajunge de la un pixel la oricare alt pixel al regiunii, trecând doar de la un pixel vecin la altul. De exemplu, obiectul gri din figura 2.13 reprezintă un obiect în sensul unei vecinătăţi V8, dar este constituit din două obiecte în sensul unei vecinătăţi V4.
Figura 2.13. Regiunea gri reprezintă un obiect (sau regiune conexă) dacă se utilizează o vecinătate V8 şi două obiecte dacă se utilizează o vecinătate V4. 29
Acelaşi lucru se poate afirma şi despre fundalul alb din figura 2.13. Pentru a se putea face o distincţie clară între fundal şi obiectele din figură se poate defini o vecinătate V4 în cazul obiectelor şi o vecinătate V8 în cazul fundalului sau invers. Aceste complicaţii nu apar numai în cazul reţelelor rectangulare. În cazul unei reţele triunghiulare se poate defini o vecinătate V3 pentru pixelii care au în comun câte o latură şi o vecinătate V12 pentru pixelii care au în comun câte un colţ (figura 2.9). În cazul unei reţele hexagonale se poate defini numai o vecinătate V6 deoarece toţi pixelii care au în comun un colţ, au în comun şi o latură, iar pixelii care au în comun o latură, au în comun şi două colţuri. În ciuda acestor dezavantaje, reţelele hexagonale sunt utilizate în mod curent în prelucrarea imaginilor deşi sistemele de achiziţie a imaginilor generează, de regulă, imagini ai căror pixeli sunt dispuşi într-o reţea rectangulară. Motivul îl reprezintă dispunerea sub formă hexagonală a senzorilor din retina ochiului uman. În cazul tri-dimensional, relaţiile de vecinătate sunt mai complexe. În acest caz există trei moduri de definire a vecinătăţilor: voxeli cu feţe comune, cu laturi comune sau cu colţuri comune. În cazul unei reţele rectangulare, aceste enunţuri permit definirea unei vecinătăţi V6, V18, respectiv V26 (figura 2.14).
k-1
l
k l-1 k
l
m-1 k+1
l
m
l+1
Figura 2.14. Cele trei tipuri de vecinătăţi posibile într-o reţea cubică 3D: (a) - V6: voxeli cu feţe comune; (b) - V18: voxeli cu laturi comune; (c) - V26: voxeli cu colţuri comune.
30
Şi în acest caz trebuie definite două tipuri de vecinătăţi pentru obiecte şi pentru fundal, pentru a putea defini în mod corect regiunile conexe. Astfel, în cazul obiectelor se poate utiliza o vecinătate V6, iar în cazul fundalului se poate utiliza o vecinătate V26 sau invers.
2.4.4. Paradoxuri de conexitate Definiţia vecinătăţii şi conexităţii pe o reţea rectangulară creează unele paradoxuri. Exemplul 1. În figura următoare sunt reprezentate trei linii digitale cu pante de 45o, respectiv -45o.
Figura 2.15. Exemplu de paradox de conexitate a liniilor. Dacă se utilizează vecinătatea V4, liniile nu sunt conexe în fiecare punct al lor. Mai mult, apar şi conflicte în raport cu înţelegerea intuitivă a proprietăţilor liniilor. Astfel, două linii perpendiculare se intersectează într-un caz (stânga-jos), dar nu se intersectează în alt caz (dreapta-sus), deoarece nu au un punct comun.
31
Exemplul 2. În figura următoare este prezentat un alt paradox.
A
D
C
B
Figura 2.15. Exemplu de paradox de conexitate a curbelor sau regiunilor. Acest paradox este cunoscut în geometria euclidiană, unde fiecare curbă (sau regiune) închisă divide planul în două regiuni neconexe. Dacă imaginea este digitizată într-o reţea pătrată, utilizând vecinătatea V8, se poate trasa o linie din partea internă a unei curbe închise până în partea externă, care nu intersectează curba. Aceasta implică faptul că părţile interne şi externe ale curbei constituie o singură regiune conexă. Exemplul 3 (paradoxul conectivităţii). Dacă se presupune vecinătatea (conexitatea) V4, figura de mai sus
conţine patru regiuni separate conexe A, B, C şi D. A ∪ B sunt neconexe, la fel ca şi C ∪ D, ceea ce reprezintă o contradicţie topologică deoarece, în mod intuitiv, dacă A ∪ B sunt neconexe, ar trebui ca C ∪ D să fie conexe. Dacă se presupune vecinătatea V8, există două regiuni A ∪ B şi C ∪ D. Cele două mulţimi conţin în întregime căile AB şi CD, dar acestea se intersectează! O soluţie de eliminare a paradoxului conexităţii este de a trata obiectele utilizând vecinătatea V4 iar fundalul utilizând vecinătatea V8 sau invers.
Problemele prezentate sunt tipice reţelelor rectangulare. Reţelele hexagonale rezolvă o mare parte a acestor probleme dar au, la rândul lor, numeroase dezavantaje. Astfel, din motive de simplitate, majoritatea 32
dispozitivelor de digitizare utilizează o reţea rectangulară, în ciuda dezavantajelor şi paradoxurilor prezentate. O alternativă pentru eliminarea problemelor de vecinătate sau conexitate este de a utiliza topologia discretă, considerând familii de mulţimi de diferite dimensiuni. De exemplu, punctele (0-dimensionale) pot fi atribuite unor mulţimi care să conţină structuri de dimensiuni mai mari (ca de exemplu, mulţimi de pixeli), care permit eliminarea paradoxurilor expuse. Liniile (1-dimensionale) permit o definiţie precisă a muchiilor şi contururilor etc.
2.4.5. Alte proprietăţi topologice şi geometrice Frontiera unei regiuni R este o mulţime de pixeli din regiune, care au unul sau mai mulţi vecini în exteriorul regiunii. Această definiţie se referă la frontiera internă, pentru a o distinge de frontiera externă, care reprezintă frontiera fundalului (complementarei) regiunii. Muchia este o proprietate a unui pixel şi a vecinătăţii sale imediate, caracterizată de o amplitudine şi o direcţie. Direcţia unei muchii este perpendiculară pe direcţia gradientului care indică direcţia de variaţie a nivelului de gri din imagine. Frontiera este un concept global relativ la o regiune, în timp ce muchia exprimă o proprietate locală a funcţiei de variaţie a nivelului de gri dintr-o imagine. Cu toate acestea, între muchii şi frontiere există o legătură. Astfel, o posibilitate de a determina frontierele este de a concatena muchiile semnificative (punctele caracterizate de un gradient mare al funcţiei de variaţie a nivelului de gri). Proprietatea de a aparţine unei muchii este caracteristică unui pixel şi vecinilor săi. Uneori este avantajos să se utilizeze proprietăţi ale unor perechi de pixeli vecini. Astfel se poate introduce noţiunea de muchie compusă, ataşată fiecărui pixel, care exprimă relaţia sa cu cei 4 vecini. Direcţia muchiilor compuse este cea de creştere a nivelului de gri şi este 33
un multiplu de 90o, în timp ce amplitudinea sa reprezintă diferenţa absolută dintre nivelurile de gri ale perechilor relevante de pixeli. Pentru descrierea proprietăţilor geometrice ale obiectelor se utilizează contururi convexe. Un contur convex este cea mai mică regiune care conţine un obiect, astfel încât oricare două puncte ale regiunii pot fi unite printr-o linie dreaptă, toate punctele liniei aparţinând regiunii. Un obiect poate fi reprezentat printr-o colecţie a componentelor sale topologice. Mulţimile de puncte din interiorul contururilor convexe, care nu aparţin unui obiect, sunt numite deficit de convexitate.
2.5. Reprezentarea spectrală a imaginilor Reprezentarea spectrală a imaginilor este utilă în analiza spectrală a acestora. Analiza spectrală oferă informaţii despre modul de variaţie a unui semnal. De exemplu, un semnal unidimensional (1D) lent variabil are un spectru concentrat în jurul originii, în timp ce un semnal rapid variabil, are un spectru mai larg. f(t)
|F(ω)|
t ω f(t)
|F(ω)|
t ω
Figura 2.16. Ilustrarea spectrului unui semnal lent şi a unuia rapid variabil.
34
În cazul imaginilor (2D) se poate determina dacă are sau nu contururi multe, prin inspecţia spectrului său, pornind de la constatarea că variaţiile rapide (frecvenţele mari) corespund contururilor. Spectrul unui semnal (sau al unei imagini) se obţine prin transformata Fourier a acestuia. Transformata Fourier unidimensională se defineşte astfel: def . ∞
ℑ{ f (t )} = ∫ f (t ) ⋅ e − j⋅ω⋅t dt = F (ω ) , f : R → C
(2.24)
−∞
not .
Transformata Fourier ℑ{ f (t )} = F (ω ) se defineşte pentru funcţiile
f(t)∈L2, unde L2 este clasa semnalelor (funcţiilor) de energie finită, pentru care există transformată Fourier directă şi inversă, adică: ∞ 2 f (t ) ∈ L2 = f (t ) | E f = ∫ f (t ) dt < ∞ , −∞ unde Ef este energia funcţiei f.
(2.25)
Dacă: f∈L2, ⇒ ∃F (ω ) ⇒ ℑ : L2 → L2 este inversabilă. Transformata Fourier inversă se defineşte ca fiind: def .
ℑ−1{F (ω )} =
1 ∞ j ⋅ω ⋅t ∫ F (ω ) ⋅ e dω = f (t ) 2π −∞
În aceste relaţii: ω este pulsaţia, iar
35
ω este frecvenţa. 2π
(2.26)
2.5.1. Transformata Fourier bidimensională Definiţie: Se consideră funcţia bidimensională f(x,y), f: R→C, unde: def . 2 f ∈ L2 = f : R 2 → C | E f = ∫∫ f ( x, y ) dxdy < ∞ R2
(2.27)
Transformata Fourier bidimensională a funcţiei f se defineşte ca fiind: def . ∞ ∞
ℑ{ f ( x, y )} = ∫ ∫ f ( x, y ) ⋅ exp[− j ⋅ (u ⋅ x + v ⋅ y )]dxdy = F (u , v)
(2.28)
− ∞ −∞
unde: x,y sunt coordonate spaţiale, iar u,v sunt frecvenţe spaţiale. Dacă: f∈L2, ⇒ ∃F (ω) ⇒ ℑ : L2 → L2 este inversabilă. Transformata Fourier bidimensională inversă se defineşte ca fiind: def .
ℑ −1 {F (u , v)} =
1 4π 2
∫∫ F (u, v) ⋅ exp[ j ⋅ (u ⋅ x + v ⋅ y )]dudv = f ( x, y ) (2.29)
R2
2.5.2. Proprietăţile transformatei Fourier bidimensionale 1. Deplasarea semnalului: Dacă funcţiei unidimensionale 1D f(t) îi corespunde transformata Fourier F(ω) atunci funcţiei f(t-t0) îi corespunde: F1D
F1D
f (t ) ↔ F (ω ) ⇒ f (t ± t 0 ) ↔ F (ω ) ⋅ e ± j⋅ω⋅t0 În cazul bidimensional (2D):
36
(2.30)
F2 D
F2 D
f ( x, y ) ↔ F (u , v) ⇒ f ( x ± x0 , y ± y 0 ) ↔ F (u , v) ⋅ e ± j⋅(u⋅x0 +v⋅ y0 ) (2.31) Demonstraţie: ℑ{ f ( x − x0 , y − y 0 )} = ∫∫ f ( x − x0 , y − y 0 ) ⋅ e − j⋅(u⋅x +v⋅ y )dxdy R
(2.32)
2
Făcând schimbările de variabile: x-x0=x’, respectiv y-y0=y’, iacobianul corespunzător este:
∂x ∂x' I = det ∂y ∂x'
∂y ∂y ' =1 ∂y ∂y '
(2.33)
ℑ{ f ( x − x 0 , y − y 0 )} = ∫∫ f ( x' , y ' ) ⋅ exp[− j ⋅ (u ⋅ ( x'+ x 0 ) + v ⋅ ( y '+ y 0 ) )]⋅ I ⋅ dx ' dy ' = R2
= exp[− j (u ⋅ x 0 + v ⋅ y 0 )] ∫∫ f ( x' , y ' ) ⋅ exp[− j ⋅ (u ⋅ x'+v ⋅ y ')] ⋅ dx' dy ' = R2
= F (u , v) ⋅ exp[− j (u ⋅ x0 + v ⋅ y 0 )] q.e.d. 2. Deplasarea spectrului: Dacă funcţiei unidimensionale 1D f(t) îi corespunde spectrul (transformata Fourier) F(ω), atunci funcţiei
f (t ) ⋅ e ± j⋅ω0 ⋅t
îi
corespunde spectrul deplasat: F1D
F1D
f (t ) ↔ F (ω ) ⇒ f (t ) ⋅ e ± j⋅ω0 ⋅t ↔ F (ω m ω 0 )
(2.34)
În cazul bidimensional: F2 D
F2 D
f ( x, y ) ↔ F (u , v) ⇒ f ( x, y ) ⋅ e ± j⋅(u⋅x0 +v⋅ y0 ) ↔ F (u m u 0 , v m v0 ) (2.35) Demonstraţia este lăsată ca exerciţiu, aceasta fiind similară demonstraţiei proprietăţii 1, de deplasare a semnalului. 37
3. Scalarea semnalului: Dacă funcţiei f(x,y) îi corespunde transformata Fourier F(u,v), atunci funcţiei scalate f (ax, by ) îi corespunde transformata Fourier: F2 D
F2 D
f ( x, y ) ↔ F (u , v) ⇒ f (a ⋅ x, b ⋅ y ) ↔
1 u v ⋅ F , a⋅b a b
(2.36)
În cazul unidimensional, acest lucru poate fi ilustrat grafic astfel: f(t)
|F(ω)|
t
ω 1 ω ⋅ F a a
f(a·t) Pt. a<1 t
ω
Figura 2.17. Ilustrarea grafică a scalării semnalului. Demonstraţie: ℑ{ f (a ⋅ x, b ⋅ y )} = ∫∫ f (a ⋅ x, b ⋅ y ) ⋅ e − j⋅(u⋅x +v⋅ y ) dxdy R
(2.37)
2
Făcând schimbările de variabile: a·x=x’, respectiv b·y=y’, iacobianul corespunzător este: ∂x ∂x' I = det ∂y ∂x'
∂y 1 ∂y ' = ∂y a ⋅ b ∂y ' 38
(2.38)
y ' 1 x' u v ⋅ F , q.e.d. ℑ{ f ( a ⋅ x, b ⋅ y )} = ∫∫ f ( x ' , y ' ) ⋅ exp − j ⋅ u ⋅ + v ⋅ ⋅ I ⋅ dx ' dy ' = a b a b ⋅ 2 a b R
4. Liniaritatea: Dacă funcţiilor f, respectiv g le corespunde transformatele Fourier F, respectiv G, atunci funcţiei compuse α ⋅ f + β ⋅ g îi corespunde transformata Fourier compusă: F2 D
f , g ↔ F,G
F2 D
⇒ α ⋅ f + β ⋅ g ↔ α ⋅ F + β ⋅ G, ∀α , β ∈ C (2.39)
Demonstraţia este trivială şi se face similar proprietăţii de deplasare a semnalului. 5. Proprietatea de simetrie: Dacă funcţiei f(x,y) îi corespunde transformata Fourier F(u,v), atunci funcţiei simetrice faţă de origine f (− x,− y ) , îi corespunde un spectru (transformata Fourier) simetric faţă de origine: F2 D
F2 D
f ( x, y ) ↔ F (u , v) ⇒ f (− x,− y ) ↔ F (−u ,−v)
(2.40)
Similar, dacă funcţiei f(x,y) îi corespunde transformata Fourier F(u,v), atunci conjugatei funcţiei simetrice faţă de origine f (− x,− y ) , îi corespunde un spectru (transformata Fourier) simetric faţă de origine, dar rotit cu 1800: F
F
2D 2D f ( x, y ) ↔ F (u , v) ⇒ f * (− x,− y ) ↔ F * (u , v)
(2.41)
Dacă funcţia f este reală (f∈R2): f ( x, y ) = f * (− x,− y ) ⇒ F (−u ,−v) = F * (u , v)
39
(2.42)
6. Teorema convoluţiei: Fie funcţiile bidimensionale f , g : R 2 → C . Produsul de convoluţie al funcţiilor f şi g se defineşte astfel: def .
f ( x, y ) ∗ g ( x, y ) = ∫∫ f (ξ ,η ) ⋅ g ( x − ξ , y − η )dξdη = R
def .
= ∫∫ f ( x − ξ , y − η ) ⋅ g (ξ ,η )dξdη = g ( x, y ) ∗ f ( x, y ) R
(2.43)
2
(2.44)
2
Enunţul teoremei convoluţiei: Dacă funcţiei f(x,y) îi corespunde transformata Fourier F(u,v), iar funcţiei g(x,y) îi corespunde transformata Fourier G(u,v), atunci produsului de convoluţie a celor două funcţii îi corespunde produsul transformatelor Fourier ale celor două funcţii, iar produsului simplu a celor două funcţii îi corespunde produsul de convoluţie a transformatelor Fourier ale celor două funcţii: F2 D
f ( x, y ) ∗ g ( x, y ) ↔ F (u , v) ⋅ G (u , v) F2 D
f ( x, y ) ⋅ g ( x, y ) ↔
1 4π 2
F (u , v) ∗ G (u , v)
(2.45)
(2.46)
Demonstraţie: def .
ℑ{ f ( x, y ) ∗ g ( x, y )} = ∫∫ f ( x, y ) ∗ g ( x, y ) ⋅ e − j⋅(u⋅x +v⋅ y ) dxdy = R2
= ∫∫ ∫∫ f (ξ ,η ) ⋅ g ( x − ξ , y − η )dξdη ⋅ e − j⋅(u⋅x +v⋅ y ) dxdy = R 2 R 2 = ∫∫ f (ξ ,η ) ⋅ ∫∫ g ( x − ξ , y − η ) ⋅ e − j⋅(u⋅ x +v⋅ y ) dxdy dξdη = R 2 R2 1444444 424444444 3 G (u ,v )⋅e − j ⋅(ξ ⋅u +η ⋅v )
40
= ∫∫ G (u , v) ⋅ f (ξ ,η ) ⋅ e − j⋅(ξ ⋅u +η ⋅v ) dξdη = F (u , v) ⋅ G (u , v) q.e.d. R2
7. Teorema lui Parceval: Fie funcţiile bidimensionale f , g : R 2 → C . Dacă funcţiei f(x,y) îi corespunde transformata Fourier F(u,v), iar funcţiei g(x,y) îi corespunde transformata Fourier G(u,v), atunci produsului scalar a celor două funcţii îi corespunde produsul scalar al transformatelor 1 . Fourier a celor două funcţii, multiplicat cu o constantă 4π 2 Deci, dacă: F2 D
F2 D
f ( x, y ) ↔ F (u , v) şi g ( x, y ) ↔ G (u , v) ,
(2.47)
atunci: ⇒ f,g =
1
F,G
4π 2
(2.48)
adică: def
⇒ f , g = ∫∫ f ( x, y ) ⋅ g * ( x, y )dxdy = R2
1
* ∫∫ F (u, v) ⋅ G (u, v)dudv =
4π 2 R 2
1 4π 2
F,G
Demonstraţie: g ( x, y ) =
1 4π
g * ( x, y ) = def
2
∫∫ G (u, v) ⋅ e
2
* − j ( u ⋅ x + v⋅ y ) dudv ∫∫ G (u, v) ⋅ e
R2
⇒ f , g = ∫∫ f ( x, y ) ⋅ R2
dudv
R2
1 4π
j (u⋅ x + v⋅ y )
1 G * (u , v) ⋅ e − j (u⋅ x +v⋅ y ) dudv dxdy = ∫∫ 4π 2 R 2
41
= ∫∫ G * (u , v) ⋅ R2
1 f ( x, y ) ⋅ e − j (u⋅ x +v⋅ y ) dxdy dudv = ∫∫ 4π 2 R 2 144444244444 3 F ( u ,v )
=
1 4π 2
* ∫∫ F (u , v) ⋅ G (u , v)dudv = R
2
1 4π 2
F,G
q.e.d.
Dacă definim energia funcţiei f ca fiind: def .
2
E f = ∫∫ f ( x, y ) dxdy R
(2.49)
2
Teorema energiei (consecinţă a teoremei lui Parceval): Energia calculată în spaţiul original (primar) este egală cu energia calculată în domeniul spectral, multiplicată cu o constantă 1 . 4π 2 Această teoremă rezultă ca un caz particular din teorema lui Parceval, pentru g=f: (2.50) 1 1 ⇒ E f = ∫∫ f ( x, y ) ⋅ f * ( x, y )dxdy = F (u , v) ⋅ F * (u , v) dudv = ⋅ EF ∫∫ 1 4 4 2 4 4 3 1 4 4 2 4 4 3 2 2 2 2 4 4 π π R R 2 2 f
F
8. Teorema simetriei: Dacă funcţiei unidimensionale f(t) îi corespunde transformata Fourier F(ω), atunci transformatei Fourier privită ca funcţie de timp F(t) îi corespunde transformata Fourier simetrică multiplicată cu constanta 2π: F1D
f (t ) ↔ F (ω )
F1D
⇒ F (t ) ↔ 2π ⋅ f (−ω )
Un exemplu este prezentat în figura următoare:
42
(2.51)
f(t)
F(ω) F1D t
ω
F(t)
f(ω) F1D t
ω
Figura 2.18. Exemple ilustrative ale teoremei simetriei. În cazul bidimensional, dacă funcţiei bidimensionale f(x,y) îi corespunde transformata Fourier F(u,v), atunci transformatei Fourier privită ca funcţie de spaţiu F(x,y) îi corespunde transformata Fourier simetrică multiplicată cu constanta 4π2: F2 D
f ( x, y ) ↔ F (u , v)
F2 D
⇒ F ( x, y ) ↔ 4π 2 ⋅ f (−u ,−v)
(2.52)
2.5.3. Proprietăţi specifice transformatei Fourier bidimensionale 9. Separabilitatea: Transformata Fourier bidimensională este separabilă: ∞ ∞
F (u , v) = ∫ ∫ f ( x, y ) ⋅ exp[− j ⋅ (u ⋅ x + v ⋅ y )]dxdy = −∞ − ∞
∞ = ∫ ∫ f ( x, y ) ⋅ exp(− j ⋅ u ⋅ x )dx ⋅ exp(− j ⋅ v ⋅ y )dy = ∞ −∞ − 14444244443 Fx (u , y ) ∞
43
∞
= ∫ Fx (u , y ) ⋅ exp(− j ⋅ v ⋅ y )dy = Fxy (u , v) = F (u , v)
(2.53)
−∞
Din această proprietate rezultă că se poate face calculul transformatei Fourier bidimensionale aplicând pe rând (pe cele două direcţii x şi y) transformata Fourier unidimensională. Cu alte cuvinte, se aplică transformata Fourier unidimensională pe direcţia x, iar asupra rezultatului se aplică transformata Fourier unidimensională pe direcţia y: F
pe x
F
pe y
1D 1D f ( x, y ) → Fx (u , y ) → Fxy (u , v) = F (u , v) sau:
F
pe y
F
pe x
1D 1D f ( x, y ) → F y ( x, v) → F yx (u , v) = F (u , v) (2.54)
Proprietatea de separabilitate are următoarele consecinţe: - Dacă se dispune de un algoritm rapid de calcul pentru cazul unidimensional (iar pentru transformata Fourier există un astfel de algoritm), atunci şi pentru transformata Fourier bidimensională există un algoritm rapid de calcul. - Dacă funcţia originală se poate scrie ca produsul a două funcţii, transformata sa Fourier este egală cu produsul transformatelor Fourier a celor două funcţii, adică: Dacă: f ( x, y ) = f1 ( x) ⋅ f 2 ( y ) ⇒ F (u , v) = F1 (u ) ⋅ F2 (v),
F1(u) = ℑ{ f1 ( x)} unde : F2(v) = ℑ{ f 2 ( y )}
(2.55)
10. Derivarea spaţială: Dacă funcţiei bidimensionale f(x,y) îi corespunde transformata Fourier F(u,v), atunci derivatei funcţiei f în raport cu cele două variabile, îi corespunde următoarele transformate Fourier:
44
∂f F2 D ∂x ←→ j ⋅ u ⋅ F (u , v) ⇒ F2 D ∂f ← → j ⋅ v ⋅ F (u , v) ∂y
F
D f ( x, y ) ←2 → F (u , v)
(2.56)
Demonstraţie: def .
f ( x, y ) = ⇒
1
∫∫ F (u , v) ⋅ exp[ j ⋅ (u ⋅ x + v ⋅ y )]dudv (2.57)
4π 2
R2
1 ∂f = ∫∫ F (u, v) ⋅ j ⋅ u ⋅ exp[ j ⋅ (u ⋅ x + v ⋅ y )]dudv = j ⋅ u ⋅ F (u, v) q.e.d. ∂x 4π 2 R 2 Această proprietate are aplicaţii în calculul diferenţial, de exemplu la calculul laplaceanului: def . ∂ 2
f
∂x
2
∆f ( x, y ) =
+
∂2 f ∂y
(2.58)
2
Deoarece: ∂f FD ←2 → j ⋅ u ⋅ F (u , v) ∂x
În mod similar: ⇒
∂2 f ∂y
2
⇒
∂2 f ∂x 2
F
D ←2 →( j ⋅ u ) 2 ⋅ F (u , v) = −u 2 ⋅ F
F
D ←2 → = −v 2 ⋅ F
F
D ⇒ ∆f ( x, y ) ←2 → = −(u 2 + v 2 ) ⋅ F (u , v)
(2.59)
11. Integrarea spaţială: Dacă funcţiei bidimensionale f(x,y) îi corespunde transformata Fourier F(u,v), atunci integralei funcţiei f în raport cu cele două variabile, îi corespunde următoarele transformate Fourier:
45
F
D → F (u , v) f ( x, y ) ←2
∞ ℑ ∫ f ( x, y )dx ←→ F (0, v) − ∞ ⇒ ∞ ℑ f ( x, y )dy ←→ F (u ,0) ∫ − ∞
(2.60)
Demonstraţie: def . ∞ ∞
F (u , v) = ∫ ∫ f ( x, y ) ⋅ exp[− j ⋅ (u ⋅ x + v ⋅ y )]dxdy −∞ −∞
∞ ∞ ∞ ∞ ⇒ F (0, v) = ∫ ∫ f ( x, y ) ⋅ exp(− j ⋅ v ⋅ y )dxdy = ∫ ∫ f ( x, y )dx ⋅ exp(− j ⋅ v ⋅ y )dy −∞ − ∞ −∞ −∞ 1 44244 3 functie de y 14o4 444424444443
transformata Fourier a functiei de y
∞ ⇒ F (0, v) = ℑ ∫ f ( x, y )dx (2.61) − ∞ 12. Teorema rotaţiei: Dacă funcţiei unidimensionale f îi corespunde transformata Fourier F, atunci funcţiei rotite cu un unghi α, fα, îi corespunde o transformată Fourier rotită în acelaşi sens şi cu acelaşi unghi α. În cazul bidimensional (deci în cazul unei imagini 2D), dacă funcţiei bidimensionale f(x,y) îi corespunde transformata Fourier F(u,v), atunci funcţiei rotite cu un unghi α, fα(x,y), îi corespunde un spectru rotit în acelaşi sens şi cu acelaşi unghi α. Rotaţia (Rotα) conservă liniaritatea şi simetriile.
46
y y'
y
α x
x'
x Figura 2.19. Rotaţia unui segment de dreaptă.
Rotaţia Rotα : R 2 → R 2 se poate scrie: x' cos α = y ' − sin α
sin α x ⋅ cos α y
(2.62)
Prin urmare, rotaţia se mai poate scrie: x' = x ⋅ cos α + y ⋅ sin α y ' = − x ⋅ sin α + y ⋅ cos α
(2.63)
⇒ fα ( x, y ) = f ( x' , y ' ) = f ( x ⋅ cos α + y ⋅ sin α ,− x ⋅ sin α + y ⋅ cos α )
Demonstraţie: ℑ{ f α ( x, y )} = ∫∫ f ( x ⋅ cos α + y ⋅ sin α,− x ⋅ sin α + y ⋅ cos α) ⋅ e − j⋅(u⋅ x +v⋅ y ) dxdy R2
x' = x ⋅ cos α + y ⋅ sin α Se face schimbarea de variabile: y ' = − x ⋅ sin α + y ⋅ cos α
47
x cos(−α ) sin( −α ) x' cos α ⋅ = ⇒ = y − sin( −α ) cos(−α ) y ' sin α
− sin α x' ⋅ cos α y '
x = x'⋅ cos α − y '⋅ sin α ⇒ y ' = x'⋅ sin α + y '⋅ cos α
(2.64)
Înlocuind:
ℑ{ f α ( x, y )} = ∫∫ f ( x' , y ' ) ⋅ exp[− j (u ⋅ x'⋅ cos α − u ⋅ y '⋅ sin α + v ⋅ x'⋅ sin α + v ⋅ y '⋅ cos α )]dx' dy ' R2
Iacobianul este: ∂x ∂x' det ∂y ∂x'
∂x cos α ∂y ' = det ∂y sin α ∂y '
⇒ ℑ{ f α ( x, y )} = ∫∫ f ( x' , y ' ) ⋅ exp− R2
− sin α =1 cos α
j (u ⋅ cos α + v ⋅ sin α) ⋅ x'+(−u ⋅ sin α + v ⋅ cos α) ⋅ y ' dx' dy ' 3 144424443 14442444 u' v'
⇒ ℑ{ fα ( x, y )} = F (u ' , v' ) = Fα (u , v) , unde: u ' cos α = v' − sin α
sin α u ⋅ cos α v
(2.65)
Prin urmare, dacă funcţia f(x,y) este cu simetrie circulară, atunci şi transformata sa Fourier F(u,v) este cu simetrie circulară. Demonstraţie: Presupunem că f(x,y) este o funcţie cu simetrie circulară. În coordonate polare (fp), această proprietate se scrie: f ( x, y ) = f ( ρ ⋅ cos ϕ , ρ ⋅ sin ϕ ) = f p ( ρ , ϕ ) , (2.66) 48
x = ρ ⋅ cos ϕ unde: y = ρ ⋅ sin ϕ Deoarece funcţia f este cu simetrie circulară: ⇒ f ( x, y )
în coordonate
=
polare
f p (ρ, ϕ)
simetrie
=
circularã
f p (ρ)
(2.67)
Pornim de la relaţia de definiţie: def . ∞ ∞
F (u , v) = ∫ ∫ f ( x, y ) ⋅ exp[− j ⋅ (u ⋅ x + v ⋅ y )]dxdy
(2.68)
−∞ −∞
Se face schimbarea de variabile carteziene în coordonate polare: x = ρ ⋅ cos ϕ y = ρ ⋅ sin ϕ
ρ = x 2 + y 2 ⇒ y ϕ = arctg x
(2.69)
Prin această schimbare de variabile, planul real se transformă în coordonate polare într-o semibandă de înălţime 2π (pentru a se acoperi tot planul φ trebuie să ia valori între 0…2π, iar ρ trebuie să ia valori între 0…∞): y
φ
R2
în coordonate polare x
2π
…
ρ
ρ=0…∞ φ=0…2π
Figura 2.20. Domeniile de valori în diferite sisteme de coordonate. 49
Iacobianul este: ∂x ∂ρ det ∂y ∂ρ
∂x cos ϕ ∂ϕ = det ∂y sin ϕ ∂ϕ
− ρ ⋅ sin ϕ =ρ ρ ⋅ cos ϕ
(2.70)
∞ ∞
⇒ F (u , v) = ∫ ∫ ρ ⋅ f p ( ρ , ϕ ) ⋅ exp[− j ⋅ (u ⋅ ρ ⋅ cos ϕ + v ⋅ ρ ⋅ sin ϕ )] ⋅ dρ ⋅ dϕ −∞ −∞
Se fac notaţiile: u = r ⋅ cos θ v = r ⋅ sin θ
not .
⇒ F p (r , θ ) = F (r ⋅ cos θ , r ⋅ sin θ )
(2.71)
not. ∞ 2π
⇒ F p (r , θ) = ∫ ∫ ρ ⋅ f p (ρ, ϕ) ⋅ exp[− j ⋅ r ⋅ ρ(cos θ ⋅ cos ϕ + sin θ ⋅ sin ϕ)] ⋅ dρ ⋅ dϕ 0 0
Dacă funcţia f este cu simetrie circulară: ⇒ f p (ρ ,ϕ ) = f p (ρ )
(2.72)
∞ 2π ⇒ F p (r ,θ ) = ∫ ρ ⋅ f p ( ρ ) ⋅ ∫ ⋅ exp[− j ⋅ r ⋅ ρ ⋅ cos(ϕ − θ )]dϕ 0 1044444 42444444 3
⋅ dρ
2π
= ∫ e − j⋅r ⋅ρ ⋅cos ϕ dϕ = J 0 ( ρ ⋅r ) 0
deoarece integrandul este periodic dupã ϕ, cu perioada 2π ∞ 2π ⇒ F p (r ,θ ) = ∫ ρ ⋅ f p ( ρ ) ⋅ ∫ e − j⋅r⋅ρ ⋅cos ϕ dϕ 0 0
dρ = F p (r )
(2.73)
Prin urmare, dacă funcţia f este cu simetrie circulară, atunci şi transformata sa Fourier F este cu simetrie circulară (q.e.d).
50
În plus, se ştie că: 2π
J 0 ( x) = ∫ exp(− j ⋅ x ⋅ cos ϕ )dϕ = funcţia Bessel de ordinul 0 0
∞
⇒ F p (r ) = ∫ ρ ⋅ f p ( ρ ) ⋅ J 0 ( ρ ⋅ r )dρ = transformata Henkel 0
def . ∞
H { f ( ρ )} = ∫ ρ ⋅ f p ( ρ ) ⋅ J 0 ( ρ ⋅ r )dρ = F (r ) = transformata 0
Henkel a unei funcţii de o singură variabilă f(ρ).
51
3. Îmbunătăţirea imaginilor La concepţia algoritmilor sau a dispozitivelor de prelucrare şi îmbunătăţire a imaginilor trebuie luat în considerare principiul percepţiei vizuale umane. Printre parametrii psiho-fizici ai percepţiei vizuale umane pot fi amintiţi: contrastul, contururile, forma, textura, culoarea etc. Percepţia umană a unei imagini poate provoca multe iluzii, înţelegerea lor furnizând explicaţii referitoare la mecanismele vederii umane şi artificiale.
3.1. Calitatea unei imagini O imagine poate fi degradată pe parcursul achiziţiei, transmisiei sau prelucrării sale. Pentru a estima această degradare se pot utiliza măsuri de calitate a imaginii. Calitatea necesară pentru o imagine depinde de scopul în care este utilizată imaginea. Metodele de apreciere a calităţii imaginii pot fi împărţite în două categorii: subiective şi obiective. Calitatea imaginii f(x,y) este estimată prin compararea cu o imagine de referinţă g(x,y). Imaginea de referinţă utilizată în acest scop este, de regulă, o imagine de sinteză. Una din clasele de metode cele mai utilizate se bazează pe diferenţa medie pătratică MSE: 2 ∑ ( g ( x, y ) − f ( x, y ) )
(3.1)
x, y
Problema acestei măsuri este că nu este posibilă distincţia între câteva diferenţe mari şi multe diferenţe mici. În locul diferenţei medii pătratice se poate utiliza eroarea medie absolută. O altă alternativă este corelaţia dintre imaginile f şi g. O măsură a degradării imaginii este reprezentată de raportul semnal-zgomot SNR. Fie f(x,y) imaginea originală şi f'(x,y)= f(x,y)+z(x,y) 52
imaginea degradată. Măsura degradării este estimată prin raportul dintre energia semnalului şi energia zgomotului, care este estimat prin relaţia:
SNR ( f ' , f ) = 10 ⋅ log10
(
∑ x , y f ( x, y ) 2 2
∑ x , y f ( x, y ) − f ' ( x , y )
)
2 2
(dB)
(3.2)
Se poate defini şi valoarea de vârf a raportului semnal-zgomot PSNR: PSNR ( f , f ' ) = 10 ⋅ log10
N ⋅ max x, y f ( x, y ) 2
(
2
∑ x , y f ( x, y ) − f ' ( x, y )
)
2 2
(dB) (3.3)
unde N este numărul de pixeli. Un raport PSNR mai mare de 32 dB corespunde unei degradări invizibile. Zgomotul care poate să apară la achiziţia, transmiterea sau prelucrarea imaginilor, poate fi dependent sau independent de conţinutul imaginii. Zgomotul este descris, în general, de caracteristicile sale probabilistice. Zgomotul alb are un spectru de putere constant, iar intensitatea sa nu se modifică odată cu frecvenţa. Acest tip de zgomot se utilizează în majoritatea cazurilor ca aproximare brută a zgomotului dintr-o imagine. Funcţia sa de auto-corelaţie este funcţia delta. Prin urmare, valorile zgomotului în doi pixeli diferiţi sunt necorelate. Avantajul acestui model de zgomot este că permite simplificarea calculelor. Un caz special de zgomot îl reprezintă zgomotul Gaussian. Zgomotul Gaussian reprezintă o aproximare foarte bună a zgomotului care intervine în majoritatea cazurilor. Probabilitatea de densitate a variabilei aleatoare ce descrie zgomotul Gaussian este dată de funcţia lui Gauss. Zgomotul Gaussian unidimensional 1D este caracterizat de media sa µ şi de deviaţia standard σ a variabilei aleatoare: p( x) =
1 2πσ
2
⋅e
−
( x− µ )2 2σ 2
53
(3.4)
Zgomotul poate fi: • zgomot aditiv, în cazul în care zgomotul η şi semnalul de imagine f sunt independente: (3.5) f’(x,y) = f(x,y) + η(x,y) În timpul transmisiei zgomotul este, în general, independent de semnalul de imagine. Prin urmare, degradarea sa poate fi modelată ca un zgomot aditiv. • zgomotul multiplicativ este o funcţie descrisă de relaţia: f'(x,y) ≈ f(x,y) · η(x,y)
(3.6)
• zgomotul impulsiv (de tip impuls) corespunde unei degradări a imaginii cu pixeli „zgomotoşi” a căror valoare diferă semnificativ de cea a pixelilor din vecinătatea lor. • zgomotul de tip „sare şi piper” este utilizat pentru a descrie zgomotul impulsiv saturat, care corespunde unei imagini degradate cu pixeli albi şi/sau negri, de exemplu. Un parametru important în aprecierea calităţii unei imagini îl constituie contrastul. Contrastul reprezintă variaţia locală a nivelului de gri şi se defineşte ca raport între nivelul mediu de gri al unui obiect şi cel al fundalului. Ochiul uman este logaritmic sensibil la iluminare şi la variaţii ale nivelurilor de gri. Acesta este motivul pentru care majoritatea monitoarelor au implementată o corecţie de tip gamma. Nivelul de gri aparent depinde foarte mult de nivelul local de gri al fundalului. Acest efect este numit contrast condiţional. Datorită acestui efect, percepţia vizuală a unor obiecte cu acelaşi nivel de gri poate fi diferită dacă acestea sunt plasate pe un fundal de culoare închisă sau deschisă.
54
3.2. Tehnici de îmbunătăţire a imaginilor Îmbunătăţirea imaginilor constă dintr-un ansamblu de tehnici de prelucrare care au ca scop scoaterea în evidenţă a anumitor caracteristici a imaginilor (de exemplu muchii sau contururi) sau eliminarea zgomotului, scopul final fiind obţinerea unei vizibilităţi superioare a componentelor imaginii. În general, termenul de îmbunătăţire este strâns legat de percepţia vizuală subiectivă a unui expert uman, considerat utilizatorul final al imaginii. Întrucât nu se pot defini standarde de calitate a imaginilor, calitatea imaginii este un criteriu subiectiv. Cei care pot face afirmaţii cu privire la calitatea unor imagini sunt experţii din domeniile din care provin imaginile. În plus, se poate afirma că îmbunătăţirea imaginilor este bine să fie interactivă şi iterativă deoarece utilizatorul poate interveni în permanenţă asupra calităţii imaginii şi fiecare utilizator o va face într-un mod caracteristic. Tehnicile de îmbunătăţire a imaginilor nu generează informaţie suplimentară despre imaginea originală, ci doar o pune pe cea existentă sub o altă formă, mai uşor de interpretat de către utilizator. Chiar şi o imagine originală, nedegradată, poate fi îmbunătăţită, obţinând o imagine modificată, dar subiectiv preferabilă. De exemplu, într-o imagine subexpusă sau supraexpusă, utilizatorul (uman sau dispozitivul tehnic) poate să nu distingă între două niveluri de luminanţă care diferă cu o cuantă; acestea sunt valori diferite în semnalul din calculator şi prin tehnici de îmbunătăţire a imaginii pot fi făcute să difere mult mai mult, astfel încât să fie depăşit pragul de sesizare a diferenţei.
55
Operatorii de îmbunătăţire a imaginilor pot fi împărţiţi în trei mari categorii: • operatori punctuali, prin care se realizează o relaţie de corespondenţă punctuală între valoarea originală a fiecărui pixel şi valoarea sa după transformare; • operatori spaţiali (locali sau de vecinătate), la care noua valoare a nivelului de gri a unui pixel se obţine din valoarea originală a pixelului respectiv şi din valorile originale ale unor pixeli din vecinătatea acestuia; • operatori integrali, în cazul cărora valoarea nouă a unui pixel depinde de valorile tuturor pixelilor din imaginea originală, obţinându-se printr-o transformare integrală a acestora. Pentru a exemplifica operaţiile de îmbunătăţire a imaginilor, se vor considera imagini de dimensiuni L×K (cu L linii şi K coloane) şi se va nota cu U imaginea iniţială şi cu V imaginea îmbunătăţită, rezultată în urma aplicării unei operaţii sau transformări de îmbunătăţire (T) asupra imaginii iniţiale: U = {U (l , k ) | 1 ≤ l ≤ L,1 ≤ k ≤ K } ,
(3.7)
V = {V (l , k ) | 1 ≤ l ≤ L,1 ≤ k ≤ K }
(3.8)
În figura 3.1 se observă că imaginea îmbunătăţită are aceleaşi dimensiuni ca şi imaginea originală.
56
K
V=T(U)
K
V
U L
L
Figura 3.1. Operaţia de îmbunătăţire a unei imagini.
3.3. Operatori punctuali de îmbunătăţire a imaginilor Operatorii punctuali de îmbunătăţire a imaginilor sunt transformări aplicate asupra nivelurilor de gri, a căror rezultat depinde doar de valoarea din pixelul considerat. Operatorii punctuali sunt definiţi prin relaţii prin care se realizează asocieri între valoarea originală a fiecărui pixel şi valoarea sa după transformare. K K v (l,k)=T(u(l,k)) v(l,k)
u(l,k) L
U
L
V
Figura 3.2. Operaţia de îmbunătăţire a unei imagini, cu operatori punctuali. Operatorii punctuali de îmbunătăţire a imaginilor pot fi împărţiţi în: 1) operatori de modificare a contrastului (engl. contrast streching); 2) transformări de decupare (engl. clipping, slicing, thresholding); 3) operatori de modificare a histogramei. 57
3.3.1. Operatori punctuali de modificare a contrastului Operaţiile de modificare a contrastului urmăresc mărirea sau micşorarea intervalului de niveluri de gri ocupat de anumite componente ale imaginii, păstrând acelaşi număr total de niveluri de gri (N).
Negativarea imaginii Cea mai simplă operaţie de modificare a contrastului este negativarea imaginii, definită de ecuaţia: def .
v(l , k ) = N − u (l , k ) ,
(3.9)
unde N este numărul de niveluri de cuantizare (de gri). Contrastul relativ perceput de un observator uman este modificat, ca urmare a diferenţei de sensibilitate între percepţia nuanţelor întunecate şi luminoase. Exemplul cel mai simplu de aplicare este percepţia unei radiografii de către un observator nespecialist: contrastul va fi apreciat ca mult mai bun pentru imaginea negativată, în care avem obiecte de interes negre pe fond alb.
Diferenţa între imagini Această operaţie poate fi definită prin relaţia: def .
v(l , k ) = u 2 (l , k ) − u1 (l , k )
58
(3.10)
Pentru obţinerea unui rezultat cât mai util, ar trebui ca imaginile U1 şi U2 să reprezinte aproximativ acelaşi lucru, dar în alte ipostaze (de exemplu un obiect în mişcare). Printre domeniile de aplicaţii se poate menţiona angiografia (grafia vaselor de sânge). În acest scop se achiziţionează o radiografie a pacientului în stare normală, după care se injectează în vasele sanguine o substanţă contrastantă în raze X şi se achiziţionează o nouă radiografie. Prin compararea şi diferenţa celor două radiografii se scot în evidenţă zonele de interes, potenţial afectate de anumite boli. Această operaţie poate fi utilizată şi pentru detecţia mişcării în secvenţe de imagini.
Modificarea liniară a contrastului Cea mai răspândită tehnică de modificare a contrastului (engl. contrast streching) este transformarea liniară pe porţiuni. Expresia analitică a acesteia este: v1 u , 0 ≤ u ≤ u1 u1 v − v1 v = f (u ) = v1 + 2 ⋅ (u − u1 ), u1 ≤ u ≤ u 2 u 2 − u1 N − 1 − v2 ⋅ (u − u 2 ), u 2 ≤ u ≤ N − 1 v 2 + N − 1 − u2
(3.11)
v N − 1 − v2 v −v , vor determina unde pantele α = 1 , β = 2 1 şi γ = u1 u 2 − u1 N − 1 − u2 variaţiile relative de contrast (figura 3.3).
59
v N-1
γ
v2
β v1 O
α
u1 u2
N-1
u
Figura 3.3. Modificarea contrastului. Astfel se vor obţine niveluri de gri cuprinse între v1 şi v2 pentru valori iniţiale cuprinse între u1 şi u2. Dacă u2-u1 < v2-v1, se va obţine o imagine cu un contrast mai mare, iar dacă u2-u1 > v2-v1 se va obţine o imagine cu un contrast mai slab pentru intervalul central al gamei de niveluri de gri. Trebuie avut în vedere faptul că prin modificarea contrastului dintr-o regiune se modifică contrastul şi în celelalte regiuni. De exemplu, prin mărirea contrastului în regiunea centrală (u1, u2), contrastul scade în celelalte regiuni (0÷u1, u2÷(N-1)). Se poate realiza şi operaţia inversă, adică dintr-o imagine cu contrast puternic să se obţină o imagine cu contrast mai slab. Din cazul general al modificării de contrast se pot obţine câteva cazuri particulare de interes. Unul dintre acestea este determinat de particularizarea v1=0 şi v2=N-1, şi constă în eliminarea extremelor şi extinderea maximă a intervalului de niveluri de gri de interes (figura 3.4).
60
v N-1
O
u1
u2
N-1
u
Figura 3.4. Extinderea nivelurilor de gri. Expresia analitică corespunzătoare este:
0, 0 ≤ u ≤ u1 N −1 (u − u1 ), u1 ≤ u ≤ u 2 v = f (u ) = u 2 − u1 N − 1, u 2 ≤ u ≤ N − 1
(3.12)
Tot din cazul general al modificării de contrast se mai pot obţine şi alte cazuri particulare, cum ar fi transformările de decupare. De exemplu, dacă intervalul central de niveluri de gri este eliminat (u1 =u2), din cazul prezentat anterior se obţine transformarea de binarizare a imaginilor sau segmentarea cu prag ("tresholding", figura 3.5) . u N-1
O
u1= u2=T
N-1
u
Figura 3.5. Binarizarea imaginilor.
61
Variantele neliniare de modificare a contrastului sunt compresia (sau compandarea) şi inversa acesteia, expandarea. Prin compresie se obţine o variaţie maximă a contrastului în zona nivelurilor de gri apropiate de 0, iar prin expandare, variaţia maximă a contrastului se obţine în zona nivelurilor de gri apropiate de N-1. Compresia logaritmică (figura 3.6) este descrisă de ecuaţia: N −1 v= ⋅ lg(1 + u ), u,v ∈ [0 ,N-1] (3.13) lg N Expandarea este descrisă de ecuaţia: v=
eu − 1 eN −1
(N − 1)
, u,v ∈ [0 ,N-1]
(3.14)
v N-1
O
N-1
u
Figura 3.6. Expandarea şi compresia.
3.3.2. Decuparea intervalelor de niveluri de gri Tehnicile de decupare a intervalelor de niveluri de gri urmăresc punerea în evidenţă numai a unei porţiuni din gama totală a nivelurilor de gri disponibile (sau ocupate efectiv de pixelii imaginii). Această punere în evidenţă este realizată în principiu prin înlocuirea tuturor celorlalte niveluri de gri cu o valoare constantă. Transformarea de "clipping" păstrează nemodificat un interval de niveluri de gri de interes (de exemplu u1…u2), restul nivelurilor de gri 62
fiind transformate într-o valoare unică, numită fundal (F). Expresia analitică a acestui operator este:
F , 0 ≤ u ≤ u1 v = f (u ) = u , u1 ≤ u ≤ u 2 F , u ≤ u ≤ N − 1 2
(3.15)
v N-1 v2 v1 F O
u1
u2
N-1
u
Figura 3.7. Decuparea nivelurilor de gri. Transformarea de decupare ("slicing") pune în evidenţă un interval de niveluri de gri prin modificarea valorilor nivelurilor de gri la 0 sau N-1, după cum acestea se situează în afara sau respectiv în interiorul intervalului considerat (figura 3.8). v N-1
O
u1
u2
N-1
u
Figura 3.8. Slicing-ul nivelurilor de gri.
63
3.3.3. Modificarea histogramei Histograma unei imagini este o funcţie ce pune în evidenţă conţinutul de niveluri de gri al acesteia. Din punct de vedere matematic, histograma se defineşte ca frecvenţa relativă de apariţie în imagine a diferitelor niveluri de gri. Dacă considerăm o imagine f, de dimensiune L×K pixeli şi notăm cu u un nivel de gri şi cu δ impulsul unitar, histograma se exprimă ca: h(u ) =
1 L −1, K −1 ∑ ∑ δ ( f (l , k ) − u ), u ∈ [0, N − 1] L ⋅ K l = 0, k = 0
(3.16)
Proprietăţile imaginii influenţează forma histogramei sale. O imagine de tip tablă de şah, formată din pătrate luminoase şi întunecate, în proporţie relativ egală, va avea o histogramă prezentând două maxime puternice, localizate în jurul valorilor de 0 şi N-1, şi valori aproape nule în zona nivelurilor de gri medii. O imagine fotografică subexpusă (şi deci foarte întunecată) are o histogramă al cărei suport (interval ce corespunde valorilor nenule) este concentrat spre valoarea 0. O imagine fotografică supraexpusă (şi deci foarte luminoasă) are o histogramă al cărei suport este situat în zona valorilor apropiat de N-1. O imagine bine contrastată, ce prezintă numeroase nuanţe, va avea o histogramă al cărei suport va acoperi aproape întreaga gamă de niveluri de gri posibile şi a cărei formă va fi neregulată. Reciproc, inspecţia formei unei histograme poate oferi informaţii despre caracteristicile imaginii, dar nu o poate individualiza, deoarece mai multe imagini pot avea aceeaşi histogramă. Această comportare corespunde faptului că histograma se comportă ca o funcţie de densitate de probabilitate a unei variabile aleatoare ale cărei realizări particulare sunt valorile nivelurilor de gri din imagine.
64
N −1
Într-adevăr, h(u)>0, ∀u, şi ∑ h(n) = 1 . Orice funcţie de densitate de n =0
probabilitate are asociată o funcţie de repartiţie. În cazul histogramei imaginilor, această funcţie de repartiţie este histograma cumulativă: u
H (u ) = ∑ h(n), u ∈ [0; N − 1] n =o
(3.17)
Imaginea ideală ar trebui să prezinte o distribuţie uniformă a nivelurilor de gri şi un contrast repartizat regulat în întreaga gamă dinamică. Pentru a obţine o asemenea imagine, operatorul de îmbunătăţire trebuie să transforme histograma originală a imaginii într-o histogramă uniformă, în care toate nivelurile de gri sunt egal probabile. Din punct de vedere matematic, problema se reduce la a transforma o funcţie de densitate de probabilitate oarecare într-o funcţie de densitate de probabilitate uniformă (constantă pe intervalul de definiţie [0,N-1]). Ţinând cont de teoria variabilelor aleatoare (funcţie de o variabilă aleatoare) şi de faptul că variabila aleatoare “nivel de gri” este discretă, formula de transformare a nivelului de gri u pentru egalizarea de histogramă este: H (u ) − H (0) (3.18) v = int ⋅ ( N − 1) + 0,5 1 − H (0) unde ”int[ ]” este operatorul parte întreagă, iar H este histograma cumulativă definită anterior.
Exemplu: Considerând o imagine de 64×64 pixeli, reprezentată cu 8 niveluri de gri (0…7), a cărei histogramă este dată în tabelul următor (căruia îi corespunde figura 3.9), să se realizeze egalizarea de histogramă.
65
Nr. de pixeli
Nivelul Nr. de pixeli de gri având acest nivel de gri 0 796 1 1023 2 850 3 650 4 329 5 245 6 122 7 81
1200 1000 800 600 400 200 0 1
2
3
4
5
6
7
8
nivelul de gri
Figura 3.9. Histograma unei imagini. Rezultatele obţinute în urma egalizării, pe baza relaţiilor de definiţie sunt cumulate în următorul tabel: Tabelul 3.1. Imaginea iniţială Imaginea transformată Nivel Nr. de h(i) Histograma Nivelul de Numărul de Transformarea nivelului de pixeli gri cumulativă de pixeli cu gri pt. nivelul transformat corespunzător gri nivelul nivelului de conform de gri i de gri i gri egalizării de histogramă transformat 0 796 0,194 0,194 0 796 0-0 1 1023 0,249 0,443 2 0 1-2 2 850 0,208 0,651 4 1023 2-4 3 650 0,159 0,81 5 0 3-5 4 329 0,08 0,89 6 850 4-6 5 245 0,06 0,95 7 650 5-7 6 122 0,031 0,981 7 329 6-7 7 81 0,019 1 7 448 7-7 66
nr de pixeli
Histograma "egalizată" este deci (figura 3.10): 1200 1000 800 600 400 200 0 1 0
21
32
43
54
65
76
78
Figura 3.10. Histograma egalizată. Din graficul prezentat se observă că histograma obţinută nu este "uniformă" şi prezintă numeroase niveluri de gri lipsă ("găuri"). Aceste efecte sunt datorate în general cuantizării nivelurilor de gri şi limitării prin trunchiere a domeniului de variaţie a valorilor (formula de transformare este dedusă pentru variabile aleatoare cu variaţie continuă). Pentru corectare au fost propuse mai multe abordări: limitarea maximelor histogramei, mutarea aleatoare a valorilor pixelilor situate pe niveluri de gri mai bine reprezentate în histogramă pe niveluri de gri absente, etc. Trebuie de asemenea remarcat faptul că egalizarea de histogramă nu asigură în toate cazurile cea mai bună calitate vizuală a imaginii transformate. Egalizarea de histogramă şi tehnicile înrudite de specificare a histogramei asigură mărirea contrastului imaginii prin redistribuirea nivelurilor de gri în cadrul gamei dinamice fixate, [0,N-1]. Sensibilitatea sistemului vizual uman este însă mult mai mare în gama color decât în cea a nivelurilor de gri. De aceea una dintre tehnicile cele mai populare de realizare a unei vizibilităţi maxime a anumitor componente dintr-o imagine este colorarea lor cu culori puternic contrastante, adică prin pseudocolorare. În cazul aplicării tehnicii de pseudocolorare, imaginea va fi afişată (vizualizată) cu o tabelă de culoare diferită de paleta originală de niveluri 67
de gri. Această nouă paletă de culoare poate fi construită după orice fel de reguli care să corespundă problemei de rezolvat: de exemplu, toţi pixelii al căror nivel de gri este 250 vor fi afişaţi cu roşu şi toţi pixelii al căror nivel de gri este cuprins între 100 şi 120 vor fi afişaţi cu verde. Se pot introduce şi condiţii relative la poziţia spaţială a pixelilor sau la alte caracteristici locale ale acestora. Schema generală a unei operaţii de pseudocolorare este detaliată în figura 3.10. index R G B Imagine iniţială
Bloc de extragere caracteristici Paletă de culori
Display
Figura 3.10. Schema unui sistem de pseudocolorare.
3.4. Operatori liniari de vecinătate pentru îmbunătăţirea imaginilor. Filtrarea liniară a imaginilor Spre deosebire de operatorii punctuali, operatorii de vecinătate (numiţi şi operatori spaţiali locali) determină valoarea nouă a unui pixel ca o funcţie de valorile pixelilor dintr-o vecinătate a sa. Dacă această funcţie (de mai multe variabile) este liniară, atunci operatorul se numeşte liniar. Se va nota cu u(l,k) imaginea iniţială (de intrare), cu l=1,…,L, k=1,…,K, cu v(l,k) imaginea rezultată (de ieşire) şi cu a(i,j), funcţia pondere a sistemului, i=1,…L, j=1,…,K. Operatorii liniari de vecinătate se implementează, în general, prin convoluţia imaginii iniţiale cu funcţia pondere a unui filtru cu răspuns finit, numit mască spaţială, adică prin aplicarea unui asemenea operator se realizează practic o filtrare bidimensională. Această filtrare este obţinută 68
prin tehnica "ferestrei glisante" (moving-window), iar fereastra ce culisează peste imagine se mai numeşte şi mască spaţială sau filtru bidimensional şi are rolul de a selecta vecinătatea pixelului curent asupra căruia operează filtrul respectiv. Expresia analitică a acestei operaţii este: v(l , k ) =
∑ ∑ a(i, j ) ⋅ u (l − i, k − k )
(i , j )∈W
(3.19)
unde u şi v sunt imaginile de intrare, respectiv de ieşire, iar a sunt coeficienţii ferestrei (măştii) de filtrare W, care are dimensiuni mai mici decât imaginile asupra cărora acţionează. După cum se observă, aceasta este de fapt o convoluţie: v=a*u, între funcţia pondere a unui filtru cu răspuns finit (a) şi imaginea de intrare, iniţială (u). Aceasta se poate exprima ca un produs punct cu punct între coeficienţii măştii şi o porţiune din imagine, de aceeaşi dimensiune. Tehnica ferestrei glisante constă în efectuarea următoarelor operaţii: • Se defineşte fereastra de filtrare adică: forma (relativ la o origine); coeficienţii din fiecare punct; • Fereastra glisează peste imaginea iniţială, adică se pune originea ferestrei în fiecare punct al imaginii. Astfel va fi selectat de către fereastră pixelul curent şi pixelii din vecinătatea acestuia. • Se face produsul punct cu punct între valoarea pixelilor din imagine selectaţi de fereastră şi coeficienţii ferestrei. • Se înlocuieşte pixelul curent cu noua valoare obţinută, ca sumă a produselor obţinute la punctul anterior. De regulă, filtrele folosite sunt de ordin impar, rectangulare şi au originea (0,0) în centrul suportului:
69
a-1-1
a-10
a-11
a0-1
a00
a01
a1-1
a10
a11
Problemele care se ridică la aplicarea acestor operatori se referă la: • Marginile imaginii, aici având 2 posibilităţi: obţinerea unei imagini de dimensiuni mai mici, atunci când glisarea începe din interiorul imaginii bordarea imaginii de intrare pentru a păstra aceleaşi dimensiuni pentru imaginea prelucrată • Numărul de operaţii necesare pentru fiecare punct. Pentru o fereastră pătrată cu latura n sunt necesare: n2 înmulţiri şi n2-1 adunări. În concluzie, din punct de vedere al volumului de calcule, este mai bine să se lucreze cu ferestre cât mai mici. Filtrele bidimensionale folosite uzual în prelucrarea imaginilor sunt nuclee cu suporturi de dimensiuni mici: 3×3, 5×5. Filtrele de dimensiuni mai mari se pot reduce adesea la aplicarea repetată asupra unei imagini a unor nuclee de dimensiuni mai mici. Nucleele pătrate de dimensiuni 3×3 sunt cele mai utilizate. Exemple de astfel de nuclee sunt: 1 1 1 1 ⋅ 1 1 1 9 1 1 1
0 1 0 1 ⋅ 1 4 1 8 0 1 0
(a)
(3.20)
(b)
Filtrul (a) realizează media între pixelul central şi vecinii săi, iar filtrul (b) realizează media ponderată între pixelul central şi vecinii săi verticali şi orizontali.
70
Forma şi coeficienţii ferestrei se aleg astfel încât să corespundă aplicaţiei concrete. Singura constrângere în ceea ce priveşte coeficienţii ferestrei, pentru filtrele de mediere (al căror efect este de FTJ) este:
∑ ∑ a(i, j ) = 1 , pentru a nu modifica regiunile uniforme. (3.21) (i , j )∈w
În cazul în care coeficienţii ferestrei îndeplinesc condiţia din relaţia (3.21) efectul filtrului (ferestrei de filtrare) este un efect de netezire. În cazul în care coeficienţii ferestrei îndeplinesc relaţia:
∑ ∑ a(i, j ) = 0
(3.22)
(i , j )∈w
efectul filtrului (ferestrei de filtrare) este un efect de reliefare (accentuare sau contrastare), respectiv de filtru trece-sus. Exemple de filtre de reliefare sunt: 0 −1 0 − 1 [4] − 1 0 − 1 0
− 1 − 1 − 1 − 1 [8] − 1 − 1 − 1 − 1
1 −2 1 − 2 [4] − 2 1 − 2 1
(3.23)
Medierea efectuată cu ajutorul operatorilor liniari de vecinătate poate fi utilă la reducerea zgomotului aditiv gaussian şi mai puţin a zgomotului impulsiv, de tip salt and pepper (engl.= "sare şi piper"). Aplicarea unui filtru de mediere asupra unei imagini afectată de zgomot impulsiv sau gaussian are ca rezultat extinderea punctelor cu zgomot (formarea de pete), respectiv apariţia efectului de bluring sau mânjeală (neclar, ceţos), ca efect special dorit sau ca o consecinţă nedorită a reducerii zgomotului. Aceste efecte au ca rezultat nedorit reducerea clarităţii imaginii filtrate, ceea ce poate produce dificultăţi suplimentare în
71
etapele ulterioare de prelucrare a imaginii (segmentare, detecţie de contur, recunoaştere de forme etc.). O primă modalitate de reducere a efectului de bluring este ponderarea pixelilor mediaţi în funcţie de distanţa faţă de centrul ferestrei:
v(l,k)= ∑ ∑ c(i,j)⋅ u(l − i,k − j) (i,j)∈W
(3.24)
unde c(i,j) sunt coeficienţii cu care se face ponderarea pixelilor din fereastra W. Valorile acestor coeficienţi c respectă în general o anumită distribuţie spaţială, cel mai adesea fiind utilizată distribuţia gaussiană. Filtrul este cunoscut în acest caz, sub denumirea de filtru gaussian. O variantă îmbunătăţită pentru reducerea efectului de bluring o constituie filtrul de netezire cu prag, care nu mai este liniar. Particularitatea sa constă în faptul că înlocuirea valorii pixelului curent cu media ponderată a vecinilor săi se face doar dacă este satisfăcută condiţia: v(l , k ) − u (l , k ) < T
(3.25)
unde T este un prag de decizie ales astfel încât să fie protejate tranziţiile din imaginea iniţială faţă de efectul de bluring. Rezultatele acestei metode sunt bune, dar apare problema selecţiei automate a pragului T.
72
3.5. Efectul în frecvenţă al operatorilor liniari de vecinătate După cum s-a arătat mai sus, aplicarea unui filtru de mediere se face prin convoluţia între imaginea iniţială u şi funcţia pondere a filtrului h: v = h*u (3.26) Coeficienţii ferestrei sunt egali şi se ştie că transformata Fourier a unei constante este un sinc (figura 3.11): y
x
Figura 3.11. Transformata Fourier a unei constante. Dacă se notează cu U, V şi H transformatele Fourier ale imaginilor de intrare, respectiv de ieşire şi a filtrului, din relaţia anterioară şi din teorema convoluţiei rezultă că: V=H⋅U
(3.27)
iar caracterizarea în frecvenţă a acţiunii filtrului se poate face pe baza lui H. Transformata Fourier a unui filtru de mediere cu coeficienţi constanţi (figura 3.12.a) este un sinc bidimensional (figura 3.12.b):
73
z
z
F x
y
x
y
(a) (b) Figura 3.12. Transformata Fourier a unui filtru de mediere cu coeficienţi constanţi. Prin urmare, medierea spaţială este echivalentă cu o filtrare trece-jos (figura 3.13): u
Mediere spaţială
v
≡
u
v FTJ
Figura 3.13. Medierea spaţială. Celelalte tipuri de filtre (în frecvenţă) se pot obţine cu un FTJ. Astfel, dacă hTJ(m,n) este funcţia de transfer a unui FTJ atunci un filtru trece-sus FTS va avea o funcţie de transfer hTS(m,n): hTS (l , k ) = δ (l , k ) − hTJ (l , k )
(3.28)
unde δ este impulsul Dirac. Deci, un FTS se poate implementa prin scăderea din imaginea iniţială a imaginii obţinute printr-un FTJ (figura 3.14):
74
u
FTS
v
≡
+
u
v -
FTJ Figura 3.14. Obţinerea unui FTS cu ajutorul unui FTJ. La fel, un filtru trece-bandă FTB poate fi caracterizat prin relaţia: hTB (l , k ) = hTJ 1 (l , k ) − hTJ 2 (l , k )
(3.29)
unde hTJ1 şi hTJ2 sunt funcţii de transfer a două FTJ. Deci un FTB se poate obţine din 2 FTJ astfel (figura 3.15): u
FTB
v
≡
FTJ hTJ1(l,k)
u
v +
-
FTJ hTJ2(l,k) Figura 3.15. Obţinerea unui FTB cu ajutorul a două FTJ. Filtrarea trece-jos se utilizează pentru atenuarea zgomotului, filtrarea trece-bandă se foloseşte pentru extragerea sau accentuarea contururilor, iar filtrarea trece-sus este utilă pentru îmbunătăţirea contururilor sau a altor caracteristici de tip trece-sus ale unei imagini, în prezenţa zgomotului. Pe baza acestor considerente, pentru îmbunătăţirea imaginilor se pot defini şi operatori integrali (numiţi si transformări integrale) în cazul cărora noua valoare într-un punct depinde de valoarea întregii imagini
75
iniţiale. Folosirea transformărilor integrale mută problema într-un plan dual planului imaginii şi anume planul frecvenţelor spaţiale. Repartiţia spaţială a spectrului diferă de la o transformare la alta: FTJ
FTJ FTB
FTB FTS
FTB
FTB
Figura 3.16. Exemplu de repartiţie spaţială a spectrului unei imagini. Exemple de transformări (Fourier, Cosinus) vor fi prezentate în capitolul referitor la transformări integrale ale imaginilor.
3.6. Filtrarea neliniară a imaginilor Filtrele liniare (de netezire sau contrastare) produc la ieşire, în fiecare punct al imaginii, o combinaţie liniară ponderată a setului de valori selectate de fereastra de filtrare plasată cu originea în acel punct. Filtrele liniare pot elimina zgomote care corespund acestui model de mediere, deci zgomote aditive şi cu distribuţie normală (gaussiană). Din această comportare de tip filtru trece-jos rezultă efecte secundare care se manifestă prin reducerea sau eliminarea din imagine a componentelor de frecvenţă înaltă (detalii, contururi). Imaginea îşi pierde claritatea şi devine mai “ceţoasă” (efect de blur). În acelaşi timp, filtrarea liniară a unor zgomote ne-aditive (de exemplu zgomotul impulsiv) produce rezultate deranjante din punct de vedere al calităţii imaginii (în particular lăţirea şi împrăştierea impulsurilor de zgomot).
76
3.6.1. Filtre neliniare de ordine Pentru eliminarea dezavantajelor filtrelor liniare apare evidentă necesitatea de a modifica structura de filtrare liniară, în care valoarea unui pixel nu este luată în considerare, tocmai printr-o triere a valorilor extrase de fereastra de filtrare, în funcţie de rangul sau importanţa lor. Evident, operaţia devine neliniară, întrucât se va baza pe compararea şi ordonarea valorilor. Acesta este modelul de filtrare de ordonare (rank-order filter), fără a fi însă singurul tip de filtrare neliniară. Dacă se notează cu: X={x1,x2,…, xN}, cele N valori extrase de o fereastră de filtrare pentru o poziţie dată, setul de valori ordonate este: X()={x(1),x(2),…,x(N)}, cu x(1) ≤x(2)≤…≤x(N). Scalarul x(k) se numeşte statistica de ordine de ordinul k (sau pe scurt, statistica de ordin k) a secvenţei X. Evident, statistica de ordinul 1 este minimul, iar statistica de ordinul N este maximul: x(1) = min X() = min X (echivalentul unei erodări morfologice) (3.30) x(N)= max X()= max X (echivalentul unei dilatări morfologice) Considerând ca exemplu setul de 5 valori: X={1, 10, 100, 10, 200}, rezultă: x1=1, x2=10, x3=100, x4=10, x5=200. Setul ordonat este: X()={1,10,10,100,200}, iar statisticile de ordine sunt: x(1)=1, x(2)=10, x(3)=10, x(4)=100, x(5)=200. Ieşirea de ordinul k a filtrului de ordine este statistica de ordin k a setului de valori selectate de fereastra de filtrare, k putând lua orice valoare între 1 şi N: rankk(X)=x(k)
(3.31)
77
Cel mai utilizat filtru de ordine este filtrul median, caracterizat de: N + 1 k= , (unde [x]= partea întreagă a lui x) adică ieşirea filtrului 2 median este statistica de ordine situată în centrul secvenţei ordonate. Considerând, de exemplu, setul de 9 valori extrase de fereastra de filtrare (N=9), X={70,200,201,75,75,198,199,255,80}, ieşirea filtrului 9 + 1 median va fi statistica de ordinul k = =[5]=5. 2 Setul de valori ordonate este: X()={70,75,75,80,198,199,200,201,255}. Ieşirea filtrului median este, în acest caz: x(5)=198. Se presupune că valorile anterioare sunt nivelurile de gri dintr-o imagine. Atunci pixelii de valori 70,75,79,80 sunt gri-închis, pixelii de valori 198,199,200,201 au o culoare de tip alb-murdar, iar pixelul de valoare 255 este alb-strălucitor (figura 3.17).
Figura 3.17. Exemplu de filtrare mediană. Dacă acest punct ar fi centrul ferestrei de filtrare, după filtrare el trebuie înlocuit cu ieşirea filtrului, deci cu 198. Rezultatul este deci eliminarea punctului a cărui valoare este extremă faţă de celelalte (sau eliminarea impulsului de zgomot). Spre comparaţie, rezultatul filtrului de mediere pentru aceleaşi date este 150 (deci o valoare ce nu corespunde nici unuia dintre valorile pixelilor imaginii). La nivelul întregii imagini, aplicarea filtrului median nu modifică în mod esenţial structura de contururi (frontiere) caracteristică obiectelor. Este însă posibil ca anumite detalii fine (sau obiecte extrem de mici, de 78
dimensiune inferioară ferestrei de filtrare) să fie eliminate. Capacitatea unui filtru de a nu modifica anumite structuri de semnal, constituie o caracteristică deterministă a acestuia şi este exprimată de semnalele rădăcină (semnalele care nu sunt modificate la trecerea prin filtru). Pentru un filtru median, semnalele rădăcină sunt compuse din paliere constante şi rampe monotone, de lungime mai mare decât dimensiunea ferestrei de filtrare. Un semnal rădăcină se poate obţine prin filtrarea repetată (până la obţinerea invarianţei) a unui semnal oarecare.
3.6.2. Filtre de ordine multi-etaj Filtrele de ordine multietaj iterează mai multe etape de filtrare de ordine, realizate eventual în ferestre de filtrare de formă şi orientare diferită. Un asemenea filtru este filtrul MIN/MAX-Median. Acesta se compune dintr-un prim etaj de 4 filtre mediane cu ferestre direcţionale, ale căror ieşiri sunt preluate de un filtru de minim sau maxim. Dacă se consideră o fereastră 3×3 în care valorile pixelilor sunt notate: x1 x4 x7
x2 x5 x8
x3 x6 x9
ieşirea filtrului multietaj va fi: max(med(x4,x5,x6),med(x2,x5,x8), med(x1,x5,x9), med(x3,x5,x7)) (3.32)
79
MAX (sau MIN)
Figura 3.18. Filtru multi-etaj. Efectul acestui filtru este mai bun decât a unui filtru median, deoarece se obţine un contrast mai bun, datorită filtrului extrem (max sau min), cu excepţia liniilor de zgomot, care trec prin filtrele mediane. Se poate arăta că filtrele de ordine pentru prelucrarea imaginilor cu niveluri de gri, pot fi dispuse sub forma unei stive, prelucrarea putându-se face independent pe fiecare nivel al stivei, suma rezultatelor fiind egală cu rezultatul filtrului compus obţinut prin însumarea filtrelor din stivă. Acest fapt are avantajul posibilităţii de implementare paralelă a filtrelor de ordine şi posibilitatea de implementare hardware a acestora, obţinându-se o viteză mare de calcul şi o paralelizare a calculelor. Dezavantajul acestor filtre îl constituie cantitatea mai mare de memorie necesară în procesul de sortare. La implementarea filtrelor de ordine, printre cele mai utilizate tehnici de sortare se numără metoda bubble-sort şi divide et impera.
80
3.6.3. Proprietăţi ale filtrelor de ordine 1. Invarianţa la translaţie Filtrarea de ordine a unui semnal x translatat cu b şi scalat cu a este: rank k (a ⋅ x + b) → a ⋅ rank k ( x) + b (3.33) 2. Filtrele de ordine păstrează caracteristicile semnalului, adică nu prezintă overshoot sau undershoot. Prin urmare, valoarea obţinută prin filtrare nu este în afara domeniului de intrare, fiind una din valorile de intrare. 3. Filtrele de ordine admit semnale rădăcină. Un semnal rădăcină este invariant la filtrarea de ordine (rămâne nemodificat). Teoremă: Orice secvenţă monotonă (crescătoare sau descrescătoare) este un semnal rădăcină al filtrelor de ordine. În practică, un semnal rădăcină se obţine prin filtrarea de ordine a unui semnal oarecare, până se obţine un semnal care nu se mai modifică. Exemplu: Prin aplicarea unui filtru median unidimensional de lungime N=3 (figura 3.19.b), cu originea în centru, asupra semnalului (figura 3.19.a): 3
2
N=3
1 (b)
(a)
Figura 3.19. Exemplu de obţinere a unui semnal rădăcină. după prima iteraţie se obţine semnalul: 81
Figura 3.20. Rezultatul după prima iteraţie a filtrării de ordine. iar după a doua iteraţie se obţine semnalul:
Figura 3.20. Rezultatul după a doua iteraţie a filtrării de ordine. După cum se observă, după a treia iteraţie se obţine un semnal monoton pe porţiuni, care dacă mai este filtrat o dată, nu se modifică, deci este un semnal rădăcină pentru un filtru median unidimensional de lungime 3. O limitare a structurilor de filtrare neliniară bazată pe ordonare este relativa lor lipsă de flexibilitate: există doar N filtre de ordine diferite, dintre care doar filtrul median este un filtru de netezire. O structură de filtrare mai flexibilă trebuie să permită reglarea gradului de netezire (sau respectiv de reliefare) între limite fixate. Astfel de filtre sunt filtrele de domeniu LUM (Lower-Upper-Middle).
82
3.6.4. Filtre de ordine de domeniu Filtrul de ordine LUM (Lower-Upper-Middle) de netezire, de ordin k, se defineşte prin: y = med x( k ) , x*, x( N −k +1) , (3.34)
{
}
N +1 , iar x* reprezintă valoarea eşantionului central 2 (din originea ferestrei): unde 1 ≤ k ≤
x( k ) , x* < x( k ) y = x( N −k +1) , x* > x( N −k +1) x*, în rest
(3.35)
Netezirea este creată prin compararea eşantionului central x* cu două statistici de ordin superior şi inferior (x(k), x(N-k+1)) şi înlocuirea sa cu un eşantion mai apropiat de mediană, dacă x* nu se încadrează în acest interval de valori „normale” (figura 3.21).
x(1)
x(N-k+1)
x(k)
x(N)
Figura 3.21. Filtru LUM de netezire. Gradul de netezire este variabil în funcţie de k. Cazurile extreme ale acestui tip de filtru sunt: N +1 • pentru: k = = filtru median 2 • pentru: k = 1 = filtru de tip trece-tot
83
Prin aplicarea unui filtru LUM de netezire se reduce contrastul imaginii. Perechea filtrului LUM de netezire este filtrul LUM de reliefare (sau conturare), pentru care ieşirea filtrului este deplasată către una din statisticile extreme (inferioară sau superioară). Filtrul LUM de reliefare se defineşte prin: x( r ) + x( N −r +1) x( r ) , x( r ) < x* ≤ 2 x( r ) + x( N −r +1) y = x( N −r +1) , < x* ≤ x( N −r +1) , 2 x*, în rest
(3.36)
N +1 . 2 Reliefarea apare ca urmare a deplasării eşantioanelor din intervalul x(r) ÷ x(N-r+1) (interpretat ca zonă de tranziţie) spre una din extremităţi (figura 3.22): unde: 1 ≤ r ≤
x(1)
x(r) x( r ) + x( N −r +1) x(N-r+1) 2
x(N)
Figura 3.21. Filtru LUM de reliefare. Gradul de reliefare poate varia în funcţie de parametrul r, de la: N +1 r= (filtrul identitate) până la reliefare maximă, pentru r=1. 2 Prin aplicarea unui filtru LUM de reliefare se obţine o mărire a contrastului.
84
3.6.5. L-filtre Un filtru a cărui ieşire este o medie ponderată a statisticilor din fereastra filtrului se mai numeşte şi L–filtru şi este de forma: N
y = ∑ a k ⋅ x( k ) ,
(3.37)
k =1
unde y este ieşirea filtrului, x(k) sunt eşantioanele din fereastra de filtrare ordonate crescător, iar ak sunt coeficienţii de ponderare. Coeficienţii de ponderare îndeplinesc aceleaşi condiţii de normalizare ca şi pentru un filtru liniar: N
N
k =1
k =1
∑ a k = 1 , pentru netezire şi ∑ a k = 0 , pentru reliefare.
(3.38)
Câteva cazuri particulare ale acestui tip de filtru sunt:
• filtrul de ordine de ordin k (inclusiv filtrul median): a k = 1 . = ≠ a 0 , j k k Elimină zgomotul impulsiv.
(3.39)
• filtrul de mediere, se obţine pentru: 1 a k = , pentru orice k = 1, 2, ..., N. N Elimină zgomotul gaussian.
(3.40)
• extractor de contur: a1 = −1 a N = 1
(3.41)
85
• filtrul de qvasi-mijloc: a i = a N +1 − i N +1 , i = 1÷ . 2 a k = 0 Elimină zgomote impulsive şi zgomotele uniforme.
(3.42)
• media α-reglabilă (α-trimmed mean): 1 , i = [α ⋅ N , N − α ⋅ N ] a i = N ⋅ (1 − 2α ) , a = 0, în rest k 1 unde α ∈ 0, . 2 Elimină zgomote de tip impulsiv şi gaussian.
86
(3.43)
4. Transformări integrale ale imaginilor În domeniul prelucrărilor de imagini, pe lângă operaţiile punctuale şi de vecinătate prezentate anterior, se folosesc adesea şi transformări integrale, al căror rezultat – o altă imagine sau o altă reprezentare a imaginii originale – are în fiecare punct o valoare ce depinde de valorile tuturor pixelilor din imagine originală. Transformata Fourier, de exemplu, transformă o imagine reprezentată prin funcţia f(x,y) de cele două variabile spaţiale x şi y, în „spectrul Fourier” al imaginii, F(u,v), care este o altă reprezentare a imaginii originale f(x,y), în planul frecvenţelor spaţiale u şi v. O filtrare (spaţială) liniară a imaginii f(x,y) cu un filtru a cărui funcţie pondere h(x,y) este un clopot Gauss, are ca rezultat o altă imagine g(x,y) obţinută prin convoluţie (bidimensională) din imaginea originală f(x,y): g ( x, y ) = ∫∫ f (ξ ,η ) ⋅ h( x − ξ , y − η )dξdη R
(4.1)
2
iar sensul fizic al variabilelor lui g este acelaşi cu cel al variabilelor lui f; f şi g sunt funcţii definite pe acelaşi plan (x,y), spre deosebire de primul exemplu cu transformata Fourier.
4.1. Transformări integrale unitare În cazul imaginilor discretizate (definite pe latici de puncte), cum e cazul imaginilor digitale, transformările integrale sunt transformări ale unor matrici (tabele 2-dimensionale de scalari, de regulă dreptunghiulare şi chiar mai mult, pătrate) în matrici de aceleaşi dimensiuni. S-a văzut că transformările punctuale sunt practic transformări neliniare, iar transformările spaţiale (de vecinătate) pot fi liniare sau neliniare. Transformările integrale sunt transformări liniare pentru a avea
87
posibilitatea de a beneficia în modelarea matematică care trebuie făcută, de aparatul algebrei liniare şi al analizei funcţionale liniare. În cele de mai jos se vor revedea noţiunile de bază privind transformările liniare pe spaţii vectoriale finit dimensionale. Se va face referire la spaţii vectoriale peste mulţimea numerelor complexe C şi transformări liniare corespunzătoare (reprezentate prin matrici pătrate cu elemente complexe) fiindcă, chiar începând cu transformata Fourier discretă, cadrul complex simplifică lucrurile. Din acest motiv se va face referire la „matrici unitare” şi nu „ortogonale”. Se va considera spaţiul CN al vectorilor N-dimensionali cu componente din C (deci numere complexe) care reprezintă secvenţe de N eşantioane ale unui semnal complex (cu o parte reală şi una imaginară). Orice transformare liniară a lui CN în el însuşi se reprezintă printr-o matrice pătrată A de dimensiune N×N: (alk ) l ,k =0, N −1 = A ∈ M N × N (C ) (4.2) iar dacă se doreşte ca transformarea să fie inversabilă, atunci matricea A trebuie să fie inversabilă, adică det A ≠ 0 . O matrice A este unitară dacă are proprietatea că este inversabilă, adică există A-1 astfel încât A·A-1= A-1·A=IN, şi că A-1=A*T, adică: A ⋅ A*T = A*T ⋅ A = I N (4.3) adică inversa ei este chiar transpusa conjugatei sale complexe. În această relaţie indicele superior T indică operaţia de transpunere a matricelor (alk→akl), iar indicele superior * indică operaţia de conjugare complexă (alk→akl*). A*T=AH se mai numeşte şi transformata Hilbert a matricii A. O matrice unitară în spaţii vectoriale peste mulţimea numerelor reale R (deoarece conjugatul unui număr real este el însuşi) se reduce la ceea ce se numeşte matrice ortogonală, adică A-1=AT; matricile ortogonale generalizează rotaţiile din plan.
88
u0 . Un vector u ∈ C N , u = (u i )i =0, N −1 = . , devine prin . u N −1 transformarea A, vectorul v, dat de relaţia v=A·u, adică: v 0 a 0, 0 . . . = . . . v N −1 a N −1,0
.
.
.
.
.
.
a 0, N −1 u 0 . ⋅ . . a N −1, N −1 u N −1
(4.4)
sau pe componente, în scriere scalară: N −1
v j = ∑ a jk ⋅ u k
(4.5)
k =0
În cazul imaginilor, dacă se notează imaginea iniţială cu U = {u (l , k ); l , k = 0,1,..., N − 1} şi imaginea transformată cu: V = {v(m, n); m, n = 0,1,..., N − 1}, expresiile matriceale ale transformării sunt mai dificil de scris, deoarece transformarea A trebuie reprezentată printr-un tablou 4-dimensional (2×2, semnalul de intrare fiind imaginea U 2-dimensională, iar ieşirea fiind V, tot 2-dimensională). În schimb, scrierea scalară (pe componente) a transformării integrale unitare bidimensionale directe, respectiv inverse, a imaginii U, se scrie imediat, prin generalizarea celei precedente: N −1 N −1
v ( m, n ) = ∑ ∑ u (l , k ) ⋅ a mn (l , k ) , respectiv l =0 k = 0
N −1 N −1
∗ (l , k ) , u (l , k ) = ∑ ∑ v ( m, n ) ⋅ a mn m = 0 n =0
unde {a mn (l , k )} este transformarea unitară.
89
(4.6)
În cazul imaginilor de dimensiuni N×N, numărul operaţiilor (multiplicări şi adunări) necesare pentru a calcula coeficienţii v(m,n) este foarte mare, fiind egal cu O(N4), adică de ordinul a N4. O reducere substanţială a complexităţii algoritmului se obţine atunci când transformarea unitară este separabilă, adică: a mn (l , k ) = a m (l ) ⋅ bn ( k ) (4.7) unde: A = {a m (l )} şi
B = {bn (k )} , sunt de asemenea matrici
unitare. În acest caz: N −1 N −1
N −1
N −1
l =0 k =0
l =0
k =0
⇒ v ( m, n) = ∑ ∑ a m (l ) ⋅ bn ( k ) ⋅ u (l , k ) = ∑ a m (l ) ∑ bn ( k ) ⋅ u (l , k )
(4.8)
După cum se observă, în cazul în care transformarea este separabilă, algoritmul este mai simplu, complexitatea sa scăzând la O(N3). În acest caz, relaţia transformării directe, respectiv inverse, se poate scrie şi sub formă matricială: V = A ⋅ U ⋅ BT
, respectiv
U = A∗T ⋅ V ⋅ B ∗
(4.9)
Dacă se defineşte energia unui semnal bidimensional u, ca fiind: N −1 N −1
E u = ∑ ∑ u (l , k )
2
(4.10)
l =0 k =0
se poate arăta că principala proprietate a transformărilor integrale unitare este aceea de conservare a energiei: Eu=Ev. Transformările integrale unitare se pot folosi pentru filtrarea şi pentru compresia imaginilor. • Filtrarea: Se presupune că se doreşte o filtrare integrală liniară a unei imagini, adică cu un filtru a cărui funcţie pondere liniară se extinde pe 90
toată imaginea. Operaţia de convoluţie spaţială a unei imagini N×N (cu funcţia pondere tot de suport N×N) necesită N4 operaţii de înmulţire. Dacă, aşa cum este cazul transformării Fourier discrete, transformarea are algoritm rapid şi este separabilă, iar echivalentul convoluţiei este o operaţie punctuală (pentru fiecare pixel avem de făcut o înmulţire), numărul de operaţii necesare este de ordinul a N2·log2N. Pentru N=1000, adică pentru imagini de 1000×1000 pixeli, raportul N4/N2·log2N devine 105 adică filtrarea care ar dura 24 ore în mod normal, prin acest procedeu se realizează într-o secundă. • Compresia: Spre deosebire de reprezentarea unei imagini prin eşantioanele sale spaţiale (în număr de N2) care, în principiu, într-o primă aproximaţie sunt variabile aleatoare independente şi uniform distribuite pe intervalul valorilor posibile [0,M-1], M-1 fiind valoarea maximă, reprezentarea imaginii prin valorile coeficienţilor transformatei (tot în număr de N2) este o reprezentare prin mărimi cu proprietăţi extrem de neuniforme: câţiva coeficienţi sunt foarte mari, iar majoritatea sunt neglijabil de mici. Transmisiei lor li se poate aplica o codare Huffman (pentru a obţine o codare fără pierderi), dar dacă unii coeficienţi sunt foarte mici, ei se pot neglija, pur şi simplu. În acest caz, etapele transmisiei sunt schiţate în figura 4.1. . . . .. U. .
T
.......
V
trunchiere
....
…Vˆ
0
canal fără pierderi
T-1 Vˆ
. . . . … Uˆ .
Figura 4.1. Compresia cu transformate. În procesul de compresie cu transformate se porneşte de la o imagine iniţială U. Acesteia i se aplică o transformare T în urma căreia se obţine o imagine V la care informaţia este concentrată în mult mai puţine componente decât în imaginea iniţială.
91
Canalul de transmisie poate fi presupus fără pierderi deoarece se pot folosi coduri cu corecţie de erori. Câştigul (reducerea numărului de eşantioane transmise) s-a obţinut prin trunchierea imaginii V prin eliminarea componentelor cu coeficienţi foarte mici (sub un anumit prag). La recepţie se obţine o imagine Vˆ care aproximează imaginea V. Dacă acestei imagini i se aplică transformarea inversă lui T (T-1) se va obţine o imagine Uˆ care este de dorit a fi cât mai apropiată de imaginea iniţială U. Un astfel de lanţ de compresie este cu atât mai bun cu cât: • eroarea (diferenţa) dintre Uˆ şi U este mai mică, relativ la un anumit criteriu, care de cele mai multe ori este eroarea medie pătratică: L −1 K −1
⇒ ∑ ∑ U (l , k ) − Uˆ (l , k ) 2 =minimă l =0 k =0
(4.11)
• factorul de compresie este mai mare. Factorul de compresie (C) se defineşte ca fiind raportul dintre numărul componentelor imaginii iniţiale (L×K) şi numărul componentelor din u, reţinute în imaginea v: L× K C= (4.12) nr. componentelor retinute din u Transformarea optimă din punct de vedere al compresiei este transformarea K-L (Karhunen-Loeve). Aceasta realizează decorelarea elementelor transformatei şi prin aceasta compactarea maximă a energiei semnalului (imaginii) în primele componente. Cu toate aceste avantaje, deoarece transformata K-L depinde de statistica imaginii şi datorită volumului mare de calcule necesare (mai ales în cazul imaginilor de dimensiuni mari), transformata K–L este dificil de implementat practic. Totuşi, pentru imagini cu corelaţie mare, transformata K-L se poate aproxima cu succes prin transformata Cosinus Discretă (DCT=Discrete Cosine Transform) mult mai rapidă şi mai uşor de implementat şi care va fi prezentată în continuare. 92
4.2. Matrici unitare Implementarea transformărilor integrale unitare a imaginilor se face utilizând matrici unitare. Fie o matrice pătrată, cu elemente complexe: A = (alk )
l = 0,l −1 k = 0 , k −1
,
unde alk ∈ C . Spunem că matricea A ∈ M N × N (C ) este unitară dacă:
A ⋅ A*T = I N , *T A ⋅ A = I N
(4.13)
not .
unde A*T = A H este transformata Hilbert a matricii A. Sistemul de mai sus este echivalent cu relaţia: A −1 = A*T . Notând A = (alk ) l ,k =0, N −1 , condiţiile (4.9) se pot scrie: N −1 1, pt. l = l ' * ∑ alk ⋅ alk ' = δ (l − l ' ) = l =0 0, pt. l ≠ l ' ⇒ N −1 a ⋅ a * = δ (k − k ' ) = 1, pt. k = k ' ∑ lk l 'k 0, pt. k ≠ k ' k =0
93
(4.14)
Proprietăţile matricilor unitare 1.) O transformare liniară dată de o matrice (A) unitară (de la CN la CN) lasă produsul scalar invariant. x1 y1 . . N ∀x, y ∈ C , adică x = . , y = . . . xN yN Se defineşte produsul scalar a doi vectori: def . N
N
x , y = ∑ xi* ⋅ y i = x *T ⋅ y = ∑ xiT ⋅ y i* = x T ⋅ y * i =1
i =1
(4.15)
(4.16)
În aceste condiţii: A ⋅ x, A ⋅ y = x, y ,
(4.17)
adică produsul scalar este invariant. Demonstraţie: A ⋅ x , A ⋅ y = ( A ⋅ x )*T ⋅ ( A ⋅ y ) . Se ştie că: (U ⋅ V )T = V T ⋅ U T .
( A ⋅ x )*T
(
⋅ ( A ⋅ y ) = x T ⋅ AT
)* ⋅ A ⋅ y = x T * ⋅ 1A4T2*4⋅3A ⋅ y = x T * ⋅ y =
x , y q.e.d.
I
Consecinţă: Transformarea dată de o matrice unitară lasă neschimbată energia. Acest lucru rezultă din proprietatea anterioară, pentru cazul particular în care x=y. x, x = A ⋅ x, A ⋅ x T ⇒1 x *2 ⋅ x = ( A ⋅ x )*T ⋅ ( A ⋅ x ) 3 1442443 x
2
A⋅ x
2
94
N
2
⇒ ∑ xi i =1
N
= ∑ a i ⋅ xi
2
⇒ E x = E A⋅ x
i =1
2.) Toate valorile proprii ale unei matrice unitare au valorile egale cu 1. Dacă x este un vector propriu al matricei A, atunci λ este valoarea (4.18) proprie a lui x dacă satisface relaţia: A ⋅ x = λ ⋅ x . Demonstraţie: S-a arătat că: ( A ⋅ x )*T ⋅ ( A ⋅ x ) = x *T ⋅ x . Înlocuind: A ⋅ x = λ ⋅ x ⇒ (λ ⋅ x )*T ⋅ (λ ⋅ x ) = x *T ⋅ x ⇒ λ* ⋅ λ ⋅ x *T ⋅ x = x *T ⋅ x 2
⇒λ ⋅ x
2
2
= x . Deoarece: x ≠ 0 (vectorul propriu este nenul)
⇒ λ = 1 q.e.d. În aceste relaţii, norma vectorului x indusă de produsul scalar este:
x
2
def .
N
= x *T ⋅ x = ∑ xi
2
i =1
= E x = energia lui x
(4.19)
3.) Vectorii proprii corespunzători unor valori proprii diferite sunt ortogonali. Fie x1, x2 doi vectori proprii cu valori proprii asociate λ1, λ2. Spunem că doi vectori sunt ortogonali dacă:
x1 ⊥ x 2 ⇔ x1 , x 2 = 0 ⇔ x1*T ⋅ x 2 = 0
(4.20)
Demonstraţie: Dacă x1, x2 sunt vectori proprii cu valorile proprii asociate λ1, λ2: ⇒ A ⋅ x1 = λ ⋅ x1 , A ⋅ x 2 = λ ⋅ x 2 Deoarece matricea A este unitară, rezultă că lasă produsul scalar invariant: ⇒ ( A ⋅ x1 )*T ⋅ ( A ⋅ x 2 ) = x1*T ⋅ x 2 (conform proprietăţii 1) ⇒ (λ1 ⋅ x1 )*T ⋅ (λ2 ⋅ x 2 ) = x1*T ⋅ x 2 95
⇒ λ1* ⋅ λ2 ⋅ x1*T ⋅ x 2 = x1*T ⋅ x 2 Dar: λ1* ⋅ λ2 ≠ 1 , deoarece λ1 ≠ λ 2 . ⇒ x1*T ⋅ x 2 ≠ 0 q.e.d. Rămâne de arătat că: λ1* ⋅ λ2 ≠ 1 . Se poate demonstra prin reducere la absurd. 4.) Teorema spectrală: Dacă A este o matrice unitară, atunci există cel puţin o matrice B care diagonalizează matricea A. Această teoremă poate fi descrisă astfel: ∀A ∈ M N × Ν (C ) o matrice unitară de dimensiune N×N, adică
A ⋅ AH = AH ⋅ A = I N ⇒ ∃B ∈ M N × Ν (C ) o matrice unitară de dimensiune N×N, adică B ⋅ B H = B H ⋅ B = I N , care diagonalizează matricea A, adică ∃Λ ∈ M N × Ν (C ) o matrice diagonală de dimensiune N×N, care satisface relaţia: B −1 ⋅ A ⋅ B = Λ , (4.21) adică matricea B diagonalizează matricea A, rezultatul fiind matricea Λ.
Demonstraţie: Dacă x1, …, xN sunt vectori proprii normaţi ai matricei A. Se va construi matricea B astfel: M B N × N = x1 M
M
x2 M
M M ... x N , unde x i M M M M
este o matrice
coloană. Dar: B −1 = B H = B T * , deoarece s-a demonstrat că vectorii proprii ai unei matrice unitare sunt ortogonali.
96
L x1*T *T T * L x 2 ⇒B = M M L x *T N
L L , unde M L
L x1*T
L
este o
matrice linie. 1 0 T* ⇒ B ⋅B = M 0
0 K 0 1 K 0 = I N (conform proprietăţii 3), M K M 0 K 1
deci matricea B astfel construită este unitară. L x1*T L x 2*T −1 ⇒ B ⋅ A⋅ B = M M L x *T N L x1*T L x 2*T = M M L x *T N
L M M L ⋅ λ1 ⋅ x1 λ2 ⋅ x 2 M M M L
L M L ⋅ A ⋅ x1 M M L
M x2 M
M ... x N = M M M
λ1 0 K 0 0 λ2 K 0 ... λ N ⋅ x N = M K M M M M 0 0 K λ N
M
M
⇒ B −1 ⋅ A ⋅ B = diag (λ1 ,..., λ N ) = Λ Prin urmare, alegând matricea B ca având pe coloane vectorii proprii normaţi ai matricea A şi matricea Λ având pe diagonală valorile proprii corespunzătoare, este verificată relaţia:
BH ⋅ A⋅ B = Λ
(4.22)
97
4.3. Transformări unitare ale unor semnale unidimensionale O transformare liniară de la C N → C N este unitară dacă este reprezentată de o matrice unitară A ∈M N × N (C ) . Dacă notăm cu u(l) un semnal unidimensional discret, cu N eşantioane ( l = 1, N ) şi cu v(m) semnalul eşantionat transformat: u (1) v(1) U = M , V = M u( N ) v( N )
În scriere matriceală: V = A⋅U ,
(4.23)
⇒ U = A −1 ⋅ V = A*T ⋅ V
(4.24)
unde A = (a m (l ) )l ,m=1, N . Indicele m semnifică indicele semnalului din baza de semnale, iar argumentul l semnifică eşantionul m din baza de semnale. Matricea A se scrie extins: a1 (1) a1 (2) K a1 ( N ) a 2 (1) a 2 (2) K a 2 ( N ) A= M M M M a (1) a (2) K a ( N ) N N N
Cu aceste notaţii, relaţiile (4.19) şi (4.20) se scriu, pe componente: N
v(m) = ∑ a m (l ) ⋅ u (l )
(4.25)
l =1 N
* u (l ) = ∑ a m (l ) ⋅ v(m)
(4.26)
m=1
N
* A ⋅ A*T = I N ⇔ ∑ a m (l ) ⋅ a m ' (l ) = δ ( m − m' ), m =1, N l =1
98
(4.27)
N
* A*T ⋅ A = I N ⇔ ∑ a m (l ) ⋅ a m (l ' ) = δ (l − l ' ), l =1, N l =1
(4.28)
1, x = 0 1, l = l ' , adică δ (l − l ' ) = unde: δ ( x) = 0, în rest 0, în rest , l ≠ l '
este impulsul Dirac. Coloanele matricei unitare A formează o bază ortonormată în spaţiul vectorial CN. Prima condiţie de mai sus (4.23) exprimă condiţia de ortonormalitate, iar cea de-a doua condiţie (4.24) exprimă condiţia de completitudine a bazei. Transformata unitară directă a unei secvenţe este: N
v(m) = ∑ a m (l ) ⋅ u (l )
(4.29)
l =1
iar transformata unitară inversă este: N
* u (l ) = ∑ a m (l ) ⋅ v(m)
(4.30)
m=1
Condiţia ca transformata să fie unitară este: N
* ∑ a m (l ) ⋅ a m ' (l ) = δ ( m − m' )
l =1
(4.32)
care exprimă condiţia de completitudine. N
* (l ' ) = δ (l − l ' ) ∑ a m (l ) ⋅ a m
m=1
exprimă condiţia de ortonormalitate.
99
(4.32)
4.4. Transformări unitare ale unor semnale bidimensionale Se spune că un semnal (v(m, n )) m,n=0, N −1 se obţine prin aplicarea unei transformări unitare asupra semnalului (u (l , k )) l ,k =0, N −1 , dacă: N −1 N −1
v(m, n) = ∑ ∑ a m,n (l , k ) ⋅ u (l , k ) este transformarea directă (4.33) l =0 k =0
N −1 N −1 * u (l , k ) = ∑ ∑ a m ,n (l , k ) ⋅ v ( m, n) este transformarea inversă (4.34) m =0 n =0
Condiţia de completitudine se scrie în acest caz: l = l ' 1, * ∑ ∑ a m,n (l , k ) ⋅ a m ,n (l ' , k ' ) = δ (l − l ' , k − k ' ) = k = k ' (4.35) m =0 n = 0 0, în rest def .
N −1 N −1
Condiţia de ortonormare se scrie: m = m ' 1, * ∑ ∑ a m,n (l , k ) ⋅ a m ',n ' (l , k ) = δ ( m − m' , n − n' ) = n = n' (4.36) l =0 k =0 0, în rest def .
N −1 N −1
În ceea ce priveşte numărul de operaţii, un algoritm are complexitatea O(NP) dacă numărul de numărul de operaţii (înmulţiri, adunări etc) este proporţional cu Np, atunci când secvenţa de date are lungimea N. În cazul transformărilor unitare unidimensionale, N −1
v(m) = ∑ a m (l ) ⋅ u (l ) , numărul de operaţii necesare Nx pentru a calcula l =0
v(m) este Nx= N2. Prin urmare, complexitatea algoritmului de calcul a lui v(m) este O(N2). 100
În cazul bidimensional, numărul de operaţii necesare pentru a N −1 N −1
calcula v(m, n) = ∑ ∑ a m,n (l , k ) ⋅ u (l , k ) este de ordinul a Nx= N4, deci k =0 l =0
complexitatea algoritmului este O(N4), însă adesea se încearcă reducerea complexităţii algoritmului. Se spune că o transformare unitară bidimensională este separabilă dacă coeficienţii transformării se pot scrie: (4.37) a mn (l , k ) = a m (l ) ⋅ bn (k ) În acest caz: N −1 N −1
v(m, n) = ∑ ∑ a m (l ) ⋅ bn (k ) ⋅ u (l , k ) l =0 k =0
(4.38)
Se poate arăta că dacă o transformare unitară bidimensională este separabilă, complexitatea algoritmului său de calcul se reduce la O(N3), adică Nx= N3, deoarece relaţia de mai sus se poate scrie şi: N −1
N −1
l =0
k =0
v(m, n) = ∑ a m (l ) ⋅ ∑ bn (k ) ⋅ u (l , k )
(4.39)
Pentru a ilustra scrierea matriceală a unei transformări unitare separabile (ultimul argument a lui bn(k) trebuie să coincidă cu primul argument a lui u(l,k)), relaţia de mai sus se poate scrie şi: N −1
N −1
l =0
k =0
v(m, n) = ∑ a m (l ) ⋅ ∑ u (l , k ) ⋅ bn (l )
(4.40)
În acest caz, transformarea unitară directă separabilă se poate scrie: v = A ⋅ u ⋅ BT ,
(4.41)
unde A = (a m (l )) m,l =0, N −1 , B = (bn (l )) n,l =0, N −1 sunt matricile unitare în care se poate separa transformarea iniţială. Transformarea inversă se poate scrie:
( )−1
u = A −1 ⋅ v ⋅ B T
(4.42)
101
Deoarece A şi B sunt, la rândul lor, matrici unitare, adică A −1 = A*T şi B −1 = B *T , deci
(BT )−1 = (B −1 )T = B* ,
transformarea
inversă se mai poate scrie: ⇒ u = A*T ⋅ v ⋅ B * (4.43) Un caz particular, îl constituie cazul în care matricile unitare A şi B sunt egale, A = B. În acest caz, transformarea directă este: N −1 N −1
v(m, n) = ∑ ∑ a m (l ) ⋅ a n (k ) ⋅ u (l , k ) l =0 k =0
(4.44)
iar transformarea inversă este: N −1 N −1
* u (l , k ) = ∑ ∑ a m (l ) ⋅ a n* (l ) ⋅ v(m, n) m =0 n =0
(4.45)
Aceste relaţii trebuie însoţite de condiţiile de ortonormalitate şi completitudine. Spaţiul matricilor pătrate M N × N (C ) are dimensiunea N2. În spaţiul matricilor pătrate de dimensiune N×N se poate defini o transformare unitară bidimensională şi prin utilizarea produsului vectorial, după cum se va arăta în continuare. Dacă matricile
X , Y ∈ M N × N (C ) , unde
X = ( xij ) i , j =0, N −1 ,
Y = ( yij ) i , j =0, N −1 , produsul scalar al celor două matrici se defineşte prin relaţia: def . N −1 N −1
X , Y = ∑ ∑ xij ⋅ yij* i =0 j =0
(4.46)
Prin urmare, produsul scalar dă spaţiului M N × N (C ) o structură de spaţiu Hilbert (un spaţiu liniar cu un produs scalar). În acest spaţiu se va fixa o matrice unitară A, care are elementele A = (a m (l )) m,l =0, N −1 , deci:
102
| A = a0 | not .
a1 ... a N −1 | | |
|
(4.47)
În acest spaţiu se poate construi o bază B în M N × N (C ) astfel:
{ }
* * T* * * B = Amn m,n =0, N −1 , unde Amn = a m ⋅ a n , vectorii a m fiind vectori
coloană (N×1), iar vectorii a n* sunt vectori linie (1×N). Pentru a verifica dacă B este o bază în spaţiul M N × N (C ) se procedează prin reducere la absurd, presupunând că B nu este o bază în spaţiul M N × N (C ) . Prin urmare rezultă că există matricile: N −1 N −1
* = 0 N ×N ⇒ ∃α mn ∈ C , nu toate nule, a.î. ∑ ∑ α mn ⋅ Amn
m = 0 n =0
⋅ al
Înmulţind relaţia de mai sus cu al, se obţine relaţia: N −1 N −1
* ⇒ ∑ ∑ α mn ⋅ a m ⋅ a n*T ⋅ al = 0
m =0 n = 0
(4.48)
1, pt. n = l în care: a n*T ⋅ al = δ (n − l ) = , 0, pt. n ≠ l deoarece matricea A este unitară. Prin urmare: N −1
* ⇒ ∑ α mnl ⋅ a m = 0,
m =0
(4.49)
adică componentele ak sunt liniar dependente, ceea ce este fals, deoarece coloanele unei matrici unitare formează o bază în spaţiul N×N, deci sunt liniar independente.
{ }
* * T* * Deoarece B = Amn m,n =0, N −1 , unde Amn = a m ⋅ a n este o bază în
spaţiul M N × N (C ) : N −1 N −1
* u N × N (l , k ) = ∑ ∑ v(m, n) ⋅ Amn
m =0 n =0
103
(4.50)
* Din algebra liniară se ştie că: v(m, n) = u , Amn . Astfel, rezultă:
N −1 N −1
(
)*
* * (l , k ) v(m, n) = u , Amn = ∑ ∑ u (l , k ) ⋅ a mn
l =0 k =0
(4.51)
* * Dar: Amn (l , k ) = a m (l ) ⋅ a n* (k ) .
N −1 N −1
⇒ v(m, n) = ∑ ∑ u (l , k ) ⋅ a m (l ) ⋅ a n (k ) l =0 k = 0
(4.52)
4.5. Transformata Fourier discretă unidimensională (DFT-1D) Deoarece transformarea Fourier este una din transformările cel mai des utilizate în domeniul prelucrărilor de imagini, în continuare va fi prezentată transformata Fourier discretă unidimensională şi proprietăţile acesteia. Transformata Fourier discretă unidimensională DFT-1D directă a unei secvenţe discrete u(k) este definită prin relaţia: N −1 2π v(m) = ∑ u (k ) ⋅ exp − j m ⋅ k , unde m = 0, N − 1 N k =0
(4.53)
Transformata Fourier discretă unidimensională DFT-1D inversă este definită prin relaţia: u (k ) =
1 N −1 2π ∑ v(m) ⋅ exp j m ⋅ k , unde k = 0, N − 1 N m =0 N
Matricea
transformării
DFT-1D
(4.54) este:
2π m ⋅ k , unde m este indicele de linie şi k indicele exp − j N m,k =0, N −1 de coloană. Se poate defini transformata DFT-1D unitară prin relaţiile: v ( m) =
1 N −1
∑ u (k ) ⋅ exp − j
N l =0
104
2π m⋅k, N
(4.55)
pentru transformata DFT-1D directă u (k ) =
1 N −1
2π m ⋅ k , N
∑ v(m) ⋅ exp j
N m =0 pentru transformata DFT-1D inversă. Se poate arăta că matricea
(4.56)
transformării,
1 2π ⋅ exp − j F = m ⋅ k , este o matrice unitară, adică N m,k =0, N −1 N F ⋅ F *T = F *T ⋅ F = I N .
4.6. Proprietăţi ale transformatei DFT-1D 1. Inversa matricei DFT-1D este egală cu conjugata matricei DFT-1D: def . 1 2π ⋅ exp − j F = m ⋅ k N m,k =0, N −1 N
(4.57)
F = F T , deoarece F * = F T * = F −1 , unde: ⇒ F −1 = F *
(4.58)
1 2π ⋅ exp j F −1 = F * = m ⋅ k N m,k =0, N −1 N
(4.59)
2. Extensia vectorului transformat este periodică, v(m + N ) = v(m) : (4.60) N −1 2π (m + N ) ⋅ k = v(m + N ) = ∑ u (k ) ⋅ exp − j N k =0 N −1 2π 2π = ∑ u (k ) ⋅ exp − j m ⋅ k ⋅ exp − j N ⋅k = N N k =0
105
N −1 2π = ∑ u (k ) ⋅ exp − j m ⋅ k = v ( m) N k =0
(4.61)
deoarece: 2π exp − j N ⋅ k = cos(2π ⋅ k ) − j ⋅ sin( 2π ⋅ k ) = 1 N 3. În vectorul transformat v, componentele sunt conjugat simetrice: N N N v ± m = v * m m , unde m = 0, − 1 2 2 2
(4.62)
Grafic: x*
x
0
N/2
N-1
Figura 4.2. Componentele vectorului transformat v sunt conjugat simetrice. 4. DFT-1D este spectrul eşantionat al semnalului discret u(k) extins cu 0 în afara intervalului [0,N-1]. Modelul matematic al semnalului eşantionat, extins cu 0 în rest este: N −1
u (k ) = ∑ u (t ) ⋅ δ (t − k )
(4.63)
k =0
A
ω
Figura 4.3. Spectrul eşantionat al semnalului discret u(k) extins cu 0. Spectrul semnalului este: 106
∞ N −1 ∞ U (ω ) = ∫ u (t ) ⋅ e − j⋅ω ⋅t dt = ∫ ∑ u (k ) ⋅ δ (t − k ) ⋅ e − j⋅ω⋅t dt = − ∞ k =0 −∞ N −1
∞
N −1
k =0
−∞
k =0
= ∑ u (k ) ⋅ ∫ δ (t − k ) ⋅ e − j⋅ω⋅t dt = ∑ u (k ) ⋅ e − j⋅ω⋅k
(4.64)
Comparând această relaţie cu transformata Fourier discretă: N −1
v ( m) = ∑ u ( k ) ⋅ e
− j⋅
2π m⋅k N
k =0
⇒ v(m) = U (ω ) ω = 2π m
(4.65)
N
A
v(1) v(2)
v(0)
-2π
2π 2π … 2 N N
-π
π
…
2π
ω
Figura 4.4. Spectrul semnalului. 5. Algoritmul rapid FFT al DFT (engl. FFT = Fast Fourier Transform) reduce complexitatea algoritmului de calcul al DFT de la N2 la N ⋅ log 2 N operaţii, unde N=2p este par:
( )FFT
O N 2 → O(log 2 N ) , dacă N=2p.
6. Transformata DFT transformă convoluţia circulară în produs. Convoluţia circulară este: not .
X (l ) ∗c Y (l ) = Z (l ) , unde l = 0, N − 1 107
(4.66)
def . N −1
Z (l ) = ∑ X ((l − i ) mod N ) ⋅ Y (i ) i =0
(4.67)
Acest mod de a defini simetria circulară limitează lungimea secvenţei Z la N. Teorema convoluţiei circulare: Dacă: X (m) = DFT {X (l )} Y (m) = DFT {Y (l )} şi Z (l ) = X (l ) ∗c Y (l ) Z (m) = DFT {Z (l )}
(4.68)
⇒ Z ( m ) = X ( m) ⋅ Y ( m )
4.7. Transformata Fourier discretă bidimensională (DFT-2D) Transformata Fourier discretă bidimensională (DFT – 2D) este una din transformările cel mai des utilizate în domeniul prelucrărilor de imagini. În mod similar cazului unidimensional, DFT – 2D este utilă la analiza spectrală a imaginilor digitale. Transformata Fourier discretă bidimensională directă a unei imagini discrete U = {u (l , k ), l,k = 0 ,1,...,N-1} se defineşte astfel: N −1 N −1 2π (m ⋅ l + n ⋅ k ) = v(m, n) = ∑ ∑ u (l , k ) ⋅ exp − j N l =0 k =0 N −1 N −1
= ∑ ∑ u (l , k ) ⋅ a mn (l , k ), 0 ≤ m,n ≤ N-1 l =0 k =0
2π (m ⋅ l + n ⋅ k ) . unde: a mn (l , k ) = exp − j N 108
(4.69)
După cum se observă, această relaţie se mai poate scrie: N −1 N −1 2π 2π v(m, n) = ∑ ∑ u (l , k ) ⋅ exp − j ⋅ m ⋅ l ⋅ exp − j ⋅n⋅k = N N l =0 k =0 N −1 N −1
= ∑ ∑ u (l , k ) ⋅ bm (l ) ⋅ c n (k ) = l =0 k =0
N −1 N −1
= ∑ ∑ u (l , k ) ⋅ W Nm⋅l ⋅ W Nn⋅k l =0 k =0
(4.70)
2π unde: bm (l ) = exp − j m⋅l N 2π c n ( k ) = exp − j n⋅k N 2π W N = exp − j N
Această rescriere a DFT – 2D directă, scoate în evidenţă proprietatea sa de separabilitate. Transformata DFT – 2D inversă este: u (l , k ) =
=
=
N −1 N −1
1 N
2
1 N
2
1 N
2
2π (m ⋅ l + n ⋅ k ) = N
∑ ∑ v(m, n) ⋅ exp j
m = 0 n =0
N −1 N −1
* * ∑ ∑ v(m, n) ⋅ bm (l ) ⋅ c n (k ) =
m =0 n = 0
N −1 N −1
− m⋅l − n⋅k ∑ ∑ v ( m, n ) ⋅ W N ⋅ W N
m =0 n = 0
109
(4.71)
unde bm(l), cn(k) şi WN au aceleaşi semnificaţii ca mai sus. Se poate arăta că transformarea astfel definită nu este o transformare unitară. Transformata DFT – 2D unitară directă se defineşte astfel: v ( m, n ) =
=
1 N −1 N −1 2π ∑ ∑ u (l , k ) ⋅ exp − j (m ⋅ l + n ⋅ k ) = N l =0 k =0 N 1 N −1 N −1 ∑ ∑ u (l , k ) ⋅ W Nm⋅l ⋅ W Nn⋅k N l =0 k = 0
(4.72)
DFT-2D unitară inversă este: u (l , k ) =
=
1 N −1 N −1 2π ∑ ∑ v(m, n) ⋅ exp − j (m ⋅ l + n ⋅ k ) = N m =0 n = 0 N 1 N −1 N −1 m⋅l n⋅k ∑ ∑ v ( m, n ) ⋅ W N ⋅ W N N m =0 n =0
(4.73)
Se poate defini matricea transformării DFT-2D unitare: 1 1 2π 2π ⋅ exp − j F = exp − j m ⋅ l n ⋅ k (4.74) N N m=0, N −1 N n =0, N −1 N l =0, N −1
k =0, N −1
(
)
care este o matrice unitară F ∗T ⋅ F = F ⋅ F ∗T = I N . Astfel, relaţiile ce definesc DFT – 2D unitară directă şi inversă se pot scrie şi sub formă matricială: • DFT – 2D unitară directă: V = F ⋅U ⋅ F T = F ⋅U ⋅ F , deoarece se observă mai sus că: F = F T 110
(4.75)
• DFT – 2D unitară inversă: U = F −1 ⋅ V ⋅ F −1 = F ∗ ⋅ V ⋅ F ∗ , deoarece: F ∗T ⋅ F = F ⋅ F ∗T = I N
(4.76)
⇒ F -1 = F ∗T = F ∗
4.8. Proprietăţi ale transformatei DFT-2D 1. Datorită separabilităţii sale, aplicarea DFT – 2D se poate face în 2 paşi: •
mai întâi se aplică DFT – 1D pe linii (sau pe coloane)
•
rezultatului i se aplică DFT – 1D pe coloane (sau pe linii). Datorită acestui fapt, se poate trage concluzia că DFT – 2D are algoritm rapid de calcul, deoarece DFT – 1D are algoritm rapid de calcul. Astfel, se poate arăta că algoritmul rapid necesar pentru
(
)
obţinerea DFT – 2D are complexitatea O N 2 ⋅ log N , în loc de
( )
O N 3 cât ar fi dacă nu ar fi separabilă. 2. Extensia DFT – 2D este periodică, adică:
v ( m + p ⋅ N , n + q ⋅ N ) = v ( m , n ) , ∀p , q ∈ N
(4.77)
3. Pentru o imagine u reală, elementele DFT – 2D sunt conjugat simetrice faţă de jumătatea imaginii transformate v (figura 4.5), adică: N N N N v ± m, ± n = v * m m, m n 2 2 2 2
111
(4.78)
N
N Figura 4.5. Elementele DFT – 2D sunt conjugat simetrice faţă de jumătatea imaginii transformate. De aici se poate trage concluzia că DFT – 2D este complet determinată de N2 numere reale şi nu de 2N2 câte ar trebui să aibă în mod normal (fiecare are parte reală şi imaginară). 4. Una dintre cele mai importante proprietăţi ale transformatei Fourier rezultă din aşa-numita ”teoremă a convoluţiei”. Conform acesteia, transformata DFT – 2D a produsului de convoluţie a două secvenţe bidimensionale x1 şi x2 este egală cu produsul simplu a transformatelor DFT – 2D a celor două secvenţe, adică: F {x1 (m, n) ⊗ x 2 (m, n)} = F {x1 (m, n)} ⋅ F {x 2 (m, n)}
(4.79)
unde cu F{x} s-a notat DFT–2D a secvenţei bidimensionale x şi cu x1 ⊗ x 2 s-a notat produsul de convoluţie a lui x1 cu x2, definit anterior. Această proprietate permite calculul convoluţiei a două secvenţe bidimensionale prin următoarea metodă: •
se calculează DFT–2D directă a celor 2 secvenţe, printr-un algoritm rapid de calcul (FFT=Fast Fourier Transform);
•
se efectuează produsul celor 2 transformate;
•
se calculează DFT–2D inversă a rezultatului. 112
Această metodă este frecvent aplicată la filtrarea digitală a imaginilor. DFT-2D este de asemenea utilă pentru determinarea unor parametri utilizaţi în procesele de analiză a imaginilor. 5. DFT-2D este spectrul eşantionat al semnalului eşantionat. Semnalul eşantionat se poate scrie: N −1 N −1
u (l , k ) = ∑ ∑ u ( x, y ) ⋅ δ ( x − l , y − k )
(4.80)
l = 0 k =0
Spectrul semnalului este: U (ξ ,η ) = ∫∫ u(x,y) ⋅ e − j (ξ ⋅ x +η ⋅ y )dxdy
(4.81)
R
Se poate arăta că: V (m, n) = U (ξ ,η ) ξ = 2π m,η = 2π n N
(4.82)
N
4.9. Transformata Cosinus discretă unidimensională Matricea transformării Cosinus discretă (DCT-1D) este C = (c m (l ) )m,l =1, N care are elementele: c m (l ) = Notând:
1
α ( m) =
1
N
, m = 0, l = 1, N
2 π (2l + 1) ⋅ m ⋅ cos , m, l = 1, N N 2N
N
unidimensională
(4.83)
,m = 0
2 , m = 1, N N
DCT-1D directă se poate rescrie: 113
(4.84)
N π (2l + 1) ⋅ m v(m) = ∑ u (l ) ⋅ α (m) ⋅ cos 2N l =1
(4.85)
DCT-1D inversă este: N π (2l + 1) ⋅ m u (l ) = ∑ v(m) ⋅ α (m) ⋅ cos 2N m=1
(4.86)
Se poate arăta că matricea C a transformării DCT-1D este unitară, adică C ⋅ C *T = C *T ⋅ C = I N . Pe de altă parte, deoarece matricea C este reală: ⇒ C −1 = C T * = C T
(4.87)
Observaţie: Transformarea cosinus discretă DCT-1D se poate obţine din DFT-1D a secvenţei simetrice de lungime 2N construită astfel: u ( N ) u ( N − 1) ... u (2) u (1) u (1) (2) ... u ( N − 1) u ( N ) (4.88) DCT-1D
DFT-1D
u(m)
u(m) …
…
u(1) u(2)
u(N-1) u(N)
…
u(N) u(N-1) u(2) u(1) u(1)u(2) u(N-1) u(N)
Figura 4.6. Obţinerea DCT-1D din DFT-1D de lungime 2N. sau: u (1) u (2) ... u ( N − 1) u ( N ) u ( N ) ( N − 1) ... u (2) u (1) (4.89)
114
u(m)
…
u(1) u(2)
…
u(N-1) u(N) u(N) u(N-1)
u(2) u(1)
Figura 4.7. Obţinerea DCT-1D din DFT-1D de lungime 2N. În acest fel, în cazul DCT-1D nu mai apar variaţii bruşte în spectru, variaţii care au ca rezultat componente importante la frecvenţe înalte, ca în cazul DFT-1D u (1) u (2) ... u ( N − 1) u ( N − 2) Consecinţă: Deoarece DCT-1D se poate construi din DFT-1D care are un algoritm de calcul rapid, rezultă că şi DCT-1D are un algoritm de calcul rapid, care reduce numărul de operaţii la 2 ⋅ ( N + 1) ⋅ log 2 ( N + 1) , adică ⇒ O( N ⋅ log 2 N ) .
4.10. Transformata Cosinus Discretă bidimensională DCT-2D Transformata Cosinus discretă bidimensională (DCT-2D) a unei imagini u(l,k) de dimensiune N×N (l,k= 1, N ) se defineşte astfel: N N (2l + 1)π (2k + 1)π ⋅ m ⋅ cos ⋅ n = v(m, n) = ∑ ∑ α (m) ⋅ α (n) ⋅ u(l , k ) ⋅ cos 2N 2N l =1k =1
N N
= ∑ ∑ u (l , k ) ⋅ a mn (l , k )
(4.90)
l =1k =1
115
unde:
m,n=1,2,...,N α ( m) =
1
, pentru m = 1
N
2 , pentru m = 1,2 ,...,N N (2l + 1)π (2k + 1)π ⋅ m ⋅ cos ⋅ n 2N 2N
amn (l , k ) = α (m) ⋅ α (n) ⋅ cos
(4.91)
(4.92)
În mod similar se defineşte şi DCT-2D inversă: N N (2l + 1)π (2k + 1)π u(l , k ) = ∑ ∑ α (m) ⋅ α (n) ⋅ v(m, n) ⋅ cos ⋅ m ⋅ cos ⋅ n = 2N 2N l =1k =1
N N
= ∑ ∑ u (l , k ) ⋅ a mn (l , k )
(4.93)
l =1k =1
iar α(m) şi α(n) au fost definite mai sus. Deoarece matricea transformării DCT-1D este: (2l + 1)π ⋅ m C (l , m) = α (m ) ⋅ cos 2N l ,m=1, N
(4.94)
DCT-2D directă şi inversă se pot scrie şi sub formă matriceală: V=C⋅U⋅CT = DCT-2D directă
(4.95)
( )−1 = C T ⋅ V ⋅ C = DCT-2D inversă, (4.96)
U = C −1 ⋅ V ⋅ C T deoarece C −1 = C T .
Deoarece DCT-2D este separabilă, rezultă că obţinerea DCT-2D se poate face în doi paşi: 116
•
întâi se aplică DCT-1D pe linii (sau pe coloane): N (2k + 1)π ⇒ v' (l , n) = ∑ α (n) ⋅ u (l , k ) ⋅ cos ⋅ n 2N k =1
•
(4.97)
rezultatului i se aplică DCT-1D pe coloane (sau pe linii): N −1 (2l + 1)π ⇒ v(m, n) = ∑ α (m) ⋅ v' (l , n) ⋅ cos ⋅ m 2N l =0
(4.98)
După cum s-a arătat, DCT-1D se poate obţine din DFT-1D a secvenţei simetrice de lungime 2N, construită astfel: u ( N ) u ( N − 1) ... u (2) u (1) u (1) (2) ... u ( N − 1) u ( N ) sau: u (1) u (2) ... u ( N − 1) u ( N ) u ( N ) ( N − 1) ... u (2) u (1) Ca şi în cazul unidimensional, se poate arăta că DCT-2D se poate obţine prin intermediul transformatei Fourier DFT. Astfel, DCT-2D se poate calcula ca fiind transformata Fourier a unei imagini u’ extinse, de dimensiune 2N×2N: 2N 2N
F (m, n) = ∑ ∑ u ' (l , k ) ⋅ W2l⋅Nm ⋅ W2kN⋅n , l =1 k =1
(4.99)
2π unde: W N = exp − j . N unde, imaginea extinsă u’ se obţine din imaginea iniţială u, astfel:
117
0 ≤ l ≤ N-1 pentru u (l , k ), 0 ≤ k ≤ N-1 N ≤ l ≤ 2 N-1 u (2 N − l − 1, k ), pentru 0 ≤ k ≤ N-1 u ' (l , k ) = 0 ≤ l ≤ N-1 u(l,2 N-k-1 ), pentru N ≤ k ≤ 2 N-1 N ≤ l ≤ 2 N-1 u (2 N − l − 1,2 N − k − 1), pentru N ≤ k ≤ 2 N-1 Relaţia care face legătura între cele două transformări este: v(m, n) = W2mN/ 2 ⋅ W2nN/ 2 ⋅ F (m, n)
(4.100)
(4.101)
În concluzie, calculul DCT-2D se reduce practic la calculul DFT-2D pe o imagine cu o arie de patru ori mai mare decât imaginea iniţială (figura 4.8): 2N N N
2N
Figura 4.8. Obţinerea DCT-2D a unei imagini, din DFT-2D pentru o imagine cu o arie de patru ori mai mare decât imaginea iniţială. Această variantă de obţinere a DCT-2D are dezavantajul că necesită un volum de calcul relativ mare datorită imaginii de dimensiuni mari (2N×2N) asupra căreia se aplică DFT-2D. O alta variantă de obţinere a DCT-2D prin intermediul DFT-2D şi care elimină dezavantajul prezentat mai sus, constă în utilizarea unei secvenţe (imagini) bidimensionale u’ obţinută din imaginea iniţială u, astfel: 118
N-1 0 ≤ l ≤ 2 u (l , k ), pentru 0 ≤ k ≤ N-1 2 N + 1 2 ≤ l ≤ N-1 pentru u (2 N − l − 1, k ), 0 ≤ k ≤ N-1 2 u ' (l , k ) = N-1 0 ≤ l ≤ 2 u(l,2 N-k-1 ), pentru N + 1 ≤ k ≤ N-1 2 N + 1 ≤ l ≤ N-1 u (2 N − l − 1,2 N − k − 1), pentru 2 N + 1 ≤ k ≤ N-1 2
(4.102)
Astfel, secvenţa (imaginea) u’(l,k) păstrează dimensiunea N×N a imaginii iniţiale. DCT –2D a imaginii iniţiale u se poate calcula în funcţie de DFT –2D a secvenţei u’(l,k) (care se va nota cu F(m,n)) astfel :
{
[
]}
{ [
]}
V (m, n) = 2 ⋅ Re W4mN ⋅ W4nN ⋅ F (m, n) + W4−Nn ⋅ F (m, N − n) = = 2 ⋅ Re W4nN W4mN ⋅ F (m, n) + W4−Nn ⋅ F (m, N − n)
(4.103)
Prin urmare, au fost prezentate 2 metode de obţinere a DCT-2D prin intermediul DFT-2D. În ambele cazuri, algoritmul de calcul al DCT–2D este următorul: •
se obţine imaginea (secvenţa) u’(l,k) din imaginea iniţială u(l,k);
•
se calculează DFT–2D a imaginii u’(l,k), care se notează cu F; 119
•
se calculează DCT–2D, utilizând relaţia de legătură corespunzătoare, între cele 2 transformate. De aici se poate trage şi concluzia că, deoarece DCT–2D se poate obţine prin intermediul DFT şi deoarece DFT are algoritm rapid de calcul, rezultă că şi DCT–2D are algoritm rapid de calcul.
4.11. Transformata Sinus discretă unidimensională DST-1D Matricea transformării este S = (s m (l ) )l ,m=1, N care are elementele: s m (l ) =
2 π ⋅ (m + 1) ⋅ (l + 1) ⋅ sin N +1 N +1
(4.104)
Pentru a arăta că transformarea S este unitară, trebuie verificate condiţiile de completitudine şi ortonormalitate. Se observă că S este o matrice reală şi simetrică, adică sk(m)=sm(k): ⇒ S −1 = S *T = S . Deoarece: V = S ⋅U este transformarea transformarea inversă este: ⇒ U = S −1 ⋅ V = S ⋅ V
(4.105) directă,
rezultă
că
(4.106)
Observaţie: Transformarea sinus unidimensională discretă DST-1D se poate obţine din transformarea Fourier unidimensională discretă DFT-1D, prin secvenţa antisimetrică de lungime 2·(N+1) construită astfel: (4.107) 0 − u ( N ) − u ( N − 1) K − u ( 2) − u (1) 0 u (1) u ( 2) K u ( N − 1) u ( N )
120
0 u(N) u(N-1) u(2) u(1) 0 …
… 0 u(1) u(2)
u(N-1) u(N) 0
Figura 4.9. Obţinerea DST-1D din DFT-1D, prin secvenţa antisimetrică de lungime 2·(N+1). Consecinţă: Deoarece DST-1D se poate construi din DFT-1D care are un algoritm de calcul rapid, rezultă că şi DST-1D are un algoritm de calcul rapid, care reduce numărul de operaţii la 2 ⋅ ( N + 1) ⋅ log 2 ( N + 1) , adică ⇒ O( N ⋅ log 2 N ) .
4.12. Transformata Sinus discretă bidimensională DST-2D Transformata Sinus discretă bidimensională DST-2D este transformarea bidimensională separabilă care are matricele A = B = S, unde S este matricea DST-1D definită anterior. Transformata Sinus directă bidimensională este definită prin relaţia: V = S ⋅ U ⋅ S T = S ⋅ U ⋅ S , deoarece S T = S
(4.108)
iar transformata Sinus inversă este definită prin relaţia: U = S −1 ⋅ V ⋅ S −1 = S ⋅ V ⋅ S , deoarece S −1 = S
(4.109)
Observaţie: Calculul DST-2D se poate face în 2 paşi: se aplică DST-1D pe linii, iar rezultatului i se aplică DST-1D pe coloane sau invers. 121
5. Restaurarea imaginilor O imagine poate fi degradată pe parcursul achiziţiei, transmisiei, prelucrării sau analizei sale. Problema restaurării unei imagini degradate se pune astfel: având o imagine digitală originală f(l,k) care a suferit un proces de degradare echivalent cu o filtrare h(l,k) se obţine o imagine degradată f’(l,k), peste care se suprapune un zgomot z(l,k), rezultatul fiind imaginea degradată f”(l,k). Prin restaurare se doreşte determinarea unei metode de estimare care pornind de la imaginea degradată f”(l,k) să conducă la un rezultat cât mai apropiat de imaginea originală f(l,k), relativ la un anumit criteriu, care poate fi, de exemplu, eroarea medie pătratică. În general, în modelul degradării (figura 5.1) se acceptă o filtrare liniară h(l,k).
f(l,k)
h(l,k)
f'(l,k)
+
f"(l,k)
z(l,k)
Figura 5.1. Modelul degradării unei imagini. De exemplu, funcţiei h(l,k) îi poate corespunde nefocalizarea obiectivului camerei sau mişcării aparatului foto în timpul expunerii. Se doreşte ca pornind de la imaginea degradată f”(l,k) să se obţină estimarea imaginii originale tot printr-o filtrare liniară g(l,k) deoarece în acest caz sunt calcule mai puţine: f"(l,k)
g(l,k)
fˆ (l , k )
fˆ (l , k ) ≅ f (l , k )
Criteriul de comparare va fi ales cel al erorii medii pătratice: 122
L K
2 ε = ∑ ∑ f (l , k ) − fˆ (l , k ) , unde f (l , k ) ∈ C
(5.1)
l =1k =1
L K
f ' (l , k ) = f (l , k ) ∗ h(l , k ) = ∑ ∑ f (l − i, k − j ) ⋅ h(i, j ) (5.2) i =1 j =1
f " (l , k ) = f ' (l , k ) + z (l , k )
(5.3)
5.1. Filtrarea inversă Dacă zgomotul ar fi inexistent (z=0) cea mai simplă metodă de restaurare ar fi ca g să fie filtrul invers lui h. Se va nota cu F(m,n) transformata Fourier discretă a funcţiei f(l,k): not .
F (m, n) = DFT { f (l , k )}
(5.4)
F " (m, n) = F (m, n) ⋅ H (m, n) + Z (m, n)
(5.5)
Dacă: z (l , k ) = 0 , ⇒ Z (m, n) = 0
(5.6)
F " (m, n) = F (m, n) ⋅ H (m, n)
(5.7)
Funcţia de transfer a filtrului de restaurare va fi, în acest caz: G (m, n) = H −1 (m, n)
(5.8)
Transformata Fourier a estimării imaginii iniţiale va fi: Fˆ (m, n) = F " (m, n) ⋅ H −1 (m, n) Pentru z=0: ⇒ Fˆ (m, n) = F (m, n) ⋅ H (m, n) ⋅ H −1 (m, n)
(5.9) (5.10)
Dacă zgomotul este nenul: z ≠ 0 ⇒ Z (m, n) ≠ 0 . În acest caz, se consideră că raportul semnal-zgomot (RSZ): • pentru semnale audio, un RSZ bun este un RSZ >20 dB • pentru imagini, un RSZ bun este un RSZ >30 dB În prezenţa zgomotului, relaţiile anterioare se scriu: 123
Fˆ (m, n) = [F (m, n) ⋅ H (m, n) + Z (m, n)] ⋅ H −1 (m, n)
(5.11)
⇒ Fˆ = F + Z ⋅ H −1
(5.12)
Dacă H are zerouri, rezultă că H-1 are poli. Rezultă că Z ⋅ H −1 are valori mari în vecinătatea polilor. Exemplu: în cazul unidimensional (1D): H-1
Fˆ = F + Z ⋅ H −1
H
F
Figura 5.2. Exemplu de filtrare inversă în cazul 1D. Dacă zgomotul z este zgomot alb, rezultă că Z(m,n) este constant. ⇒ Z ⋅ H −1 are valori mari în vecinătatea zerourilor lui H. ⇒ Fˆ este mult diferit de F ⇒ fˆ este mult diferit de f.
fˆ f
Figura 5.3. Semnalul original şi semnalul restaurat prin filtrare inversă, din semnalul original degradat cu zgomot alb.
124
În figura 5.3 fˆ reprezintă rezultatul restaurării, adică a estimării semnalului original f.
5.2. Filtrul invers cu constrângeri Lanţul complet al procesului de degradare şi restaurare este reprezentat în figura 5.4. f
f'
h
+
degradare
f"
z
g
fˆ
restaurare
Figura 5.4. Lanţul degradare-restaurare al unei imagini. Dacă s-ar filtra fˆ cu h s-ar obţine ceva apropiat de f’. Dar diferenţa dintre fˆ filtrat cu h ( fˆ ∗ h ) şi f” este zgomotul z. Energia zgomotului este presupusă a fi cunoscută: L
K
2
E z = ∑ ∑ Z (m, n) = cunoscută m =1 n =1
(5.13)
Prin urmare, energia diferenţei dintre fˆ ∗ h şi f” se doreşte a fi egală cu energia zgomotului: E f "− fˆ ∗h = E z
(5.14) L K
Dar: fˆ (l , k ) ∗ h(l , k ) = ∑ ∑ f (l − i, k − j ) ⋅ h(i, j ) i =1 j =1
L K
2 ⇒ ∑ ∑ f " (l , k ) − fˆ (l , k ) ∗ h(l , k ) = E z l =1K =1
125
(5.15) (5.16)
Se presupune cunoscut filtrul de degradare h. În acest caz rezultă un sistem de ecuaţii cu L×K necunoscute, pornind de la f(l,k), unde l = 1, L , k = 1, K .
Se impune ca energia derivatei lui fˆ să fie minimă, pentru ca abaterile faţă de f să fie cât mai mici. Pentru aceasta se consideră un nucleu de filtrare c=c(l,k) care să reprezinte o măsură a derivatei ( de ex. un laplacean sau un gaussian). Operaţia de filtrare cu nucleul c este echivalentă cu o convoluţie discretă, proporţională cu derivata: fˆ (l , k ) ∗ c(l , k ) , în domeniul spaţial (5.17) Fˆ (m, n) ∗ C (m, n) , în domeniul spectral
(5.18)
Se doreşte ca E fˆ ∗c să fie minimă, adică: L K
2 ∑ ∑ fˆ (l , k ) ∗ c(l , k ) să fie minimă.
(5.19)
l =1k =1
Prin urmare, trebuie minimizată relaţia (5.19) cu constrângerea (5.16). Din teorema lui Parceval, se ştie că: L K
1D
2 E f = ∑ ∑ fˆ (k , l ) → l =1k =1
L
K
∑ ∑ F (m, n)
2
m=1 n =1
Aplicând teorema lui Parceval relaţiei (5), după aplicarea prealabilă a transformatei Fourier, se obţine: (5.16) ⇔
L K 2 1 ∑ ∑ F " (m, n) − Fˆ (m, n) ⋅ H (m, n) = E z (5.20) L ⋅ K m=1, n=1
(5.19) ⇔
L K 2 1 ∑ ∑ Fˆ (m, n) ⋅ C (m, n) = minimă (5.21) L ⋅ K m=1, n=1
sau: L
K
2 ⇔ ∑ ∑ F " (m, n) − Fˆ (m, n) ⋅ H (m, n) = L ⋅ K ⋅ E z = T m=1, n =1
126
(5.22)
L
K
2 ⇔ ∑ ∑ Fˆ (m, n) ⋅ C (m, n) = minimă m=1, n =1
(5.23)
Se ştie că atunci când trebuie minimizată o funcţie f(x1, …, xn) cu constrângerile: g1 ( x1 ,..., x n ) = 0 M g ( x ,..., x ) = 0 n m 1
(5.24)
trebuie construită funcţia Lagrange: Ψ ( x1 ,..., x n ) = f ( x1 ,..., x n ) − λ1 ⋅ g1 ( x1 ,..., x n ) − ... − λ m ⋅ g m ( x1 ,..., x n ) după care se minimizează funcţia Lagrange, pornind de la relaţia uzuală: ∂Ψ = 0, (5.25) ∂xi unde i = 1, n , iar λi sunt coeficienţii Lagrange. În cazul de faţă, constrângerea este constituită de relaţia (5.16) iar funcţia de minimizat este dată de relaţia (5.19). Se construieşte funcţia Lagrange: L K 2 2 Ψ = ∑ ∑ Fˆ (m, n) ⋅ C (m, n) − λ ⋅ F " (m, n) − Fˆ (m, n) ⋅ H (m, n) − T m=1, n =1 (5.26) L K 2 2 ⇒ Ψ = ∑ ∑ Fˆ (m, n) ⋅ C (m, n) − λ ⋅ F " (m, n) − Fˆ (m, n) ⋅ H ( m, n) + λ ⋅ T m =1, n =1
Funcţiile Fˆ (m, n) sunt argumentele funcţiei Ψ, în număr de L×K. Fˆ (m, n) = A(m, n) + j ⋅ B(m, n)
(5.27)
Rezultă că sunt 2×L×K necunoscute, adică argumentele funcţiei Ψ. 127
Se obţine sistemul: ∂Ψ ∂A(m, n) = 0 , m = 1, L, n = 1, K ∂Ψ = 0 ∂B (m, n)
(5.28)
Sistemul are 2×L×K ecuaţii cu 2×L×K necunoscute. F "− Fˆ ⋅ H
2
)(
(
)
= F "− Fˆ ⋅ H ⋅ F "* − Fˆ * ⋅ H * = 2
= F " − F "⋅Fˆ * ⋅ H * − Fˆ ⋅ H ⋅ F "* + Fˆ
⇒
[
∂Ψ 2 = 2 ⋅ A(m, n) ⋅ C (m, n) − ∂A(m, n)
2
⋅H
2
(5.29)
(5.30)
]
2 − λ ⋅ F " (m, n) ⋅ H * (m, n) − Fˆ "* (m, n) ⋅ H (m, n) − 2 ⋅ A(m, n) ⋅ H (m, n) = 0
⇒
[
∂Ψ 2 = 2 ⋅ B ( m, n ) ⋅ C ( m, n ) − ∂B ( m, n)
(5.31)
− λ ⋅ − j ⋅ F " (m, n) ⋅ H * (m, n) + j ⋅ F "* (m, n) ⋅ H (m, n) − 2 ⋅ B(m, n) ⋅ H (m, n)
Rezultă: A(m, n) =
B ( m, n ) =
{
C (m, n) + λ ⋅ H (m, n)
(5.32)
}
(5.33)
2
{
λ ⋅ Im F " (m, n) ⋅ H * (m, n) 2
C (m, n) + λ ⋅ H (m, n)
2
⇒ Fˆ (m , n ) = A(m, n) + j ⋅ B(m, n) =
128
]= 0
}
λ ⋅ Re F " (m, n) ⋅ H * (m, n) 2
2
=
λ ⋅ F " (m, n) ⋅ H * (m, n) 2
C (m, n) + λ ⋅ H (m, n)
⇒ G (m, n) =
2
= G ( m, n ) ⋅ F " ( m, n )
λ ⋅ H * ( m, n ) 2
C (m, n) + λ ⋅ H (m, n)
2
(5.34)
(5.35)
G este funcţia de transfer a filtrului de restaurare, cu necunoscuta λ. S-a arătat că: L
K
2
∑ ∑ F " (m, n) − Fˆ (m, n) ⋅ H (m, n) = L ⋅ K ⋅ E z = T
m=1, n =1
(5.36)
Relaţiile (5.35) şi (5.36) formează un sistem cu L×K+1 ecuaţii cu L×K+1 necunoscute, constituite de λ plus necunoscutele lui Fˆ (m, n) , în număr de L×K. Sistemul format din relaţiile (5.35) şi (5.36) este neliniar, iar rezolvarea sa se poate face prin metode numerice, acesta fiind dezavantajul filtrului invers cu constrângeri.
129
6. Morfologie matematică
Morfologia matematică (în limba greacă morphos = formă, logos = ştiinţă, deci ştiinţa formelor) constă într-o abordare bazată pe formă, a prelucrării imaginilor. Ideea de bază a unei prelucrări morfologice este considerarea imaginii ca un ansamblu (mulţime, reuniune de puncte sau obiecte) asupra căruia se aplică transformări a căror esenţă este comparaţia cu mulţimi mai simple, numite elemente structurante. Prin urmare, transformările morfologice (care sunt neliniare şi neinversabile) se bazează pe compararea imaginii (sau a unui obiect conţinut în imagine) cu un obiect mai mic, de formă cunoscută, numit element structurant. În urma acestei comparaţii sunt extrase din imaginea iniţială zonele ce corespund proprietăţilor (de formă şi dimensiune) specifice elementului structurant folosit. De exemplu, recunoaşterea unei forme implică identificarea locală a părţilor sale componente, deci o simplă operaţie de potrivire de măşti ("pattern matching").
6.1. Transformarea Hit or Miss Transformarea morfologică de bază este transformarea “Hit or Miss”, ce ar putea fi numită şi “Totul sau Nimic” (sau “Ochit sau Ratat”, într-o traducere cuvânt cu cuvânt, din limba engleză). Efectul aplicării acestei transformări de identificare este extragerea din imagine a punctelor a căror vecinătate este identică cu elementul structurant folosit. Transformarea “Hit or Miss” a mulţimii A prin elementul structurant B se defineşte ca fiind:
{
A ∗ B = x | ( B1 ) x ⊂ A & ( B2 ) x ⊂ A c
unde: A c este complementara mulţimii A; 130
}
(6.1)
B1 şi B2 formează o partiţie netă a
lui B, adică:
B1 ∪ B2 = B şi B1 ∩ B2 = Φ ;
(6.2)
B x = {b + x b ∈ B}= translaţia mulţimii B cu vectorul x, sau
translaţia mulţimii B cu originea, în punctul x (figura 6.1). y
Bx
x B
x
Figura 6.1. Translaţia mulţimii B cu vectorul x. Trebuie specificat faptul că oricărui element structurant trebuie să i se ataşeze o origine. Exemplu: Să se efectueze transformarea “Hit or Miss” a mulţimii A prin elementul structurant B, unde:
• •
x
x
• • A=• x
x x
x x
• x x x • • • •
originea lui B originea nu aparţine lui B2 x x • • x • x x x ⊗ • x • x , B2 = x • • , B = x x • , B1 = x x • • • x x x ⊗ • •
•
•
•
A este reprezentată de punctele marcate cu x, iar cu “•” au fost marcate punctele care aparţin fundalului.
131
Rezultatul transformării este dat de cele 2 puncte marcate (încercuite) pe mulţimea A (punctele peste care se suprapun perfect B1 şi B2 , adică B).
Transformarea “Hit or Miss” prezintă un interes mai mult teoretic, dar datorită structurii sale stă la baza construcţiei teoretice a morfologiei matematice. Pe baza acestei transformări se pot defini operaţiile morfologice fundamentale (erodarea şi dilatarea).
6.2. Erodarea Erodarea mulţimii A prin elementul structurant B se defineşte ca: AΘB = x (B ) x ⊂ A = A ∗ B B2 =∅ B1 = B
{
}
(6.3)
După cum se observă din definiţie, erodarea se poate obţine ca un caz particular al transformării “Hit or Miss” şi anume pentru B2 = ∅ şi B1 = B .
Erodarea unei mulţimi A prin elementul structurant B este mulţimea punctelor, pentru care elementul structurant translatat cu originea în punctul respectiv este inclus în mulţimea ce se erodează. Rezultatul transformării se numeşte erodata (sau eroziunea) mulţimii A prin elementul structurant B. Efectul general al erodării este acela de subţiere a corpurilor, subţiere care depinde de structura elementului structurant. Mulţimile considerate pot fi continue sau discrete. Exemplu: - în cazul mulţimilor continue: Dacă: 132
B=
A=
⇒
AΘB=
A=
⇒
AΘB=
- în cazul mulţimilor discrete: • • AΘB = x •
• • x x ⊗ ⊗
• x x
x
•
x
• • x • • ⊗ ⊗ x AΘB = ⊗ ⊗ ⊗ x •
⊗
x
x B=x
x x
B= x
x
x
•
A este reprezentată de punctele marcate cu x, iar cu “•” au fost marcate punctele care aparţin fundalului. Punctele încercuite sunt
rezultatele erodării AΘB.
133
6.3. Dilatarea Dilatarea unei mulţimi A prin elementul structurant B se defineşte prin relaţia: A ⊕ B = {x | ( B ) x ∩ A ≠ ∅} = A * B | B1 =∅ B2 = B
C
(6.4)
Dilatarea unei mulţimi A prin elementul structurant B este mulţimea punctelor pentru care elementul structurant deplasat cu originea în punctul respectiv are puncte comune (cel puţin unul) cu mulţimea A ce se dilată. Rezultatul transformării se numeşte dilatată (sau dilatarea) mulţimii A prin elementul structurant B. Efectul general al operaţiei de dilatare este acela de îngroşare a obiectelor. Transformarea poate fi aplicată atât pe mulţimi continue, cât şi pe mulţimi discrete. Exemplu: - în cazul mulţimilor continue:
B=
A=
A ⊕ B=
A=
A ⊕ B=
134
- în cazul mulţimilor discrete: • •
x B=x
x x
A⊕ B =
x
• •
• A⊕ B = •
x
• Θ • Θ ⊗ Θ
• •
• Θ ⊗ ⊗ ⊗ Θ Θ ⊗ ⊗ ⊗ ⊗ Θ
• B= x
• •
Θ ⊗ Θ Θ • Θ • • •
•
•
•
• • •
• Θ ⊗ • • Θ ⊗ ⊗ ⊗ •
Θ ⊗ ⊗ ⊗ ⊗ • • Θ ⊗ • • • A este reprezentată de punctele marcate cu x, iar cu “•” au fost marcate punctele care aparţin fundalului. Punctele încercuite reprezintă rezultatele dilatării A ⊕ B . Dilatarea şi erodarea nu sunt transformări inverse una alteia (şi nici nu admit inversă): ( A ⊕ B )ΘB ≠ A (6.5)
( AΘB ) ⊕ B ≠ A Dilatarea este o operaţie extensivă
(A ⊆ A ⊕ B),
în timp ce
eroziunea este o operaţie antiextensivă (A ⊇ AΘB) numai în cazul folosirii elementelor structurante ce îşi conţin originea. Elementele structurante clasice sunt variantele discretizate ale discului unitar, vecinătatea V4 şi vecinătatea V8 : x V4 = x
x x
x
135
x
x
x
V8 = x x
x x
x x
Dilatarea şi erodarea sunt operaţii intensive din punct de vedere matematic, deoarece este evident că aplicarea transformărilor morfologice implică verificarea condiţiilor de definiţie pentru fiecare punct al imaginii (evitând evidentele efecte de margine), deci complexitatea algoritmică este comparabilă cu a unei operaţii de filtrare în domeniul spaţial. Pentru evitarea sau micşorarea complexităţii calculelor ar trebui găsită o metodă care să nu implice verificarea fiecărui punct al imaginii, ci eventual fiecare punct al structurantului. Aceasta se poate obţine prin rescrierea operaţiei de erodare: AΘB = {x ∀b ∈ B, ∃a ∈ A, a.î . b + x = a}=
{x ∀b ∈ B,
∃a ∈ A, a.î . x = a − b} = I A−b = I Ab b∈B
b∈B S
(6.6)
unde B-S=-B este simetrica mulţimii B faţă de origine, iar Ab este translaţia mulţimii A în punctul b. Astfel, se fac operaţii mult mai puţine şi în plus se permite implementarea paralelă deoarece fiecare translaţie se poate face pe unităţi diferite, iar la sfârşit se face intersecţia rezultatelor acestora. În mod similar, pentru dilatare se obţine: A ⊕ B = {x | B x I A ≠ ∅} = U Ab (6.7) b∈B S
6.4. Proprietăţile operaţiilor morfologice 1. Invarianţa la translaţie: At ⊕ B = ( A ⊕ B )t
(6.8)
At ΘB = ( AΘB )t
136
A ⊕ Bt = ( A ⊕ B )−t În acest caz intervine semnul ”-” deoarece în
definiţia erodării şi dilatării intervine mulţimea simetrică AΘBt = ( AΘB )−t 2. Invarianţa la scalare: 1 A ⊕ B = (λ ⋅ A ⊕ B )
(6.9)
λ
AΘB =
1
λ
(λ ⋅ AΘB ) , λ ∈ [0,1]
3. Erodarea şi dilatarea sunt transformări duale una alteia:
(
)C C AΘB = (AC ⊕ B )
A ⊕ B = AC ΘB
(6.10)
Această dualitate se manifestă faţă de operaţia de complementare a mulţimilor. 4. Dilatarea şi erodarea nu sunt inverse una alteia şi nici nu admit inversă: ( AΘB ) ⊕ B ≠ A (6.11)
( A ⊕ B )ΘB ≠ A 5. Descompunerea elementului structurant: ( A ⊕ B ) U ( A ⊕ C ) = A ⊕ (B U C )
(6.12)
( AΘB ) I ( AΘC ) = AΘ(B U C )
6. Proprietatea de iterare:
( A ⊕ B ) ⊕ C S = A ⊕ (B ⊕ C ) ( AΘB )ΘC S = AΘ(BΘC ) Caz particular: B=C ; B=BS
⇒ ( A ⊕ B ) ⊕ B = A ⊕ 2B , 137
(6.13)
( AΘB )ΘB = AΘ2 B
6.5. Transformări morfologice derivate Prin iterarea unor operaţii morfologice de bază se obţin transformări morfologice mai complexe, numite şi transformări morfologice derivate sau filtre morfologice.
6.5.1. Operatori de extragere a conturului Dintre extractoarele morfologice de contur, cele mai utilizate sunt: • conturul exterior: ∆A = ( A ⊕ B ) − A ,
(6.14)
unde ”-” reprezintă diferenţa între mulţimi. • conturul interior: δA = A − ( AΘB )
(6.15)
• gradientul morfologic: gradA = ( A ⊕ B ) − ( AΘB )
(6.16)
Exemplul 1:
A= AΘB
• •
•
•
• •
• •
•
x
• •
• •
x
⊗ x •
x B=x
• x ⊗ ⊗ x • • • • •
x •
x •
• • • •
δA = conturul interior
∆A = conturul exterior 138
x x
x
În acest caz (structurant simetric), gradientul morfologic va fi reuniunea celor două contururi. Cu cât structurantul este mai mic, conturul va fi mai subţire. Exemplul 2: • • • • • •
A=
• • •
x • •
• •
x
x
x •
• x x
x
x •
• • x x • • • • • • • •
B=x
x
x
∆A = conturul exterior
δA = conturul interior În acest caz, structurantul este simetric după o singură direcţie. Exemplul 3:
• • • A= • •
• • • • • • • x • • • x x x • B= x x x x x x • • x x • • δA = conturul interior
• • • • • • ∆A = conturul exterior În acest caz, structurantul este asimetric şi se pierde conturul exterior, corespunzător direcţiei de simetrie, din cazul anterior. Prin urmare, pentru obţinerea de contururi direcţionale, trebuie utilizate elemente structurante direcţionale.
139
6.5.2. Deschiderea şi închiderea Deschiderea mulţimii A prin elementul structurant B se defineşte ca fiind:
A o B = ( AΘB ) ⊕ B S unde
(6.17)
B S = − B = {x − x ∈ B}, reprezintă mulţimea simetrică a
mulţimii B faţă de origine (se mai numeşte şi transpusa mulţimii B). Mulţimea rezultată după o deschidere este diferită de mulţimea originală. Din punct de vedere al acţiunii asupra obiectelor (figura 6.2), în urma unei deschideri rezultă (efectul deschiderii este reprezentat cu linie punctată) o lărgire a golurilor înglobate în obiect, eliminarea obiectelor mici (mai mici decât elementul structurant folosit), netezirea contururilor prin teşirea convexităţilor şi separarea obiectelor unite prin “istmuri” (zone) mai mici decât dimensiunea elementului structurant.
deschidere închidere
Figura 6.2. Ilustrarea efectului deschiderii şi închiderii morfologice. Închiderea mulţimii A prin elementul structurant B este: A • B = ( A ⊕ B )ΘB S
(6.18)
Închiderea este operaţia duală deschiderii faţă de complementarea mulţimilor: 140
(AC • B )C = ((AC ⊕ B )ΘB S )C = (( AΘB )C ΘB S ) = ( AΘB ) ⊕ B S = A o B (6.19)
(AC o B)C = ((AC ΘB)⊕ B S )C = (( A ⊕ B)C ⊕ B S ) = ( A ⊕ B)ΘB S = A • B
Datorită acestei proprietăţi de dualitate, închiderea va avea asupra fundalului aceleaşi efecte pe care le are deschiderea asupra obiectelor. Prin urmare, efectul închiderii asupra obiectelor (reprezentat cu linie îngroşată în figura 6.2) poate fi găsit prin complementare, adică închiderea va umple golurile înglobate în obiecte, va netezi contururile prin umplerea concavităţilor şi va fuziona obiectele foarte apropiate (umplerea strâmtorilor de dimensiuni mai mici decât a elementului structurant). Deschiderea este o operaţie antiextensivă ( A o B ⊆ A) , în timp ce închiderea este o transformare extensivă
( A • B ⊇ A) .
În acelaşi timp,
deschiderea şi închiderea sunt transformări idempotente, adică iterarea deschiderilor sau închiderilor succesive, cu acelaşi element structurant, nu mai produc modificări: (A o B) o B = A o B
(A • B) • B = A • B
(6.20)
După cum s-a arătat, închiderea şi deschiderea sunt de fapt nişte filtre deoarece au un efect de netezire a formelor. Pentru o netezire mai puternică se foloseşte un structurant mai mare. Mărirea gradată a elementului structurant folosit şi aplicarea alternativă a deschiderii şi închiderii (pentru a beneficia de efectele lor complementare) a dus la definirea filtrelor alternate secvenţial (FAS): FAS ( A) = ((((( A o B ) • B ) o 2 B ) • 2 B ) o 3B ) • 3B... (6.21) unde: k ⋅ B = B ⊕ B ⊕ B... ⊕ B (de k ori). Filtrele alternate secvenţial constau în deschideri şi închideri alternate succesiv, cu elemente structurante de dimensiune crescătoare. 141
Aplicarea FAS poate fi oprită în orice moment, obţinându-se astfel o netezire gradată, interactivă a imaginii.
6.6. Trierea dimensională a obiectelor Cu ajutorul transformărilor morfologice se poate realiza şi trierea dimensională a obiectelor. Aceasta se obţine prin utilizarea transformării “Top Hat”. Transformarea de tip “Top Hat” (TH) a mulţimii A prin elementul structurant B se defineşte ca: TH ( A) = A − ( A o B )
(6.22)
Rezultatul aplicării acestei transformări este o imagine ce conţine toate punctele eliminate de deschiderea imaginii prin elementul structurant folosit. După cum se poate observa, transformarea “Top Hat” are efect invers celui al unei site: păstrează obiectele mai mici. Dar, prin deschidere, obiectele nu îşi mai păstrează forma iniţială, deci, pentru extragerea formei exacte a obiectului trebuie realizată o operaţie de reconstrucţie. Există şi varianta transformării de tip “Top Hat” generalizată: TH g ( A) = ( A o n1 B ) − ( A o n2 B ) , n 2 > n1
(6.23)
Transformarea “Top Hat” simplă se obţine ca un caz particular: TH ( A) = THg ( A) n1 =0
n2 =1
142
(6.24)
Exemplu:
~4B ~2B A=
~B (obiecte cu dimensiuni mai mici decât B) TH(A)= A − ( A o B)
Ao B =
A o 2B =
TH(A)= A − ( A o 2 B)
Se observă că transformarea “Top Hat” simplă se comportă ca un filtru de tip trece-sus, iar “Top Hat” generalizat se comportă ca un filtru de tip trece-bandă (ne permite să obţinem obiecte cu dimensiune cuprinsă între n1 B şi n2 B ).
6.7. Caracterizarea morfologică a formelor Cu ajutorul transformărilor morfologice se poate realiza şi caracterizarea morfologică a formelor, care poate implica două aspecte: extragerea de informaţii asupra formei date sau compararea formei respective cu o altă formă (etalon). 143
6.7.1. Reconstrucţia după marker Reconstrucţia (φ) unei imagini (I) pe baza unor markeri (M), poate fi descrisă prin următoarea relaţie: ϕ M (I ) = U I j (6.25) I j I M ≠∅
Reconstrucţia imaginii I după markerul M este egală cu intersecţia elementelor Ij ale elementelor conexe ale imaginii I cu proprietatea că intersecţia lor cu markerul este nevidă. I1
I1 Ik
IN
IN Imaginea I
Markerii M
Figura 6.3. Reconstrucţia după markeri. O modalitate de determinare a unor markeri este de a utiliza transformarea Top Hat. O problemă care se pune în cazul reconstrucţiei după markeri este etichetarea imaginii, pentru a putea fi comparată cu markerii. Metoda de etichetare trebuie să fie cât mai rapidă.
144
6.7.2. Distanţa Haussdorf Pentru a se compara sau calcula asemănarea dintre două obiecte, se compară obiectul dat cu un obiect martor, se defineşte o distanţă între obiecte şi se impune ca distanţa să fie cât mai mică. Distanţa Hausdorff dintre mulţimile K1 şi K2, în domeniul continuu se defineşte ca fiind: d ( K1 , K 2 ) = inf {ε K1 ⊆ K 2 ⊕ ε ⋅ B & K 2 ⊆ K1 ⊕ ε ⋅ B}
(6.26)
În domeniul discret, distanţa Hausdorff dintre mulţimile K1 şi K2 este:
d ( K1 , K 2 ) = inf {n K1 ⊆ K 2 ⊕ n ⋅ B & K 2 ⊆ K1 ⊕ n ⋅ B} (6.27)
Distanţa Hausdorff este o măsură bună a asemănării obiectelor doar dacă acestea sunt centrate. În plus, este şi greu de calculat şi din acest motiv se foloseşte mai rar în practică.
6.7.3. Extragerea skeletonului morfologic Un concept important în prelucrarea şi analiza imaginilor, îndeosebi la recunoaşterea formelor, îl constituie skeletonul morfologic, cu aplicaţii la compresia imaginilor binare, recunoaşterea, aproximarea şi reconstrucţia formelor. Extragerea skeletonului morfologic se bazează pe conceptul de disc maximal inclus într-o mulţime (A). Astfel, pentru o mulţime binară plană închisă A, se defineşte discul maximal Br(x) de centru x şi rază r, prin relaţia:
x = x' B r ( x ) ⊆ Br ' ( x ' ) ⊆ A ⇔ r = r ' 145
(6.28)
Aceasta însemnă că discul maximal este inclus în mulţimea A şi nu există nici un alt disc inclus în A care să-l conţină Skeletonul morfologic al unei forme este egal cu reuniunea centrelor discurilor maximale incluse în forma respectivă. În practică se foloseşte formularea echivalentă: N max
SK ( A) = U S n ( A) ,
(6.29)
n =∅
unde SK(A) este skeletonul mulţimii A, iar Sn(A) se numesc seturi skeleton de ordinul n: S n ( A) = ( AΘnB ) − [( AΘ(n + 1)B ) ⊕ B ] = ( AΘnB ) − [( AΘnB )o B ] (6.30) Discul unitate, B, poate fi V4 sau V8, în funcţie de metrica utilizată în spaţiul discret, astfel:
• dacă se utilizează distanţa euclidiană, discul unitar este V4; • dacă se utilizează distanţa inter-bloc, discul unitar este V8. Distanţa inter-bloc, pentru punctul de coordonate (i, j) este:
(
d i ( x, y ) = max i x − i y , j x − j y
)
(6.31)
Exemplu: Se dă mulţimea A şi elementul structurant B:
A=
x x
x x
x x
x x
x x
x x
x x
x x
x , B=x
x x
x x
x x
x x
x x
x x
x x
x x
x
x
x
x
x
x
x
x
x x
x
Figura 6.4. Exemplul unei mulţimi A şi a uni element structurant B.
146
Setul skeleton de ordinul 0 este: S 0 ( A) = A − [( AΘ B ) ⊕ B ] :
S 0 ( A) =
x
x
x
x
x x
⊗ ⊗ ⊗ ⊗
x x
⊗ ⊗
x
(6.32)
- punctele încercuite = AΘB - punctele încadrate = ( AΘB ) ⊕ B
- punctele îngroşate: S 0 ( A) = A − ( AΘB ) ⊕ B x x ⊗ ⊗ ⊗ ⊗ x x x ⊗ ⊗ ⊗ ⊗ ⊗ ⊗ x x x x x x x x x x
Figura 6.5. Setul skeleton de ordinul 0 al mulţimii A din figura 6.4. În mod similar: S1 ( A) = ( AΘ B ) − [( AΘ2 B ) ⊕ B ], unde AΘ2 B = ( AΘB )ΘB
S 2 ( A) = ( AΘ2 B ) − [( AΘ3B ) ⊕ B ]
(6.33) (6.34)
S 3 ( A) = ( AΘ3B ) − [( AΘ4 B ) ⊕ B ] = ∅ ,
(6.35)
⇒ N max = 2 , pentru exemplul prezentat. N max
⇒ SK ( A) = U S n ( A)
(6.36)
n =∅
x
SK ( A) = x x
x
x
x x
x
x
x
x
x
x x
x x
Figura 6.6. Setul skeleton de ordinul 3 al mulţimii A din figura 6.4.
147
Se observă că skeletonul unui obiect are mai puţine puncte decât obiectul. Skeletonul morfologic este o transformare reversibilă, adică se poate obţine forma iniţială A, cunoscând skeletonul acesteia: N max
A = U (S n ( A) ⊕ nB )
(6.37)
n =0
Se observă că pentru a putea reconstitui forma iniţială, pentru fiecare punct al skeletonului, trebuie menţionat setul skeleton căruia îi aparţine:
0
0 1 1 1 1
SK ( A) = 0
2 2 1 1
0 1 1
0
0
Figura 6.7. Seturile skeleton al mulţimii A din figura 6.4. Dacă 2 obiecte au skeletoanele identice, atunci au aceeaşi formă. Detaliile fine sunt conţinute în seturile skeleton de ordin mic, iar cele grosiere în seturi skeleton de ordin mare. Deci, o reconstituire aproximativă Ã a formei A se poate face eliminând detaliile fine, adică seturile skeleton de ordin mic: N max
Ã= U (S n ( A) ⊕ nB )
(6.38)
n=k
În cazul continuu, skeletonul unui disc este centrul său, iar skeletonul unui pătrat este format din diagonalele sale (figura 6.8).
148
x
Figura 6.8. Skeletoanele unor forme continue cunoscute. Utilizarea skeletonului morfologic pentru recunoaşterea formelor este restricţionată de puternica sa sensibilitate la zgomote, deoarece o mică schimbare a formei duce la o modificare puternică a skeletonului său. De exemplu, skeletonul unui disc lipsit de centru se transformă dintrun punct într-un cerc (figura 6.9).
Figura 6.9. Skeletonul unui disc lipsit de centru său.
6.7.4. Skeletonul generalizat Skeletonul generalizat (GSK) este definit prin elemente structurante generalizate. Fie {E n } un set de mulţimi având perioada T, adică: E n+ kT = E n , ∀n, k ∈ Ζ . Pe baza acestui set generator (constructor) se poate construi un set de elemente structurante generalizate, pe baza relaţiilor: = elementul structurant de ordin 0 (originea) B0 = {0 n } (6.39) Bn = Bn −1 ⊕ E n = elementul structurant de ordin n
149
Exemplu: Pentru: T=1 şi E1 =
⇒ B0 = x
B1 =
x x
x x x x x x x x x x B3 = ... (6.40) x x x x
x x x
x x x x
B2 = x x x x x x
x x x x
Pentru extragerea skeletonului morfologic trebuie construită o hartă de distanţe a obiectului, adică fiecărui punct al obiectului (mulţimii) i se ataşează ordinul elementului structurant generalizat maximal centrat în punctul respectiv: 0 , dacã x ∉ A D(x ) = n , dacã (Bn-1 ) x ⊂ A si (Bn ) ⊄ A Exemplu: x x x x x x
A=
x x
4 3 2 1 4 3 2 1 4 3 2 1 4 3 2 1
x x
x x
x x
x x
D=
x x
x x
x x
x x
x x
x x
x x
x x
3 3 3 3 3 3 2 1 2 2 2 2 2 2 2 1
x
x
x
x
x
x
x
x
1 1 1 1 1 1 1 1
(6.41)
(6.42)
Skeletonul generalizat al unei forme este mulţimea centrelor elementelor structurante maximale incluse în formă, adică: GSK ( A) = x | B D ( x )−1 ⊄ B D ( y )−1 , ∀x ≠ y (6.43)
{
(
)x (
150
)y
}
4 3 2 1
GSK( A) =
4 3 2 1
4 3 2 1 4 3 2 1
4 3 2 1
4 3 2 1
=
4 3 2 1
4 3 2 1
3
3
3 3 3 3 2 1
3 3 3 3 3 3 2 1
2
2
2 2 2 2 2 1
2 2 2 2 2 2 2 1
1
1
1 1 1 1 1 1
1 1 1 1 1 1 1 1
(6.44)
După cum se observă, skeletonul generalizat conduce la obţinerea unei rate de compresie mai bune decât în cazul skeletonului morfologic (exemplele prezentate pun în evidenţă acest fapt: skeletonul morfologic e alcătuit din 16 puncte în timp ce skeletonul generalizat conţine doar 7 puncte). În ultima figură punctele încercuite reprezintă punctele suficiente pe baza cărora se poate reconstitui forma iniţială, exemplificându-se astfel că skeletonul generalizat nu este minimal din punct de vedere al ratei de compresie.
6.8. Extinderea morfologiei matematice la imagini cu niveluri de gri Prin extinderea morfologiei matematice la imagini cu niveluri de gri se realizează trecerea de la mulţimi la funcţii şi invers. Pentru aceasta se consideră o mulţime A inclusă în mulţimea părţilor lui Zn. Un element al acesteia este de forma: x = ( x1 ,x 2 ,...,x n−1 , x n ) (6.45) { 14 4244 3 domeniu spatial sup rafata A( n , n ) A( 1,n −1 )
151
6.7.1. Trecerea de la mulţime la funcţie Transformarea prin care se realizează trecerea de la o mulţime A la o funcţie f se numeşte vârf (în engleză Top): T [ A] = f ,
unde f : A(1, n − 1) → A(n, n)
f ( z ) = max{y ( z , y ) ∈ A} ,
unde z ∈ A(1, n − 1)
Pentru n=2 avem o mulţime planară: ( x1 , x 2 ) = (i, j )
(6.46) (6.47)
(6.48)
Dacă alegem i = domeniul spaţial, iar j = suprafaţa, topul s-ar obţine prin fixarea z ∈ A(1, n − 1) şi ar arăta ca în figura: j
x x x xxxx xxxxxx
i
Se fixează z∈A(1,n-1)
Figura 6.10. Topul unei mulţimi.
152
6.8.2. Trecerea de la funcţie la mulţime Transformarea prin care se realizează trecerea de la o funcţie f la o mulţime A se numeşte umbră: U[ f ] = A,
unde f : A(1, n − 1) ⊂ Z n−1 → A(n, n) ⊂ Z
z ∈ A(1, n − 1) U [ f ] = ( z , y ) , f ( z) ≥ y , A ( n , n ) Exemplu: j
x x x xxxx xxxxxx xxxxxx xxxxxx xxxxxx xxxxxx xxxxxx
y
(6.49) (6.50)
i
umbra unei funcţii este semiinfinită
Figura 6.11. Umbra unei funcţii. Deoarece umbra unei funcţii este o mulţime semiinfinită, se poate introduce o limitare, la un nivel λ: z ∈ A(1, n − 1) U λ [ f ] = ( z , y ) , f ( z) ≥ y ≥ λ y, A(n, n)
(6.51)
⇒ U [ f ] = lim U λ [ f ]
(6.52)
λ →∞
153
Proprietăţi: 1. T [U [ f ]] = f
(6.53)
2. U [T [ A]] ⊇ A ,
(6.54)
deoarece umbra unei funcţii este o mulţime semiinfinită
6.8.3. Lucrul cu funcţii Dilatarea unei funcţii f cu un structurant (sau funcţie structurant) g: f ⊕ g = T [U [ f ] ⊕ U [ g ]] (6.55) În mod similar, se poate defini erodarea unei funcţii f cu un structurant (sau funcţie structurant) g: fΘg = T [U [ f ]ΘU [ g ]] (6.56) Acestea sunt: f ⊕ g = sup( f ( x − y ) + g ( y )) , y∈supp ( g )
fΘg = inf( f ( x − y ) − g ( y )) y∈supp ( g )
(6.57) (6.58)
unde supp ( g ) = {y g ( y ) > −∞}, reprezintă suportul lui g. Exemplu:
suppV4: toate punctele lui V4 dar şi cele de sub V4 (până la -∞) aparţin suportului.
154
x x x x x
Figura 6.12. Suportul lui V4. Suportul plat (engl. flat) a funcţiei g se defineşte astfel: g ( y ) = 0, ∀y ∈ supp ( g ) ,
(6.59)
adică nu mai avem punctele de sub V4, în cazul exemplului de mai sus. În acest caz: f ⊕g = fΘg =
max ( f ( x − y ))
y∈supp ( g )
min
y∈supp ( g )
( f ( x − y ))
(6.60) (6.61)
În cazul în care se utilizează un structurant non-flat, pot să apară situaţii când valoarea minimă (min) şi maximă (max) nu sunt cuprinse în gama de niveluri de gri, adică sunt mai mici decât 0 sau mai mari decât N-1. În acest caz, pentru niveluri de gri cuprinse în intervalul [0,N-1], există două alternative: 1. se rescalează domeniul, păstrându-se nivelul relativ de gri:
155
N-1
0
N-1+k1
N-1+k2
Figura 6.13. Rescalarea domeniului. 2. se limitează valorile negative la 0 şi cele mai mari ca N-1 la N-1, în acest fel pierzându-se informaţie. Astfel de elemente structurante pot fi utilizate la îmbunătăţirea contrastului imaginilor. Pentru f1 ≤ f ≤ f 2 : f1 (m, n), f (m, n) ≤ kf1 + (1 − k ) f 2 , k∈[0,1] f ' (m, n) = f 2 (m, n), f (m, n) > kf1 + (1 − k ) f 2 f(t) f1
(6.62)
pentru k=0 sau 1 f
f2
pentru k = t
1 2
Figura 6.14. Îmbunătăţirea contrastului imaginilor. O îmbunătăţire a performanţelor se obţine luând: f1 = f o g şi f 2 = f • g sau f1 = fΘg şi f 2 = f ⊕ g
156
(6.63) (6.64)
7. Segmentarea imaginilor Segmentarea reprezintă o categorie de tehnici de prelucrare care permite descompunerea unei scene în componentele sale sau extragerea din imagini a unor elemente constituente de interes (obiecte, fundal etc) în scopul analizei lor ulterioare şi eventual al clasificării lor. În urma procesului de segmentare, din imagine se extrag regiuni omogene închise de puncte de frontieră (contur), obiecte distincte sau regiuni omogene care satisfac anumite criterii de uniformitate. Metodele de segmentare a imaginilor se pot clasifica în două mari clase: •
segmentare orientată pe regiuni
•
segmentare orientată pe contururi
7.1. Segmentarea orientată pe regiuni Prin detecţia regiunilor omogene se înţelege gruparea pixelilor din imagine în categorii distincte în funcţie de proprietăţile lor (de exemplu nivelul de gri), astfel încât să se pună în evidenţă regiuni caracterizate de o relativă uniformitate. Zonele astfel determinate permit în ultimă instanţă o separare a obiectului ce trebuie analizat de fundalul imaginii şi de eventuale alte obiecte aflate în scenă. Tehnicile de segmentare orientate pe detecţia regiunilor omogene se pot clasifica în: •
metode bazate pe etichetarea grupurilor conexe de pixeli cu caracteristici similare;
•
metode de segmentare bazate pe histograma imaginii;
•
tehnici de creştere şi fuziune a regiunilor.
157
7.1.1. Etichetarea componentelor Această metodă de segmentare a imaginilor binare constă în asocierea unui acelaşi număr (numit etichetă), tuturor punctelor unui obiect conex. Metoda bazată pe etichetarea componentelor se implementează prin succesiuni de baleieri normale (sus-jos, stânga-dreapta) şi inverse a imaginii de segmentat. 1 1 1 1 1
1 1 1 1 1
2 2 2 2
3
2 2 2 2
2 2
2
Figura 7.1. Etichetarea componentelor. La întâlnirea unui punct al unui obiect care nu are nici un punct vecin deja etichetat, i se atribuie o nouă etichetă. Se continuă baleierea până la întâlnirea unui alt punct, căruia i se atribuie eticheta vecinului (sus, jos, stânga, dreapta, diagonală), dacă are. Dacă are doi (sau mai mulţi) vecini cu etichete diferite i se atribuie eticheta cu valoarea cea mai mică. Decizia cărui obiect îi aparţine se face la baleierea inversă, respectând algoritmul de mai sus. Baleierea se repetă până nu se mai schimbă nimic. Dezavantajul acestei tehnici este că nu asigură obţinerea etichetelor în ordine (figura 7.2), iar numărul de baleieri depinde de conţinutul imaginii.
158
1111 4444 2222222
Figura 7.2. Exemplu de rezultat al segmentării prin etichetarea componentelor. O altă tehnică de a implementa metoda de etichetare a componentelor se bazează pe analiza conexităţii pe secvenţe. Pentru a implementa această metodă, se analizează imaginea, linie cu linie şi se etichetează punctele obiectelor, pe linii. În etapa următoare se analizează adiacenţa secvenţelor pentru a defini obictele. X X X X X X
I (2,3), (6,2) – pe linia I există un obiect pe poziţia 2, de lungime 3 X X X X II (1,3), (7,2) – pe linia II există un obiect pe poziţia 1, de lungime 3 X X X III (5,3) – pe linia III există un obiect pe poziţia 5, de lungime 3 X X IV (5,2) – pe linia IV există un obiect pe poziţia 5, de lungime 2
Figura 7.3. Exemplu de segmentare prin analiza conexităţii pe secvenţe.
7.1.2. Metoda aborelui cuaternar (quad-tree) Această metodă, care se mai numeşte şi "Divide şi contopeşte", se bazează pe împărţirea recursivă a imaginii în câte 4 regiuni sau sferturi de imagine până la obţinerea de regiuni uniforme sau regiuni formate dintr-un singur pixel. Acestei împărţiri i se poate asocia un arbore cuaternar, în care fiecare nod terminal are patru descendenţi.
159
Figura 7.4. Principiul segmentării prin metoda arborelui cuaternar. Nodul principal este constituit de întreaga imagine, iar nodurile secundare reprezintă câte un sfert de imagine. Împărţirea imaginii se repetă până când se obţin numai careuri uniforme (care conţin aceeaşi valoare).
Figura 7.5. Împărţirea imaginii la segmentarea prin metoda arborelui cuaternar. Pentru a implementa segmentarea prin metoda arborelui cuaternar, imaginea trebuie să fie pătrată şi de dimensiune putere a lui 2. Segmentarea se face prin împărţirea succesivă până se obţine zone uniforme sau s-a ajuns la nivel de pixel. În continuare se concatenează zonele ce conţin „1” logic, rezultând obiectele.
160
X X X X X
X X X X XX XX XX
X X X X
Figura 7.6. Împărţirea imaginii la segmentarea prin metoda arborelui cuaternar şi concatenarea zonelor cu aceeaşi valoare.
7.2. Segmentarea imaginilor cu niveluri de gri 7.2.1. Segmentarea bazată pe histogramă Aceste metode se bazează pe histograma imaginilor, adică pe numărul de apariţii a nivelurilor de gri. Tehnicile de segmentare bazate pe prăguirea (“tresholding”) histogramelor sunt utile şi eficiente atunci când există o separare relativ clară a nivelurilor de gri între obiectele analizate. Valorile caracteristice de amplitudine corespunzătoare obiectelor sunt alese astfel încât un interval dat de niveluri de gri să reprezinte o clasă unică de obiecte. Cea mai generală metodă de tresholding (multi-nivel) constă în alegerea unui număr N de praguri T1, T2, … TN şi crearea unei imagini de etichete v, pe baza imaginii iniţiale u, astfel: eticheta 1, daca T1 ≤ u(i,j) ≤ T2 eticheta 2 , daca T ≤ u(i,j) ≤ T 2 3 v(i, j ) = ... eticheta N, daca TN-1 ≤ u(i,j) ≤ TN
161
(7.1)
Numărul de apariţii
α T1 β T2
γ
T3
Nivelul de gri
Figura 7.7. Segmentarea bazată pe prăguirea histogramei. În exemplul prezentat, punctelor cu valori cuprinse între [0,T1] li se va atribui eticheta α, punctelor cu valori cuprinse între [T1,T2] li se va atribui eticheta β etc. Etichetele alocate sunt, de regulă, numere întregi. Pragurile de segmentare se aleg, în general ca fiind minimele histogramei, deci nivelurile de gri cele mai slab reprezentate în imagine. Se porneşte de la ideea că obiectele identice cu acelaşi nivel de gri, dacă sunt bine reprezentate, au cam aceleaşi valori maxime în histogramă. Prin urmare, pragurile filtrării (Ti, unde i reprezintă numărul pragurilor luate în considerare) vor corespunde minimelor histogramei (unde avem un număr mic de pixeli cu valoarea respectivă), care trebuie detectate. Dacă aceste minime nu sunt bine reliefate, se poate aplica fie o filtrare a imaginii (înlăturarea zgomotului, mărirea contrastului), fie se poate construi o histogramă ponderată, care să ia în considerare doar pixelii care se află pe platouri de intensitate (şi nu în regiuni de tranziţie). Această distincţie poate fi făcută pe baza laplacianului, iar histograma ponderată sumează pentru fiecare pixel cu nivelul de gri, nu L, ci valoarea 1/(1+|L|), unde L este valoarea laplacianului calculat pe o vecinătate a punctului respectiv. Pe histograma ponderată minimele apar mai bine reliefate, deci pragurile de segmentare se pot alege mai uşor. 162
Dacă pentru imagine (şi obiectele conţinute în aceasta) se dispune de informaţie suplimentară (de tipul unei caracterizări statistice a conţinutului imaginii şi a modului de degradare a acesteia) este posibilă o abordare derivată din teoria deciziilor optimale. În acest caz, pragul optim de segmentare pentru o imagine cu două tipuri de obiecte (corpuri şi fundal) depinde de caracteristicile statistice ale zgomotului şi ale obiectelor, după formula: Pob nob + n f 1 − Pob T= − d 2 ⋅ ln (7.2) 2 nob − n f unde nob este nivelul mediu de gri a obiectelor, nf este nivelul mediu de gri al fundalului, Pob este suprafaţa relativă din imagine ocupată de pixeli obiect, iar d2 este varianţa zgomotului aditiv Gaussian de medie nulă, aplicat imaginii. Ca domeniu de aplicaţii s-ar putea aminti medicina, de exemplu pentru numărarea nucleelor de un anumit tip dintr-o imagine, prin numărarea etichetelor care ne interesează.
Figura 7.8. Exemplu din medicină, la numărarea nucleelor. Dezavantajul metodei este că vor fi luate în considerare şi zgomotele care au acelaşi nivel de gri.
163
Pentru segmentarea imaginilor binare se poate utiliza şi histograma L×K
1
i =1
i =0
cumulativă, H i = ∑ hi sau cu varianta sa normată H i = ∑ hi .
(7.3)
Număr de apariţii 1 P
α
T
β
N-1
Nivel de gri
Figura 7.9. Histograma cumulativă. În acest caz, se stabileşte un procent de pixeli (p) sau un prag de nivel de gri (T) cărora li se atribuie o anumită etichetă (α), iar la restul li se atribuie o altă etichetă (β). Metoda se aplică cu succes la imagini cu histograme bimodale. Pentru aceste metode există şi varianta de prăguire adaptivă, adică de a împărţi imaginea în zone în care se calculează histograma locală şi se aplică una din metodele anterioare, obţinându-se un rezultat cu caracteristici locale. Există şi metode de segmentare bazate pe potrivirea sau căutarea de măşti. Acestea se pot implementa de exemplu utilizând transformarea Hit or Miss din morfologia matematică sau cu filtre adaptive la care funcţia de transfer (funcţia de intercorelaţie dintre intrare şi ieşire) se poate transforma în funcţie de autocorelaţie. De exemplu, în cazul unui timbru cu ştampila poştei se poate identifica numai ştampila şi caracteristicile ei:
164
timbru cu ştampila poştei
masca
Figura 7.10. Exemplu în cazul căutării timbrului poştal.
7.2.2. Segmentarea bazată pe creşterea şi fuziunea regiunilor Aplicarea tehnicilor de segmentare pe histogramă este condiţionată de posibilitatea reprezentării diferitelor clase de obiecte din imagine pe intervale de niveluri de gri diferite care nu se suprapun (sau se suprapun parţial pe porţiuni foarte mici). În plus este necesară cunoaşterea numărului de tipuri de obiecte diferite. În fine, se presupune că valorile prag corespunzătoare se pot determina cu o precizie corespunzătoare. Chiar în cazurile în care toate aceste condiţii enunţate sunt îndeplinite, nu se poate garanta condiţia de conexitate a regiunilor obţinute în urma segmentării. Acest lucru este evident, atât timp cât la două obiecte de acelaşi tip, neconexe, li se atribuie prin segmentarea pe histogramă o aceeaşi etichetă, şi formează în imaginea de etichete o regiune neconexă. O metodă care respectă toate condiţiile impuse prin definiţia matematică a segmentării, este creşterea regiunilor. Principiul pe care se bazează creşterea regiunilor este simplu: se aleg în imagine puncte reprezentative pentru fiecare obiect individual şi categorie de obiecte, pe baza cărora are loc un proces de aglomerare a pixelilor vecini acestora, care au aceleaşi proprietăţi (în particular acelaşi nivel de gri). În urma acestui proces de aglomerare sau adăugare de puncte, se obţin zone (sau regiuni) de pixeli cu aceleaşi caracteristici, deci obiecte individuale. Procesul se opreşte în momentul în care fiecare punct 165
al imaginii a fost alocat unei regiuni. Evident, metoda astfel descrisă pe scurt, are două etape esenţiale: •
iniţializarea sau alegerea punctelor de start (puncte iniţiale), numite germeni sau seminţe (engl. seed);
•
creşterea propriu-zisă a regiunilor. Numărul final de regiuni rezultate este egal cu numărul de germeni aleşi iniţial pentru creştere, deci alegerea, respectiv granulaţia (densitatea) acestor puncte este foarte importantă. În principiu, este de dorit ca fiecare obiect individual aflat în imagine să fie marcat cu câte un germene. Dacă în interiorul unui aceluiaşi obiect se găsesc mai mulţi germeni, pentru fiecare dintre ei va fi crescută o regiune. Aceasta face ca obiectul iniţial să fie împărţit artificial prin segmentare în mai multe regiuni. Parţial, acest neajuns se poate corecta printr-o etapă ce urmează creşterii regiunilor, şi anume fuziunea regiunilor adiacente ce au proprietăţi asemănătoare. Dacă în interiorul unui obiect nu este ales nici un germene, obiectul respectiv va fi înglobat de regiunile ce cresc pornind de la germeni din vecinătatea spaţială. Astfel, respectivul obiect nu apare ca o regiune distinctă şi este pierdut, rezultând o eroare gravă de segmentare. Pornind de la germenii aleşi, regiunile sunt obţinute printr-un proces de creştere aproape simultană, început de la aceştia, până când toţi pixelii imaginii sunt repartizaţi unei regiuni (figura 7.11).
Figura 7.11. Principiul creşterii regiunilor.
166
Cvasi-simultaneitatea creşterii poate fi realizată cu un algoritm serial, prin alocarea pixelilor ce sunt adiacenţi (vecini) zonelor deja segmentate. Această alocare trebuie să ţină seama de criteriul ca regiunile crescute să fie uniforme, adică nivelul de gri al pixelului ce se adaugă nu trebuie să difere cu mai mult de un prag prestabilit faţă de nivelul de gri al germenului regiunii la care se alocă. În acelaşi timp, la o singură trecere, numărul de puncte ce se adaugă unei regiuni nu poate depăşi un număr prestabilit (condiţia încearcă să asigure creşterea relativ uniformă şi izotropă a tuturor regiunilor). Dacă adăugarea de noi pixeli se blochează (criteriul de uniformitate nu mai este respectat), diferenţa maximă admisă pentru nivelul de gri poate fi crescută în etape, până la epuizarea pixelilor imaginii. Avantajele pe care le are o asemenea tehnică de creştere a regiunilor sunt acelea că nu mai este necesară nici o informaţie privind conţinutul imaginii, regiunile crescute sunt conexe şi nu există puncte neetichetate (nealocate vreunei regiuni), iar poziţia frontierelor percepute subiectiv în imagine se conservă. Fuziunea regiunilor deja determinate în etapa de creştere, are drept scop reducerea numărului de regiuni în care a fost împărţită iniţial imaginea, pentru a evita fenomenul de supra-segmentare. Regiunile candidate la fuzionare trebuie să fie învecinate, iar decizia de fuzionare se ia în funcţie de pixelii aflaţi pe frontiera comună. Astfel, punctele slabe (în număr de ns) sunt punctele pentru care diferenţa nivelurilor de gri între vecinii din regiunile adiacente este foarte mică (mai mică decât un anumit prag fixat). Punctele tari (în număr de nt) sunt acele puncte pentru care diferenţa nivelurilor de gri între vecinii din regiunile adiacente este foarte mare (mai mare ca un anumit prag fixat). Cu aceste definiţii, se poate afirma că regiunile Ri şi Rj vor fuziona dacă: •
numărul de puncte slabe (ns) raportat la perimetrul minim (Pm) este n mare: s > θ1 , (7.4) Pm unde Pm=min(Perimetrul(Ri), Perimetrul(Rj)); 167
•
•
numărul de puncte slabe de pe frontiera comună e mare: n s > θ 2 , (7.5) P unde P este numărul de puncte aflate pe frontiera comună a regiunilor Ri şi Rj; n numărul de puncte tari de pe frontiera comună este mic: t < θ 3 , (7.6) P θ1 = 0,5 Valori tipice ale pragurilor θ1 ,θ 2 ,θ 3 sunt: θ 2 = 0,75 θ = 0,2 3
•
(7.7)
distanţa de similaritate dintre regiuni este mare. Distanţa de similaritate este o măsură a asemănării. Aceasta poate fi o distanţa euclidiană ponderată sau un produs vectorial sau scalar între vectorii corespunzători.
7.3. Segmentarea orientată pe contururi După cum s-a arătat, segmentarea reprezintă o categorie de tehnici de prelucrare care permite extragerea din imagini a unor elemente de interes în scopul analizei lor ulterioare şi eventual a clasificării lor, prin două metode principale: detecţia muchiilor (contururilor) şi detecţia regiunilor omogene. Prin urmare, pentru analiza imaginilor (îndeosebi la segmentare) o etapă esenţială constă în detecţia muchiilor şi liniilor (frontierelor) care reprezintă grupuri de pixeli aflaţi în zona de tranziţie (de variaţie bruscă a nivelului de gri) dintre două regiuni relativ uniforme ale imaginii iniţiale. Muchiile sunt graniţa de separaţie fie între obiecte şi fundal, fie între două zone omogene ale aceluiaşi obiect. O asemenea muchie poate fi formată din puncte şi segmente de dreaptă.
168
7.3.1. Operatori de tip gradient Din punct de vedere practic, un punct de contur dintr-o imagine reprezintă un pixel (sau un grup mic şi omogen de pixeli conecşi) având o valoare net diferită faţă de cea a vecinilor săi. Detecţia acestora se poate face prin comparaţie (diferenţă) între nivelul de gri propriu şi media nivelurilor de gri dintr-o vecinătate a sa. Măsura diferenţei de intensitate (de nivel de gri) poate fi calculată printr-o filtrare liniară cu o mască de forma: − 1 − 1 − 1 − 1 [8] − 1 − 1 − 1 − 1
(7.8)
Deoarece pixelii aparţinători muchiilor se caracterizează prin faptul că se află la graniţa dintre două regiuni omogene între care există diferenţe mari de niveluri de gri, detecţia lor se poate face cu operatori de tip gradient (operatori de derivare locală). Pentru o imagine iniţială u(l,k), putem scrie: •
gradientul pe direcţia orizontală: g x (l , k ) = u (l + 1, k ) − u (l − 1, k ) (7.10)
•
gradientul pe direcţia verticală: g y (l , k ) = u (l , k + 1) − u (l , k − 1) (7.11)
•
gradientul după o direcţie oarecare r (care face unghiul θ cu orizontala): ∂u ∂u ∂x ∂u ∂y = ⋅ + ⋅ = g x ⋅ cos θ + g y ⋅ sin θ (7.12) ∂θ ∂x ∂r ∂y ∂r
Direcţia în care se află muchia va fi dată de direcţia după care gradientul are valoarea maximă: ∂ ∂u =0 ∂θ ∂r
(7.13)
Rezultă vectorul gradient g(l,k) cu: •
modulul (amplitudinea): g (l , k ) = g x2 (l , k ) + g 2y (l , k )
169
(7.14)
•
g y (l , k ) g x (l , k )
α (l , k ) = tg
direcţia:
(7.15)
De multe ori, pentru simplificarea calculelor, amplitudinea (modulul) gradientului se defineşte ca: g (l , k ) = g x (l , k ) + g y (l , k )
(7.16)
Implementarea operatorilor de derivare (tip gradient) se face prin convoluţia imaginilor cu măşti (ferestre), care sunt deci, filtrări liniare. Operatorii de tip gradient sunt reprezentaţi de o pereche de filtre (H1,H2) care măsoară gradientul imaginii după două direcţii ortogonale. Spre exemplu, gradienţii după direcţia orizontală respectiv verticală, definiţi mai sus (gx, gy), se obţin prin convoluţia cu măştile:
[1 [0]
1 [0] − 1
− 1]
orizontal
(7.17)
vertical
Un pixel este declarat ca “punct de frontieră” dacă g(l,k) depăşeşte un prag dat, t. Pixelilor care depăşesc acest prag li se alocă în imaginea finală eticheta C (contur), iar celorlalţi, eticheta F (fundal). Un bloc care implementează un astfel de extractor de contur poate fi reprezentat ca în figura 7.12. h1(-l,-k)
gx
g = g12 + g 22
Amplitudine g(l,k) prag
u(l,k) Direcţie αg(l,k)
h2(-l,-k)
gy
t O
g2 g1
α g = arctg
Figura 7.12. Principiul extragerii contururilor. 170
Harta muchiilor
Prin binarizarea imaginii “amplitudine g(l,k)” cu pragul “t” se obţine harta de muchii. Harta muchiilor oferă informaţiile necesare pentru trasarea contururilor din imagine. În general, pragul “t” se alege folosindu-se histograma cumulativă pentru g(l,k), 5-10% din pixelii imaginii g(l,k) fiind declaraţi muchii. În continuare sunt prezentaţi operatori de derivare clasici, care determină diferenţe orizontale şi verticale şi realizează însumări locale, reducându-se efectul zgomotului (elementul încadrat indică originea): − 1 0 1 Prewitt: H 1 = − 1 [0] 1 − 1 0 1
− 1 − 1 − 1 H 2 = 0 [0] 0 1 1 1
(7.18)
− 1 0 1 H 1 = − 2 [0] 2 − 1 0 1
− 1 − 2 − 1 H 2 = 0 [0] 0 1 2 1
(7.19)
3 3 3 H 2 = 3 [0] 3 − 5 − 5 − 5
(7.20)
Sobel:
− 5 3 3 Kirsch: H 1 = − 5 [0] 3 − 5 3 3 0 −1 Izotrop: H 1 = − 2 [0] − 1 0
− 1 − 2 1 [0] 2 H 2 = 0 1 1 2
− 1 0 1
(7.21)
Aceşti operatori au proprietatea de a avea un efect nul în cazul regiunilor uniforme. Operatorii definiţi anterior sunt performanţi în cazul unei tranziţii bruşte a nivelurilor de gri din imagine. Când aceste tranziţii sunt mai lente, este indicat să se folosească derivatele de ordin doi, combinate într-un operator laplacian:
171
∇ 2 f ( x, y ) =
∂ 2 f ( x, y ) ∂x 2
+
∂ 2 f ( x, y )
(7.22)
∂y 2
Operatorul de tip laplacian poate fi implementat în discret prin convoluţia imaginii de prelucrat cu una dintre măştile: 0 −1 0 − 1 [4] − 1 0 − 1 0
− 1 − 1 − 1 − 1 [8] − 1 − 1 − 1 − 1
1 −2 1 − 2 [4] − 2 1 − 2 1
(7.23)
Datorită derivatei de ordinul doi, operatorul laplacian este mai sensibil la zgomot decât cei definiţi anterior. Amplitudinea binarizată pentru ∇ 2 f produce muchii duble, ceea ce duce la apariţia în harta de muchii, a contururilor îngroşate. Acest dezavantaj se elimină dacă se consideră ca fiind punct de contur nu cel corespunzător maximelor sau minimelor laplacianului, ci trecerilor sale prin zero (figura 7.13). f'(x)
f(x)
O
O
x
x
f''(x) O
x
Figura 7.13. Graficele corespunzătoare laplaceanului.
172
7.3.2. Operatori de tip compas Operatorii prezentaţi, cât şi topologia planului discret, favorizează numai câteva direcţii: axele şi cele două bisectoare, fiecare cu câte două sensuri. Pornind de la această observaţie, se pot construi măşti pentru determinarea gradientului pe fiecare dintre aceste 3 direcţii. Operatorii rezultaţi se numesc operatori compas. Aceştia se pot obţine din oricare din operatorii prezentaţi anterior. De exemplu, operatorii compas derivaţi din operatorii Prewitt, corespunzători celor 8 orientări sunt: 1 1 1 0 [0] 0 − 1 − 1 − 1
N
0 1 1 1 [0] − 1 0 − 1 − 1
1 1 1
NV
− 1 − 1 − 1 0 [0] 0 1 1 1
− 1 − 1 0 − 1 [0] 1 0 1 1
S
SE
− 1 − 1 − 1
0
[0] 0
V − 1 − 1 − 1
0 − 1 − 1 1 [0] − 1 1 1 0
SV 0
[0] 0
E
1 1 1
1 1 0 − 1 [0] 1 − 1 − 1 0
NE
Măştile corespunzătoare diferitelor orientări se obţin prin rotirea unei măşti de bază în jurul originii. Pentru a obţine rezoluţii unghiulare mai bune (mai multe direcţii) se pot utiliza operatori de tip compas de dimensiuni mai mari de 3×3. Filtrele prezentate se pretează la obţinerea prin permutare a unui set de patru, opt sau mai multe măşti, sensibile la orientarea muchiei, în funcţie de dimensiunea operatorului. Pentru determinarea punctelor de contur se consideră direcţia după care gradientul este maxim. Operatorul compas este evident neliniar, datorită operaţiei de maxim.
173
Hk
gk
max{gk}
Figura 7.14. Determinarea punctelor de contur ca puncte de gradient maxim. Se poate utiliza şi varianta adaptivă a acestui operator (compas adaptiv), utilizând aceleaşi măşti ca şi în cazul anterior, dar pentru fiecare mască se va verifica şi situarea punctului curent pe mijlocul tranziţiei. În acest scop, pentru nivelurile de gri din vecinătatea punctului curent se pot utiliza notaţiile: 1 1 1 P1 P2 P3 P4 P7
[P5 ] P8
P6 P9
−2 1
[− 2] 1
−2 1
(7.24)
Cu aceste notaţii, se verifică condiţiile: (P4 > P1 & P4 > P7 ) & (P5 > P2 & P5 > P8 ) & (P6 > P3 & P6 > P9 ) (7.25) adică se verifică faptul că valorile centrale sunt mai mari decât cele laterale, deci dacă prin aceste puncte trece o linie, care determină un profil de forma:
Gradientul după direcţia respectivă, care rezultă prin aplicarea unei astfel de măşti, se ia în considerare doar dacă sunt îndeplinite condiţiile (7.25). Dezavantajul operatorilor de tip compas este că zgomotele suprapuse, de exemplu peste o linie, au ca rezultat o tranziţie, respectiv un contur fals. Acest dezavantaj se poate elimina printr-o prefiltrare a imaginii. Astfel rezultă operatorii compas cu netezire. Netezirea imaginii se face adaptiv, după o anumită direcţie, rezultând o netezire adaptivă direcţional. În acest scop se alege un număr 174
de măşti de netezire direcţionale (S1, S2,..., Sk). Pentru fiecare punct al imaginii calculăm rezultatul netezirii (filtrării) cu aceste măşti. Rezultă valorile netezirii, v1, v2,..., vk, care sunt nişte valori medii ponderate de coeficienţii măştilor. Pentru fiecare din aceste măşti se calculează şi o dispersie locală, σ1, σ2,..., σk, care reprezintă dispersia valorilor de gri ponderate cu coeficienţii măştii, faţă de media ponderată vk. Se va folosi masca care are ca rezultat dispersia minimă σ j = min (σ i ) . Prin urmare, valoarea pixelului curent se va înlocui cu i =1,...k
valoarea vj, care a avut ca rezultat σ j = min (σ i ) . i =1,...k
Masca poate fi aleasă şi după alt criteriu cum ar fi, de exemplu, obţinerea unei valori minime pentru diferenţa între valoarea punctului curent (în care a fost centrată masca) şi valoarea rezultată în urma filtrării,
respectiv medierii: min ( u ( x, y ) − vi ) . i =1,...k
O altă variantă de operator extractor de contururi este gradientul morfologic, definit ca diferenţă între imaginea dilatată şi imaginea erodată, cu un acelaşi element structurant, aşa cum s-a arătat în capitolul referitor la morfologie matematică. Dacă elementul structurant este V8flat (plat), această operaţie poate fi asimilată cu un operator compas. Toate aceste tehnici se dovedesc a fi sensibile la zgomot şi prin urmare, pentru imagini afectate de zgomot sunt necesare soluţii diferite. O asemenea prelucrare este diferenţa de gaussieni (DOG) care reprezintă diferenţa între două imagini netezite cu filtre gaussiene de dimensiuni diferite.
175
8. Compresia imaginilor O clasificare a metodelor de compresie a imaginilor este prezentată în figura de mai jos: codarea conturului nivel înalt nivel bloc
skeleton quad-tree
morfologic generalizat
WBS (White Block Skipping)
binare
RLE (Run Length Encoding) - cu cuvânt fix
nivel bit (pixel)
Compresia imaginilor
- cu cuvânt trunchiat - cu cuvânt modificat Ziv-Lempel entropică (Huffman)
cu niveluri de gri
pe plane (ca la imagini binare) DCT cu transformate (cu pierderi) KL (Karhunen Loeve) cu predicţie, bazate pe DPCM wavelet (Differential Pulse Code Modulation) cu fractali
etc.
VQ (Vector Quantization)
Există şi metode combinate: RLE+Hufmann, skeleton+RLE (sau WBS), DCE+Hufmann (utilizată în standardul JPEG Joint Picture Expert Group) etc.
176
8.1. Compresia imaginilor binare 8.1.1. Codarea de nivel înalt Pentru compresia de nivel înalt a imaginilor binare, pe lângă compresia cu skeleton morfologic, se mai poate utiliza codarea contururilor. Descriptorii de contur sunt coduri care înmagazinează şi compactează într-un volum mic de memorie informaţiile esenţiale cu privire la conturul unui obiect. Determinarea acestor descriptori trebuie precedată de detecţia contururilor şi eventual de netezirea acestora.
8.1.1.1. Aproximări poligonale
Aproximarea poligonală a unui contur este un descriptor care descinde direct din metoda de detecţie a conturului, cu acelaşi nume. Procedura de determinare a aproximării poligonale a unui contur constă în următoarele etape (figura 8.1): •
determinarea diagonalei principale (AB) a conturului;
•
marcarea punctelor de pe contur (C şi D) având distanţa maximă faţă de diagonala principală;
•
rezultă prima aproximare: poligonul ABCD;
•
marcarea punctelor maxim depărtate de laturile patrulaterului (EFGH);
•
rezultă a doua aproximare: poligonul ABCDEFGH;
•
se repetă ultimele etape până când distanţa maximă este mai mică decât o valoare prestabilită.
177
A
A
A E
C
D
H
C
D
F B
B
A
G B
E
H
C
D
F
G B
Figura 8.1. Procedura de determinare a aproximării poligonale a unui contur. Codul rezultat conţine coordonatele colţurilor poligonului, ceea ce practic defineşte laturile şi unghiurile acestuia şi implicit dă o informaţie suficient de corectă asupra conturului obiectului studiat. Aplicaţiile acestei metode sunt limitate, deoarece modelul obţinut nu este invariant la transformările geometrice şi deci, o simplă scalare sau rotire duce la schimbarea codului. Pentru a înlătura într-o oarecare măsură acest dezavantaj, pentru codarea poligonului se poate utiliza codul Freeman.
8.1.1.2. Codul Freeman
Codul Freeman se mai numeşte şi “cod lanţ” şi este folosit pentru codarea contururilor. Principiul acestuia constă în codificarea (numerotarea) vectorilor de “direcţie” dintre pixelii succesivi de pe contur prin numerotarea într-o anumită ordine a direcţiilor posibile de deplasare de-a lungul unui contur. Această metodă de compresie se bazează pe ideea că pentru un obiect „plin” este suficientă informaţia de contur, pentru a putea reface obiectul, iar conturul este codificat reţinând poziţia relativă a pixelilor de pe contur. Codul Freeman de bază foloseşte 8 direcţii de mişcare (figura 8.2), vectorii fiind codificaţi prin cuvinte de 3 biţi:
178
2
C
B
5
1
3
4
3
D 0
4
A 2
6
H
E 7
5
1
7 F
6
0
G
Figura 8.2. Codul Freeman. Tipic, codul Freeman conţine adresa pixelului de start (în exemplul prezentat, A), urmată de un şir de cuvinte cod, care codifică direcţia vectorului spre pixelul următor şi adresa acestuia (A34567012). Codul obţinut este invariant la translaţie şi este dependent de alegerea punctului de start pentru parcurgerea conturului. Din acest motiv, ca alternative la acest cod Freeman primar, se pot folosi şi alte variante ale acestuia, cum ar fi: •
codul Freeman diferenţial, obţinut prin scrierea diferenţei modulo-8 între cifrele succesive. Practic se numără întotdeauna în acelaşi sens numărul de direcţii care separă două orientări consecutive ale conturului, adică două cifre din codul primar.
•
“numărul de formă” este cifra minimă obţinută din permutarea circulară a codului diferenţial. Acest descriptor de contur este constant, independent de poziţia punctului de start.
179
Exemplu: •
codul Freeman primar = 60026442;
•
codul Freeman diferenţial = 20246064 6-4=2 4-4=0 4-2=2
=se începe cu diferenţa ultimelor două cifre;
•
numărul de formă: 02460642 (cel mai mic). Această metodă de compresie se utilizează la comanda plotterelor, pentru realizarea de hărţi, desene sau amprente digitale. Compresia cu coduri Freeman permite şi o serie de prelucrări pe contururi, fără a fi nevoie să fie memorată sau analizată întreaga imagine. Astfel, se pot calcula sau realiza: 1. perimetrul unui obiect, pe baza informaţiei referitoare la conturul său: perimetrul = n par + 2 ⋅ nimpar ,
(8.1)
unde npar este numărul legăturilor pare, iar nimpar este numărul legăturilor impare din codul Freeman. 2. ariile închise de contururi. 3. netezirea frontierelor obiectelor. Codul Freeman este folosit în mod frecvent în caracterizarea formelor, putând fi folosit şi pentru determinarea altor parametri de formă. Parametrii de formă reprezintă scalari sau funcţii asociate unei forme, pe care o caracterizează. Astfel, formele asemănătoare sunt caracterizate de parametri de formă de valori apropiate. Parametrii de formă compun un fel de “fişă de identitate” a formei respective, pe baza căreia această formă poate fi recunoscută în mod unic. În mod ideal aceşti parametri trebuie să fie invarianţi la translaţie, rotaţie şi scalare. Printre cei mai cunoscuţi parametri (descriptori) de formă sunt: 180
•
descriptori primari (geometrici): lungimi, perimetre, arii etc.;
•
descriptori de contur: Freeman, Fourier, aproximări poligonale:
•
momente statistice (invariante): Hu, Zernike, afine etc. Parametrii geometrici se bazează pe măsurarea unor atribute geometrice simple (sau combinaţii ale acestora) cum ar fi: •
perimetrul (P): P = ∫ x 2 (t ) + y 2 (t )dt
(8.2)
unde t este parametrul de contur, dar nu în mod necesar lungimea acestuia. •
aria (A): A = ∫∫ dxdy = ∫ y (t ) ∂R
R
dx(t ) dy (t ) dt − ∫ y (t ) dt dt dt ∂R
(8.3)
unde R şi ∂R reprezintă regiunea obiectului şi respectiv, conturul acesteia. •
•
excentricitatea formei (sau circularitatea sa, adică măsura în care forma se deosebeşte de un disc) se defineşte ca raport între raza cercului circumscris (R) şi raza cercului înscris (r) în forma studiată: R c= (8.4) r raportul de compactizare (sau rotunjimea formei) reprezintă raportul dintre pătratul perimetrului şi suprafaţa formei: K=
P2
(8.5)
4πA 2
8.1.1.3. Descriptori Fourier
Descriptorii Fourier reprezintă o metodă utilă pentru reprezentarea şi descrierea conturilor. Conturul poate fi descris cu o pereche de funcţii x(n) şi y(n) (periodice, cu perioada N, egală cu numărul punctelor de pe 181
contur sau cu un submultiplu al acestuia), în continuare putând fi utilizate toate tehnicile de reprezentare unidimensională pentru semnale. Astfel, pentru orice contur eşantionat, care conţine N puncte, se poate scrie: u (n) = x(n) + j ⋅ y (n) , unde n = 0,1,..., N − 1 (8.6) În cazul unui contur închis, semnalul u(n) va fi un semnal periodic cu perioada N. x
α
a(k)=x(k)+j⋅y(k)
y
Figura 8.3. Descriptori Fourier. Descriptorii Fourier sunt chiar coeficienţii a(k) ai transformatei Fourier ai funcţiei u(n), obţinuţi prin aplicarea transformatei Fourier unidimensionale. Astfel, a(0) reprezintă centrul de greutate al curbei, iar ceilalţi coeficienţi a(k) conţin informaţiile privind variaţiile locale ale conturului descris: N −1 2π ⋅ k ⋅ n a (k ) = ∑ u (n) ⋅ exp − j , unde 0 ≤ k ≤ N − 1 N n =0
(8.7)
Descriptorii Fourier pot fi utilizaţi în recunoaşterea de forme de acelaşi tip, chiar dacă acestea au dimensiuni şi orientări diferite. Se poate arăta că vectorul descriptor Fourier al conturului unei forme: V F = {a (0) , a (1) , a (2) ,..., a ( N − 1) }
(8.8)
este invariant la scalare şi translaţie. În plus, faza arg(a (k ) ) , unde k=0,1,2,…,N-1, este invariantă la rotaţie. 182
Deci, anumite transformări geometrice ale unui contur sau forme se reflectă, în cazul descriptorilor Fourier corespondenţi, în transformări (operaţii) simple ale acestora. De exemplu, dacă un contur este translatat cu: u 0 = x0 + j ⋅ v0 , ⇒ u ' ( n) = u ( n) + u 0 (8.9) În acest caz, noii descriptori Fourier rămân identici, cu excepţia celui pentru k=0, a ' (k ) = a (k ) + u 0 ⋅ δ (k ) . Efectul de scalare (mărire sau micşorare a conturului) are ca efect o scalare a coeficienţilor a(k), deci dacă: u ' ( n) = α ⋅ u ( n) ⇒ a' (k ) = α ⋅ a(k ) (8.10) Modificarea punctului de referinţă (de start) a conturului duce la o modulare a coeficienţilor a(k). Deci, dacă: u ' ( n) = u ( n) ⋅ e
− jΦ 0
⇒ a' (k ) = a(k ) ⋅ e
−j
2π ⋅k ⋅n0 N
(8.11)
Rotaţia conturului cu un unghi φ0 produce un defazaj suplimentar constant φ0 al descriptorilor: u ' ( n ) = u ( n ) ⋅ e − jΦ 0
⇒ a ' ( k ) = a ( k ) ⋅ e − jΦ 0
(8.12)
8.1.2. Codarea la nivel de bloc 8.1.2.1. Metoda arborelui cuaternar (Quad-tree)
Compresia bazată pe metoda arborelui cuaternar (quad-tree) se implementează prin segmentarea succesivă a imaginii în sferturi succesive, până se obţin doar zone omogene, ca şi în cazul segmentării prin metoda quad-tree. Codarea se face apoi prin atribuirea pentru fiecare „frunză” a arborelui a unei valori corespunzătoare valorii din zona (pătratul) adresată de frunză.
183
1 1 0 0 1
1 2
0
0
1
0
0
1
2
3
0
3
0
1
3
2
convenţie
Figura 8.4. Compresia prin metoda arborelui cuaternar. Conform acestei metode, pentru a reface imaginea din codul acesteia, pentru fiecare „frunză” trebuie ştiut (memorat):
•
valoarea: este conţinută în fiecare frunză;
•
dimensiunea: poate fi dedusă din nivelul pe care se află frunza (de exemplu, frunzele, deci pătratele de pe primul nivel au latura egală cu ½ din latura imaginii);
•
poziţia: rezultă din calea pe care se ajunge la frunză. Prin urmare codul obţinut prin metoda quad-tree trebuie să conţină toate informaţiile de mai sus. Din aceste informaţii se poate elimina informaţia referitoare la dimensiunea frunzei, deoarece aceasta rezultă în mod implicit, deoarece din lungimea codului rezultă şi nivelul pe care se află frunza şi deci dimensiunea acesteia. De exemplu, codul corespunzător valorii pătratului negru din partea dreapta-jos a imaginii de mai sus este 111 (cel al poziţiei sale este 22), iar cel al valorii pătratului alb din stângajos este 10 (cel al poziţiei sale este 3).
8.1.2.2. Metoda WBS (White Block Skipping) Această metodă de codare la nivel de bloc se bazează pe omiterea zonelor albe din imagine. Principiul implementării acestor metode constă în: 184
•
împărţirea şirului de biţi ce defineşte imaginea în blocuri de câte m biţi;
•
un bloc nul (toţi biţii sunt 0-logic) se înlocuieşte cu un singur bit de 0;
•
blocurile nenule (în care cel puţin un bit este diferit de 0-logic) sunt înlocuite cu un bit de 1-logic, urmat blocul respectiv. De exemplu, pentru m=3 şi pentru secvenţa următoare: 000 110 000 101 000 011 codul corespunzător este: 0 1110 0 1101 0 1011 După cum se observă, acest tip de compresie este eficient doar în cazul imaginilor cu mult alb, în caz contrar putându-se obţine chiar un efect contrar compresiei. Se poate arăta că: 1 R = 1 + − p nul , unde m
(8.13)
R este inversul ratei de compresie (C), pnul este probabilitatea ca blocul de m biţi să fie nul. 1 numarul de biti din imaginea initiala (8.14) C= = R numarul de biti din imaginea finala
8.1.3. Codarea la nivel de bit 8.1.3.1. Codarea RLE (Run Length Encoding) Principiul codării RLE (Run Length Encoding) constă în reţinerea numărului de biţi succesivi care au aceeaşi valoare în secvenţa de biţi ce descrie imaginea. De exemplu, pentru secvenţa iniţială: 000111010111110000111111 care are 24 biţi prin codarea conform principiului RLE se obţine secvenţa codată: 033111546 185
Adică primul bit este 0-logic, care se repetă de 3 ori, după care se schimbă simbolul (în 1-logic, dar acesta nu se mai memorează, deoarece în imaginile binare tranziţia din 0-logic se face doar în 1-logic), acesta se repetă de 3 ori, după care se schimbă din nou simbolul, acesta apare o singură dată, după care avem 2 schimbări succesive de simbol, acestea apărând câte o singură dată, după care avem simbolul 1-logic care se repetă de 5 ori, se schimbă simbolul (în 0-logic) care apare de 4 ori, iar apoi ultimul simbol apare de 6 ori. Se observă că în secvenţa codată cel mai mare număr care apare este 6, deci pentru codarea binară a acestei secvenţe este nevoie de 8 coduri (apar 8 numere în secvenţa codată 33111546) de câte 3 biţi (pentru codarea celui mai mare număr este nevoie de 3 biţi) plus bitul ce reprezintă valoarea primului bit din secvenţa iniţială. Prin urmare, se poate deduce că metoda este eficientă în codarea imaginilor care conţin zone uniforme de dimensiuni mari. Metoda RLE este utilizată la faxuri.
8.1.3.2. Codarea entropică (Huffman) Deoarece codarea entropică se studiază în detaliu în teoria informaţiei, nu se va insista pe detalii ale metodei, ci doar pe descrierea principiului metodei, prin intermediul unui exemplu. Pentru implementarea codării entropice se porneşte de la un număr de simboluri disponibile, care au anumite probabilităţi de apariţie. Principiul codării entropice (Huffman) constă în a coda cu un număr mic de biţi simbolurile cu probabilitate de apariţie mare, iar simbolurile sau secvenţele rare (cu probabilitate de apariţie mică) să fie codate cu un număr mai mare de biţi. Pentru exemplificarea principiului metodei, se presupun şapte simboluri cu probabilităţile de apariţie din tabelul 8.1.
186
Tabelul 8.1. Simbolul PA
S1
6
S2 S3
3 3
S4 S5 S6 S7
2 2 1 1
Iniţial CA
Etapa 1 PA CA
Etapa 2 PA CA
6
6
4 (G2) 3 3 3 3 010 2 (G1) 2 011 2 100 2 101 0110
Etapa 3 PA CA
Etapa 4 PA CA
7 (G4) 6 6 00 01 5 (G3) 5 4 10 3 11
Etapa 5 PA CA 11 (G5) 0 7 1
Etapa 6 CA
00 11 010 100 101 0110 0111
PA = probabilitatea de apariţie CA = codul alocat În faza iniţială, se ordonează simbolurile în ordinea descrescătoare a probabilităţilor, după care se grupează simbolurile două câte două. Simbolurile cele mai puţin probabile sunt grupate şi alcătuiesc o grupă restrânsă de probabilităţi (G1, în tabelul 8.1). Procesul se repetă (G2, ..., G5) până când avem o grupă restrânsă cu doar două elemente (G5 şi G4), cărora le alocăm câte un simbol diferit (0,1), după care se refac în sens invers grupele, respectiv simbolurile iniţiale, prin alocarea unor biţi corespunzători care se adaugă celor din grupele superioare (anterioare), din care provin grupele sau simbolurile curente. Prin codarea entropică se obţine codul cu lungimea medie minimă, în raport cu alte coduri.
187
8.2. Compresia imaginilor cu niveluri de gri 8.2.1. Codarea pe plane În cazul codării pe plane a imaginilor cu niveluri de gri se pot aplica metode similare metodelor de codare a imaginilor binare, pe fiecare plan ce codează informaţia din imagine.
…
…
b7
b1 b0
Figura 8.5. Principiul codării pe plane. În acest caz trebuie prevăzute metode de detecţie a erorilor la nivelurile înalte pentru a nu avea salturi bruşte de valori deoarece, de exemplu, o eroare a bitului cel mai semnificativ MSB (Most Significant Bit) se poate traduce într-o eroare de ½ din valoarea iniţială. Pentru o imagine codată pe 8 biţi, aceasta înseamnă un salt de la 255 la 127 (11111111 -> 01111111). În cazul imaginilor color se poate reduce numărul de culori, înlocuind componentele RGB cu coeficienţii cromatici r,g,b: R G B r= , g= , b= ⇒ r + g + b = 1 (8.15) R+G+B R+G+ B R+G+B Prin, această înlocuire se obţine un spaţiu bidimensional, deoarece componentele r,g,b nu sunt independente (dacă cunoaştem două componente, cea de-a treia rezultă dintr-o combinaţie liniară a primelor 188
două). Spaţiul bidimensional obţinut determină un spaţiu fizic reprezentabil în planul celor două componente cunoscute:
r
b
Figura 8.6. Spaţiul bidimensional r-b. Pentru a se obţine un spaţiu fizic reprezentabil extins, se poate utiliza spaţiul cromatic modificat sau o reprezentare de tip JNC (Just Noticeable Color) prin intermediul elipselor lui McAdam, care au proprietatea că nu se poate deosebi culoarea din centrul elipselor de restul culorilor din interiorul elipsei:
…
r
…
b Figura 8.7. Elipsele lui McAdam.
189
8.2.2. Metode predictive de compresie Transmisia digitală clasică a imaginilor constă în transmisia secvenţială a eşantioanelor obţinute prin baleierea clasică, linie după linie. Prin urmare, avem de transmis o secvenţă u(k) unidimensională, în care fiecare eşantion u(k) este un eşantion cuantizat pe n biţi. Transmisia a n biţi este însă redundantă: eşantioanele vecine sunt puternic corelate, nefiind nici pe departe variabile aleatoare independente. Dacă se transmite în loc de eşantioane diferenţele dintre ele, se obţine o decorelare importantă, iar diferenţele pot fi transmise cu un număr mai mic de biţi fiindcă nu sunt variabile aleatoare cu repartiţie uniformă pe domeniul de valori [0, M-1], ci cu un maxim pronunţat în jurul lui 0 şi cu lungimea de cod semnificativ mai mică decât n. De fapt, folosind acest procedeu nu înseamnă decât să se facă ceea ce se numeşte predicţie de ordinul 0: se estimează: uˆ (k ) = u (k − 1) (8.16) şi se transmite eroarea: ⇒ e(k ) = u (k ) − uˆ (k ) = u (k ) − u (k − 1)
(8.17)
Evident procedeul se poate generaliza: dacă uˆ (k ) = ψ (u (k − 1),..., u (k − l ) ) ≅ u (k )
(8.18)
În acest exemplu, e este eroarea de predicţie, iar ψ este predictorul sau regula de predicţie de ordinul l-1, ψ : R l → R (dacă secvenţa este reală, u(k) ∈ R ), care se deduce pe baza unor considerente statistice. Un caz particular este cazul predicţiei liniare, în care regula de predicţie este o funcţie liniară: l
uˆ (k ) = ∑ ai ⋅ u (k − i )
(8.19)
i =1
Coeficienţii ai se determină pe baza unor considerente statistice.
190
Transmiterea eşantioanelor cuantizate şi folosind un cod (de exemplu 8 biţi pentru 256 de niveluri, folosind convenţia de scriere binară) se numeşte modulaţie de impulsuri în cod PCM (Pulse Code Modulation, în engleză). Transmiterea diferenţelor între eşantioane se numeşte DPCM (Differential Pulse Code Modulation). Codarea diferenţei se poate face în diverse moduri. Cel mai simplu caz din punct de vedere al biţilor de transmis este ca eroarea e(n) să fie codată pe un singur bit (bitul de semn). În acest caz, se fixează o cuantă q şi la transmisie:
•
dacă e(k)>0 se transmite 1-logic;
•
dacă e(k)<0 se transmite 0-logic. La recepţie:
•
dacă s-a recepţionat un 1-logic se adaugă cuanta q la eşantionul anterior;
•
dacă s-a recepţionat un 0-logic se scade cuanta q la eşantionul anterior. Această codare DPCM particulară se numeşte modulaţie delta. În figura 8.8 se observă care este dezavantajul esenţial al modulaţiei delta: eroarea de neurmărire, în momentele când semnalul se modifică rapid. u(k)
uˆ (k )
eroare de neurmărire
zgomot granulat
q
t
Figura 8.8. Principiul modulaţiei delta.
191
În zonele cu pantă mare, deci cu variaţii bruşte, ale u(k) apar distorsiuni (erori) de neurmărire, iar în zonele cu pantă mică, deci cu variaţii lente, ale u(k) apare un zgomot granulat. Modulaţia delta a fost îmbunătăţită mult în raport cu această variantă iniţială prin procedee de modificare adaptivă a pantei de creştere în funcţie de secvenţa de cifre binare transmise. Eliminarea distorsiunii de urmărire se face prin modificarea adaptivă a cuantei de la un pas la altul, adică (figura 8.9):
•
creşterea cuantei pentru două creşteri succesive ale semnalului;
•
scăderea cuantei pentru două scăderi consecutive;
•
păstrarea valorii cuantei de la un pas la altul, pentru o creştere şi o descreştere succesivă (sau o descreştere şi o creştere succesive). u(k)
uˆ (k )
2q q t
Figura 8.9. Principiul modulaţiei delta adaptive. Pentru implementarea compresiei cu predicţie, uˆ (k ) = ψ (u (k − 1),..., u (k − l ) ) ≅ u (k ) , se poate utiliza schema bloc din figura 8.10.
192
u(k)
+
+
e(k)
cuantizor
eq(k)
-
Predictor cu întârzieri
uˆ (k )
Figura 8.10. Principiul compresiei cu predicţie. Această metodă nu se foloseşte în practică deoarece în e(k) se vor acumula erorile de cuantizare a lui u(k). Pentru a elimina acest dezavantaj se foloseşte metoda: uˆ (k ) = ψ (uˆ (k − 1),..., uˆ (k − l ) ) (8.20) e(k) u(k) + + -
eq(k)= e(k)+∆e(k) cuantizor
codor
canal de comunicaţie
decodor
eq(k)
+
+
u'(k)
+
uˆ ( k ) Predictor (cu întârzieri)
Predictor (cu întârzieri) uˆ (k )
-
Figura 8.11. Principiul predicţiei cu compresie modificat. În această figură: u ' (k ) ≅ u (k ) +erori de cuantizare. În cazul 2D al compresiei cu predicţie a imaginilor, în mod similar, pentru o imagine u(m,n), imaginea prezisă va fi:
})
({
uˆ (m, n) = ψ u (i, j ) (i , j )∈Wm,n ,
(8.21)
unde Wm,n este fereastra de predicţie, adică vecinătatea punctului curent, în care se face predicţia. Fereastra de predicţie trebuie aleasă astfel încât să fie cauzală. De exemplu, pentru o baleiere normală a imaginii (sus-jos, stânga-dreapta), 193
pentru punctul curent, fereastra de predicţie trebuie aleasă ca în figura de mai jos (Wm,n şi nu W’m,n), pentru ca valorile ce intervin în predicţie să fie deja calculate:
Wm,n W’m,n
Figura 8.12. Constrângerea de cauzalitate a ferestrei de predicţie.
8.2.3. Compresia cu transformate După cum s-a arătat şi în capitolul de transformări ale imaginilor, transformările integrale unitare se pot aplica cu succes şi la compresia imaginilor. Într-un proces de compresie cu transformate se porneşte de la o imagine iniţială I în care energia este repartizată relativ uniform. Acesteia i se aplică o transformare T în urma căreia se obţine o imagine în care energia este concentrată în mult mai puţine componente decât în imaginea iniţială (de dorit în componentele cu indici mici). . . . .. I. . .
T
....
....... T(I)
canal
codare
T-1
Tˆ ( I )
. . . . Iˆ .
decodare
Figura 8.13. Principiul compresiei cu transformate. La codare se reţin (respectiv se transmit) numai componentele în care este concentrată majoritatea energiei. La decodare (la recepţie), 194
imaginii recepţionate, ce conţine informaţie doar în unele componente, i se aplică transformarea inversă, obţinându-se imaginea iniţială I (în practică se obţine de fapt o aproximaţie a acesteia). Dacă transformarea este liniară şi unitară, deoarece transformările unitare conservă energia ( E ( I ) = E (T ( I ) ) ), se poate controla cât anume se pierde din energia imaginii iniţiale prin netransmiterea (sau nerecepţionarea) componentelor nesemnificative din imaginea transformată. Mărimile cantitative care definesc eficienţa compresiei sunt: •
raportul de compresie. Raportul de compresie se defineşte ca fiind raportul dintre cantitatea de informaţie din imaginea iniţială şi cantiatea de informaţie din imaginea transformată care se transmite. În prezent, factorii de compresie uzuali sunt de 10...20.
•
raportul semnal-zgomot. Raportul semnal-zgomot exprimă o măsură a calităţii compresiei. Raportul semnal-zgomot al compresiei se defineşte ca fiind raportul dintre energia imaginii iniţiale şi energia erorii de compresie: RSZ =
EI = E I − Iˆ
∑ I (i, j )
2
i, j
2 ∑ I (i, j ) − Iˆ(i, j )
(8.22)
i, j
În decibeli: RSZ dB = 10 ⋅ log10 RSZ . Valori acceptabile ale RSZdB sunt cele de peste 30 dB: RSZdB >30. Pentru valori peste 30 dB, ochiul uman nu mai poate distinge diferenţele între imaginea originală şi cea transformată. Pot să apară şi zgomote cu RSZdB >30 dB, de exemplu pentru imagini cu multe contururi, prin aplicarea transformatei Fourier discrete DFT. În acest caz, la aplicarea transformatei inverse, contururile vor fi puternic afectate, deoarece componentele de frecvenţă înaltă (contururile) vor fi eliminate.
195
Dintre transformările prezentate anterior, cea mai utilizată transformare pentru compresie este transformata Cosinus discretă DCT, care elimină dezavantajele transformatei Fourier DFT. numărul de operaţii sau complexitatea algoritmului de compresie. În capitolul de transformări ale imaginilor, s-a arătat că pentru transformările unitare bidimensionale, ale unor imagini de dimensiuni N×N, care admit algoritm rapid de calcul, numărul de operaţii este: n x ≈ N 2 ⋅ log 2 N
(8.23)
Complexitatea algoritmului se poate reduce dacă nu se aplică transformarea pe întreaga imagine, ci se împarte imaginea în blocuri şi se face codarea fiecărui bloc transformat. T T n
...
n
...
•
...
...
Figura 8.14. Principiu de reducere a complexităţii compresiei cu transformate. Pentru blocuri de dimensiunea n×n, complexitatea algoritmului este: n x'
(
)
N ≈ n ⋅ log 2 n ⋅ n 2
2
(8.24)
Prin urmare, câştigul în ce priveşte complexitatea algortimului de compresie este:
196
nx n x'
=
N 2 ⋅ n 2 ⋅ log 2 N N 2 ⋅ n 2 ⋅ log 2 n
= log n N
(8.25)
Pentru N=512 şi n=8 (n fiind standardizat în cadrul JPEG): log 2 N 9 = =3 (8.26) log 2 n 3 ⇒ se obţine un algoritm de 3 ori mai rapid.
Prin urmare, rămâne de stabilit transformarea optimă pentru compresie, care să prezinte caracteristici optime, relativ la cele prezentate anterior. Transformarea optimă va fi stabilită relativ la criteriul erorii pătratice a mediei statistice. Dacă se consideră o variabilă aleatoare, ξ, reală şi se fac multe experimente, în care se măsoară ξ, se vor obţine diferite valori care pot fi ordonate pe o axă, ca în figura 8.15. x x xx xxxx x x
Figura 8.15. Ordonarea valorilor unei variabile aleatoare pe o axă. Dacă aceste valori ar fi infinit de multe, densitatea lor pe axă ar arăta ca funcţia de densitate de probabilitate ce caracterizează variabila ξ:
Figura 8.16. Densitatea valorilor unei variabile aleatoare. Dacă ξ nu este o variabilă aleatoare reală scalară ci un vector cu N componente(reale), la fiecare realizare se pune un punct în spaţiul n-dimensional. ξ va putea fi scris: 197
ξ0 ξ1 . . ξ = . sau ξ = . (8.27) . . ξ ξ N −1 N Prima numerotarea componentelor poate fi avantajoasă uneori, din considerente practice (realizările lui ξ sunt notate cu majuscule): X0 X1 . . (8.28) X = . , respectiv X = . . . X X N −1 N Cel mai simplu caz este N=2, în care rezultatele pot fi reprezentate ca nişte puncte în plan. Rezultatele unor experimente repetate poate arăta ca în figura 8.17.a. Un alt exemplu este prezentat în figura 8.17.b. X2
X2
X1
X1
(a) (b) Figura 8.17. Exemple de reprezentare a realizărilor unei variabile aleatoare. În cazul din figura 8.17.b legătura între valoarea lui X1 şi valoarea lui X2 este mai strânsă. În acest caz, dacă s-ar roti axele ca în figura 8.18, 198
adică dacă în loc de valorile calculate pentru X1 şi X2 s-ar face calcule cu nişte combinaţii ale acestora, care se pot obţine ca noi coordonate ale punctelor într-un sistem de axe rotit, pentru caracterizarea experimentului ar fi aproape suficient y1, adică o singură variabilă şi nu două, deoarece y2 este mic în comparaţie cu y1 şi adesea neglijabil, în asemenea situaţii. X2 y1 y2
X1
Figura 8.18. Exemplu de reprezentare a realizărilor unei variabile aleatoare, într-un sistem de axe rotit. Din astfel de motive, prezintă interes „transformări” ale vectorilor aleatori (în general N-dimensionali) şi fiindcă cele mai simple sunt transformările liniare, vor fi studiate doar acestea. O transformare liniară de la RN se scrie cu ajutorul unei matrici A: η = A ⋅ξ ,
a 00 . unde: A = . . a 0, N −1
.
.
.
.
.
.
199
a 0 N −1 . , . . a N −1, N −1
(8.29)
ξ0 η0 . . ξ = . şi η = . . (8.30) . . ξ η N −1 N −1 Scopul este de a găsi o transformare care să fie optimă dintr-un anumit punct de vedere, şi anume să facă cât mai mici cât mai multe componente ale lui η astfel încât ξ să se poată aproxima cât mai bine prin cât mai puţine numere. Bineînţeles, fiind vorba de variabile aleatoare, minimizarea acestora trebuie înţeleasă în sensul de eroare medie. Mai exact, instrumentul matematic adecvat este eroarea medie pătratică fiindcă este o funcţie derivabilă: la eroarea medie ar trebui vorbit de modul, dar modulul nu este o funcţie derivabilă. Pentru calculele ce urmează se presupune, de asemenea, că toate componentele lui ξ (ξ0, ξ1,…, ξN-1) sunt variabile aleatoare de medie nulă. Se presupune că se cunosc nişte informaţii minime despre aceste variabile şi anume, se presupun cunoscute mediile şi corelaţia între componente, respectiv dispersia fiecăreia, adică momentele de ordinul 1 ( ξ k ) şi de ordinul 2 ( ξ i ξ j ). Dacă nu ar fi variabile aleatoare de medie nulă, s-ar putea scade media (cunoscută), iar rezultatul ar fi variabile aleatoare de medie nulă. Se va limita căutarea optimului printre matricile unitare (care, de fapt, sunt generalizări ale rotaţiilor din plan), deoarece ele păstrează produsul scalar şi norma, adică dacă: η ' = A ⋅ ξ ' şi η " = A ⋅ ξ " (8.31) atunci:
ξ 'T ⋅ξ " = (ξ ' , ξ " ) = (η ' ,η " ) = η 'T ⋅η "
(8.32)
şi, în particular:
ξ'
2
= (ξ ' , ξ ' ) = (η ' ,η ' ) = η '
200
2
(8.33)
În aceste relaţii, vectorii sunt consideraţi ca fiind vectori coloană de dimensiune N×1, ca în relaţia (8.30). Produsul scalar a doi vectori ξ şi η a fost notat cu (ξ, η). În notaţie matriceală, acest produs se scrie ca un produs ξ T·η, între o matrice (vector) 1×N şi o matrice N×1 şi are ca N −1 rezultat o matrice (scalar) 1×1: ξ T ⋅ η = ∑ ξ i ⋅ η i . i =0
Aceste matrici mai au proprietatea că A −1 = AT (vezi capitolul 4, paragrafele 4.1 şi 4.2), adică transformarea inversă se obţine foarte simplu, şi că dacă se notează cu (a 0T , a1T ,..., a TN −1 ) liniile lui A, respectiv cu (a0 , a1 ,..., a N −1 ) coloanele lui AT:
[
]
AT = a0 , a1 ,..., a N −1 , (8.34) atunci vectorii ai sunt ortogonali doi câte doi şi de normă 1. Prin urmare, se va căuta printre matricile unitare o matrice L cu proprietăţile de optimalitate descrise. Cu alte cuvinte, se porneşte de la relaţiile:
η0 . η = L ⋅ ξ şi ξ = LT ⋅ η , unde η = . . (8.35) . η N −1 Primul pas ar fi să se determine cum ar trebui să fie L pentru ca eroarea medie pătratică (asupra lui ξ) să fie minimă, dacă se înlocuieşte ηN-1 cu zero, adică dacă se neglijează o componentă.
η0 . ~ Fie: η~ = . şi ξ = LT ⋅ η~ η N −2 0
201
(8.36)
~2 2 Este evident că: ξ − ξ = η − η~ , deoarece: ~2
ξ −ξ
(
) (
) [
] [
]
T ~T ~ = ξ − ξ ⋅ ξ − ξ = LT ⋅ (η − η~ ) ⋅ LT ⋅ (η − η~ ) =
= (η − η~ )T ⋅ L ⋅ LT ⋅ (η − η~ ) = (η − η~ )T ⋅ (η − η~ ) = η − η~
2
(8.37)
Prin urmare:
ε 2 = η − η~
2
N −1
= ∑ (η k − η~k )2 = η N2 −1 = η N −1 ⋅ η N −1 = k =0
= l TN −1 ⋅ ξ ⋅ ξ T ⋅ l N −1 = l TN −1 ⋅ ξ ⋅ ξ T ⋅ l N −1 = l TN −1 ⋅ K ξ ⋅ l N −1 (8.38) unde Kξ este matricea de covariaţie (şi de corelaţie, deoarece este de medie nulă) a vectorului ξ. Alegerea lui lN-1 (lk fiind coloanele matricii LT) trebuie făcută astfel încât acest ε2 să fie minim, cu constrângerea că lN-1 este un vector unitar: l TN-1 ⋅ l N −1 = 1
(8.39)
Folosind metoda multiplicatorilor lui Lagrange, aceasta înseamnă să minimizăm liber funcţia:
(
)
Ψ = l TN −1 ⋅ K ξ ⋅ l N −1 − λ ⋅ l TN −1 ⋅ l N −1 − 1
(8.40)
Derivarea lui Ψ în raport cu componentele lui lN-1 şi anularea ei, se scrie mai compact ca anularea gradientului lui Ψ, iar:
(
)
∇ v Ψ vT ⋅ A ⋅ v = 2 ⋅ A ⋅ v
(8.41)
pentru orice vector v şi orice matrice A, astfel că: ∇ lN −1 Ψ = 2 ⋅ K ξ ⋅ l N −1 − 2 ⋅ λ ⋅ l N −1 = 0
(8.42)
adică lN-1 trebuie să fie vector propriu al matricii de autocovariaţie a lui ξ: K ξ ⋅ l N −1 = λ ⋅ l N −1
(8.43)
Cum în acest caz eroarea devine:
ε = l TN −1 ⋅ K ξ ⋅ l N −1 = l TN −1 ⋅ (λ ⋅ l N −1 ) = λ ⋅ l TN −1 ⋅ l N −1 = λ ⋅ l N −1 202
2
=λ
(deoarece lN-1 este de modul 1) este clar că, pentru a avea un minim, trebuie ca λ să fie cea mai mică valoare proprie a lui Kξ (se ştie că Kξ este o matrice simetrică pozitiv definită şi prin urmare are N valori proprii pozitive, iar vectorii proprii sunt ortogonali, de normă 1). η0 . Vectorul η = . la care se doreşte ca ultima componentă să . η N −1 fie minimă, are proprietatea interesantă că:
(
η k ⋅ η N −1 = l kT ⋅ ξ ⋅ l TN −1 ⋅ ξ
)T = l KT ⋅ ξ ⋅ ξ T ⋅ l N −1 = l KT ⋅ ξ ⋅ ξ T ⋅ l N −1 =
= l kT ⋅ K ξ ⋅ l N −1 = l kT ⋅ (λ min ⋅ l N −1 ) = λ min ⋅ l kT ⋅ l N −1 = 0
(8.44)
pentru ∀k ≠ N − 1 , deoarece L este o matrice ortogonală (unitară, reală). Cu alte cuvinte ηN-1 este necorelată cu toate celelalte componente ale lui η. Adică Kη are forma:
η 0 ⋅ η 0 η 0 ⋅ η1 η1 ⋅ η 0 . Kη = . . 0 0
.
.
. η 0 ⋅ η N −2
η1 ⋅ η N −2
.
.
.
0
0 0 λ min
(8.45)
Dacă ηN-1 este necorelată cu toate celelalte componente ηk ale lui η, ea este necorelată şi cu orice combinaţie liniară a acestor componente. Aceste lucruri permit să se spună că se poate repeta raţionamentul precedent pentru vectorul:
203
η0 . η'= . . η N −2 adică să se caute o transformare unitară care să minimizeze ultima componentă (media ei pătratică) şi apoi tot aşa. De fapt, deoarece ηN-1 este necorelată cu toate celelalte componente ηk, transformarea căutată se poate considera şi ea N-dimensională, cu ultima linie l TN −1 cunoscută. Se ajunge la concluzia că matricea L care compactează cel mai bine energia lui ξ adică pe ξ
2
în câteva componente, este cea pentru care
coloanele lui LT sunt vectorii proprii ai matricii Kξ în ordinea descrescătoare a valorilor proprii, adică dacă: LT = [l 0 ,..., l N −1 ], atunci lk satisface relaţia: K ξ ⋅ l k = λk ⋅ l k
(8.46)
unde: λ0 ≥ λ1 ≥ λ 2 ≥ ... ≥ λ N −1 . Transformarea unitară L = (l1 ...l n )T , în care K ξ ⋅ l k = λ k ⋅ l k se numeşte transformata Karhunen-Loeve (KL). Eroarea ε dată de compresia cu transformata KL când se înlocuiesc cu constante ultimele „N-m” componente ale lui Φ este: N
N
N
k =m +1
k = m+1
ε = ∑ Φ Tk ⋅ K x ⋅ Φ k = ∑ λ k ⋅ Φ Tk ⋅ Φ k = ∑ λ k k = m +1
1 424 3 λk ⋅Φ k
Construcţia matricii transformării KL se poate face astfel: • se determină valorile proprii (λ1,..., λN) ale lui Kξ; • se ordonează descrescător: λ(1)>...> λ(N) 204
(8.55)
• vectorii proprii corespunzători sunt ordonaţi în aceeaşi ordine:
(
)
L = l (1) l ( 2) ...l ( N ) T
(8.56)
În acest mod, energia lui η = L ⋅ ξ va fi concentrată în primele componente. Exemplu:
În cazul unei transformări separabile, aplicată unei imagini IN×N, se determină matricea de autocovariaţie a liniilor Kl, matricea de autocovariaţie a coloanelor Kc (figura 8.19):
(
)(
Kc = c − c ⋅ c − c
)T , K l = (l − l )⋅ (l − l )T ,
(8.57)
unde mediile statistice pot fi estimate ca medii aritmetice pe mai multe imagini. c1 c2 cN l1 l2 ...
... I
lN
Figura 8.19. Calculul matricilor de autocovariaţie pe linii şi pe coloane. Media statistică pe o singură imagine este: l≅ ⇒ Kl ≅
1 N 1 N , l c ≅ ∑i ∑ ci N i =1 N i =1
(
1 N ∑ li − l N i =1
(8.58)
)T ⋅ (li − l ), K c ≅ N1 ∑ (ci − c )T ⋅ (ci − c ) N
i =1
205
(8.59)
După determinarea matricelor Kl şi Kc, se determină vectorii proprii şi valorile proprii: K l → Ll → λl , K c → Lc → λc (8.60) Cu aceste mărimi calculate, pentru o imagine dată U, imaginea V obţinută prin transformarea KL bidimensională a imaginii iniţiale U, este: V = Lc ⋅ U ⋅ LTl . . . .. U. .
T
....... V
Figura 8.16. Principiul compresiei cu transformata T=KL. În concluzie, transformarea optimă din punct de vedere al compresiei este transformarea K-L (Karhunen-Loeve). Aceasta realizează decorelarea elementelor imaginii şi transformatei şi prin aceasta compactarea maximă a energiei semnalului (imaginii) în primele componente. Cu toate aceste avantaje, deoarece transformata K-L depinde de statistica imaginii şi datorită volumului mare de calcule necesare (mai ales în cazul imaginilor de dimensiuni mari), transformata K–L este dificil de implementat practic. Totuşi, pentru imagini cu corelaţie mare, transformata K-L se poate aproxima cu succes prin transformata Cosinus Discretă (DCT=Discrete Cosine Transform) mult mai rapidă şi mai uşor de implementat.
206
Bibliografie
1. Buzuloiu V., Note de curs – Prelucrarea şi analiza imaginilor, Universitatea din Oradea, 1996-2000. 2. Castleman K. R., “Digital Image Processing”, Editura Prentice-Hall, 1996. 3. Gonzales R.C., Woods R.E., “Digital Image Processing”, Editura Prentice-Hall, 2002. 4. Jahne B., “Digital Image Processing”, Editura Springer, 2002. 5. Jain A. K., “Fundamentals of Digital Image Processing”, Editura Prentice-Hall Inc., 1989. 6. Jiang M., “Digital Image Processing”, Curs, Peking University, 2001. 7. Pitas I., “Digital Image Processing. Algorithms and applications”, Editura Wiley-Interscience, 2000. 8. Pratt W.K., “Digital Image Processing”, Editura John Wiley & Sons 2001. 9. Spătaru A., “Teoria transmiterii informaţiei”, Editura Didactică şi Pedagogică, Bucureşti, 1984. 10. Vertan C., Ciuc M., “Tehnici fundamentale de prelucrarea şi analiza imaginilor”, Editura MatrixRom, Bucureşti, 2007. 11. Wahl F., “Digital Image Signal Processing”, Editura Artech House, 1987.
207
