clear;clc; %%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%% %基于各向异性热扩散方程的SAR图像分割 %王子路 %2007.2.2 %%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%% Img = imread('hanguangsar.bmp'); % Img = imread('d:\周琳图像分割\wangzi\1212121.bmp'); Img = rgb2gray(Img); % Img = imnoise(Img,'speckle',0.08); figure; imshow(Img); Img = double(Img); Img = exp(Img/17); %11 6 % Img = imadjust(Img,[0.45,0.7],[],1.2322); % imshow(uint8(Img)); SIZE = size(Img); SIZEX = SIZE(1,1); SIZEY = SIZE(1,2); % P_Class = 1/3*ones(SIZEX,SIZEY,3); % P_Class_tem = P_Class; Total = 0;time = 40; %设定初始阈值 % Threhold1 = 56;Threhold2 = 254;%不加后续处理16%加后续处理处理MSTAR为36%处理coast用56 Vector_1 = Img(Img<3);Vector_1 = Vector_1.^2;%2,180 7 3 Vector_2 = Img(Img>=3 & Img<=254);Vector_2 = Vector_2.^2; Vector_3 = Img(Img>254);Vector_3 = Vector_3.^2; Mark_1 = size(Vector_1);Mark_2 = size(Vector_2);Mark_3 = size(Vector_3); P_Class(:,:,1) = Mark_1(1,2)/(SIZEX*SIZEY)*ones(SIZEX,SIZEY); P_Class(:,:,2) = Mark_2(1,2)/(SIZEX*SIZEY)*ones(SIZEX,SIZEY); P_Class(:,:,3) = Mark_3(1,2)/(SIZEX*SIZEY)*ones(SIZEX,SIZEY); % Mark_1 = 1;Mark_2 = 1;Mark_3 = 1; % for i = 1:SIZEX % for j = 1:SIZEY % if Img(i,j) < Threhold1 % Flag(i,j) = 1; % Vector_1(Mark_1) = Img(i,j)^2; % Mark_1 = Mark_1+1; % elseif Img(i,j) >= Threhold1 & Img(i,j) <= Threhold2 % Flag(i,j) = 2; % Vector_2(Mark_2) = Img(i,j)^2; % Mark_2 = Mark_2+1; % elseif Img(i,j) > Threhold2; % Flag(i,j) = 3; % Vector_3(Mark_3) = Img(i,j)^2; % Mark_3 = Mark_3+1; % end % end % end % P_Class(:,:,1) = (Mark_1-1)/(SIZEX*SIZEY)*ones(SIZEX,SIZEY); % P_Class(:,:,2) = (Mark_2-1)/(SIZEX*SIZEY)*ones(SIZEX,SIZEY); % P_Class(:,:,3) = (Mark_3-1)/(SIZEX*SIZEY)*ones(SIZEX,SIZEY); % 进行迭代 Sigma(1) = mean(Vector_1); Sigma(2) = mean(Vector_2); Sigma(3) = mean(Vector_3); Mark_1 = 1;Mark_2 = 1;Mark_3 = 1; disp('准备开始初始分割......'); for m = 1:time Vector_1 = 0;Vector_2 = 0;Vector_3 = 0;k = 1:3; str = '初始分割已经完成'; str = strcat(str,num2str(m/time*100.0),'%'); disp(str); for i =1:SIZEX for j = 1:SIZEY Total = 0; for k = 1:3 Total = Total+1/(Sigma(k)+0.00000001)*exp(-1*Img(i,j)^2/(Sigma(k)+0.00000001))*P_Class(i,j,k); end % Total = 1./(Sigma(k)+0.00000001).*exp(-1*Img(i,j)^2./(Sigma(k)+0.00000001)).*reshape(P_Class(i,j,k),1,3); % Total = sum(Total); % P_Class_tem(i,j,k) = 1./(Sigma(k)+0.00000001).*exp(-1*Img(i,j)^2./(Sigma(k)+0.00000001)).*reshape(P_Class(i,j,k),1,3)/(Total+1/inf); for k = 1:3 % P_Class_tem(i,j,k) = 1/Sigma(k)*exp(-1*Img(i,j)^2/Sigma(k))*P_Class(:,:,k)/(Total+0.00000001); P_Class_tem(i,j,k) = 1/(Sigma(k)+0.00000001)*exp(-1*Img(i,j)^2/(Sigma(k)+0.00000001))*P_Class(i,j,k)/(Total+1/inf); end end end P_Class = P_Class_tem; [Value,Flag] = max(P_Class,[],3); Img1 = Img.^2; Vector_1 = Img1(Flag==1); Vector_2 = Img1(Flag==2); Vector_3 = Img1(Flag==3); % for i = 1:SIZEX % for j = 1:SIZEY % [Value,Flag(i,j)] = max([P_Class(i,j,1),P_Class(i,j,2),P_Class(i,j,3)]); % if Flag(i,j) == 1 % Vector_1(Mark_1) = Img(i,j)^2; % Mark_1 = Mark_1+1; % elseif Flag(i,j) == 2 % Vector_2(Mark_2) = Img(i,j)^2; % Mark_2 = Mark_2+1; % elseif Flag(i,j) == 3 % Vector_3(Mark_3) = Img(i,j)^2; % Mark_3 = Mark_3+1; % end % end % end Sigma(1) = std(Vector_1,1); Sigma(2) = std(Vector_2,1); Sigma(3) = std(Vector_3,1); Mark_1 = 1;Mark_2 = 1;Mark_3 = 1; end figure, imshow(Flag,[]); P_Class = abs(1000*P_Class); disp('初始分割已经完成'); disp('开始平滑后验概率矩阵......'); for i = 1:3 if i == 3 P_Class3(:,:,i) = RSRADFilter(P_Class(:,:,i),1.06); else P_Class3(:,:,i) = RSRADFilter(P_Class(:,:,i),2.36); end end for i = 1:SIZEX for j = 1:SIZEY [Value,Flag(i,j)] = max([P_Class3(i,j,1),P_Class3(i,j,2),P_Class3(i,j,3)]); end end Flag = WipeSpeckles(Flag); % map = [0,0,0.5;0,1,0;1,0,0]; figure; % image(Flag),colormap(map); imshow(Flag,[]);