clc
clear all
%% 输入图像，这种输入方式可以同时实现横向与纵向的梯度检测与叠加，无需二次输入，但一次操作只能做两组数据的差分
I=imread('4.1.06.tiff');
I=rgb2gray(I);
I=im2double(I);
subplot(1,4,1);
imshow(I);
title('Input');
N=4%倍增倍数
% syms A B C D E F G H I J K L
% I=[A B C D E F; G H I J K L]%输入图像
II=zeros(length(I),length(I)*N*2);
II2=zeros(length(I),length(I)*N*2);
%% 对输入图像矩阵做顺序变换，使其适配梯度矩阵的顺序
count=0;%计数
 for i=1:1:(length(I)-1)*N*2 %按列遍历
    if mod(i,2*N)==1
        II(:,i)=I(:,mod(i,8)+floor(i/8));
    elseif mod(i,2*N)==2
        II(:,i)=I(:,mod(i,8)+floor(i/8));
    elseif mod(i,2*N)==3
        II(:,i)=I(:,mod(i,8)+floor(i/8)-2);
    elseif mod(i,2*N)==4
%                II(:,i)=I(:,mod(i,8)+count-3);
        II(1:end-1,i)=I(2:end,mod(i,8)+floor(i/8)-3);
        II(end,i)=I(1,floor(i/8)+1);
    elseif mod(i,2*N)==5
        II(:,i)=I(:,mod(i,8)+floor(i/8)-4);
    elseif mod(i,2*N)==6
        II(:,i)=I(:,mod(i,8)+floor(i/8)-4);
    elseif mod(i,2*N)==7
        II(:,i)=I(:,mod(i,8)+floor(i/8)-6);
    elseif mod(i,2*N)==0
        II(1:end-1,i)=-I(2:end,mod(i,8)+floor(i/8));
        II(end,i)=-I(1,floor(i/8));

    end
    end
    II(:,(length(I)-1)*N*2+1)=I(:,length(I));
    II(:,(length(I)-1)*N*2+2)=I(:,1);
    II(:,(length(I)-1)*N*2+3)=I(:,length(I));
    II(1:end-1,(length(I)-1)*N*2+4)=I(2:end,length(I));
    II(end,(length(I)-1)*N*2+4)=I(1,length(I));
    II(:,(length(I)-1)*N*2+5)=I(:,length(I));
    II(:,(length(I)-1)*N*2+6)=I(:,1);
    II(:,(length(I)-1)*N*2+7)=I(:,length(I));
    II(1:end-1,(length(I)-1)*N*2+8)=-I(2:end,length(I));
    II(end,(length(I)-1)*N*2+8)=-I(1,length(I));



%% 对输入图像做倍增变换，使其适配梯度矩阵的衰减和相移
% for i=1:1:length(I)*N*2 %遍历变换后的输入图像II
%     % for i=1:1:6*N
%     if mod(i,2*N)==1
%         II2(:,i)=II(:,i)/0.54/exp(-j*96/180*pi);
%     elseif mod(i,2*N)==2
%         II2(:,i)=II(:,i)/0.54/exp(-j*96/180*pi);
%     elseif mod(i,2*N)==3
%         II2(:,i)=II(:,i)/0.54/exp(-j*96/180*pi);
%     elseif mod(i,2*N)==4
%         II2(:,i)=II(:,i)/0.54/exp(-j*96/180*pi);
%     elseif mod(i,2*N)==5
%         II2(:,i)=II(:,i)/0.54/exp(-j*96/180*pi);
%     elseif mod(i,2*N)==6
%         II2(:,i)=II(:,i)/0.54/exp(-j*96/180*pi);
%     elseif mod(i,2*N)==7
%         II2(:,i)=-II(:,i)/0.54/exp(-j*96/180*pi);
%     elseif mod(i,2*N)==0
%         II2(:,i)=-II(:,i)/0.54/exp(-j*96/180*pi);
%     end
% end
for i=1:1:length(I)*N*2 %遍历变换后的输入图像II
    % for i=1:1:6*N
    if mod(i,2*N)==1
        II2(:,i)=II(:,i)/0.54/exp(-j*96/180*pi);
    elseif mod(i,2*N)==2
        II2(:,i)=II(:,i)/0.54/exp(-j*96/180*pi);
    elseif mod(i,2*N)==3
        II2(:,i)=II(:,i)/0.54/exp(-j*96/180*pi);
    elseif mod(i,2*N)==4
        II2(:,i)=II(:,i)/0.54/exp(-j*96/180*pi);
    elseif mod(i,2*N)==5
        II2(:,i)=II(:,i)/0.54/exp(-j*96/180*pi);
    elseif mod(i,2*N)==6
        II2(:,i)=II(:,i)/0.54/exp(-j*96/180*pi);
    elseif mod(i,2*N)==7
        II2(:,i)=-II(:,i)/0.54/exp(-j*96/180*pi);
    elseif mod(i,2*N)==0
        II2(:,i)=-II(:,i)/0.54/exp(-j*96/180*pi);
    end
