圆斑检测

一 简介

本次实验利用拉普拉斯边缘信息改进全局阈值处理的方法，对原图进行自适应阈值分割;利用函数regionprops函数计算圆斑区域描绘子，进而得到符合要求的圆斑。此次试验在上一次试验的基础上改进了两点，1）采用了自适应阈值分割，提高了自动化程度；2）在处理过程中，只使用了一次二值化，避免了圈定圆斑时出现重复的问题。

二 算法的步骤

第一步 预处理
1.将原图转化为二值图像。
2.高斯滤波。

第二步 自适应阈值分割
自适应阈值分割采用拉普拉斯边缘信息改进全局阈值处理的方法，具体算法步骤如下：
1.f为原图。求图像拉普拉斯的绝对值，并将其归一化到[0, 1]范围内。
2.画出归一化后的直方图，使用99.5%估计拉普拉斯绝对值的阈值，记为Q。再用Q对拉普拉斯的绝对值做阈值处理，形成标记图像。原图与标记图相乘，得到结果的记为fp，并画出其直方图，记为hp 。
3.图像fp包含围绕背景和物体边界的f的像素，所以，fp的直方图受0的控制。但是，我们的兴趣在于物体边界周围的分割值，所以需要消除0对直方图的贡献。因此，把fp的第一个元素排除在外。
4.用最终得到的直方图hp ，计算Otsu阈值。

第三步 得到符合要求的圆斑
利用函数regionprops函数计算圆斑区域描绘子，通过计算候选圆斑的偏心率，面积大小，确定符合要求的圆斑。再通过计算目标圆斑的质点坐标，半径，确定圆斑的位置和大小，并将其画出。

三 处理结果

图1 pos71_run3_img230原图

图2 pos71_run3_img230处理后

四 源码

源码第1/6 speckle_main.m

close all;
clear;
clc;

%读取图片
ori=imread('pos71_run7_img190.jpg'); 
% pos71_run3_img230  两张图的函数名，方便复制
% pos71_run7_img190

%转换为灰度图像
grayImg=rgb2gray(ori);

%高斯滤波
sigma = 1.6;      
gausFilter = fspecial('gaussian',[5 5],sigma);
src=imfilter(grayImg,gausFilter,'replicate');
figure(1);
imshow(src);
title('高斯滤波后');

% 拉普拉斯边缘信息改进全局阈值处理
img = L_edge_threshold( src );
%disp(strcat('阈值:',num2str(??)));%在command window里显示出 :阈值
figure(3)
imshow(img)
title('拉普拉斯边缘信息改进全局阈值处理后结果');

% 计算偏心率，偏心率为0的椭圆是圆，偏心率为1的椭圆是直线
imLabel = bwlabel(img);               
stats = regionprops(imLabel,'Eccentricity');    
eccentricity = cat(1,stats.Eccentricity); 
index = find(eccentricity < 0.75);%偏心率小于0.75被认为是候选圆斑        
img = ismember(imLabel,index);
figure(4);
imshow(img);
title('偏心率限制后结果');

% 计算圆斑的面积，质心坐标，半径
imLabel = bwlabel(img);               
stats = regionprops(imLabel,'Area','Centroid','EquivDiameter');    
area = cat(1,stats.Area);  
index = find(area > 10); %面积小于10的候选圆斑，被认为是噪声      
img = ismember(imLabel,index);
figure(5);
imshow(img);
title('面积限制后结果');

figure(6);
imshow(ori);
title('目标圆斑显示');
hold on;
for i=1:length(index)
    x = stats(index(i)).Centroid(1);  
    y = stats(index(i)).Centroid(2);  
    ar = stats(index(i)).EquivDiameter / 2;  
    circle(x,y,ar);
end

%保存
F=getframe(gcf); % 获取整个窗口内容的图像
imwrite(F.cdata,'G:\工作汇报\假期1\圆斑检测 周立刚\实验结果\pos71_run7_img190.jpg')

源码2/6 L_edge_threshold.m

function [ g ] = L_edge_threshold( src )
% 拉普拉斯边缘信息改进全局阈值处理

%输入图像的直方图
f = tofloat(src);
imhist(f)           
hf = imhist(f);

%用函数graythresh对直方图进行分割
[Tf SMf] = graythresh(f);
gf = im2bw(f, Tf);      
figure, imshow(gf)  

