Wykład nr 7
Dwuwymiarowa transformata Fouriera - zastosowanie w przetwarzaniu obrazów
Przygotowała Klaudia Wrzask

Grafika w Mathematice - import, pobieranie danych, export

Poniżej: 
import grafiki z podaną ścieżką Import[]
zamiana grafiki na trojwymiarowa tablicę danych ImageData[]
ImageDimensions[] - na wyjściu otrzymujemy szerokość i wysokość grafiki
Dimensions[] - na wyjściu otrzymujemy rzędy i kolumny z tablicy danych; Uwaga! na odwrócenie wartości w porównaniu z ImageDimensions[]
Min[] Max[] - na wyjściu otrzymujemy minimalną i maksymalną wartość danych

fingerprint=Import["c:/Users/wrzask/Desktop/PO/fingerprint.jpg"]
datafinger=ImageData[fingerprint];
ImageDimensions[fingerprint]
Dimensions[datafinger]
{Min[datafinger],Max[datafinger]}

{Show[fingerprint],
Image[datafinger,ColorSpace->"RGB"]}

W tym miejscu wybieramy 128x128 pikseli z pasma R:

pix=datafinger[[Range[10,137],Range[5,132],1]];
Dimensions[pix]

fingerplot=ListDensityPlot[pix,Mesh->False]

Maski dolno (low-pass mask) i górno przepustowe (high-pass mask)

Tworzymy filtry (maski): ten przepuści niskie częstotliwości (zostawiamy “ogóły”)
Na wykresie kolor biały to 1, niebieski 0. Zakładamy także, że środek to niskie częstotliwości.

lpmask[r_]:=Table[If[(x-64)^2+(y-64)^2≤ r^2,1.0,0.0],{x,1,128},{y,1,128}];
ListDensityPlot[lpmask[20],Mesh->False]

ten przepuści wysokie częstotliwości (zostawiamy szczegóły)

hpmask[r_]:=Table[If[(x-64)^2+(y-64)^2≥r^2,1.0,0.0],{x,1,128},{y,1,128}];
ListDensityPlot[hpmask[20],Mesh->False]

Interpretacja dwuwymiarowej transformaty Fouriera

Przedstawienie dwuwymiarowej transformaty Fouriera naszego odcisku palca: niskie częstotliwości mamy tutaj w pobliżu 0 i 128, natomiast wysokie w środku czyli 64. Pamiętamy o własnościach symetrii transformat Fouriera z poprzednich ćwiczeń.
Jasny kolor na wykresie to duży wkład częstotliwości, ciemny (niebieski) kolor to mały wkład:

Fourier[pix];
ListDensityPlot[Abs[Fourier[pix]],Mesh->False]

Przedstawienie tej samej transformaty Fouriera, tym razem wykres został przesuniety o 64 pikseli. Taki zabieg może się wydawać dziwny, jednak w ten sposób możemy myśleć o niskich częstotliowściach jako tych w środku, a wysokich jako tych na krańcach.

ListDensityPlot[RotateRight[Abs[Fourier[pix]],{64,64}],Mesh->False]

Jeszcze przedstawienie tej samej transformaty, tym razem w 3D:

ListPlot3D[RotateRight[Abs[Fourier[pix]],{64,64}],Mesh->False,PlotRange->{0,2}]

Inne obrazki

obrazek1=Table[If[Sin[2 Pi x/4]≤ 0.5,1.0,0.0],{x,1,128},{y,1,128}];
plot1=ListDensityPlot[obrazek1,Mesh->False];
fplot1=ListDensityPlot[RotateRight[Abs[Fourier[obrazek1]],{64,64}],Mesh->False];
obrazek2=Table[If[Sin[2 Pi x/16]≤ 0.5,1.0,0.0],{x,1,128},{y,1,128}];
plot2=ListDensityPlot[obrazek2,Mesh->False];
fplot2=ListDensityPlot[RotateRight[Abs[Fourier[obrazek2]],{64,64}],Mesh->False];
obrazek3=Table[If[Sin[2 Pi x/32]≤ 0.5,1.0,0.0],{x,1,128},{y,1,128}];
plot3=ListDensityPlot[obrazek3,Mesh->False];
fplot3=ListDensityPlot[RotateRight[Abs[Fourier[obrazek3]],{64,64}],Mesh->False];

GraphicsGrid[{{plot1,plot2, plot3},{fplot1,fplot2,fplot3}}]

