jtc.m

在matlab中

%JTC
%% 读入图像
clc;clear;close all;
N=1024;
ob_1=imread('xp.bmp');
ob_2=imread('xp.bmp');
ob_1=double(rgb2gray(ob_1));
ob_2=double(rgb2gray(ob_2));
% figure,subplot(121),imshow(ob_1,[]);
% subplot(122),imshow(ob_2,[]);
%% 把输入图像置于同一平面；
[m10,n10]=size(ob_1);
[m20,n20]=size(ob_2);
m1=floor(m10/2);
m2=floor(m20/2);
n1=floor(n10/2);
n2=floor(n20/2);
input=zeros(N);
input(N/2-m1+1:N/2+m10-m1,N/4-n1+1:N/4+n10-n1)=ob_1;
input(N/2-m2+1:N/2+m20-m2,3*N/4-n2+1:3*N/4+n20-n2)=ob_2;
figure,imshow(input,[]);
%%
input_edge=edge(input,'canny');%滤波
% figure,imshow(input_edge,[]);
%% %相关计算
input_fft=fft2(input_edge);
input_fftshift=fftshift(input_fft);
% figure,imshow(log(1+abs(input_fftshift)),[]);
input_abs=input_fftshift.*conj(input_fftshift);
output_ifft=ifft2(input_abs);
output=fftshift(output_ifft);
figure;
mesh(abs(output));
figure,imshow(output,[]);
%%
%高通滤波
r=N/4;
for i=1:N;
    for j=1:N;
       if((i-N/2).^2+(j-N/2).^2<=r)
           h(i,j)=0;
       else
           h(i,j)=1;
       end
    end
end
output_filt=output.*h;
figure,mesh(output_filt);