%%
close all force
clear all
clc
%%
im1=imread('1.png');
ReIm=im1;
im2=rgb2gray(im1);
%% 转化为灰度 大津算法计算阈值分割图像
level = graythresh(im2);
level=level*1.2;
bw = im2bw(im2, level);
imtool(bw,[]);
bw1=bw;
for ii=1+1:size(bw1,1)-1
for jj=1+1:size(bw1,2)-1
if bw1(ii,jj)==1
if bw(ii+1,jj)+bw(ii,jj+1)+bw(ii-1,jj)+bw(ii,jj-1)+bw(ii+1,jj+1)+bw(ii-1,jj-1)<6
bw1(ii,jj)=0;
end
end
end
end
imtool(bw1,[]);
im2=double(im2);
im2(293,338)=im2(293,338)+30;
im2(272,521)=100;
imtool(im2,[]);
mask=~bw1;
imtool(mask,[]);
mask1=zeros(size(mask));
stt=120;
mask1(find(im2>stt))=1;
mask2=mask.*mask1;
%% 根据高斯分布的对称性进一步计算mask
mask3=zeros(size(mask));
th=2;
for ii=1+th:size(bw1,1)-th
for jj=1+th:size(bw1,2)-th
temp1=mask(ii,jj-th:jj+th);
temp2=mask(ii-th:ii+th,jj);
if sum(temp1)==length(temp1)
left1=im2(ii,jj-th:jj);
right1=im2(ii,jj:jj+th);
left2=diff(left1);
right2=diff(right1);
flag1=find(left2<0);
flag2=find(right2>0);
if isempty(flag1)&&isempty(flag2)
mask3(ii,jj)=1;
end
end
end
end
%%
mask4=zeros(size(mask));
for ii=1+th:size(bw1,1)-th
for jj=1+th:size(bw1,2)-th
temp1=mask(ii,jj-th:jj+th);
temp2=mask(ii-th:ii+th,jj);
if sum(temp2)==length(temp2)
left1=im2(ii-th:ii,jj);
right1=im2(ii:ii+th,jj);
left2=diff(left1);
right2=diff(right1);
flag1=find(left2<0);
flag2=find(right2>0);
if isempty(flag1)&&isempty(flag2)
mask4(ii,jj)=1;
end
end
end
end
%%
mask5=(mask2.*mask3).*mask4;
th=10;
PosLis=zeros(0,2);
VaLis=zeros(0,11);
for ii=1+th:size(mask,1)-th
for jj=1+th:size(mask,2)-th
if mask5(ii,jj)>0
lis=[];
for kk=1:th
temp1=im2(ii-kk:ii+kk,jj-kk:jj+kk);
temp2=im2(ii-kk+1:ii+kk-1,jj-kk+1:jj+kk-1);
temp3=sum(sum(temp1))-sum(sum(temp2));
temp4=temp3/(numel(temp1)-numel(temp2));
lis=[lis, temp4];
end
lis=[im2(ii,jj), lis];
if lis(end)<91
PosLis=[PosLis; [ii,jj]];
VaLis=[VaLis; lis];
end
end
end
end
%% 高斯函数拟合
ReLis=zeros(0,2);
for ii=1:size(VaLis,1)
VaLis1=VaLis(ii,:);
VaLis1=VaLis1(1:end);
% VaLis1=VaLis1-VaLis(end)+1;
xx=1:length(VaLis1);
% xx=xx-1;
ff=@(a,x)a(1)*exp(-x.^2/a(2))+a(3);
[aa, res]=lsqcurvefit(ff, [1,1,0], xx, VaLis1);
if aa(2)>2
ReLis=[ReLis; PosLis(ii,:)];
figure;
plot(xx,VaLis1,'--o');
title(num2str([PosLis(ii,1),PosLis(ii,2)]));
hold on;
plot(xx,ff(aa,xx),'-r*');
res=res/sum(VaLis1);
end
end
%% plot
th=5;
for kk=1:size(ReLis,1)
ii=ReLis(kk,1);
jj=ReLis(kk,2);
for mm=ii-th:ii+th
ReIm(mm,jj-th,1)=255;
ReIm(mm,jj-th,2)=0;
ReIm(mm,jj-th,3)=0;
ReIm(mm,jj+th,1)=255;
ReIm(mm,jj+th,2)=0;
ReIm(mm,jj+th,3)=0;
end
for nn=jj-th:jj+th
ReIm(ii-th,nn,1)=255;
ReIm(ii-th,nn,2)=0;
ReIm(ii-th,nn,3)=0;
ReIm(ii+th,nn,1)=255;
ReIm(ii+th,nn,2)=0;
ReIm(ii+th,nn,3)=0;
end
end
th=th+1;
for kk=1:size(ReLis,1)
ii=ReLis(kk,1);
jj=ReLis(kk,2);
for mm=ii-th:ii+th
ReIm(mm,jj-th,1)=255;
ReIm(mm,jj-th,2)=0;
ReIm(mm,jj-th,3)=0;
ReIm(mm,jj+th,1)=255;
ReIm(mm,jj+th,2)=0;
ReIm(mm,jj+th,3)=0;
end
for nn=jj-th:jj+th
ReIm(ii-th,nn,1)=255;
ReIm(ii-th,nn,2)=0;
ReIm(ii-th,nn,3)=0;
ReIm(ii+th,nn,1)=255;
ReIm(ii+th,nn,2)=0;
ReIm(ii+th,nn,3)=0;
end
end
imshow(ReIm);
