clc;
clear all;
close all;
%点检测测试图像（Detection of Isolated Points）
I=im2double(imread('D:\Gray Files\10-4.tif'));
[M,N]=size(I);

%%
%===========================孤立点检测=====================================
%输出图像
g=zeros(M,N);
%拉普拉斯内核（Laplacian kernel），各向同性
L=[1 1 1;1 -8 1;1 1 1];
n_l=floor(size(L,1)/2);
%对原图进行扩展，方便处理边界
I_pad=padarray(I,[n_l,n_l],'symmetric');
%卷积翻转
L=-L;
%门限
T=0.89903;
for i=1:M
    for j=1:N
        %获得图像子块区域
        Block=I_pad(i:i+2*n_l,j:j+2*n_l);
        %用拉氏内核对子区域卷积
        if abs(sum(sum(Block.*L)))>T
            g(i,j)=1;
        end
    end
end
%对孤立点做开操作，删除小点
g=Opening(g,3);
imshow(g)

开操作函数，Opening如下：

%开操作
% I 输入图像
% n 结构元素长度
function J_Opening=Opening(I,n)
    B=ones(n,n);
    n_B=length(find(B==1));
    %这里需要B对其原点进行翻转，因为B是对称的，所以翻转后的结果与其本身相同
    l_m=floor(n/2);
    l_n=floor(n/2);
    [M,N]=size(I);
    %将原始图像进行扩展，这里采用了镜像扩展，以进行图像边缘计算    
    J_Erosion=zeros(M,N);
    J_Opening=zeros(M,N);
    I_pad=padarray(I,[l_m,l_n],'symmetric');
    %腐蚀操作
    for x=1:M
        for y=1:N
            %从扩展图像中取出子图像
             Block=I_pad(x:x+2*l_m,y:y+2*l_n);
             %将结构元素与子图像点乘,即逻辑“与”操作
             c=B.*Block;
             %比较结构元素与c中的1的数量，如果一样多，则该点的值为1
             ind=find(c==1);
             if length(ind)==n_B
                 J_Erosion(x,y)=1;
             end
        end
    end
    %膨胀操作
    J_Erosion_pad=padarray(J_Erosion,[l_m,l_n],'symmetric');
    for x=1:M
        for y=1:N
            %从扩展图像中取出子图像
             Block=J_Erosion_pad(x:x+2*l_m,y:y+2*l_n);
             %将结构元素与子图像点乘,即逻辑“与”操作
             c=B.*Block;
             %比较结构元素与c中的0的数量，如都为0，则该点的值为0，否则为1
             ind=find(c==0);
             if length(ind)~=n_B
                 J_Opening(x,y)=1;
             end
        end
    end
end