obrazek4=Table[If[Sin[2 Pi y/4]≤ 0.5,1.0,0.0],{x,1,128},{y,1,128}];
plot4=ListDensityPlot[obrazek4,Mesh->False];
fplot4=ListDensityPlot[RotateRight[Abs[Fourier[obrazek4]],{64,64}],Mesh->False];
obrazek5=Table[If[Sin[2 Pi y/16]≤ 0.5,1.0,0.0],{x,1,128},{y,1,128}];
plot5=ListDensityPlot[obrazek5,Mesh->False];
fplot5=ListDensityPlot[RotateRight[Abs[Fourier[obrazek5]],{64,64}],Mesh->False];
obrazek6=Table[If[Sin[2 Pi y/32]≤ 0.5,1.0,0.0],{x,1,128},{y,1,128}];
plot6=ListDensityPlot[obrazek6,Mesh->False];
fplot6=ListDensityPlot[RotateRight[Abs[Fourier[obrazek6]],{64,64}],Mesh->False];

GraphicsGrid[{{plot4,plot5,plot6},{fplot4,fplot5,fplot6}}]

obrazek7=Table[If[Sin[2 Pi (x+y)/8]≤ 0.5,1.0,0.0],{x,1,128},{y,1,128}];
plot7=ListDensityPlot[obrazek7,Mesh->False];
fplot7=ListDensityPlot[RotateRight[Abs[Fourier[obrazek7]],{64,64}],Mesh->False];
obrazek8=Table[If[Sin[2 Pi (2 x+y)/8]≤ 0.5,1.0,0.0],{x,1,128},{y,1,128}];
plot8=ListDensityPlot[obrazek8,Mesh->False];
fplot8=ListDensityPlot[RotateRight[Abs[Fourier[obrazek8]],{64,64}],Mesh->False];
obrazek9=Table[If[Sin[2 Pi (x+2 y)/8]≤ 0.5,1.0,0.0],{x,1,128},{y,1,128}];
plot9=ListDensityPlot[obrazek9,Mesh->False];
fplot9=ListDensityPlot[RotateRight[Abs[Fourier[obrazek9]],{64,64}],Mesh->False];

GraphicsGrid[{{plot7,plot8,plot9},{fplot7,fplot8,fplot9}}]

Zastosowanie masek (filtrów)

Idealne maski low-pass i high-pass

A teraz sprawdźmy co otrzymamy poprzez zastosowanie maski (filtra dwuwymiarowego) dolno-przepustowej:
Postępujemy analogicznie jak w przypadku filtrów jednowymiarowych z poprzednich wykładów, czyli:
Krok 1 Wykonujemy dwuwymiarową transformatę Fouriera grafiki i ustawiamy dane w taki sposób aby niskie częstotliwości były w środku RotateRight[]
Krok 2 Nakładamy maskę (przemnażamy każdy element macierzy transformaty przez element macierzy maski) 
Krok 3 Wykonujemy odwrotną transformatę Fouriera z kroku 2; nie ma potrzeby przestawiania macierzy z kroku 2, ponieważ wykonując odwrotną transformate wykonujemy sumę, czyli operację przemienna.

mod1= RotateRight[Fourier[pix],{64,64}] lpmask[20];
ListDensityPlot[Abs[mod1],Mesh->False]
low=ListDensityPlot[Abs[InverseFourier[mod1]],Mesh->False]
ListDensityPlot[Abs[InverseFourier[RotateRight[mod1,{64,64}]]],Mesh->False]

I górno-przepustowej....

mod2= RotateRight[Fourier[pix],{64,64}] hpmask[20];
ListDensityPlot[Abs[mod2],Mesh->False]
high=ListDensityPlot[Abs[InverseFourier[mod2]],Mesh->False]

Poniżej mamy przedstawienie grafiki odcisku palca z zastosowaniem filtrów dolnoprzepustowych (zostawiamy niskie częstotliwości i otrzymujemy ogólny wygląd obrazka) oraz z zastosowaniem filtra górnoprzepustowego (zostawiamy wysokie częstotliwości i otrzymujemy szczegółowy wygląd obrazka )

GraphicsRow[{fingerplot,low,high}]

Zastosowanie masek dolnoprzepustowych:
- blurrowanie, wygładzanie obrazu 
Zastosowanie masek górnoprzepustowych:
- wyostrzanie obrazu, wykrywanie krawędzi

Idealne maski dolno i górno przepustowe mają wadę tzw. ringing effect. Korzystając np. z filtrów Butterworth lub Gausianowskiego możemy załagodzić ten efekt.

