INF2310 - V?r 2017

INF2310 v?r 2017 - L?sningshint 9

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 -1

Konvolusjonsfilteret 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 H₂ vil utjevne mindre fordi det bruker et mindre naboskap til ? utjevne enn det filteret H₁ gj?r (og filtrene ellers har lik relativ nedgang fra sitt senterpunkt). Filteret H₂ vil derimot for?rsake mindre markant ringing fordi ringene vil forekomme n?rmere intensitetskantene og hver ring vil utbre seg kortere enn for H₁.
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.

Bidraget langs aksene er normalt sett prim?rt for?rsaket av bilderanddiskontinuiteten og vil derfor ofte bli tydelig p?virket av ? bruke en vindusfunksjon.

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];
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;
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:

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;

Publisert 24. apr. 2017 22:02 - Sist endret 24. apr. 2017 22:02