%求拉普拉斯绝对值边缘
w = [-1 -1 -1; -1 8 -1; -1 -1 -1];
lap = abs(imfilter(f, w, 'replicate')); 
lap = lap/max(lap(:));%lap的值归一化到[0, 1]范围内
h = imhist(lap);

%用99.5%估计梯度的阈值，保留绝对值边缘图像中较大的值
%这个值应该发生在接近物体和背景的边界处：
Q = percentile2i(h, 0.995);
markerImage = lap > Q;%标记图像
fp = f.*markerImage;%标记图和输入图的乘积
figure;
imshow(fp); 
hp = imhist(fp);

% 图像fp包含围绕背景和物体边界的f的像素，所以，fp的直方图受0的控制。
% 因为我们的兴趣在于物体边界周围的分割值，所以需要消除0对直方图的贡献。
% 因此，把fp的第一个元素排除在外

hp(1) = 0;
figure, bar(hp, 0) 

% 用结果的直方图得到Otsu阈值
T = otsuthresh(hp);
g = im2bw(f, T);   

end

源码第3/6 tofloat.m

% [out,revertclass] = tofloat(inputimage)
% inputimage 输入图像
% out            输出浮点图像,
% revertclass 输入图像的类型
% 另外，tofloat为冈萨雷斯数中的IPT函数，附上源码：
function [out,revertclass] = tofloat(inputimage)
%Copy the book of Gonzales
identify = @(x) x;
tosingle = @im2single;
table = {'uint8',tosingle,@im2uint8 
         'uint16',tosingle,@im2uint16 
         'logical',tosingle,@logical
         'double',identify,identify
         'single',identify,identify};
classIndex = find(strcmp(class(inputimage),table(:,1)));
if isempty(classIndex)
    error('不支持的图像类型');
end
out = table{classIndex,2}(inputimage);
revertclass = table{classIndex,3};

源码4/6 otsuthresh.m

function [T, SM] = otsuthresh(h)
%OTSUTHRESH Otsu's optimum threshold given a histogram.
% [T, SM] = OTSUTHRESH(H) computes an optimum threshold, T, in the
% range [0 1] using Otsu's method for a given a histogram, H.
% Normalize the histogram to unit area. If h is already normalized,
% the following operation has no effect.
h = h/sum(h);
h = h(:); % h must be a column vector for processing below.
% All the possible intensities represented in the histogram (256 for
% 8 bits). (i must be a column vector for processing below.)
i = (1:numel(h))';
% Values of P1 for all values of k.
P1 = cumsum(h);
% Values of the mean for all values of k.
m = cumsum(i.*h);
% The image mean.
mG = m(end);
% The between-class variance.
sigSquared = ((mG*P1 - m).^2)./(P1.*(1 - P1) + eps);
% Find the maximum of sigSquared. The index where the max occurs is
% the optimum threshold. There may be several contiguous max values.
% Average them to obtain the final threshold.
maxSigsq = max(sigSquared);
T = mean(find(sigSquared == maxSigsq));
% Normalized to range [0 1]. 1 is subtracted because MATLAB indexing
% starts at 1, but image intensities start at 0.
T = (T - 1)/(numel(h) - 1);
% Separability measure.
SM = maxSigsq / (sum(((i - mG).^2) .* h) + eps);

源码5/6 percentile2i.m

function I=percentile2i(h,P)
%PERCENTILE2I Computes an intensity value given a percentile.
% I=PERCENTILE2I(H,P) Given a percentile ,p,and a histogram,
% H, this function computes an intensity ,I representing 
% the Pth percentile and returns the value  in I. P must be in the 
% range [0,1] and I is returned as a value in the range [0,1] also.
% Example:
% I=percentile2i(h,0.5)

% Check value of P
if P<0 || P>1
    error('The percentile must be in the range [0,1].');
end

% Normalized the histogram to unit area.if it is already normalized
% the following computation has no effect.
h=h/sum(h);

% Camulative distribution
C=cumsum(h);
% Calculations.
idx=find(C >= P,1,'first');
% Subtract 1 from idx because indexing starts at 1,but intensities
% start at 0. Also ,normalize to the range [0,1].
I=(idx-1)/(numel(h)-1);

源码6/6 circle.m

function [] = circle( x,y,r )  
%画圆
theta=0:0.1:2*pi;  
Circle1=x+r*cos(theta);  
Circle2=y+r*sin(theta);  

plot(Circle1,Circle2,'b','linewidth',2);  
axis equal  
end