Butterworth filtry górno i dolno przepustowe 2D

Filtr Butterworth dolno przepustowy (low pass)
BLP=(1/(1+(((((f(x,y))/fc)))^(2n)))),   f(x,y)=(Sqrt[((((x-Dim(x)/2)/Dim(x)))^2)+((((y-Dim(y)/2)/Dim(y)))^2)]) to poprostu odległość punktu (x,y) od środka

blpmask[p_,fc_] := Table[1/(1+(Sqrt[((x-64)/128)^2+((y-64)/128)^2]/fc)^(2 p)),{x,1,128},{y,1,128}];
ListDensityPlot[blpmask[1,0.1], Mesh -> False]
ListPlot3D[blpmask[1,0.1], Mesh -> False,PlotRange->All]

Filtr Butterworth górno przepustowy (high pass)
BHP=1-BLP

bhpmask[p_,fc_] := Table[1-(1/(1+(Sqrt[((x-64)/128)^2+((y-64)/128)^2]/fc)^(2 p))),{x,1,128},{y,1,128}];
ListDensityPlot[bhpmask[1,0.3], Mesh -> False]
ListPlot3D[bhpmask[4,0.3], Mesh -> False,PlotRange->All]

Poniżej porównanie masek górnoprzepustowych: po lewej Butterwortha po prawej przy zastosowaniu maski idealne widać “ringing effect”

mod3= RotateRight[Fourier[pix],{64,64}] bhpmask[2,0.3];
ListDensityPlot[Abs[mod3],Mesh->False]
highbuterworth=ListDensityPlot[Abs[InverseFourier[mod3]],Mesh->False];
GraphicsRow[{highbuterworth,high}]

Poniżej porównanie masek dolnoprzepustowych: po lewej Butterwortha po prawej przy zastosowaniu maski idealnej także widać “ringing effect”

mod4= RotateRight[Fourier[pix],{64,64}] blpmask[10,0.1];
ListDensityPlot[Abs[mod4],Mesh->False]
lowbuterworth=ListDensityPlot[Abs[InverseFourier[mod4]],Mesh->False];
GraphicsRow[{lowbuterworth,low}]

Gaussianowskie filtry górno i dolno przepustowe 2D

Filtr Gaussianowski dolno przepustowy (low pass)
GLP=(e^(-((f(x,y))^2)/2(σ^2))),   f(x,y)=(Sqrt[((((x-Dim(x)/2)/Dim(x)))^2)+((((y-Dim(y)/2)/Dim(y)))^2)]) to poprostu odległość punktu (x,y) od środka

glpmask[σ_] := Table[Exp[-(((x-64)/128)^2+((y-64)/128)^2)/(2 σ^2)],{x,1,128},{y,1,128}];
ListDensityPlot[glpmask[0.1], Mesh -> False]
ListPlot3D[glpmask[0.1], Mesh -> False,PlotRange->All]

Filtr Gaussianowski górno przepustowy (high pass)
GHP=1-GLP,

ghpmask[σ_] := Table[1-Exp[-(((x-64)/128)^2+((y-64)/128)^2)/(2 σ^2)],{x,1,128},{y,1,128}];
ListDensityPlot[ghpmask[0.1], Mesh -> False]
ListPlot3D[ghpmask[0.1], Mesh -> False,PlotRange->All]

Maski band-pass

Poniżej przedstawiam maski pozostawiąjące zakres częstotliwości

bpmask[r1_,r2_] := Table[If[r1^2≤(x-64)^2+(y-64)^2≤ r2^2,1.0,0.0],{x,1,128},{y,1,128}];
ListDensityPlot[bpmask[20,40], Mesh -> False]

Maski orientacyjne

Poniżej przedstawiam cztery maski orientacyjne

orientmask1 = Table[If[(x >= 64 && y >= 64 && y >= x) || (x <= 64 && y <= 64 && y <= x), 1.0, 0.0], {x, 1, 128}, {y, 1, 128}];
ListDensityPlot[orientmask1, Mesh -> False]
orientmask2 = Table[If[(x >= 64 && y >= 64 && y <= x) || (x <= 64 && y <= 64 && y >= x), 1.0, 0.0], {x, 1, 128}, {y, 1, 128}];
ListDensityPlot[orientmask2, Mesh -> False]
orientmask3 = Table[If[(x <= 64 && y >= 64 && y <= -x + 128) || (x >= 64 && y <= 64 && y >= -x + 128), 1.0, 0.0], {x, 1, 128}, {y, 1, 128}];
ListDensityPlot[orientmask3, Mesh -> False]
orientmask4 = Table[If[(x <= 64 && y >= 64 && y >= -x + 128) || (x >= 64 && y <= 64 && y <= -x + 128), 1.0, 0.0], {x, 1, 128}, {y, 1, 128}];
ListDensityPlot[orientmask4, Mesh -> False]