end
%% 8入8出板子所实现梯度矩阵，只涉及横纵向，不涉及交叉梯度
D=zeros(length(I)*N*2,length(I));
for i=1:1:length(I) %矩阵按照列遍历
    D((i-1)*N*2+1,i)=0.54*exp(-j*96/180*pi);
    D((i-1)*N*2+2,i)=-0.57*exp(-j*80/180*pi);
    D((i-1)*N*2+3,i)=0.52*exp(-j*90/180*pi);
    D((i-1)*N*2+4,i)=-0.52*exp(-j*78/180*pi);
    D((i-1)*N*2+5,i)=0.54*exp(-j*96/180*pi);
    D((i-1)*N*2+6,i)=-0.57*exp(-j*80/180*pi);
    D((i-1)*N*2+7,i)=-0.67*exp(-j*80/180*pi);
    D((i-1)*N*2+8,i)=-0.67*exp(-j*80/180*pi);
end
%% 输出图像
Iout=II2*D;
I_edge_1=Iout;%是横向纵向处理之和，边缘检测
I_Sharpen=Iout+I;%边缘检测叠加原始图像，是图像锐化
subplot(1,4,2);
imshow(I_edge_1);
title('Gradient');
% subplot(1,3,3);
% imshow(I_Sharpen);
% title('图像锐化');
% sgtitle('预处理幅度，简单梯度算子边缘检测/增强与图像锐化')
% suptitle('简单梯度算子边缘检测/增强与图像锐化')

%% 输入图像，这种输入方式可以同时实现横向与纵向的梯度检测与叠加，无需二次输入，但一次操作只能做两组数据的差分
I=imread('4.1.06.tiff');
I=rgb2gray(I);
I=im2double(I);
% figure
% subplot(1,4,1);
% imshow(I);
% title('原图像');
N=4%倍增倍数
% syms A B C D E F G H I J K L
% I=[A B C D E F; G H I J K L]%输入图像
II=zeros(length(I),length(I)*N*2);
II2=zeros(length(I),length(I)*N*2);
%% 对输入图像矩阵做顺序变换，使其适配梯度矩阵的顺序
count=0;%计数
for i=1:1:(length(I)-1)*N*2 %按列遍历
    if mod(i,2*N)==1
        II(:,i)=I(:,mod(i,8)+floor(i/8));
    elseif mod(i,2*N)==2
        %         II(:,i)=I(:,mod(i,8)+floor(i/8));
        II(1:end-1,i)=I(2:end,mod(i,8)+floor(i/8));
        II(end,i)=I(1,floor(i/8)+1+1);
    elseif mod(i,2*N)==3
        %         II(:,i)=I(:,mod(i,8)+floor(i/8)-2);
        II(:,i)=I(:,mod(i,8)+floor(i/8)-1);
    elseif mod(i,2*N)==4
        %                II(:,i)=I(:,mod(i,8)+count-3);
        II(1:end-1,i)=I(2:end,mod(i,8)+floor(i/8)-3);
        II(end,i)=I(1,floor(i/8)+1);
    elseif mod(i,2*N)==5
        II(:,i)=I(:,mod(i,8)+floor(i/8)-4);
    elseif mod(i,2*N)==6
        II(1:end-1,i)=I(2:end,mod(i,8)+floor(i/8)-4);
        II(end,i)=I(1,floor(i/8)+1+1);
    elseif mod(i,2*N)==7
        II(:,i)=I(:,mod(i,8)+floor(i/8)-5);
    elseif mod(i,2*N)==0
        II(1:end-1,i)=-I(2:end,mod(i,8)+floor(i/8));
        II(end,i)=-I(1,floor(i/8));
        
    end
end
II(:,(length(I)-1)*N*2+1)=I(:,length(I));

II(1:end-1,(length(I)-1)*N*2+2)=I(2:end,1);
II(end,(length(I)-1)*N*2+2)=I(1,1);

II(:,(length(I)-1)*N*2+3)=I(:,1);

