clear all;clc;close all;
X=imread('Prediction2.png');
K = Bianhuan(X);
image_size = size(X);
demen = numel(image_size);
if demen == 3
X = rgb2gray(X);
end
[M, N] = size(X);
X = im2double(X);
X = X*255;
figure(1)
imshow(X);
%X = medfilt2(X, [3, 3]); %图像预处理
C0 = max(max(X)); %求矩阵的最大值
delta = 0.45;
theta = edge(X, 'canny', [0.3, 0.6]);%返回一个与X相同大小的二值化图像BW
theta = im2double(theta);
sigma2 = 0.3;
rho = 1; %曲线饱满值
Rho = [rho^2 0; 0 rho^(-2)];
M1 = zeros(3, 3); %输入权重矩阵
H = zeros(M+2, N+2);
Y = zeros(M,N);
Y1 = ones(M,N);
T = zeros(M,N);
%求频谱残差
sigma1 = 8;
A = abs(fft2(X));%二维快速傅里叶变换 频谱残差法
P = angle(fft2(X)); %求相位角
C = log(A);
C1 = imfilter(C,fspecial('average',3)); %用局部平均滤波器对其进行平滑获得log谱的大致形状
R = C - C1; %谱残差
g = fspecial('gaussian',[M,N],sigma1);
S = double(g).*(ifft2(exp(P+R))).^2; %谱残差能够描述一幅图像中的异常部分。将谱残差和相位进行二维离散傅里叶反变换。得到重构图
S = abs(S);
%通过频谱残差来计算连接强度系数
Beta = S/max(S(:));
Beta = double(Beta);
W = [2^(-0.5) 1 2^(-0.5);
1 0 1;
2^(-0.5) 1 2^(-0.5)]; %网络中各神经元的连接权系数矩阵
n1 = 1;
E = zeros(1, 100);
pi = 3.14;
for i = 2:M+1
for j = 2:N+1
Z = [K(i-1, j-1) K(i-1, j) K(i-1, j+1);
K(i, j-1) K(i, j) K(i, j+1);
K(i+1, j-1) K(i+1, j) K(i+1, j+1)];
Htheta = [cos(theta(i-1, j-1)) sin(theta(i-1, j-1));
-sin(theta(i-1, j-1)) cos(theta(i-1, j-1))]; %旋转矩阵
M2 = 2*C0*pi*sigma2^2;
for k = 1:3
for l = 1:3
if k ~= 2 && l ~= 2
x = [1;1];
elseif k ~= 2 && l == 2
x = [1;0];
elseif k == 2 && l ~= 2
x = [0;1];
elseif k == 2 && l == 2
x = [0;0];
end
M4 = x'*Htheta'*Rho*Htheta*x;
M3 = exp(-M4/2*sigma2^2);
M1(k,l) = M3/M2;%%%%%%高斯内核输入权重矩阵
end
end
F(i-1, j-1) = double(sum(sum(Z.*M1))); %由权重矩阵计算得出
end
end
while(sum(sum(Y1 == Y)) ~= M*N) %改进的脉冲耦合神经网络部分关键部分
%检验分割结果是否变化，没有变化则输出
Y1 = Y;
for i = 2:M+1
for j = 2:N+1
V = [H(i-1, j-1) H(i-1, j) H(i-1, j+1);
H(i, j-1) H(i, j) H(i, j+1);
H(i+1, j-1) H(i+1, j) H(i+1, j+1)];
%利用公式依次计算出耦合输出L,反馈输出F，调制域输出U,动态阈值E以及脉冲输出序列Y
L(i-1, j-1) = double(sum(sum(V.*W)));
U(i-1, j-1) = F(i-1, j-1)*(1 + Beta(i-1, j-1)*L(i-1, j-1));
if U(i-1, j-1) >= delta*E(n1) %判断输出脉冲Y
Y(i-1,j-1) = 1;
else
Y(i-1,j-1) = 0;
end
end
end
if n1 ~= 1
H = Bianhuan(Y);%进行边缘检测
end
n1 = n1 + 1;
e(n1) = sum(sum(X.*Y));
y(n1) = sum(sum(Y));
E(n1) = e(n1)/y(n1);
figure(n1)
imshow(Y);
end