Dwuwymiarowa Transformata Fouriera - czyszczenie obrazu metodą progową

Procedura odszumiania przy użyciu transformaty Fouriera w 2D;
parametr h, bedę dalej nazywać parametrem odcięcia:
1. Wykonujemy transformatę Fouriera w 2D
2. Zerujemy wszystkie współczynniki tej transformaty poniżej progu h
3.  Wykonujemy odwrotną transformatę Fouriera

Odszumianie[danep_List,h_]:=
Module[
{danek},
danek=Chop[Abs[Fourier[danep]],h];
Chop[InverseFourier[ReplacePart[Fourier[danep],0,Position[danek,0]]]]
]

ReplacePart[expr, new, n] yields an expression in which the n(^th) part of expr is replaced by new. 
Position[expr, pattern] gives a list of the positions at which objects matching pattern appear in expr.

Dowolna funkcja w 2D i jej przedstawienie graficzne w kolorze:

dane=Table[Sin[x]+Cos[y],{x,0,8π,0.1},{y,0,8π,0.1}];

ListDensityPlot[dane,Mesh->False,ColorFunction->Hue]

ListDensityPlot[Re[Fourier[dane]],Mesh->False]

A teraz do naszej funkcji dodamy szumy:

danezszumem=Table[Sin[x]+Cos[y]+Random[] ,{x,0,8π,0.1},{y,0,8π,0.1}];

ListDensityPlot[danezszumem,Mesh->False,ColorFunction->Hue]

Sprawdźmy w jaki sposób "zadziała" procedura odszumiania dla parametru odcięcia h=0.3:

ListDensityPlot[Odszumianie[danezszumem,0.3],Mesh->False,ColorFunction->Hue]

dla parametru odcięcia h=0.5:

ListDensityPlot[Odszumianie[danezszumem,0.5],Mesh->False,ColorFunction->Hue]

dla parametru odcięcia h=3 wykres wygląda już "ładnie":

ListDensityPlot[Odszumianie[danezszumem,3],Mesh->False,ColorFunction->Hue]

Ten parametr odcięcia jest zdecydowanie za duży i wyciął wszystkie składowe czestotliwości horyzontalne:

ListDensityPlot[Odszumianie[danezszumem,126],Mesh->False,ColorFunction->Hue]

Ten wyciana wszystko...

ListDensityPlot[Odszumianie[danezszumem,130],Mesh->False,ColorFunction->Hue]

Transformata Fouriera a dyfrakcja światła

Dyfrakcja na pojedynczej szczelinie:
Slit Diffraction:

Poniżej przedstawiona jest szczelina i jej wykres:
Below we have a slit and its plot:

slit1[m_]:=If[126≤m≤130,1.0,0.0]
tablslit=Table[slit1[m],{m,1,256}];
ListPlot[tablslit,Joined->True]

Wykres Transformaty Fouriera możemy interpretować jako obraz dyfrakcyjny.
The Fourier Transform can be interpreted as the diffraction pattern.

Fourier[tablslit];
ListPlot[Abs[Fourier[tablslit]],Joined->True,PlotRange->All]

Dokonamy jeszcze przesunięcia obrazu dyfrakcyjnego.
We can rearrange this list:

ListPlot[Abs[RotateLeft[Fourier[tablslit],128]],Joined->True,PlotRange->All]

Dyfrakcja na otworze kwadratowym: 
Diffraction on a square hole:

Poniżej przedstawiony jest w postaci tablicy otwór kwadratowy i jego wykres:
Below we have square hole and its plot:

square=Table[If[30≤x≤35 &&30≤y≤35,0.0,1.0 ],{x,1,64},{y,1,64}];
ListDensityPlot[square,Mesh->False]

Wykres Transformaty Fouriera możemy interpretować jako obraz dyfrakcyjny.
The Fourier Transform can be interpreted as the diffraction pattern.

Fourier[square];
ListDensityPlot[Abs[Fourier[square]],Mesh->False]

Dokonamy jeszcze przesunięcia obrazu dyfrakcyjnego.
We can rearrange this list:

