% ========== 步骤7: 对每条裂隙轮廓提取并计算 JRC ==========
fprintf('📐 开始计算每条裂隙的 JRC 粗糙度...\n');
jrcResults = struct(); % 存储每条裂隙的不同采样方法下的 JRC
for fissureIdx = 1:length(fissureResults)
fprintf('🔍 处理第 %d 条裂隙...\n', fissureIdx);
% 获取该裂隙内点
% 获取该裂隙内点
x_inlier = fissureResults(fissureIdx).inliers_x;
y_inlier = fissureResults(fissureIdx).inliers_y;
% === 步骤1: 对每个 x 位置的多个 y 取平均 ===
[unique_x, ~, idx] = unique(round(x_inlier));
mean_y = accumarray(idx, y_inlier, [], @mean);
% 排序
[unique_x, sort_idx] = sort(unique_x);
mean_y = mean_y(sort_idx);
X_raw = unique_x;
Y_raw = mean_y;
% === 步骤2: 高密度插值生成平滑轮廓 ===
if length(X_raw) < 2
fprintf('⚠️ 第%d条裂隙有效点不足，跳过。\n', fissureIdx);
continue;
end
% 生成等间距的查询点
xq = linspace(min(X_raw), max(X_raw), 1000);
% 使用 spline 插值（或 pchip 防止过冲）
try
Y_interp = interp1(X_raw, Y_raw, xq, 'pchip', 'extrap');
catch ME
fprintf('❌ 插值失败：%s\n', ME.message);
continue;
end
% 去除 NaN
valid = ~isnan(Y_interp);
X_dense = xq(valid);
Y_dense = Y_interp(valid);
% 提取轮廓：使用 Canny 边缘检测从当前裂隙掩膜中获取清晰边界
mask_current = false(height, width);
R_fit = fissureResults(fissureIdx).params(1);
P_fit = fissureResults(fissureIdx).params(2);
beta_fit = fissureResults(fissureIdx).params(3);
C_fit = fissureResults(fissureIdx).params(4);
x_all = 1:width;
y_curve = R_fit * sin(2*pi/P_fit*(x_all - beta_fit)) + C_fit;
half_width = 30;
for x = 1:width
yc = y_curve(x);
if isnan(yc) || yc < 1 || yc > height, continue; end
y_low = max(1, floor(yc - half_width));
y_high = min(height, ceil(yc + half_width));
mask_current(y_low:y_high, x) = true;
end
se = strel('disk', 3);
mask_current = imclose(double(mask_current), se) > 0;
% 使用 canny 边缘检测提取上下两个边界
edge_map = edge(uint8(mask_current), 'canny', [], 0.1);
[y_edge, x_edge] = find(edge_map);
if isempty(x_edge)
fprintf('⚠️ 第%d条裂隙未提取到边缘！跳过。\n', fissureIdx);
continue;
end
% 将边缘点聚类为上边缘和下边缘（根据纵坐标中位数划分）
midline_y = median([fissureResults(fissureIdx).inliers_y]);
upper_edge = (y_edge < midline_y);
lower_edge = ~upper_edge;
% 取其中一个边缘作为代表（比如下边缘更完整）
x_contour = x_edge(lower_edge);
y_contour = y_edge(lower_edge);
if length(x_contour) < 10
fprintf('⚠️ 边缘点太少，跳过第%d条裂隙。\n', fissureIdx);
continue;
end
% 按 x 排序
[~, idx_sort] = sort(x_contour);
x_contour = x_contour(idx_sort);
y_contour = y_contour(idx_sort);
% 平滑轮廓（可选）
y_contour = smoothdata(y_contour, 'gaussian', 5);
% 准备三种采样方式
num_samples = 100; % 统一采样点数量便于比较
methods = {'uniform', 'adaptive', 'random'};
jrcValues = zeros(1, 3);
for m = 1:3
method = methods{m};
x_sampled = [];
y_sampled = [];
switch method
case 'uniform'
% 等间隔采样 x
idx_uniform = round(linspace(1, length(x_contour), num_samples));
idx_uniform = unique(idx_uniform);
x_sampled = x_contour(idx_uniform);
y_sampled = y_contour(idx_uniform);
case 'adaptive'
% 自适应采样：曲率大处多采样
dx = diff(x_contour);
dy = diff(y_contour);
ddx = diff([0, dx]); % 近似二阶导
ddy = diff([0, dy]);
% 计算每个点的曲率近似 |d²y/dx²| / (1 + (dy/dx)^2)^(3/2)
slopes = dy ./ (dx + eps);
curvature = abs([0, diff(slopes)]) ./ (1 + [slopes(1:end-1).^2, 0]).^1.5;
% 根据曲率加权概率分布采样
prob = curvature + 0.01; % 加小常数避免零权重
cumprob = cumsum(prob);
total = cumprob(end);
rand_vals = rand(1, num_samples) * total;
idx_adaptive = interp1(cumprob, 1:length(cumprob), rand_vals, 'linear', 'nearest');
idx_adaptive = unique(round(idx_adaptive));
x_sampled = x_contour(idx_adaptive);
y_sampled = y_contour(idx_adaptive);
case 'random'
% 随机采样
idx_rand = randperm(length(x_contour), min(num_samples, length(x_contour)));
idx_rand = sort(idx_rand); % 保持顺序
x_sampled = x_contour(idx_rand);
y_sampled = y_contour(idx_rand);
end
% 计算 Z 参数
dx_samp = diff(x_sampled);
dy_samp = diff(y_sampled);
N = length(dx_samp);
if N == 0
Z = 0;
else
Z = sqrt(sum((dy_samp ./ (dx_samp + eps)).^2) / N);
end
% 计算 JRC
JRC = 51.85 * (Z^0.6) - 10.37;
JRC = max(JRC, 0); % 防止负值
jrcValues(m) = JRC;
fprintf('📊 [%s 采样] Z=%.4f, JRC=%.2f\n', method, Z, JRC);
end
% 保存结果
jrcResults(fissureIdx).x_contour = x_contour;
jrcResults(fissureIdx).y_contour = y_contour;
jrcResults(fissureIdx).JRC_uniform = jrcValues(1);
jrcResults(fissureIdx).JRC_adaptive = jrcValues(2);
jrcResults(fissureIdx).JRC_random = jrcValues(3);
jrcResults(fissureIdx).mean_JRC = mean(jrcValues);
end
% 输出总体统计
all_jrc = [jrcResults.JRC_uniform, jrcResults.JRC_adaptive, jrcResults.JRC_random];
fprintf('✅ 所有裂隙 JRC 计算完成！平均 JRC ≈ %.2f\n', mean(all_jrc));
修改语法错误
扫一扫
390
被折叠的 条评论
为什么被折叠?
到【灌水乐园】发言