II(1:end-1,(length(I)-1)*N*2+4)=I(2:end,length(I));
II(end,(length(I)-1)*N*2+4)=I(1,length(I));

II(:,(length(I)-1)*N*2+5)=I(:,length(I));

II(1:end-1,(length(I)-1)*N*2+6)=I(2:end,1);
II(end,(length(I)-1)*N*2+6)=I(1,1);

II(:,(length(I)-1)*N*2+7)=I(:,1);



%% 对输入图像做倍增变换，使其适配梯度矩阵的衰减和相移
for i=1:1:length(I)*N*2 %遍历变换后的输入图像II
    % for i=1:1:6*N
    if mod(i,2*N)==1
        II2(:,i)=II(:,i)/0.54/exp(-j*96/180*pi);
    elseif mod(i,2*N)==2
        II2(:,i)=II(:,i)/0.54/exp(-1i*96/180*pi);
    elseif mod(i,2*N)==3
        II2(:,i)=II(:,i)/0.54/exp(-j*96/180*pi);
    elseif mod(i,2*N)==4
        II2(:,i)=II(:,i)/0.54/exp(-j*96/180*pi);
    elseif mod(i,2*N)==5
        II2(:,i)=II(:,i)/0.54/exp(-j*96/180*pi);
    elseif mod(i,2*N)==6
        II2(:,i)=II(:,i)/0.54/exp(-j*96/180*pi);
    elseif mod(i,2*N)==7
        II2(:,i)=-II(:,i)/0.54/exp(-j*96/180*pi);
    elseif mod(i,2*N)==0
        II2(:,i)=-II(:,i)/0.54/exp(-j*96/180*pi);
    end
end
%     II(:,length(I)*2)=II(:,length(I)*2)/0.93/exp(j*109/180*pi);
%     II(:,1)=II(:,1)/0.75/exp(j*95/180*pi);
%% 8入8出板子所实现梯度矩阵，交叉梯度
D=zeros(length(I)*N*2,length(I));
for i=1:1:length(I) %矩阵按照列遍历
    D((i-1)*N*2+1,i)=0.54*exp(-j*96/180*pi);
    D((i-1)*N*2+2,i)=-0.57*exp(-j*80/180*pi);
    D((i-1)*N*2+3,i)=0.52*exp(-j*90/180*pi);
    D((i-1)*N*2+4,i)=-0.52*exp(-j*78/180*pi);
    D((i-1)*N*2+5,i)=0.54*exp(-j*96/180*pi);
    D((i-1)*N*2+6,i)=-0.57*exp(-j*80/180*pi);
    D((i-1)*N*2+7,i)=-0.67*exp(-j*80/180*pi);
    D((i-1)*N*2+8,i)=-0.67*exp(-j*80/180*pi);
end
%% 输出图像
Iout=II2*D;
I_edge_2=Iout;%同时是横向纵向处理之和，边缘检测
I_Sharpen=Iout+I;%边缘检测叠加原始图像，是图像锐化
subplot(1,4,3);
imshow(I_edge_2);
title('Robert');
% subplot(1,3,3);
% imshow(I_Sharpen);
% title('图像锐化');
% sgtitle('预处理幅度，Robert梯度算子边缘检测/增强与图像锐化')
% suptitle('Robert算子边缘检测/增强与图像锐化')

%% 输入图像，首先对图像做横向的梯度处理
% I=imread('ccc.png');
I=imread('4.1.06.tiff');
I=rgb2gray(I);
I=im2double(I);
% subplot(1,3,1);
% imshow(I);
% title('原图像');
N=3
I2=I';%竖向输入
II=zeros(length(I),length(I)*N*2);
II2=zeros(length(I),length(I)*N*2);
%% 对输入图像矩阵做顺序变换，使其适配梯度矩阵的顺序

for i=1:1:length(I)-(N-1) %按列遍历输入图像I
    II(:,(2*N)*(i-1)+1)=I(:,i);
    II(:,(2*N)*(i-1)+2)=I(:,i+1);
    II(:,(2*N)*(i-1)+3)=I(:,i+2);

    II(1:end-(N-1),(2*N)*(i-1)+4)=I(N:end,i);
    II(end-1:end,(2*N)*(i-1)+4)=I(1:2,i);

    II(1:end-(N-1),(2*N)*(i-1)+5)=I(N:end,i+1);
    II(end-1:end,(2*N)*(i-1)+5)=I(1:2,i+1);

    II(1:end-(N-1),(2*N)*(i-1)+6)=I(N:end,i+2);
    II(end-1:end,(2*N)*(i-1)+6)=I(1:2,i+2);
