Oppgave 1 - Konvolusjonsteoremet
-
Sirkelkonvolusjon i bildedomenet tilsvarer punktmultiplikasjon i Fourier-domenet, og omvendt.
-
Siden konvolusjonsteoremet sier at effekten av å sirkelkonvolvere et konvolusjonsfilteret og et bilde i bildedomenet er det samme som å punktmultiplisere 2D DFT-en til filteret med 2D DFT-en til bildet i Fourier-domenet, vil verdiene i Fourier-spekteret til filteret direkte reflektere hvilke frekvenser som dempes og hvilke som bevares.
-
Ved å designe filtre i Fourier-domenet kan vi lage vi konvolusjonsfiltre med bestemte frekvensegenskaper.
Filtrene vi designer i Fourier-domenet bør være konjugert symmetriske fordi dette er ekvivalent med at representasjonen av filtret er reelt i bildedomenet.
Det er naturlig å la realdelverdiene ligger i intervallet [0,1] og sette imaginærdelverdiene til 0. Førstnevnte gjør at vi ikke forsterker noen frekvenser og sammen gjør de to kriteriene at det designede konvolusjonsfilteret kun påvirker Fourier-spekteret til bildet det anvendes på, og er dermed lettere å kontrollere og tolke ((en mild versjon av) det første kriteriet trengs fordi realdelene må være ikke-negative for å ikke påvirke fasen). -
Det første spekteret er et vertikalt lavpassfilter, dvs. et filter som demper høyere vertikale frekvenser og bevarer lavere vertikale frekvenser. Dette passer med det første konvolusjonsfilteret ettersom det filteret vil vertikalt smøre ut bildet det konvolveres med. Det andre spekteret er et horisonalt høypassfilter, dvs. et filter som demper lavere horisontale frekvenser og (stort sett) bevarer høyere horisontale frekvenser. Dette passer med det andre konvolusjonsfilteret ettersom det filteret tilnærmer den deriverte i horisontal retning.
Alternativt kunne man resonnert at summen av det første konvolusjonsfilteret er større enn 0, mens summen av det andre konvolusjonsfilteret er lik 0. Dette gjør at DC av det første filteret må være positivt, mens DC av det andre filteret må være 0, og dermed får man assosieringene mellom konvolusjonsfiltrene og Fourier-spektrene. -
Fra konvolusjonsteoremet vet vi at punktmultiplikasjonen av 2D DFT-ene tilsvarer å sirkelkonvolvere nullutvidelsen av de to opprinnelige konvolusjonsfiltrene i bildedomenet. Siden konvolusjonsfiltrene, som har størrelse 3x1 og 1x3, har blitt nullutvidet til 200x200 før sirkelkonvolusjon, så er vi garantert at ikke-null koeffisienter av det ene filteret kun overlapper ikke-null koeffisienter av det andre filteret i ett 3x3-område (når vi skal beregne responsen i pikselposisjonene til 200x200-filtrene, dvs. bevare "inn-bildets" størrelse). I dette 3x3-området vil responsen være konvolusjonen av de opprinnelige filtrene i bildedomenet når vi antar nullutvidelse, som er:
1 0 -1
2 0 -2
1 0 -1Konvolusjonsfilteret som har det viste Fourier-spekteret er nullutvidelsen av dette 3x3-filteret til en størrelse på 200x200. (Vi gjenkjenner forøvrig 3x3-filteret som den horisontale estimatoren i Sobel-operatoren.)
Oppgave 2 - Design av filtre i Fourier-domenet: Ideelle filtre
-
P = 512; Q = 512;
D0 = 0.2;
H1 = zeros(P,Q);
for i = 1:P
for j = 1:Q
if sqrt( ((i-floor(P/2+1))/(P/2))^2 + ...
((j-floor(Q/2+1))/(Q/2))^2 ) <= D0
H1(i,j) = 1;
end
end
end
figure(1); imshow(H1); axis on; -
figure(2); imshow(fftshift(imag(ifft2(ifftshift(H1)))), []); axis on; colorbar;
-
figure(1); imshow(fftshift(real(ifft2(ifftshift(H1)))), []); axis on;
Ringingen i den romlige representasjonen er forårsaket av de maksimalt raske overgangene mellom 0 og 1 i Fourier-domenet. -
D0 = 0.4;
H2 = zeros(P,Q);
for i = 1:P
for j = 1:Q
if sqrt( ((i-floor(P/2+1))/(P/2))^2 + ...
((j-floor(Q/2+1))/(Q/2))^2 ) <= D0
H2(i,j) = 1;
end
end
end
figure(2); imshow(fftshift(real(ifft2(ifftshift(H2)))), []); axis on;Filteret H2 vil utjevne mindre fordi det bruker et mindre naboskap til å utjevne enn det filteret H1 gjør (og filtrene ellers har lik relativ nedgang fra sitt senterpunkt). Filteret H2 vil derimot forårsake mindre markant ringing fordi ringene vil forekomme nærmere intensitetskantene og hver ring vil utbre seg kortere enn for H1. -
f = double(imread('car.png'));
[M, N] = size(f);
F = fftshift( fft2(f, P, Q) );
gLP1 = real( ifft2( ifftshift( F.*H1 ) ) );
gLP1 = gLP1(1:M, 1:N);
gLP2 = real( ifft2( ifftshift( F.*H2 ) ) );
gLP2 = gLP2(1:M, 1:N);
figure(1); imshow(f, []); axis on;
figure(2); imshow(gLP1, []); axis on;
figure(3); imshow(gLP2, []); axis on; -
gHP1 = real( ifft2( ifftshift( F.*(1-H1) ) ) );
gHP1 = gHP1(1:M, 1:N);
gHP2 = real( ifft2( ifftshift( F.*(1-H2) ) ) );
gHP2 = gHP2(1:M, 1:N);
figure(1); imshow(f, []); axis on;
figure(2); imshow(gHP1, []); axis on;
figure(3); imshow(gHP2, []); axis on; -
all(abs(gLP1(:) + gHP1(:) - f(:)) < 10^-10)
all(abs(gLP2(:) + gHP2(:) - f(:)) < 10^-10)Ringingen i bildedomenet er bølgelignende intensitetsvariasjoner som visuelt ser ut til å utgå fra kantene i det filtrerte bildet. Disse variasjonene vil alternere mellom å påvirke intensitetene negativt og positivt i forhold til det opprinnelige bildet. Dersom et filtrert bilde har markant ringing og summen av det bildet og et annet bilde skal bli lik det opprinnelige (altså ufiltrerte) bildet, må det andre bildet ha tilsvarende markant ringing som "demmer opp mot" det filtrerte bildets intensitetsavvik. Derfor må et lavpassfilter og det tilhørende høypassfilteret forårsake tilsvarende grad av ringing, men med omvendt påvirkningsretning av intensitetene (positivt og negativt eller negativt og positivt) - ettersom summen av de filtrerte bildene skal være lik det opprinnelige bildet. (En detalj er at det lavpassfiltrerte resultatet ringer rundt den opprinnelige intensitetsverdien, mens det høypassfiltrerte resultatet ringer rundt 0.)
Oppgave 4 - Aliasing
f = double(imread('forresampling.png'));
[M N] = size(f);
for d = 1:4
figure(d); imshow(f(1:d:end,1:d:end), [0 255]); axis on;
end
for d = 1:4
lnFdS = log( abs(fftshift(fft2(f(1:d:end,1:d:end)))) + 1 );
figure(d); imshow(lnFdS, [0 max(lnFdS(:))]); axis on;
end
Oppgave 5 - Vindusfunksjoner og Fourier-spektre
f = double(imread('lena.png'));
[M,N] = size(f);
alpha = 0.5;
w = (tukeywin(N,alpha)*tukeywin(M,alpha)')';
figure(1); imshow(w); axis on;
lnFS = log( abs(fftshift(fft2(f))) + 1);
figure(2); imshow(lnFS, [0 max(lnFS(:))]); axis on;
lnFwS = log( abs(fftshift(fft2(f.*w))) + 1);
figure(3); imshow(lnFwS, [0 max(lnFwS(:))]); axis on;
To forklaringer på hva vindusfunksjoner gjør (to sider av samme sak):
- Reduserer bidraget i Fourier-spekteret som er forårsaket av diskontinuiteten langs bilderanden.
- Glatter Fourier-spekteret.
Oppgave 7 (ekstraoppgave) - Design av filtre i Fourier-domenet: Praktisk eksempel
f = imread('tekna133.png');
[M,N] = size(f);
P = 2*M; Q = 2*N;
F = fft2(f);
Fp = fft2(f, P, Q);
-
FS = abs(fftshift(F));
lnFS = log( FS + 1 );
figure(1); imshow(lnFS, [0 max(lnFS(:))]); axis on;
FpS = abs(fftshift(Fp));
lnFpS = log( FpS + 1 );
figure(2); imshow(lnFpS, [0 max(lnFpS(:))]); axis on;Periodisitetsantagelsen fører "kun" til markant diskontinuitet i vertikal retning i bildet og dermed tilhørende markant bidrag langs den vertikale aksen i Fourier-spekteret, u-aksen, mens nullutvidelsen fører til markant diskontinuitet langs hele bilderanden og dermed markant bidrag langs både u- og v-aksene. Bidraget fra diskontinuiteten langs bilderanden er klart størst som følge av nullutvidelsen, også langs u-aksen. -
u = [ 38 116 136 157];
v = [ 85 37 38 38];
u_sym = M + 2 - rem(M,2) - u;
v_sym = N + 2 - rem(N,2) - v;
u = [u u_sym];
v = [v v_sym];
up = [ 76 231 273 315];
vp = [170 75 75 75];
up_sym = P + 2 - rem(P,2) - up;
vp_sym = Q + 2 - rem(Q,2) - vp;
up = [up up_sym];
vp = [vp vp_sym]; -
D0 = 0.1;
% Alternativt vba av den tilhørende cut-off-frekvensen:
%D0 = 11.05;
n = 2;
D = zeros(M,N,length(u));
H = ones (M,N);
for x = 1:M
for y = 1:N
for k = 1:length(u)
D(x,y,k) = sqrt(((x-u(k))/(M/2))^2 + ...
((y-v(k))/(N/2))^2);
% Alternativt vba av cut-off-frekvens:
%D(x,y,k) = sqrt((x-u(k))^2 + (y-v(k))^2);
H(x,y) = H(x,y) / (1 + (D0/D(x,y,k))^(2*n));
end
end
end
G = F.*ifftshift(H);
g = real(ifft2(G));
figure(1); imshow(f, [0 255]); axis on;
figure(2); imshow(g, [0 255]); axis on;Støyen er hovedsakelig fjernet. Wraparound-feil er visuelt synlig langs øvre del av bildet. Ringing finnes. -
D0 = 0.1;
% Alternativt vba av den tilhørende cut-off-frekvensen:
%D0 = 22.1;
n = 2;
Dp = zeros(P,Q,length(up));
Hp = ones (P,Q);
for x = 1:P
for y = 1:Q
for k = 1:length(up)
Dp(x,y,k) = sqrt(((x-up(k))/(P/2))^2 + ...
((y-vp(k))/(Q/2))^2);
% Alternativt vba av cut-off-frekvens:
%Dp(x,y,k) = sqrt((x-up(k))^2 + (y-vp(k))^2);
Hp(x,y) = Hp(x,y) / (1+(D0/Dp(x,y,k))^(2*n));
end
end
end
Gp = Fp.*ifftshift(Hp);
gp = real(ifft2(Gp));
gp = gp(1:M,1:N);
figure(3); imshow(gp, [0 255]); axis on;Den nye feilen, som er de vertikale stripene langs hver vertikal bildekant, er forårsaket av at notch-stoppfilteret delvis demper noe av det horisontale aksebidraget. Siden dette aksebidraget er forårsaket av diskontinuiteten langs de vertikale bildekantene som følge av nullutvidelsen, vet vi at det er med å lage den skarpe overgangen fra 0 til gråtoneverdiene langs de vertikale bildekantene, og når vi demper noe av dette bidraget er det dermed naturlig at denne skarpe overgangen ikke lenger representeres perfekt. -
Alternativ 1: Vi har sett at wraparound-feilen forekom vertikalt, mens feilen relatert til bilderanddiskontinuiteten forekom horisontalt. Derfor kan vi (omtrentlig) unngå begge disse feilene samtidig ved å bruke sirkulær indeksering i horisontal retning, mens nullutvidelse i vertikal retning. For å enkelt kunne definere et tilsvarende Butterworth-filtre som i de to forrige punktene kan du sirkulært duplikere horisontalt (i stedet for å kun vertikalt doble størrelse ved nullutvidelse) - da kan vi bruke up og vp direkte. En implementering av dette er:
fp2 = [f f];
Alternativ 2: Alternativet over fungerer veldig bra for akkurat dette bildet, men er litt spesialisert for det aktuelle inn-bildet. Mer generelt kan vi utvide inn-bildet ved bruke refleksiv indeksering. Dette vil forårsake at hvert tidligere frekvenspar nå gir opphav til fire punkter i Fourier-spekteret gjennom å speile det opprinnelige frekvensparet om u-aksen (eller om v-aksen, resultatet er det samme). Likevel er det slik at siden disse nye punktene tilhører utvidelsen vil vi ikke være avhengig av å fjerne de for å fjerne støyen i inn-bildet. Dermed blir en implementering slik:
Fp2 = fft2(fp2, P, Q);
Gp2 = Fp2.*ifftshift(Hp);
gp2 = real(ifft2(Gp2));
gp2 = gp2(1:M,1:N);
figure(4); imshow(gp2, [0 255]); axis on;Fps = fft2(padarray(f, [M N], 'symmetric', 'post'));
FpsS = abs(fftshift(Fps));
lnFpsS = log( FpsS + 1 );
figure(5); imshow(lnFpsS, [0 max(lnFpsS(:))]); axis on;
ups = up;
vps = vp;
% Ved å kommentere ut linjene under kan du se hvordan
% resultatet blir hvis man også fjerner de nye punktene
% i Fourier-spekteret. Du vil isåfall oppdage at
% filtreringsresultatet er noe dårligere pga. at
% vi da fjerner/demper flere frekvenser, og disse
% frekvensene representerte litt "sann" bildeinformasjon
% (altså, disse frekvensene tilstedet i bildet, selv om
% de er svake, og de representerer ikke støy).
% ups = [up(1:4) up up(5:end)];
% vps = [vp vp];
Dps = zeros(P,Q,length(ups));
Hps = ones (P,Q);
for x = 1:P
for y = 1:Q
for k = 1:length(ups)
Dps(x,y,k) = sqrt(((x-ups(k))/(P/2))^2 + ...
((y-vps(k))/(Q/2))^2);
% Alternativt vba av cut-off-frekvens:
%Dps(x,y,k)= sqrt((x-ups(k))^2+(y-vps(k))^2);
Hps(x,y)= Hps(x,y)/(1+(D0/Dps(x,y,k))^(2*n));
end
end
end
Gps = Fps.*ifftshift(Hps);
gps = real(ifft2(Gps));
gps = gps(1:M,1:N);
figure(6); imshow(gps, [0 255]); axis on;