Re: [問題] 影像處理DFT已回收
function DFTg3(filename)
% [版本一] 內建code的寫法
% f1 = imread('t.jpg');
% af = fftshift(fftn(double(f1))); fftshow1(af);
% bf = ifftn(double(af)); fftshow1(bf);
% [版本二]，跑完所需時間大約是內建code的 2~15 倍
% fft2a.m 跑完所需時間大約是 fftn.m 的 27~55 倍
% 讀取影像
% af = fftn(double(f1)); ff = af - F; max(ff(:)); fb = bf1 - bf; max(fb(:));
% format short g; format compact; f1 = imread('t.jpg');
% clear all; tic; DFTg3('t.jpg'); toc;
if ischar(filename)
f1 = imread(filename);
else
f1 = filename;
end
F = fft2a(double(f1),1);
% 平移轉換，參考 fftshift.m 寫出來的。
% 平移轉換並不是單純的 F(x,y) * (-1)^(x+y)，
%"平移轉換"是指改變 F(x,y) 的位置。
% FFTSHIFT(X) swaps the first and third quadrants
% and the second and fourth quadrants.
nd = ndims(F); idx = cell(1,nd);
for k = 1 : nd
m = size(F, k); p = ceil(m/2); idx{k} = [p+1:m 1:p];
end
Fi = F(idx{:}); % fftshow1(Fi);
bf = fft2a(double(Fi),-1); fftshow1(bf);
function F = fft2a(f,id)
% id = 1 or -1; for DFT transform, id = 1; for IDFT transform, id = -1;
[M,N] = size(f); f = double(f); F = zeros(M,N); c = sqrt(-1);
if M ~= N
a = exp(-id*2*pi*c/M); b = exp(-id*2*pi*c/N);
d = zeros(M,M-1); d(1,:) = 1;
for i = 2 : M; d(i,1) = d(i-1,1) * a; end;
for i = 2 : M-1; d(:,i) = d(:,1) .* d(:,i-1); end;
D = zeros(N,N-1); D(1,:) = 1;
for i = 2 : N; D(i,1) = D(i-1,1) * b; end;
for i = 2 : N-1; D(:,i) = D(:,1) .* D(:,i-1); end;
else
% 使用 * 的次數 = 3 + M-1 + M*(M-2) = M^2 - M + 1;
a = exp(-id*2*pi*c/M);
d = zeros(M,M-1); d(1,:) = 1;
for i = 2 : M; d(i,1) = d(i-1,1) * a; end;
for i = 2 : M-1; d(:,i) = d(:,1) .* d(:,i-1); end;
D = d;
end
% 使用 * 的次數 = (M-1)*M + (N-1)*(M*N + (M-1)*M);
F(1,1) = sum(f(:)); s1 = sum(f,2)';
for u = 2 : M; F(u,1) = s1 * d(:,u-1); end
for v = 2 : N
s = [f * D(:,v-1)]'; F(1,v) = sum(s);
for u = 2 : M; F(u,v) = s * d(:,u-1); end
end
if id == -1;
MN = M*N; F = F/MN;
end
A = zeros(M,N); A(1,:) = F(1,:);
for i = 2 : M; A(i,:) = F(M-i+2,:); end;
F = conj(A);
function fftshow1(f,type)
if nargin < 2; type = 'log'; end;
if type == 'log'
f1 = log(1+abs(f)); fm = max(f1(:)); imshow(double(f1/fm));
elseif type == 'abs'
fa = abs(f); fm = max(fa(:)); imshow(fa/fm);
else
errow('TYPE must be abs or log.');
end
※ 編輯: buddin 來自: 140.112.7.59 (05/09 03:26)
※ 編輯: buddin 來自: 140.112.7.59 (05/09 03:46)
討論串 (同標題文章)
完整討論串 (本文為第 3 之 3 篇)：
0
3
MATLAB 近期熱門文章
PTT數位生活區 即時熱門文章
