function zernike7(I)
%% 1:读图
I=imread('cropped_image.bmp'); %读图
subplot(231) %创建 2x2 布局,选择第 1 个子图
imshow(I)
title('原图');
%% 2:7*7Zernik模板,定义一个模板矩阵
%M00 0阶0次模板 常数项 突出中心像素
M00=...
[
0 0.0287 0.0686 0.0807 0.0686 0.0287 0
0.0287 0.0815 0.0816 0.0816 0.0816 0.0815 0.0287
0.0686 0.0816 0.0816 0.0816 0.0816 0.0816 0.0686
0.0807 0.0816 0.0816 0.0816 0.0816 0.0816 0.0807
0.0686 0.0816 0.0816 0.0816 0.0816 0.0816 0.0686
0.0287 0.0815 0.0816 0.0816 0.0816 0.0815 0.0287
0 0.0287 0.0686 0.0807 0.0686 0.0287 0
];
%M11R Zernike矩1阶1次 实部模板
M11R=...
[
0 -0.015 -0.019 0 0.019 0.015 0
-0.0224 -0.0466 -0.0233 0 0.0233 0.0466 0.0224
-0.0573 -0.0466 -0.0233 0 0.0233 0.0466 0.0573
-0.069 -0.0466 -0.0233 0 0.0233 0.0466 0.069
-0.0573 -0.0466 -0.0233 0 0.0233 0.0466 0.0573
-0.0224 -0.0466 -0.0233 0 0.0233 0.0466 0.0224
0 -0.015 -0.019 0 0.019 0.015 0
];
%与 M11R 配合,覆盖全方向边缘检测
M11I=...
[
0 -0.0224 -0.0573 -0.069 -0.0573 -0.0224 0
-0.015 -0.0466 -0.0466 -0.0466 -0.0466 -0.0466 -0.015
-0.019 -0.0233 -0.0233 -0.0233 -0.0233 -0.0233 -0.019
0 0 0 0 0 0 0
0.019 0.0233 0.0233 0.0233 0.0233 0.0233 0.019
0.015 0.0466 0.0466 0.0466 0.0466 0.0466 0.015
0 0.0224 0.0573 0.069 0.0573 0.0224 0
];
%M20核心作用:通过模板与图像卷积结果,计算边缘到像素中心的距离
M20=...
[
0 0.0225 0.0394 0.0396 0.0394 0.0225 0
0.0225 0.0271 -0.0128 -0.0261 -0.0128 0.0271 0.0225
0.0394 -0.0128 -0.0528 -0.0661 -0.0528 -0.0128 0.0394
0.0396 -0.0261 -0.0661 -0.0794 -0.0661 -0.0261 0.0396
0.0394 -0.0128 -0.0528 -0.0661 -0.0528 -0.0128 0.0394
0.0225 0.0271 -0.0128 -0.0261 -0.0128 0.0271 0.0225
0 0.0225 0.0394 0.0396 0.0394 0.0225 0
];
%% 3:图像预处理
if length(size(I))==3 I=rgb2gray(I); end %转为灰度图像
%二值化
% 设定阈值为0.6(0-1范围),图像中凡是超过255*0.6=153的灰度值均变为255,否则都为0,即二值化
I=im2bw(I,0.2);
%UNIT变为double型
K=double(I);
subplot(232)
imshow(I)
title('二值化图');
%% 4:图像K与模版卷积
A11I=conv2(M11I,K); % M11I(虚部模板)与K卷积,结果A11I为边缘虚部特征图(对应对角方向梯度)
A11R=conv2(M11R,K); % M11R(实部模板)与K卷积,结果A11R为边缘实部特征图(对应水平/垂直方向梯度)
A20=conv2(M20,K); % M20(2阶模板)与K卷积,结果A20为边缘曲率/距离特征图
% 去除边缘多余部分7×7模板卷积后,图像会向四周各扩展3个像素(共扩展6个像素),需裁剪回原尺寸m×n。
% 行:从第4行到倒数第3行;列:从第4列到倒数第3列,恢复原尺寸
A11I=A11I(4:end-3,4:end-3);
A11R=A11R(4:end-3,4:end-3);
A20=A20(4:end-3,4:end-3);
subplot(234)
imshow(A11I)
title('A11I实部-对角梯度图');
subplot(235)
imshow(A11R)
title('A11R虚部-水平/垂直梯度图');
subplot(236)
imshow(A20)
title('A20-边缘曲率/距离特征图');
%% 5:亚像素边缘参数计算(突破像素级精度的核心)
%初始化J变量,boundar也初始化,K为数组1944x2592
J=zeros(size(K)); %初始化一个与K同大小的零矩阵J
%计算θ l k 三个参数
theta=atan2(A11I,A11R); % 计算theta 边缘法线方向角θ
A11C=A11R.*cos(theta)+A11I.*sin(theta); % 计算边缘幅值A11C
l=A20./A11C; % 计算亚像素距离l:边缘到当前像素中心的距离
k=1.5*A11C./((1-l.^2).^1.5); % 计算边缘曲率k:描述边缘的弯曲程度,k值越大,边缘越陡峭
e=abs(l)>1/3.5; % 筛选异常值:l的绝对值>1/3.5≈0.286时,视为无效边缘(超出合理亚像素距离范围)
k(e)=0; % 在条件 e 为真的情况下,将 k 的对应元素置为 0 无效位置的曲率设为0,避免后续边缘筛选时误判为有效边缘
%边缘判断条件 陡峭(k 大)且连续(l 稳定)
a=abs(l)<1/sqrt(2)*2/7; % 创建条件 a,l的绝对值<0.202(筛选距离在合理亚像素范围内的边缘,避免极端值)
b=abs(k)>max(I(:))/10; % 创建条件 b,曲率绝对值>图像最大值的10%(筛选强度足够的边缘,去除弱噪声边缘)
% a,b分别为距离和边缘强度判断结果
J(a&b)=1;
figure(2)
subplot(221) %创建 2x2 布局,选择第 1 个子图
imshow(I)
title('预处理后的图像');
% [x1sub y1sub],[x2sub,y2sub]
subplot(222)
imshow(J)
title('满足条件后的提取边界图');
subplot(223);
I41=imfill(J,'holes'); % 使用 imfill 函数填充二值图像 J 中的孔洞
imshow(I41)
title('填充孔洞');
subplot(224);
I4=bwperim(I41); % 提取二值图像 I41 的最外围边界
imshow(I4);
title('提取外边界');
%% 6:8连通图像处理
[L,num]=bwlabel(J,8);%对二值图像8邻域连通性分析,L为与J相同大小的矩阵,但为以组件序号的矩阵,num为连通区域数
%自动搜索28个连通区域组件,如查找L==i的组件
s=zeros(1,num);%设S为1*num(28)的全0矩阵1*28,初始化s为0
for i=1:num
%查看L==i的,即第标签为i的组件的大小为行大小,列为1
% 统计第i个轮廓的像素数:find(L==i)找到该轮廓的所有像素坐标,size(...,1)返回坐标行数(即像素数)
s(i)=size(find(L==i),1);
end
[bwL,label]=sort(s,'descend'); %对S中数据按照降序排列,并返回在原sz中的位置在label中
%降序排列的s数组bwL为最大为8943,最小为1,共28个数据,并返回原s中对应的标签位置
if label(1)<label(2)
index1=label(1);
index2=label(2);
else
index1=label(2);
index2=label(1);
end
%% 7:特征点亚像素定位 功能:计算第二个特征点的亚像素坐标,原理同第一个点
%第一个主要轮廓的特征点提取
%计算左边探针的最前端坐标,找到L==index1=1的连通组件1的坐标位置r1和c1
[r1,c1]=find(L==index1); % 找到 L 中等于 index1 的所有行和列索引
A1=[r1 c1]; % 将行列索引组合为一个坐标矩阵
y1=max(A1(:,2)); %该连通域中y最大为连通域,行任意,A1中所有第2列的最大值2592
x1=max(A1(find(A1(:,2)==y1),1)); % 找到对应于最大 y 值的 x1 值(行索引)
%% 亚像素坐标点计算
% 亚像素原理:
%l(x1,y1):Zernike 矩计算的距离参数(边缘到像素中心的距离,单位:像素)
%theta(x1,y1):边缘的法线方向角(单位:弧度)
%3.5:7×7 模板的半径系数(模板尺寸的一半)
%三角函数转换:将极坐标形式的偏移量转换为直角坐标系下的 x、y 偏移
%精度提升:将传统像素级定位(±0.5 像素误差)提升到亚像素级(通常 ±0.1 像素误差)
x1sub=x1+3.5*l(x1,y1)*cos(theta(x1,y1)); % 计算 x1 的偏移
y1sub=y1+3.5*l(x1,y1)*sin(theta(x1,y1)); % 计算 y1 的偏移
% 找到标签为 index2 的连通组件(通常是次大的连通区域)
[r2,c2]=find(L==index2); % 找到 L 中等于 index2 的所有行和列索引
A2=[r2 c2]; % 将行列索引组合为一个坐标矩阵
y2=min(A2(:,2)); % 在连通区域中第 2 列中的最小值 y2
x2=max(A2(find(A2(:,2)==y2),1));
x2sub=x2+3.5*l(x2,y2)*cos(theta(x2,y2));
y2sub=y2+3.5*l(x2,y2)*sin(theta(x2,y2));
%% 亚像素边缘精确定位
% 计算边缘强度(用于阈值筛选)
edge_strength = sqrt(A11R.^2 + A11I.^2);
%设置为边缘强度的 20%(可调参数) 阈值设置:采用自适应阈值(最大值的 20%),确保在不同光照条件下都能筛选出有效边缘
strength_threshold = 0.2 * max(edge_strength(:)); % 设置强度阈值可调参数
% 获取候选边缘点位置(满足强度阈值且l在有效范围内)
[rows, cols] = size(K);
[y, x] = find(edge_strength > strength_threshold & abs(l) < 1/10 & ~isnan(l));
%% 亚像素坐标计算
subpixel_points = zeros(length(x), 2); %存储亚像素坐标的矩阵(N×2)
valid_idx = true(length(x), 1); %标记有效点的逻辑数组(初始全为 true)
for i = 1:length(x)
% 获取当前点的Zernike矩参数 获取当前点的方向角(current_theta)和距离参数(current_l)
current_theta = theta(y(i), x(i));
current_l = l(y(i), x(i));
% 计算亚像素偏移量(注意坐标方向)
x_offset = current_l * cos(current_theta);
y_offset = current_l * sin(current_theta);
% 计算亚像素坐标
x_sub = x(i) + x_offset;
y_sub = y(i) + y_offset;
% 验证坐标是否在图像范围内
if x_sub < 1 || x_sub > cols || y_sub < 1 || y_sub > rows
valid_idx(i) = false; % 如果超出范围,标记为无效
continue; % 跳过当前循环
end
subpixel_points(i, :) = [x_sub, y_sub]; % 存储有效的亚像素坐标
end
% 去除无效点
subpixel_points = subpixel_points(valid_idx, :);
%% 边缘点后处理
% % 1. 去除重复点(四舍五入到整数像素后去重)
% rounded_points = round(subpixel_points); % 对亚像素点进行四舍五入
% [~, unique_idx] = unique(rounded_points, 'rows', 'stable'); % 找到唯一的点索引
% subpixel_points = subpixel_points(unique_idx, :); % 只保留唯一的亚像素点
% 边缘点后处理
% 1. 去除重复点
rounded_points = round(subpixel_points);
[~, unique_idx] = unique(rounded_points, 'rows', 'stable');
subpixel_points = subpixel_points(unique_idx, :);
% 2. 新增:过滤偏离整体趋势的异常点(基于局部距离)
if size(subpixel_points,1) >= 5
% 计算每个点到前后2个点的平均距离
dists = zeros(size(subpixel_points,1),1);
for i = 3:size(subpixel_points,1)-2
prev2 = subpixel_points(i-2,:);
prev1 = subpixel_points(i-1,:);
next1 = subpixel_points(i+1,:);
next2 = subpixel_points(i+2,:);
% 计算当前点与前后点的距离偏差
d_prev = mean(sqrt(sum(([prev2; prev1] - subpixel_points(i,:)).^2,2)));
d_next = mean(sqrt(sum(([next1; next2] - subpixel_points(i,:)).^2,2)));
dists(i) = mean([d_prev, d_next]);
end
% 过滤距离过大的异常点(超过均值2倍)
if ~isempty(dists)
thresh = 2 * mean(dists(dists>0));
valid = dists < thresh | dists == 0; % 保留合理点和边界点
subpixel_points = subpixel_points(valid,:);
end
end
%% 可视化结果
figure;
subplot(222);
imshow(I);
hold on; % 保持当前图像,使后续绘图能够叠加在此图上
plot(subpixel_points(:,1), subpixel_points(:,2), 'y.', 'MarkerSize', 6); % 绘制亚像素点
%{
参数逐解析:
subpixel_points(:,1):取亚像素坐标矩阵的第 1 列(x 轴坐标,对应图像的列索引);
subpixel_points(:,2):取亚像素坐标矩阵的第 2 列(y 轴坐标,对应图像的行索引);
'r.':设置绘图样式为 "红色圆点"(r= 红色,.= 点型标记),红色在黑白 / 灰度原图中对比度高,便于观察;
'MarkerSize', 6:设置圆点标记的大小为 6(数值越大,点越明显,避免因点太小被原图细节掩盖)。
功能:将计算出的所有亚像素边缘点叠加在原图上,可直接验证边缘点是否连续、是否准确贴合零件边缘(如机械零件的外圆、孔洞轮廓)。
%}
title('亚像素边缘检测结果'); % 设置子图的标题
hold off; % 设置子图的标题
% 显示边缘强度图
subplot(223);
imagesc(edge_strength);% 将边缘强度矩阵以彩色图的形式可视化
axis image; % 设置坐标轴的比例,使得每个像素在图中显示为正方形
colormap jet; % 使用 'jet' 色图来映射数据值
colorbar; % 添加颜色条,显示色图和数据值的对应关系
title('边缘强度图'); % 设置子图的标题
% 显示距离参数l
subplot(224);
imagesc(l);
axis image;
colormap jet;
colorbar;
title('距离参数 l');
figure(7);
imshow(I);
hold on;
h_original_points = plot(subpixel_points(:,1), subpixel_points(:,2), 'b.', 'MarkerSize', 3);% 绘制亚像素点
title('亚像素边缘检测结果');
% 新增代码:椭圆分割与分段拟合
if ~isempty(subpixel_points) && size(subpixel_points,1)>=5
% 1. 先按角度排序亚像素点(确保点按椭圆轮廓顺序排列)
centroid = mean(subpixel_points); % 椭圆质心
angles = atan2(subpixel_points(:,2)-centroid(2), subpixel_points(:,1)-centroid(1));
[~, idx] = sort(angles);
sorted_points = subpixel_points(idx, :); % 按角度排序后的点(形成闭合轮廓)
% 2. 分割为圆弧段和直线段(调整曲率阈值控制分割)
curvature_thresh = 0.05; % 阈值越小,越容易识别为圆弧(可根据实际数据调整)
[arc_segs, line_segs] = segment_contour(sorted_points, curvature_thresh);
% 3. 拟合圆弧段(蓝色虚线)
for i = 1:length(arc_segs)
arc_points = arc_segs{i};
[xc, yc, r] = circle_fit_segment(arc_points);
% 绘制圆弧(只画对应段的角度范围)
arc_angles = atan2(arc_points(:,2)-yc, arc_points(:,1)-xc);
theta_start = min(arc_angles);
theta_end = max(arc_angles);
theta = linspace(theta_start, theta_end, 100);
x_arc = xc + r*cos(theta);
y_arc = yc + r*sin(theta);
plot(x_arc, y_arc, 'b--', 'LineWidth', 1.5); % 蓝色虚线表示圆弧
end
% 4. 拟合直线段(绿色实线)
for i = 1:length(line_segs)
line_points = line_segs{i};
[k, b] = line_fit_segment(line_points);
% 绘制直线(只画对应段的x范围)
x_min = min(line_points(:,1));
x_max = max(line_points(:,1));
x_line = linspace(x_min, x_max, 100);
y_line = k*x_line + b;
plot(x_line, y_line, 'g-', 'LineWidth', 1.5); % 绿色实线表示直线
end
% 添加图例
legend(['红色亚像素点 (', num2str(size(subpixel_points,1)), '个)'], ...
'蓝色虚线:拟合圆弧段', '绿色实线:拟合直线段', 'Location', 'best');
end
hold off;
%% 输出结果
disp(['检测到亚像素边缘点数量: ', num2str(size(subpixel_points,1))]);
disp(['使用的强度阈值: ', num2str(strength_threshold)]);
fprintf('发现 %d 个独立的边缘轮廓\n', num);
curve_colors = lines(num); %lines(num)是 MATLAB 内置的颜色生成函数,生成一个num×3的 RGB 颜色矩阵(num为连通域数量,即轮廓数量),每一行对应一种独特的颜色
% 原始点颜色(红色)
original_point_color = [1 0 0]; % 红色
% 圆拟合点颜色(绿色)
circle_fit_color = [0 1 0]; % 绿色
% 存储每类的亚像素点坐标 cell(细胞数组),每个单元格可存储任意类型 / 尺寸的数据(
classified_points = cell(num, 1);
% 存储每个轮廓的圆度
circularity = zeros(num, 1);
% 存储圆拟合结果
circle_results = struct('xc', cell(num, 1), 'yc', cell(num, 1), ...
'R', cell(num, 1), 'diameter', cell(num, 1), ...
'fit_error', cell(num, 1), 'circularity', cell(num, 1));
% 新增代码(添加在circle_results之后)
ellipse_results = struct('A', cell(num, 1), 'B', cell(num, 1), 'C', cell(num, 1), ...
'D', cell(num, 1), 'E', cell(num, 1), 'F', cell(num, 1), ...
'center', cell(num, 1), 'major_axis', cell(num, 1), ...
'minor_axis', cell(num, 1), 'orientation', cell(num, 1), ...
'fit_error', cell(num, 1));
% 圆度阈值设置
% CIRCULARITY_THRESHOLD = 0.74; % 圆度大于此值视为圆形轮廓
CIRCULARITY_THRESHOLD = 0.6; % 降低圆阈值(圆形轮廓需更高圆度)
ELLIPTICITY_THRESHOLD = 0.3; % 新增椭圆阈值(圆度在0.3~0.6之间视为椭圆)
% 将亚像素点转换为整数坐标,确保在图像范围内
rounded_coords = round(subpixel_points);
valid_indices = rounded_coords(:, 2) >= 1 & rounded_coords(:, 2) <= size(J, 1) & ...
rounded_coords(:, 1) >= 1 & rounded_coords(:, 1) <= size(J, 2);
% 处理每个连通区域
for k = 1:num
% 收集属于当前连通域的亚像素点
cluster_points = [];
for i = 1:size(subpixel_points, 1)
if valid_indices(i) && L(rounded_coords(i, 2), rounded_coords(i, 1)) == k
cluster_points = [cluster_points; subpixel_points(i, :)];
end
end
% 保存分类结果
classified_points{k} = cluster_points; % 将当前轮廓的亚像素点存入细胞数组对应位置
n_points = size(cluster_points, 1); % 获取当前轮廓的亚像素点数量
fprintf('轮廓 %d: 包含 %d 个点\n', k, n_points);
%{
圆度是衡量 "轮廓与理想圆形接近程度" 的指标,公式为4π×面积/周长?,取值范围 0-1:
1:完美圆形;
0.7-1:接近圆形(如机械零件的孔、轴);
<0.7:非圆形(如矩形、条形孔)。
%}
% 计算轮廓圆度(只有点数足够时)
if n_points >= 3 % 圆度计算的前提:至少3个点(3点确定一个圆,点数越多拟合越准确)
circularity(k) = calculate_circularity(cluster_points); % 调用圆度计算函数
fprintf('轮廓 %d 圆度: %.4f\n', k, circularity(k)); % 输出圆度值(保留4位小数)
else
circularity(k) = 0; % 点数不足时圆度为0
end
end
%% 圆拟合
% 创建图例对象数组
legend_handles = gobjects(num + 1, 1); % +1 用于原始点
legend_labels = cell(num + 1, 1);
% 第一个图例项是原始点
legend_handles(1) = h_original_points;
legend_labels{1} = '原始亚像素点';
% 对每类点进行曲线连接和圆拟合(仅对圆形轮廓)
h_circle_fit = gobjects(num, 1); % 存储圆拟合的句柄
for k = 1:num
cluster_points = classified_points{k};
n = size(cluster_points, 1);
if n < 2
fprintf('轮廓 %d 点数不足,无法连接\n', k);
continue;
end
% 1. 计算点集的质心 计算当前轮廓所有点的质心(平均值),作为排序参考点
centroid = mean(cluster_points, 1);
% 2. 按角度排序点,形成轮廓顺序
%计算每个点相对于质心的角度(使用 atan2 函数)
%将负角度转换为 0-2π 范围内的正角度
%按角度排序并重新排列点的顺序
% angles = atan2(cluster_points(:, 2) - centroid(2), cluster_points(:, 1) - centroid(1));
% angles(angles < 0) = angles(angles < 0) + 2*pi;
% [~, idx] = sort(angles);
% sorted_points = cluster_points(idx, :);
% 2. 优化排序:先过滤异常点,再按角度排序
% 计算点到质心的距离,过滤过远点(异常值)
dists_to_centroid = sqrt( (cluster_points(:,1)-centroid(1)).^2 + (cluster_points(:,2)-centroid(2)).^2 );
valid = dists_to_centroid < 1.5*mean(dists_to_centroid); % 保留1.5倍均值内的点
cluster_points = cluster_points(valid, :);
if size(cluster_points,1) < 3
continue; % 点数不足则跳过
end
% 重新计算质心和角度
centroid = mean(cluster_points, 1);
angles = atan2(cluster_points(:, 2) - centroid(2), cluster_points(:, 1) - centroid(1));
angles(angles < 0) = angles(angles < 0) + 2*pi;
% 按角度排序,并确保连续(避免因角度突变导致的顺序错误)
[~, idx] = sort(angles);
sorted_points = cluster_points(idx, :);
% 检查相邻点角度差,修正跳变(如从350度到10度应为连续,而非跳变)
angle_diff = diff(angles(idx));
jump_idx = find(angle_diff > pi); % 角度差超过180度视为跳变
if ~isempty(jump_idx)
% 将跳变点后的点移到前面,保证连续性
sorted_points = [sorted_points(jump_idx+1:end, :); sorted_points(1:jump_idx, :)];
end
% 3. 添加第一个点到最后,形成闭合轮廓
closed_points = [sorted_points; sorted_points(1, :)];
% 4. 计算累积弦长作为参数
dx = diff(closed_points(:, 1));
dy = diff(closed_points(:, 2));
distances = sqrt(dx.^2 + dy.^2);
t = [0; cumsum(distances)];
% 5. 样条插值生成光滑轮廓线 相比直接连接原始点,插值后的轮廓线更加平滑自然
tq = linspace(0, t(end), 300);
xx = spline(t, closed_points(:, 1), tq);
yy = spline(t, closed_points(:, 2), tq);
% 6. 绘制光滑轮廓线(使用曲线颜色)
h_curve = plot(xx, yy, 'Color', curve_colors(k, :), 'LineWidth', 2);
% 7. 添加图例项(轮廓线)
legend_handles(k+1) = h_curve;
legend_labels{k+1} = sprintf('轮廓 %d 轮廓线 (圆度: %.2f)', k, circularity(k));
fprintf('轮廓 %d 轮廓线连接完成\n', k);
% 8. 对圆形轮廓进行圆拟合
if circularity(k) > CIRCULARITY_THRESHOLD
if numel(xx) >= 3
try
% 圆拟合函数
[xc, yc, R] = circle_fit(xx, yy);
diameter = 2 * R*1000*0.0000722705;
% 保存圆拟合结果
circle_results(k).xc = xc;
circle_results(k).yc = yc;
circle_results(k).R = R;
circle_results(k).diameter = diameter;
circle_results(k).fit_error = 0;
circle_results(k).circularity = circularity(k);
% 绘制拟合的圆(绿色虚线)
theta = linspace(0, 2*pi, 100);
x_circle = xc + R*cos(theta);
y_circle = yc + R*sin(theta);
% h_circle = plot(x_circle, y_circle, '--', 'Color', circle_fit_color, 'LineWidth', 1.5);
h_circle = plot(x_circle, y_circle, 'Color', circle_fit_color, ...
'LineStyle', '--', 'LineWidth', 1.5);
% 标记圆心
plot(xc, yc, '+', 'Color', circle_fit_color, 'MarkerSize', 10, 'LineWidth', 1.5);
% 显示直径文本(带背景色)
text(xc, yc, sprintf('直径: %.2f', diameter), ...
'Color', 'white', 'BackgroundColor', circle_fit_color, ...
'FontSize', 8, 'Margin', 0.5, 'HorizontalAlignment', 'center');
h_circle_fit(k) = h_circle;
fprintf('轮廓 %d (圆度 %.3f) 进行圆拟合 - 直径: %.2f\n', k, circularity(k), diameter);
catch ME
circle_results(k).fit_error = ME.message;
warning('轮廓 %d 圆拟合失败: %s', k, ME.message);
end
end
elseif circularity(k) > ELLIPTICITY_THRESHOLD && circularity(k) <= CIRCULARITY_THRESHOLD
% 1. 确保xx是一维数组(处理可能的二维情况)
xx = xx(:); % 强制转换为列向量(关键修正)
yy = yy(:); % 同步处理yy
% 确保xx和yy是插值后的有效点集(检查是否为空或点数不足)
if numel(xx) < 5
fprintf('轮廓 %d 点数量不足5个,跳过椭圆拟合\n', k);
continue;
end
% 3. 计算方差(过滤过于集中的点集)
x_var = var(xx);
y_var = var(yy);
if x_var < 1e-3 && y_var < 1e-3
fprintf('轮廓 %d 点集过于集中,跳过椭圆拟合\n', k);
continue;
end
% 4. 构造线性拟合矩阵(此时xx已为列向量,与ones维度匹配)
A_line = [xx, ones(length(xx), 1)]; % 现在可正常拼接
b_line = yy;
params_line = A_line \ b_line;
y_pred = A_line * params_line;
line_residual = mean(abs(y_pred - yy));
if line_residual < 1e-2
fprintf('轮廓 %d 接近直线,跳过椭圆拟合\n', k);
continue;
end
try
% 调用椭圆拟合函数
[A,B,C,D,E,F, center, major_axis, minor_axis, orientation] = ellipse_fit(xx, yy);
% 保存椭圆拟合结果
ellipse_results(k).A = A;
ellipse_results(k).B = B;
ellipse_results(k).C = C;
ellipse_results(k).D = D;
ellipse_results(k).E = E;
ellipse_results(k).F = F;
ellipse_results(k).center = center;
ellipse_results(k).major_axis = major_axis;
ellipse_results(k).minor_axis = minor_axis;
ellipse_results(k).orientation = orientation;
ellipse_results(k).fit_error = ''; % 无错误
% 绘制椭圆(参数方程)
t = linspace(0, 2*pi, 100);
x_ellipse = center(1) + major_axis*cos(orientation).*cos(t) - minor_axis*sin(orientation).*sin(t);
y_ellipse = center(2) + major_axis*sin(orientation).*cos(t) + minor_axis*cos(orientation).*sin(t);
plot(x_ellipse, y_ellipse, 'm--', 'LineWidth', 1.5); % 品红色虚线
plot(center(1), center(2), 'm+', 'MarkerSize', 10); % 标记椭圆中心
fprintf('轮廓 %d (圆度 %.3f) 椭圆拟合成功\n', k, circularity(k));
catch ME
ellipse_results(k).fit_error = ME.message;
warning('轮廓 %d 椭圆拟合失败: %s', k, ME.message);
end
else
fprintf('轮廓 %d (圆度 %.3f) 非圆形轮廓,跳过圆拟合\n', k, circularity(k));
end
end
% 添加圆拟合到图例
if any(ishandle(h_circle_fit))
% 添加圆拟合图例句柄(如果存在成功的圆拟合)
valid_circles = ishandle(h_circle_fit);
circle_handle = h_circle_fit(find(valid_circles, 1)); % 使用第一个有效的圆
if ~isempty(circle_handle)
legend_handles(end+1) = circle_handle;
legend_labels{end+1} = '圆拟合';
end
end
% 添加图例(原始点 + 所有轮廓线 + 圆拟合)
valid_handles = ishandle(legend_handles);
legend(legend_handles(valid_handles), legend_labels(valid_handles), ...
'Location', 'bestoutside');
% 单独显示圆拟合结果(仅圆形轮廓)
figure(9);
imshow(I);
hold on;
title(sprintf('圆形轮廓圆拟合结果 (圆度>%.2f)', CIRCULARITY_THRESHOLD));
% 绘制原始点
plot(subpixel_points(:, 1), subpixel_points(:, 2), 'r.', 'MarkerSize', 8);
for k = 1:num
cluster_points = classified_points{k};
n = size(cluster_points, 1);
if n < 2 || circularity(k) <= CIRCULARITY_THRESHOLD
continue;
end
% 只绘制成功的圆拟合 *1000*0.0000722705
if ~isempty(circle_results(k).xc) && circle_results(k).fit_error == 0
xc = circle_results(k).xc;
yc = circle_results(k).yc;
R = circle_results(k).R;
diameter = circle_results(k).diameter;
% 绘制轮廓线(原始曲线颜色)
centroid = mean(cluster_points, 1);
angles = atan2(cluster_points(:, 2) - centroid(2), cluster_points(:, 1) - centroid(1));
angles(angles < 0) = angles(angles < 0) + 2*pi;
[~, idx] = sort(angles);
sorted_points = cluster_points(idx, :);
closed_points = [sorted_points; sorted_points(1, :)];
dx = diff(closed_points(:, 1));
dy = diff(closed_points(:, 2));
distances = sqrt(dx.^2 + dy.^2);
t = [0; cumsum(distances)];
% 5. 改进插值:先平滑点集,再用分段三次Hermite插值(避免过拟合)
% 对排序后的点进行平滑(移动平均)
window = 3; % 平滑窗口大小(奇数)
if size(sorted_points,1) >= window
smoothed_x = movmean(sorted_points(:,1), window, 'Endpoints', 'shrink');
smoothed_y = movmean(sorted_points(:,2), window, 'Endpoints', 'shrink');
else
smoothed_x = sorted_points(:,1);
smoothed_y = sorted_points(:,2);
end
closed_points = [smoothed_x, smoothed_y; smoothed_x(1), smoothed_y(1)]; % 闭合轮廓
% 重新计算累积距离
dx = diff(closed_points(:, 1));
dy = diff(closed_points(:, 2));
distances = sqrt(dx.^2 + dy.^2);
t = [0; cumsum(distances)];
% 用PCHIP插值(保形分段三次插值,避免超调凸起)
tq = linspace(0, t(end), 300);
xx = pchip(t, closed_points(:,1), tq); % 替换spline为pchip
yy = pchip(t, closed_points(:,2), tq);
% tq = linspace(0, t(end), 300);
%
%
%
% xx = spline(t, closed_points(:, 1), tq);
% yy = spline(t, closed_points(:, 2), tq);
plot(xx, yy, 'Color', curve_colors(k, :), 'LineWidth', 1, 'LineStyle', '-.');
% 绘制圆
theta = linspace(0, 2*pi, 100);
x_circle = xc + R*cos(theta);
y_circle = yc + R*sin(theta);
plot(x_circle, y_circle, 'g-', 'LineWidth', 2);
% 标记圆心
plot(xc, yc, 'r+', 'MarkerSize', 15, 'LineWidth', 2);
% 显示直径文本(带背景色)
text(xc, yc+100, sprintf('轮廓%d 直径: %.4f\n', k, diameter), ...
'Color', 'white', 'BackgroundColor', 'blue', ...
'FontSize', 10, 'Margin', 0.1, 'HorizontalAlignment', 'center');
end
end
% 新增椭圆拟合结果图
figure(10);
imshow(I);
hold on;
title(sprintf('椭圆轮廓拟合结果 (圆度: %.2f~%.2f)', ELLIPTICITY_THRESHOLD, CIRCULARITY_THRESHOLD));
% 绘制原始亚像素点
plot(subpixel_points(:, 1), subpixel_points(:, 2), 'r.', 'MarkerSize', 5);
% 遍历椭圆轮廓并绘制
for k = 1:num
if circularity(k) > ELLIPTICITY_THRESHOLD && circularity(k) <= CIRCULARITY_THRESHOLD && isempty(ellipse_results(k).fit_error)
% 提取椭圆参数
center = ellipse_results(k).center;
major = ellipse_results(k).major_axis;
minor = ellipse_results(k).minor_axis;
theta = ellipse_results(k).orientation;
% 绘制椭圆
t = linspace(0, 2*pi, 100);
x_ellipse = center(1) + major*cos(theta).*cos(t) - minor*sin(theta).*sin(t);
y_ellipse = center(2) + major*sin(theta).*cos(t) + minor*cos(theta).*sin(t);
plot(x_ellipse, y_ellipse, 'm-', 'LineWidth', 2); % 品红色实线
% 标记中心和参数
plot(center(1), center(2), 'm+', 'MarkerSize', 15);
text(center(1), center(2)+50, sprintf('轮廓%d: 长轴%.2f, 短轴%.2f', ...
k, major, minor), 'Color', 'white', 'BackgroundColor', 'purple');
end
end
hold off;
hold off;
%% 输出圆拟合结果到控制台
fprintf('\n=== 圆拟合结果 (仅圆形轮廓) ===\n');
for k = 1:num
if circularity(k) > CIRCULARITY_THRESHOLD && ~isempty(circle_results(k).xc) && circle_results(k).fit_error == 0
fprintf('轮廓 %d (圆度 %.2f): 圆心(%.2f, %.2f), 直径 %.4f\n', ...
k, circularity(k), circle_results(k).xc, circle_results(k).yc, circle_results(k).diameter);
end
end
% 保存圆拟合结果到文件
save('circle_fit_results.mat', 'circle_results', 'circularity');
fprintf('已保存圆拟合结果到 circle_fit_results.mat\n');
%% 圆拟合函数
function [xc, yc, R] = circle_fit(x, y)
% 最小二乘法圆拟合 (Pratt方法)
% 输入: x, y - 要拟合的点
% 输出: xc, yc - 圆心坐标, R - 半径
% 1. 转换为列向量
x = x(:);
y = y(:);
% 2. 确保点数量足够
if numel(x) < 3
error('圆拟合需要至少3个点');
end
% 3. 构建矩阵
% 方程为: (x - xc)^2 + (y - yc)^2 = R^2
% 展开: x^2 - 2xc*x + xc^2 + y^2 - 2yc*y + yc^2 - R^2 = 0
% 重写为: 2xc*x + 2yc*y + d = x^2 + y^2, 其中d = R^2 - xc^2 - yc^2
% 4. 构造线性方程组
A = [x, y, ones(size(x))];
b = x.^2 + y.^2;
% 5. 解超定方程组
params = A \ b;
% 6. 提取圆心和半径
xc = params(1)/2;
yc = params(2)/2;
R = sqrt(params(3) + xc^2 + yc^2);
% 7. 检查结果有效性
if ~isreal(R) || R <= 0
error('圆拟合失败: 无效的半径 %.2f', R);
end
end
function circ = calculate_circularity(points)
% 计算轮廓圆度:(4π×面积) / (周长?),圆形为1,椭圆<1,越扁越接近0
centroid = mean(points, 1); % 质心
n = size(points, 1);
if n < 3
circ = 0;
return;
end
% 1. 按角度排序点(确保点按轮廓顺序排列,避免面积/周长计算错误)
angles = atan2(points(:,2) - centroid(2), points(:,1) - centroid(1));
[~, idx] = sort(angles);
sorted_points = points(idx, :);
% 闭合轮廓(首尾相连)
closed_points = [sorted_points; sorted_points(1, :)];
% 2. 计算面积(使用多边形面积公式,更准确)
x = closed_points(:,1);
y = closed_points(:,2);
area = polyarea(x, y);
if area <= 0
circ = 0;
return;
end
% 3. 计算周长(累加相邻点距离,更准确)
perimeter = 0;
for i = 1:n
dx = closed_points(i+1,1) - closed_points(i,1);
dy = closed_points(i+1,2) - closed_points(i,2);
perimeter = perimeter + sqrt(dx^2 + dy^2);
end
if perimeter <= 0
circ = 0;
return;
end
% 4. 计算圆度:(4π×面积) / (周长?)
circ = (4 * pi * area) / (perimeter^2);
% 限制范围在 [0,1]
circ = max(0, min(1, circ));
end
% % 圆度计算函数
% function circ = calculate_circularity(points)
% % 计算轮廓圆度 (4π*面积/周长?0?5)
%
% % 1. 排序点以形成多边形
% centroid = mean(points, 1);
% angles = atan2(points(:,2) - centroid(2), points(:,1) - centroid(1));
% [~, idx] = sort(angles);
% sorted_points = points(idx, :);
%
% % 2. 计算多边形面积
% x = sorted_points(:,1);
% y = sorted_points(:,2);
% area = polyarea(x, y);
%
% % 3. 计算多边形周长
% perim = 0;
% n = size(sorted_points, 1);
% for i = 1:n-1
% perim = perim + norm(sorted_points(i,:) - sorted_points(i+1,:));
% end
% perim = perim + norm(sorted_points(end,:) - sorted_points(1,:));
%
% % 4. 计算圆度
% if perim > 0
% %圆度 = (4π × 面积) / (周长 ?)
% circ = (4 * pi * area) / (perim * perim);
% else
% circ = 0; % 避免除零错误
% end
%
% % 5. 处理特殊情况
% if circ > 1 || isnan(circ) || isinf(circ)
% circ = 0;
% end
% end
% function [A,B,C,D,E,F, center, major_axis, minor_axis, orientation] = ellipse_fit(x, y)
% % 椭圆拟合:输入点坐标x,y,输出椭圆参数
% n = length(x);
% if n < 5
% error('椭圆拟合需要至少5个点');
% end
%
% % 构造最小二乘矩阵
% D = zeros(n, 6);
% for i = 1:n
% D(i, :) = [x(i)^2, x(i)*y(i), y(i)^2, x(i), y(i), 1];
% end
%
% % 求解最小二乘问题 (D'*D)*params = 0,取最小特征值对应的特征向量
% [V, ~] = eig(D'*D);
% [~, idx] = min(diag(V));
% params = V(:, idx);
% A = params(1);
% B = params(2);
% C = params(3);
% D = params(4);
% E = params(5);
% F = params(6);
%
% % 计算椭圆中心、长短轴、方向
% % 中心 (xc, yc)
% denominator = 2*(A*C - (B/2)^2);
% xc = (B*E - 2*C*D) / denominator;
% yc = (B*D - 2*A*E) / denominator;
% center = [xc, yc];
%
% % 标准化椭圆方程,计算长短轴
% x_shifted = x - xc;
% y_shifted = y - yc;
% a = A*xc^2 + B*xc*yc + C*yc^2 + D*xc + E*yc + F;
% b = 2*A*xc + B*yc + D;
% c = 2*C*yc + B*xc + E;
% d = A;
% e = B/2;
% f = C;
%
% % 计算长短轴长度
% theta = 0.5*atan2(B, A - C);
% cos_theta = cos(theta);
% sin_theta = sin(theta);
% a_rot = a / (d*cos_theta^2 + 2*e*cos_theta*sin_theta + f*sin_theta^2);
% b_rot = a / (d*sin_theta^2 - 2*e*cos_theta*sin_theta + f*cos_theta^2);
% major_axis = max(sqrt(a_rot), sqrt(b_rot));
% minor_axis = min(sqrt(a_rot), sqrt(b_rot));
% orientation = theta;
% end
function [A,B,C,D,E,F, center, major_axis, minor_axis, orientation] = ellipse_fit(x, y)
% 椭圆拟合(改进版):基于最小二乘法,提高数值稳定性
% 输入:x, y - 椭圆边缘点的坐标(列向量或行向量)
% 输出:
% A,B,C,D,E,F - 椭圆方程参数 (Ax?+Bxy+Cy?+Dx+Ey+F=0)
% center - 椭圆中心 [xc, yc]
% major_axis - 长轴长度
% minor_axis - 短轴长度
% orientation - 椭圆旋转角(弧度,范围 [-π/2, π/2])
% 1. 输入处理:确保为列向量,去除重复点
x = x(:); % 转为列向量
y = y(:);
if length(x) ~= length(y)
error('x和y的长度必须一致');
end
if length(x) < 5
error('椭圆拟合至少需要5个点');
end
% 去除重复点(避免矩阵奇异)
points = unique([x, y], 'rows');
x = points(:,1);
y = points(:,2);
n = length(x);
if n < 5
error('去重后点数量不足5个,无法拟合椭圆');
end
% 2. 构造椭圆方程的最小二乘矩阵
% 椭圆方程:Ax? + Bxy + Cy? + Dx + Ey + F = 0
% 约束:A*C - (B/2)? > 0(确保是椭圆而非其他二次曲线)
M = [x.^2, x.*y, y.^2, x, y, ones(n,1)]; % n×6矩阵
% 3. 求解最小二乘问题(避免矩阵奇异)
% 使用SVD分解求解,比直接 eig 更稳定
[U, S, V] = svd(M, 0); % 经济型SVD分解
params = V(:, end); % 最小奇异值对应的特征向量(解)
A = params(1);
B = params(2);
C = params(3);
D = params(4);
E = params(5);
F = params(6);
% 4. 强制约束为椭圆(修正非椭圆情况)
ellipse_check = A*C - (B/2)^2; % 椭圆判别式
if ellipse_check <= 0
% 若为抛物线/双曲线,添加微小扰动修正为椭圆
perturb = 1e-8;
A = A + perturb;
C = C + perturb;
ellipse_check = A*C - (B/2)^2;
if ellipse_check <= 0
error('无法将点集修正为椭圆,可能点集非椭圆');
end
end
% 5. 计算椭圆中心 (xc, yc)
% 计算椭圆中心(修改这部分)
denominator = 2*(A*C - (B/2)^2); % 原分母计算
% 新增:检查分母是否接近零(增加保护机制)
if abs(denominator) < 1e-10 % 阈值可根据数据调整
% 方案1:尝试修正点集(去除离群点后重新拟合)
points = [x, y];
% 计算点到质心的距离,去除离群点(超过均值2倍的点)
centroid = mean(points);
dists = sqrt(sum((points - centroid).^2, 2));
valid = dists < 2*mean(dists); % 保留合理点
if sum(valid) >=5 % 修正后点数仍足够
x = x(valid);
y = y(valid);
% 重新构造矩阵求解
M = [x.^2, x.*y, y.^2, x, y, ones(length(x),1)];
[U, S, V] = svd(M, 0);
params = V(:, end);
A = params(1); B = params(2); C = params(3);
D = params(4); E = params(5); F = params(6);
denominator = 2*(A*C - (B/2)^2); % 重新计算分母
else
% 方案2:若无法修正,直接报错(被try-catch捕获)
error('点集不符合椭圆特征,无法计算中心');
end
end
% 计算中心(确保分母安全)
xc = (B*E - 2*C*D) / denominator;
yc = (B*D - 2*A*E) / denominator;
center = [xc, yc];
% 6. 计算旋转角(椭圆长轴方向)
if A ~= C
orientation = 0.5 * atan2(B, A - C); % 弧度
else
orientation = 0.5 * atan2(B, 1e-12); % 避免A=C时分母为零
end
% 限制角度在 [-π/2, π/2]
if orientation > pi/2
orientation = orientation - pi;
elseif orientation < -pi/2
orientation = orientation + pi;
end
% 7. 计算长短轴长度
% 平移坐标到中心
x_shift = x - xc;
y_shift = y - yc;
% 旋转坐标(消除旋转影响)
cos_theta = cos(orientation);
sin_theta = sin(orientation);
x_rot = x_shift .* cos_theta + y_shift .* sin_theta; % 旋转后的x
y_rot = -x_shift .* sin_theta + y_shift .* cos_theta; % 旋转后的y
% 拟合轴对齐椭圆(x'^2/a? + y'^2/b? = 1)
% 构造方程:x'^2*b? + y'^2*a? = a?*b?
M_rot = [x_rot.^2, y_rot.^2, ones(n,1)];
b_rot = zeros(n,1); % 右边项为0(移项后)
params_rot = M_rot \ b_rot; % 求解 [b?, a?, -a?b?]^T
b_sq = params_rot(1);
a_sq = params_rot(2);
ab_sq = -params_rot(3);
% 计算长短轴(确保为正值)
a = sqrt(abs(a_sq)); % 半长轴
b = sqrt(abs(b_sq)); % 半短轴
% 确保长轴 >= 短轴
if a < b
[a, b] = swap(a, b);
end
major_axis = 2*a; % 长轴长度(直径)
minor_axis = 2*b; % 短轴长度(直径)
% 8. 最终检查(避免异常值)
if isnan(major_axis) || isnan(minor_axis) || major_axis <= 0 || minor_axis <= 0
error('椭圆轴长计算异常,可能点集不符合椭圆特征');
end
end
% 辅助函数:交换两个变量
function [a, b] = swap(a, b)
temp = a;
a = b;
b = temp;
end
% 1. 计算轮廓点曲率(用于区分圆弧和直线)
function k = curvature(points)
n = size(points, 1);
k = zeros(n, 1); % 曲率存储每个点的曲率
if n < 3
return; % 点数不足,返回0
end
for i = 2:n-1
% 取当前点和前后点计算曲率
p_prev = points(i-1, :);
p_curr = points(i, :);
p_next = points(i+1, :);
% 向量计算
v1 = p_curr - p_prev;
v2 = p_next - p_curr;
% 曲率近似公式:|v1 × v2| / (|v1| * |v2| * 平均距离)
cross_prod = v1(1)*v2(2) - v1(2)*v2(1); % 叉积
norm_v1 = norm(v1);
norm_v2 = norm(v2);
if norm_v1 > 1e-6 && norm_v2 > 1e-6 % 避免除以0
avg_dist = (norm_v1 + norm_v2) / 2;
k(i) = abs(cross_prod) / (norm_v1 * norm_v2 * avg_dist);
end
end
% 首尾点曲率用相邻点近似
k(1) = k(2);
k(end) = k(end-1);
end
% 2. 分割轮廓为圆弧段和直线段
function [arc_segments, line_segments] = segment_contour(points, curvature_thresh)
arc_segments = {}; % 存储圆弧段(cell数组,每个元素是一段点)
line_segments = {}; % 存储直线段
if size(points, 1) < 3
return; % 点数不足,不分割
end
% 计算每个点的曲率
k = curvature(points);
% 标记:曲率 > 阈值 → 圆弧;否则 → 直线
is_arc = k > curvature_thresh;
% 遍历点,分割连续的圆弧/直线区域
current_segment = points(1, :); % 当前段的起点
current_type = is_arc(1); % 当前段类型(圆弧/直线)
for i = 2:length(is_arc)
if is_arc(i) == current_type
% 同类型,继续添加到当前段
current_segment = [current_segment; points(i, :)];
else
% 类型切换,保存当前段
if current_type
arc_segments{end+1} = current_segment;
else
line_segments{end+1} = current_segment;
end
% 开始新段
current_segment = points(i, :);
current_type = is_arc(i);
end
end
% 保存最后一段
if current_type
arc_segments{end+1} = current_segment;
else
line_segments{end+1} = current_segment;
end
% 过滤过短的段(少于3个点的段不拟合)
arc_segments = arc_segments(cellfun(@(x) size(x,1)>=3, arc_segments));
line_segments = line_segments(cellfun(@(x) size(x,1)>=2, line_segments));
end
% 3. 圆弧段拟合为圆
function [xc, yc, r] = circle_fit_segment(points)
x = points(:, 1);
y = points(:, 2);
n = length(x);
% 最小二乘圆拟合
A = [2*x, 2*y, ones(n, 1)];
b = x.^2 + y.^2;
params = A \ b;
xc = params(1); % 圆心x
yc = params(2); % 圆心y
r = sqrt(params(3) + xc^2 + yc^2); % 半径
end
% 4. 直线段拟合为直线
function [k, b] = line_fit_segment(points)
x = points(:, 1);
y = points(:, 2);
n = length(x);
% 最小二乘直线拟合(y = kx + b)
A = [x, ones(n, 1)];
params = A \ y;
k = params(1); % 斜率
b = params(2); % 截距
end
end 警告: 轮廓 8 椭圆拟合失败: 椭圆轴长计算异常,可能点集不符合椭圆特征
> 位置:zernike_uint5zuizhong (第 617 行)