Wykład nr 8
Transformata Fouriera - zastosowanie w przetwarzaniu obrazów cd.
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 136x136 pikseli z pasma R:

pix=datafinger[[Range[2,137],Range[1,136],1]];
Dimensions[pix]

fingerplot=ListDensityPlot[pix,Mesh->False]

Rotacja grafiki
Na podstawie: 
Fast Fourier method for the accurate rotation of sampled images
Kieran G. Larkin , Michael A. Oldfield , Hanno Klemm
http://www.nontrivialzeros.net/KGL_Papers/27_Rotation_Paper_1997_qualityscan_OCR.pdf

Przesunięcia

Rozważmy pewien sygnał

sygnal[t_]:=Exp[-t/2]*(Sin[2 Pi t] UnitStep[t (3-t)]+Sin[4 Pi t] UnitStep[(t-3) (6-t)])
samplesig=Table[sygnal[n/20],{n,1,20*6}];
dim=Dimensions[samplesig][[1]]
syg1=ListPlot[samplesig,Joined->True,PlotRange->All]

Poniżej:
1. Wykonanie transformaty Fouriera
2. Przemnożenie transformaty przez pewną fazę
3. Wykonanie odwrotnej transformaty
W wyniku otrzynujeny przesunięty sygnał

shift=55;
FT=Fourier[samplesig];
ftprzes=Flatten[Table[FT[[j]]*Exp[2 Pi I shift (j-1)/dim ],{j,1,dim}]];
przes=Re[InverseFourier[ftprzes]];
ListPlot[przes,Joined->True,PlotRange->All]

Obroty

Poniżej macierz obrotu w dwóch wymiarach o kąt θ :

Clear[θ,M]
M[θ_]:={{Cos[θ], Sin[θ]},{-Sin[θ],Cos[θ]}};

Taki obrót można także zrealizować za pomocą dwóch macierzy translacyjnych w kierunku powiedzmy x i y :

(Subscript[M, x])={{1,Tan[θ/2]},{0,1}};
(Subscript[M, y])={{1,0},{-Sin[θ],1}};

Poniżej macierz obrotu i złożenie obrotu z macierzy translacyjnych

M[θ]//MatrixForm
Simplify[(Subscript[M, x]).(Subscript[M, y]).(Subscript[M, x])]//MatrixForm

Poniżej działanie macierzy translacyjnych na pewnien wektor (możemy przyjąć, że x i y to współrzędne pewnego punktu)

Clear[θ]
(Subscript[M, x]).{x,y}//MatrixForm
(Subscript[M, y]).{x,y}//MatrixForm

Poniżej pokazano działanie macierzy obrotu na pewien punkt w kartezjańskiem układzie współrzędnych

θ=0.2Pi;
p1={10,3};
p2=N[M[θ].p1];
p3=N[(Subscript[M, x]).p1];
p4=N[(Subscript[M, y]).p3];
p5=N[(Subscript[M, x]).p4];
ListPlot[{p1,p3,p4,p5},PlotRange->All,PlotStyle->PointSize[0.02]]

Obroty grafiką

Krok 1 - Ustalamy kąt obrotu

θ=0.25Pi;
xobr=68;
yobr=68;
a=Tan[θ/2];
b=-Sin[θ];

Krok 2 - Wykonujemy transformatę Fouriera grafiki na wierszach
Krok 3 - Mnożymy transformatę z kroku 2 przez czynnik fazowy
Krok 4 - Wykonujemy odwrotną transformatę z kroku 3 w wyniku której dane zostają przesunięte
Krok 5 - Całość wykonujemy na wszystkich wierszach (w tym przypadku 136)
Krok 6 - Potrzebujemy transponować otrzymaną macierz

rotx=Transpose[Table[
InverseFourier[
Table[
Exp[ -2 Pi I Round[ a  (y-yobr)] (x2-1)/136]*Fourier[
Table[
pix[[x,y]],{x,1,136}
]
][[x2]],{x2,1,136}
]
],{y,1,136}]];
rotxplot=ListDensityPlot[rotx,Mesh->False]

Krok 7 - Na macierzy z kroku 6 wykonujemy transformatę Fouriera tym razem na kolumnach
Krok 8 - Mnożymy transformatę z kroku 7 przez czynnik fazowy
Krok 9 - Wykonujemy odwrotną transformaę z kroku 8 w wyniku której dane zostają przesunięte
Krok 10 - Całość wykonujemy na wszystkich kolumnach (w tym przypadku 136)

rotyx=Table[
InverseFourier[
Table[
Exp[ -2 Pi I Round[b  (x-xobr)] (y2-1)/136]*Fourier[
Table[
rotx[[x,y]],{y,1,136}
]
][[y2]],{y2,1,136}
]
],{x,1,136}];

rotyxplot=ListDensityPlot[rotyx,Mesh->False]

