-
Notifications
You must be signed in to change notification settings - Fork 1
/
Copy pathspatial_frequency_DOG.m
102 lines (77 loc) · 2.42 KB
/
spatial_frequency_DOG.m
1
2
3
4
5
6
7
8
9
10
11
12
13
14
15
16
17
18
19
20
21
22
23
24
25
26
27
28
29
30
31
32
33
34
35
36
37
38
39
40
41
42
43
44
45
46
47
48
49
50
51
52
53
54
55
56
57
58
59
60
61
62
63
64
65
66
67
68
69
70
71
72
73
74
75
76
77
78
79
80
81
82
83
84
85
86
87
88
89
90
91
92
93
94
95
96
97
98
clear all
close all
clc
% params of filter mentioned in the paper
sj=0.03;
q=2;
% create random image for test n size of 401x401
I=randn(401,401);
% calculate centered frequcy of given image
F=fftshift(fft2(I));
% spatial domain filter (filter size nxn)
n=51;
[x y]=meshgrid(-(n-1)/2:(n-1)/2);
r=sqrt(x.^2+y.^2);
% calculate first part of filter
p1=2*pi*sj^2*(q*exp(-2*pi^2*q^2*sj^2*r.^2));
% then second part
p2=2*pi*sj^2*(exp(-2*pi^2*sj^2*r.^2));
% make sure both parts are normalized
p1=p1/sum(p1(:));
p2=p2/sum(p2(:));
% and set the spatial filter
sDoG=p1-p2;
% frequency domain filter
% create meshgrid in the same resolution of given image
[U V]=meshgrid(-fix(size(F,2)/2):fix(size(F,2)/2),-fix(size(F,1)/2):fix(size(F,1)/2));
% note the last index must be -0.5 to 0.5 hence we must use normalized
% frequency
U=0.5*U/max(U(:));
V=0.5*V/max(V(:));
% calculate ro value
p=sqrt(U.^2+V.^2);
fDoG=exp(-0.5*(p/(q*sj)).^2)-exp(-0.5*(p/sj).^2);
% filter image in frequency domain
FD=F.*fDoG;
% transform the filtered image in to spatial domain
fI=ifft2(ifftshift(FD));
% filter image in spatial domain
oI=conv2(I,sDoG,'same');
figure;colormap(gray);imagesc(I);
figure;mesh(U,V,abs(F));
xlabel('u');
ylabel('v');
zlabel('Frequency component''s magnitude');
figure;
subplot(1,3,1);imagesc(fI);title('frequency filter result');
subplot(1,3,2);imagesc(oI);title('spatial filter result');
subplot(1,3,3);imagesc(fI-oI);title('differences');
figure;
subplot(1,3,1);mesh(U,V,abs(fftshift(fft2(fI))));title('frequency of frequency filter result');
xlabel('u');
ylabel('v');
zlabel('Frequency component''s magnitude');
subplot(1,3,2);mesh(U,V,abs(fftshift(fft2(oI))));title('frequency of spatial filter result');
xlabel('u');
ylabel('v');
zlabel('Frequency component''s magnitude');
zlim([0 600]);
subplot(1,3,3);mesh(U,V,abs(fftshift(fft2(fI)))-abs(fftshift(fft2(oI))));title('differences');
% impulse response
[H,f1,f2] = freqz2(sDoG);
figure;
subplot(1,3,1);mesh(f1,f2,imresize(fDoG,[size(H,1) size(H,2)]));
title('frequency filter''s impulse response');
xlabel('u');
ylabel('v');
zlabel('Frequency component''s magnitude');
subplot(1,3,2);mesh(f1,f2,H);title('spatial filter''s impulse response');
zlim([0 0.5]);
xlabel('u');
ylabel('v');
zlabel('Frequency component''s magnitude');
subplot(1,3,3);mesh(f1,f2,H-imresize(fDoG,[size(H,1) size(H,2)]));
title('impulse responses differences');
xlabel('u');
ylabel('v');
zlabel('Frequency component''s magnitude');