end
II(:,(2*N)*(length(I)-2)+1)=I(:,(length(I)-1));
II(:,(2*N)*(length(I)-2)+2)=I(:,(length(I)));
II(:,(2*N)*(length(I)-2)+3)=I(:,1);

II(1:end-(N-1),(2*N)*(length(I)-2)+4)=I(N:end,length(I)-1);
II(end-1:end,(2*N)*(length(I)-2)+4)=I(1:2,length(I)-1);

II(1:end-(N-1),(2*N)*(length(I)-2)+5)=I(N:end,length(I));
II(end-1:end,(2*N)*(length(I)-2)+5)=I(1:2,length(I));

II(1:end-(N-1),(2*N)*(length(I)-2)+6)=I(N:end,1);
II(end-1:end,(2*N)*(length(I)-2)+6)=I(1:2,1);

II(:,(2*N)*(length(I)-1)+1)=I(:,(length(I)));
II(:,(2*N)*(length(I)-1)+2)=I(:,1);
II(:,(2*N)*(length(I)-1)+3)=I(:,2);

II(1:end-(N-1),(2*N)*(length(I)-1)+4)=I(N:end,length(I));
II(end-1:end,(2*N)*(length(I)-1)+4)=I(1:2,length(I));

II(1:end-(N-1),(2*N)*(length(I)-1)+5)=I(N:end,1);
II(end-1:end,(2*N)*(length(I)-1)+5)=I(1:2,1);

II(1:end-(N-1),(2*N)*(length(I)-1)+6)=I(N:end,2);
II(end-1:end,(2*N)*(length(I)-1)+6)=I(1:2,2);
%% 对输入图像做倍增变换，使其适配梯度矩阵的衰减和相移
for i=1:1:length(I2)*N*2 %遍历变换后的输入图像II
    if mod(i,2*N)==1
        II1(:,i)=II(:,i)/0.54/exp(-j*96/180*pi);
    elseif mod(i,2*N)==2
        II1(:,i)=II(:,i)/0.54/exp(-j*96/180*pi);
    elseif mod(i,2*N)==3
        II1(:,i)=II(:,i)/0.54/exp(-j*96/180*pi);
    elseif mod(i,2*N)==4
        II1(:,i)=II(:,i)/0.54/exp(-j*96/180*pi);
    elseif mod(i,2*N)==5
        II1(:,i)=II(:,i)/0.54/exp(-j*96/180*pi);
    elseif mod(i,2*N)==0
        II1(:,i)=II(:,i)/0.54/exp(-j*96/180*pi);

    end
end
%% 竖向输入图像做相同处理
%% 首先是顺序变换
for i=1:1:length(I2)-(N-1) %按列遍历输入图像I
    II2(:,(2*N)*(i-1)+1)=I2(:,i);
    II2(:,(2*N)*(i-1)+2)=I2(:,i+1);
    II2(:,(2*N)*(i-1)+3)=I2(:,i+2);

    II2(1:end-(N-1),(2*N)*(i-1)+4)=I2(N:end,i);
    II2(end-1:end,(2*N)*(i-1)+4)=I2(1:2,i);

    II2(1:end-(N-1),(2*N)*(i-1)+5)=I2(N:end,i+1);
    II2(end-1:end,(2*N)*(i-1)+5)=I2(1:2,i+1);

    II2(1:end-(N-1),(2*N)*(i-1)+6)=I2(N:end,i+2);
    II2(end-1:end,(2*N)*(i-1)+6)=I2(1:2,i+2);
end
II2(:,(2*N)*(length(I2)-2)+1)=I2(:,(length(I2)-1));
II2(:,(2*N)*(length(I2)-2)+2)=I2(:,(length(I2)));
II2(:,(2*N)*(length(I2)-2)+3)=I2(:,1);

II2(1:end-(N-1),(2*N)*(length(I2)-2)+4)=I2(N:end,length(I2)-1);
II2(end-1:end,(2*N)*(length(I2)-2)+4)=I2(1:2,length(I2)-1);

II2(1:end-(N-1),(2*N)*(length(I2)-2)+5)=I2(N:end,length(I2));
II2(end-1:end,(2*N)*(length(I2)-2)+5)=I2(1:2,length(I2));

II2(1:end-(N-1),(2*N)*(length(I2)-2)+6)=I2(N:end,1);
II2(end-1:end,(2*N)*(length(I2)-2)+6)=I2(1:2,1);