Krok 11 - Wykonujemy transformatę Fouriera macierzy z kroku 10 na wierszach
Krok 12 - Mnożymy transformatę z kroku 11 przez czynnik fazowy
Krok 13 - Wykonujemy odwrotną transformatę z kroku 12 w wyniku której dane zostają przesunięte
Krok 14 - Całość wykonujemy na wszystkich wierszach (w tym przypadku 136)
Krok 15 - Potrzebujemy transponować otrzymaną macierz

rotθ=Transpose[Table[
InverseFourier[
Table[
Exp[ -2 Pi I Round[ a  (y-yobr)] (x2-1)/136]*Fourier[
Table[
rotyx[[x,y]],{x,1,136}
]
][[x2]],{x2,1,136}
]
],{y,1,136}]];
rotplot=ListDensityPlot[Abs[rotθ],Mesh->False]

GraphicsRow[{fingerplot,rotxplot,rotyxplot,rotplot}]

Ciekawostka

pix4=Re[InverseFourier[Fourier[pix],FourierParameters->{0,4}]];
finger4plot=ListDensityPlot[pix4,Mesh->False]

Rozpoznanie obrazów twarzy

Analiza transformaty Fouriera

Import obrazka:
ImageDimensions[] - zwraca szerokość i wysokość grafiki
Dimensions[] - zwraca rzędy i kolumny danych, a także ilość pasm

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

{Show[woman],
Image[datawoman,ColorSpace->"RGB"]}

W tym miejscu wybieramy 400x320 pikseli z pasma R:

pixwoman=datawoman[[All,All,1]];
Dimensions[pixwoman]

{ListDensityPlot[pixwoman,Mesh->False,AspectRatio->400/320],
womanplot=Image[pixwoman]}

ftwoman=Fourier[pixwoman];

W tym miejscu przedstawiamy transformatę Fouriera twarzy: w środku mamy niskie częstotliwości na zewnątrz wysokie, jasny kolor oznacza duży wkład, ciemny zerowy. Poniżej przedstawiono moduł i argument transformaty:

absftw=ListDensityPlot[RotateLeft[Abs[ftwoman],{200,160}],Mesh->False]
argftw=ListDensityPlot[RotateLeft[Arg[ftwoman],{200,160}],Mesh->False]

Image[RotateLeft[Abs[ftwoman],{200,160}]]
Image[RotateLeft[Arg[ftwoman],{200,160}]]

Pytanie: gdzie zapisane są najważniejsze informacje o twarzy w module czy argumencie transformaty Fouriera?

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

{Show[man],
Image[dataman,ColorSpace->"RGB"]}

pixman=dataman[[All,All,1]];
Dimensions[pixman]

{ListDensityPlot[pixman,Mesh->False,AspectRatio->400/320],
manplot=Image[pixman]}

ftman=Fourier[pixman];
absftm=ListDensityPlot[RotateLeft[Abs[ftman],{200,160}],Mesh->False]
argftm=ListDensityPlot[RotateLeft[Arg[ftman],{200,160}],Mesh->False]

Poniżej mamy przedstawione: w pierwszym wierszu transformatę Fouriera twarzy kobiety, w drugim wierszu twarzy mężczyzny

GraphicsGrid[{{absftw,argftw},{absftm,argftm}}]

A teraz wykonamy pewien eksperyment: tzw. “hybrydy”

mixmw=Abs[ftman] Exp[I Arg[ftwoman]];
mw=InverseFourier[mixmw];
absmargw=Image[Chop[mw]]

mixwm=Abs[ftwoman] Exp[I Arg[ftman]];
wm=InverseFourier[mixwm];
abswargm=Image[Chop[wm]]

Poniżej mamy przedstawione: w pierwszym wierszu oryginalne twarze kobiety i mężczyzny, w drugim wierszu po lewej: obraz uzyskany z modułu TF twarzy mężczyzny i argumentu FT twarzy kobiety, po prawej: obraz uzyskany z modułu TF twarzy kobiety i argumentu FT twarzy mężczyzny

GraphicsGrid[{{womanplot,manplot},{absmargw,abswargm}}]

Powyżej został pokazany eksperyment Oppenheima, który pokazuje, że większość informacji o twarzy zawarta jest w argumencie transformaty Fouriera.

I jeszcze jedna twarz do analizy widma:

woman2=Import["c:/Users/wrzask/Desktop/PO/woman2.jpg"];
datawoman2=ImageData[woman2];
pixwoman2=datawoman2[[Range[21,420],Range[25,344],1]];
Dimensions[pixwoman2]
{ListDensityPlot[pixwoman2,Mesh->False,AspectRatio->400/320],
woman2plot=Image[pixwoman2]}
ftwoman2=Fourier[pixwoman2];

