IN2070 - V?r 2021

IN2070 v?r 2021 - 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. Om filteret er lite i utstrekning, vil en nullutvidelse typisk v?re n?dvendig for ? f? visualisert effekten filteret vil ha p? "alle" frekvenser.
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 amplitude-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 nullfrekvensen av det f?rste filteret m? v?re positivt, mens nullfrekvensen 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).
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.

Publisert 6. apr. 2021 08:19 - Sist endret 6. apr. 2021 08:19