II2(:,(2*N)*(length(I2)-1)+1)=I2(:,(length(I2)));
II2(:,(2*N)*(length(I2)-1)+2)=I2(:,1);
II2(:,(2*N)*(length(I2)-1)+3)=I2(:,2);

II2(1:end-(N-1),(2*N)*(length(I2)-1)+4)=I2(N:end,length(I2));
II2(end-1:end,(2*N)*(length(I2)-1)+4)=I2(1:2,length(I2));

II2(1:end-(N-1),(2*N)*(length(I2)-1)+5)=I2(N:end,1);
II2(end-1:end,(2*N)*(length(I2)-1)+5)=I2(1:2,1);

II2(1:end-(N-1),(2*N)*(length(I2)-1)+6)=I2(N:end,2);
II2(end-1:end,(2*N)*(length(I2)-1)+6)=I2(1:2,2);
%% 对输入图像做倍增变换
for i=1:1:length(I2)*N*2 %遍历变换后的输入图像II
    if mod(i,2*N)==1
        II22(:,i)=II2(:,i)/0.54/exp(-j*96/180*pi);
    elseif mod(i,2*N)==2
        II22(:,i)=II2(:,i)/0.54/exp(-j*96/180*pi);
    elseif mod(i,2*N)==3
        II22(:,i)=II2(:,i)/0.54/exp(-j*96/180*pi);
    elseif mod(i,2*N)==4
        II22(:,i)=II2(:,i)/0.54/exp(-j*96/180*pi);
    elseif mod(i,2*N)==5
        II22(:,i)=II2(:,i)/0.54/exp(-j*96/180*pi);
    elseif mod(i,2*N)==0
        II22(:,i)=II2(:,i)/0.54/exp(-j*96/180*pi);

    end
end
%% 8入8出板子所实现梯度矩阵，只涉及横纵向，不涉及交叉梯度
D=zeros(length(I)*N*2,length(I));
for i=1:1:length(I) %鐭╅樀鎸夌収鍒楅亶鍘?
D((i-1)*N*2+1,i)=0.54*exp(-j*96/180*pi);
    D((i-1)*N*2+2,i)=0.52*exp(-j*90/180*pi);
    D((i-1)*N*2+3,i)=0.54*exp(-j*96/180*pi);
    D((i-1)*N*2+4,i)=-0.57*exp(-j*80/180*pi);
    D((i-1)*N*2+5,i)=-0.52*exp(-j*78/180*pi);
    D((i-1)*N*2+6,i)=-0.57*exp(-j*80/180*pi);
end
%% 输出图像
Iout=II1*D;%横向处理
Iout2=II22*D;%纵向处理
Iout2=Iout2';
I_edge_3=Iout+Iout2;%横向纵向处理之和，边缘检测
I_Sharpen=Iout+I+Iout2;%边缘检测叠加原始图像，是图像锐化
subplot(1,4,4);
imshow(I_edge_3);
title('Prewitt');
% subplot(1,3,3);
% imshow(I_Sharpen);
% title('图像锐化');
% % suptitle('Prewitt梯度算子边缘检测/增强与图像锐化')
% sgtitle('只预处理传输损耗，而不处理相位，Prewitt梯度算子边缘检测/增强与图像锐化')
kkk1=abs( I_edge_1);%basic灰度
kkk2=abs( I_edge_2);%robert灰度
kkk3=abs( I_edge_3);%prewitt灰度
% T1 = graythresh( I_edge_1)
% T2 = graythresh( I_edge_2)
% T3= graythresh( I_edge_3)
% F1=find( I_edge_1<T1)
% F2=find( I_edge_2<T2)
% F3=find( I_edge_3<T3)
B=abs((I_edge_1(23,1)+I_edge_1(24,1))/(I_edge_1(22,1)+I_edge_1(21,1)))
R=abs((I_edge_2(23,1)+I_edge_2(22,1))/(I_edge_2(20,1)+I_edge_2(21,1)))
P=abs((I_edge_3(23,1)+I_edge_3(22,1))/(I_edge_3(20,1)+I_edge_3(21,1)))
% B=(kkk1(23,1)+kkk1(24,1))/(kkk1(22,1)+kkk1(21,1))
% R=(kkk2(23,1)+kkk2(22,1))/(kkk2(20,1)+kkk2(21,1))
% P=(kkk3(23,1)+kkk3(22,1))/(kkk3(20,1)+kkk3(21,1))
I_edge_basic=I_edge_1;
I_edge_Robert=I_edge_2;
I_edge_Prewitt=I_edge_3;