ListDensityPlot[RotateRight[Abs[Fourier[square]],{32,32}],Mesh->False]

Warto także przedstawić taki odraz dyfrakcyjny za pomoca wykresu 3D.
Lets also see the diffraction pattern in 3D:

ListPlot3D[RotateRight[Abs[Fourier[square]],{32,32}],PlotRange->{0,0.6}]

Pytania na kolokwium

1. Jakie może być zastosowanie masek górno i dolno przepustowych w grafice? p.p. 
2. Na czym polega algorytm czyszczenia obrazów metoda progową? p.p 
3. Filtry Butterwortha i Gausianowskie. p.r.

Zadania na laboratorium

Zadanie 10 - 3p.
Należy:
1. Przedstawić dowolna funkcję 2D graficznie w kolorze
2. Zaszumić funkcję z punktu 1.
3. Przetestować procedurę odszumiania dla paru współczynników odcięcia

Dodatkowe zadania na laboratorium (projekty)

Zadanie A - 5p.
Poniżej znajdują się obrazki i ich transformaty Fouriera (w środku znajdują się najniższe częstotliwości). Obrazki przedstawiono od lewiej: obrazek1, obrazek2, obrazek3, transformaty są przemieszane i trzeba je dopasować  p.p i p.r (można także oddać jako zadanie dodatkowe)

obrazek1=Table[Sin[2 Pi x/8],{x,1,128},{y,1,128}];
obrazek2=Table[Sin[2 Pi y/8],{x,1,128},{y,1,128}];
obrazek3=Table[Sin[2 Pi (x+y)/8],{x,1,128},{y,1,128}];

obrazek1=Table[Sin[2 Pi x/8],{x,1,128},{y,1,128}];
obrazek2=Table[Sin[2 Pi x/8]+Sin[2 Pi x/16],{x,1,128},{y,1,128}];
obrazek3=Table[Sin[2 Pi x/8]+Sin[2 Pi x/16]+Sin[2 Pi x/32],{x,1,128},{y,1,128}];

obrazek1=Table[(x-64)^2+(y-64)^2,{x,1,128},{y,1,128}];
obrazek2=Table[(x-64)^2+(y-64)^2/64,{x,1,128},{y,1,128}];
obrazek3=Table[(x-64)^2/64+(y-64)^2,{x,1,128},{y,1,128}];

obrazek1=Table[If[200≤(x-64)^2+(y-64)^2≤ 250  ,1.0,0.0],{x,1,128},{y,1,128}];
obrazek2=Table[If[200≤(x-64)^2+(y-64)^2≤ 800,1.0,0.0],{x,1,128},{y,1,128}];
obrazek3=Table[If[200≤(x-64)^2+(y-64)^2≤ 1600,1.0,0.0],{x,1,128},{y,1,128}];

obrazek1=Table[If[200≤(x-64)^2+(y-64)^2≤ 250  ,1.0,0.0],{x,1,128},{y,1,128}];
obrazek2=Table[If[200≤(x-64)^2+(y-64)^2≤ 250 || 800≤(x-64)^2+(y-64)^2≤ 870 ,1.0,0.0],{x,1,128},{y,1,128}];
obrazek3=Table[If[200≤(x-64)^2+(y-64)^2≤ 250 || 800≤(x-64)^2+(y-64)^2≤ 870 || 1600≤(x-64)^2+(y-64)^2≤ 1700 ,1.0,0.0],{x,1,128},{y,1,128}];

obrazek1=Table[If[200≤(x-64)^2+(y-64)^2≤ 250 ,1.0,0.0],{x,1,128},{y,1,128}];
obrazek2=Table[If[200≤(x-64)^2+(y-64)^2≤ 250 ,1.0,0.0]+0.1 Random[],{x,1,128},{y,1,128}];
obrazek3=Table[If[200≤(x-64)^2+(y-64)^2≤ 250 ,1.0,0.0]+Random[],{x,1,128},{y,1,128}];

Zadanie B - 3p.
Należy:
Przetestować maski orientacyjne na dowolnej grafice, może być odcisk palca, gdyż układ kierunkowy linii papilarnych, sprzyja temu ćwiczeniu. Opisać wnioski z ćwiczenia.

Zadanie C 
Zaproponować projekt przetestowania masek górno i dolno przepustowych oraz ich zastosowań.

Źródła/literatura

http : // www.cs.tut.fi/~moncef/SGN - 3016 - DIP/Chap04.pdf