absftw2=ListDensityPlot[RotateLeft[Abs[ftwoman2],{200,160}],Mesh->False]
argftw2=ListDensityPlot[RotateLeft[Arg[ftwoman2],{200,160}],Mesh->False]

GraphicsGrid[{{absftw,absftm,absftw2},{womanplot,manplot,woman2plot}}]

Testy - city landscape

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

Show[city1]
pixcity1=datacity1[[All,All,1]];
ftcity1=Fourier[pixcity1];
{city1plot=Image[pixcity1],
absftcity1=ListDensityPlot[RotateLeft[Abs[ftcity1],{84,150}],Mesh->False]}

city2=Import["c:/Users/wrzask/Desktop/PO/city2.jpg"];
datacity2=ImageData[city2];
Dimensions[datacity2]

Show[city2]
pixcity2=datacity2[[All,All,1]];
ftcity2=Fourier[pixcity2];
{city2plot=Image[pixcity2],
absftcity2=ListDensityPlot[RotateLeft[Abs[ftcity2],{97,130}],Mesh->False]}

city3=Import["c:/Users/wrzask/Desktop/PO/city3.jpg"];
datacity3=ImageData[city3];
Dimensions[datacity3]

Show[city3]
pixcity3=datacity3[[All,All,1]];
ftcity3=Fourier[pixcity3];
{city3plot=Image[pixcity3],
absftcity3=ListDensityPlot[RotateLeft[Abs[ftcity3],{84,150}],Mesh->False]}

city4=Import["c:/Users/wrzask/Desktop/PO/city4.jpg"];
datacity4=ImageData[city4];
Dimensions[datacity4]

Show[city4]
pixcity4=datacity4[[All,All,1]];
ftcity4=Fourier[pixcity4];
{city4plot=Image[pixcity4],
absftcity4=ListDensityPlot[RotateLeft[Abs[ftcity4],{137,91}],Mesh->False]}

GraphicsGrid[{{absftcity1,absftcity2,absftcity3,absftcity4},{city1plot,city2plot,city3plot,city4plot}}]

Testy - beach landscape

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

Show[beach1]
pixbeach1=databeach1[[All,All,1]];
ftbeach1=Fourier[pixbeach1];
{beach1plot=Image[pixbeach1],
absftbeach1=ListDensityPlot[RotateLeft[Abs[ftbeach1],{97,129}],Mesh->False]}

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

Show[beach2]
pixbeach2=databeach2[[All,All,1]];
ftbeach2=Fourier[pixbeach2];
{beach2plot=Image[pixbeach2],
absftbeach2=ListDensityPlot[RotateLeft[Abs[ftbeach2],{91,137}],Mesh->False]}

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

Show[beach3]
pixbeach3=databeach3[[All,All,1]];
ftbeach3=Fourier[pixbeach3];
{beach3plot=Image[pixbeach3],
absftbeach3=ListDensityPlot[RotateLeft[Abs[ftbeach3],{112,112}],Mesh->False]}

GraphicsGrid[{{absftbeach1,absftbeach2,absftbeach3},{beach1plot,beach2plot,beach3plot}}]

Pytania na kolokwium

nie ma pytań

Zadania na laboratorium

Zadanie 11 - 4p.
Należy:
1. Przygotować dwa zdjęcia najlepiej twarzy o tych samych wymiarach (można je “przyciąć”).
2. Wykonać dwuwymiarową transformatę Fouriera (TF) obu grafik. W obu przypadkach na wykresach narysować: Re[], Im[], Abs[] i Arg[] z TF.
3. Wykonać i przedstawić na wykresach dla obu grafik (twarzy) odwrotną transformatę Fouriera z:
    TF
    Re[] z TF
    I Im[] z TF (z przodu I-jednostka urojona)
    -I Im[] z TF (z przodu I-jednostka urojona)
    Abs[] z TF
    Exp[I Arg[]] z TF (I-jednostka urojona)
    Exp[-I Arg[]] z TF (I-jednostka urojona)
4. Podac wnioski z ćwiczenia
W razie problemów z przesłaniem ćwiczenia można wykasować rysunki.
5. Wykonać “hybrydy” z obu grafiki, jak w materałach z wykładu lub wg własnego pomysłu.
5. Opisać zastosowaną w ćwiczeniu funkcję Mathematici, podać trzy własne przykłady zastosowania.

Dodatkowe zadania na laboratorium (projekty)

Można zaproponować własny projekt dotyczący aktykułów poniżej. Np. analiza widm martwej i żywej natury lub twarzy męskich i damkskich itp.

Źródła/literatura

http://www.nontrivialzeros.net/KGL_Papers/27_Rotation_Paper_1997_qualityscan_OCR.pdf
http://www.naun.org/journals/computers/19-632.pdf
http://web.mit.edu/torralba/www/ne3302.pdf