“function main_q3_nopad()
% Q3：复杂裂隙定量建模 + JRC（零 padarray / 零 imfilter 版）
% 输入：C:\Users\Administrator\Desktop\C题\附件3\
% 输出：<base>\out_q3\ q3_results.csv 及关键图像叠加
%% ---------- 固定路径 ----------
baseDir = 'C:\Users\Administrator\Desktop\C题\附件3';
outDir = fullfile(baseDir,'out_q3'); if ~exist(outDir,'dir'), mkdir(outDir); end
fprintf('[Q3] baseDir = %s\n', baseDir);
%% ---------- 常量 ----------
P0 = pi*30; % 钻孔周长(mm)，直径30 mm
segDepth_mm = 500; % 单张图对应孔深(mm)（题面给定）
% 预处理/分割参数（轻量、不用形态学开闭）
cfg.gamma = 0.85; % γ增强
cfg.bgLen = 61; % 纵向均值背景（conv2）
cfg.minArea = 120; % 最小连通域像素
cfg.minRows = 40; % 至少跨越的行数（抑制小噪声）
%
拟合&JRC
fft_peak_guard = [0.6 1.6]; % 允许的 P/P0 范围，用于 FFT 初值（复杂裂隙允许偏离）
jrc_func = @(Z2) 32.2 + 32.47*log10(max(Z2,1e-8)); % 经验式
%% ---------- 遍历图像 ----------
IM = dir(fullfile(baseDir, '*.*'));
IM = IM(~[IM.isdir]);
IM = IM( endsWith({IM.name},{'.jpg','.png','.jpeg','.bmp','.tif','.tiff'},'IgnoreCase',true) );
% 按文件名中的“图3-数字”排序（若无则字典序）
ord = regexp({IM.name},'图?3[-_－—]?\s*(\d+)','tokens','once');
tryNum = nan(numel(IM),1);
for i=1:numel(ord), if ~isempty(ord{i}), tryNum(i)=str2double(ord{i}{1}); end, end
[~,sidx] = sortrows([isnan(tryNum) tryNum(:)]);
IM = IM(sidx);
rows_all = {};
for k = 1:numel(IM)
fpath = fullfile(IM(k).folder, IM(k).name);
I0 = im2double(imread(fpath));
if size(I0,3)==3, I0 = rgb2gray(I0); end
[H,W] = size(I0);
sx = P0 / max(W,1); % x: mm/px（周向）
sy = segDepth_mm / max(H,1); % y: mm/px（轴向）
% -------- 预处理（全部无形态学） --------
I = I0.^cfg.gamma; % γ增强
kerV = ones(cfg.bgLen,1) / cfg.bgLen; % 纵向均值背景
bg = conv2(I, kerV, 'same');
Ienh = clip01(I - bg);
% --------
边缘/对比度响应 + Otsu --------
Gx = [diff(Ienh,1,2), zeros(H,1)];
Gy = [diff(Ienh,1,1); zeros(1,W)];
R = clip01( sqrt(Gx.^2 + Gy.^2) );
tt = safe_otsu(R) * 0.95;
BW = R >= max(tt, 0.15);
% 去小块（不依赖 bwareaopen）
BW = area_open(BW, cfg.minArea);
% -------- 取主要裂隙连通域（可能多条） --------
CC = bwconncomp(BW,8);
if CC.NumObjects==0, continue; end
areas = cellfun(@numel, CC.PixelIdxList);
% 只保留面积前 6 条，避免过多噪声
[~,aord] = sort(areas,'descend');
keepIdx = aord(1:min(6,numel(aord)));
for ci = keepIdx(:).'
pix = CC.PixelIdxList{ci};
[rr,cc] = ind2sub([H W], pix);
% 行聚类：要求跨越的行数达到阈值
if numel(unique(rr)) < cfg.minRows, continue; end
% 行中心线：每一行的列均值（再平滑）
x_center = nan(H,1);
urows = unique(rr).';
for r0 = urows
cols = cc(rr==r0);
x_center(r0) = mean(cols);
end
sel = ~isnan(x_center);
y_px = find(sel).'; x_px = x_center(sel).';
if numel(x_px) < 40, continue; end
% 平滑中心线
x_px = smooth(x_px, 9, '
lowess').';
% ---- 像素->物理量 ----
x_mm = (x_px-1) * sx; % 0..P0
y_mm = (y_px-1) * sy; % 0..segDepth
% ---- 复杂裂隙的“等效正弦”
拟合 ----
% 初值：用 FFT 找主频，约束 P 在 [0.6,1.6]P0 内
[P_init, beta_init] = guess_period_phase_fft(x_mm, y_mm, P0, fft_peak_guard);
C_init = mean(y_mm);
R_init = 0.5*(max(y_mm)-min(y_mm));
p0 = [R_init, P_init, beta_init, C_init];
lb = [0.0, 0.6*P0, -2*pi, min(y_mm)-50];
ub = [segDepth_mm, 1.6*P0, 2*pi, max(y_mm)+50];
[Rfit, Pfit, betafit, Cfit] = fit_sine_firstharm(x_mm, y_mm, p0, lb, ub);
% ---- 生成
离散轮廓 & 计算 JRC（多采样方案对比） ----
% 1) 等间距采样（沿 x）
N1 = 300; [x1,y1] = resample_profile(x_mm, y_mm, N1);
Z2_1 = z2_metric(x1,y1,P0);
JRC1 = jrc_func(Z2_1);
% 2) 曲率自适应（在 |Δ²y| 大处加密）
[x2,y2] = adaptive_resample_curvature(x_mm,y_mm, N1);
Z2_2 = z2_metric(x2,y2,P0);
JRC2 = jrc_func(Z2_2);
% 3) 面积加权（按局部面积密度分配点数；更关注大幅度段）
[x3,y3] = area_weighted_resample(x_mm,y_mm, N1);
Z2_3 = z2_metric(x3,y3,P0);
JRC3 = jrc_func(Z2_3);
% 推荐报告值（默认用等间距）
JRC = JRC1;
% ---- 记录结果 ----
rows_all(end+1,:) = {IM(k).name, ci, Rfit, Pfit, betafit, Cfit, ...
JRC, JRC1, JRC2, JRC3, Z2_1, Z2_2, Z2_3}; %#ok<AGROW>
% ---- 关键图像（只对“图3-1/2/3”保存叠加） ----
if any(contains(IM(k).name, {'图3-1','图3-2','图3-3'}))
hfig = figure('Color','w','Visible','off');
subplot(1,2,1);
imshow(I0,[]); hold on;
plot(x_px, y_px, 'r.', 'MarkerSize', 6); title([IM(k).name '：中心线']);
subplot(1,2,2);
plot(x_mm, y_mm,'k.'); hold on;
xf = linspace(min(x_mm), max(x_mm), 800);
yf = Cfit + Rfit * sin(2*pi*xf/Pfit + betafit);
plot(xf, yf, 'r-', 'LineWidth', 1.5);
legend('数据','等效正弦','Location','best'); grid on;
xlabel('x (mm)'); ylabel('y (mm)');
title(sprintf('R=%.2f, P=%.2f, \\beta=%.2f, C=%.2f | JRC=%.2f', ...
Rfit,Pfit,betafit,Cfit,JRC));
saveas(hfig, fullfile(outDir, [IM(k).name '_fit.png']));
close(hfig);
end
end
end
%% ---------- 汇总表 ----------
varNames = {'图像编号','裂隙编号','振幅R_mm','周期P_mm','相位beta_rad','中心线C_mm', ...
'JRC推荐','JRC_等间距','JRC_曲率自适应','JRC_面积加权', ...
'Z2_等间距','Z2_曲率自适应','Z2_面积加权'};
if isempty(rows_all)
T = table('Size',[0 numel(varNames)], ...
'VariableTypes', {'string','double','double','double','double','double', ...
'double','double','double','double','double','double','double'}, ...
'VariableNames', varNames);
else
T = cell2table(rows_all, 'VariableNames', varNames);
end
writetable(T, fullfile(outDir,'q3_results.csv'),'WriteMode','overwrite');
fprintf('[Q3] Done. Results -> %s\n', outDir);
end
%% ============== 工具函数（全部无 padarray / 无 imfilter） =================
function A = clip01(A), A = max(0,min(1,A)); end
function BW = area_open(BW, minA)
CC = bwconncomp(BW,8);
if CC.NumObjects==0, return; end
rm = false(CC.NumObjects,1);
for i=1:CC.NumObjects
if numel(CC.PixelIdxList{i}) < minA, rm(i)=true; end
end
if any(rm)
idx = cat(1, CC.PixelIdxList{rm});
BW(idx) = 0;
end
end
function t = safe_otsu(I)
try
t = graythresh(I);
catch
% 简单 Otsu 实现（直方图 256 级）
nbins = 256; h = imhist(I, nbins); p = h/sum(h);
omega = cumsum(p); mu = cumsum((0:nbins-1)'.*p);
mu_t = mu(end);
sigma2 = (mu_t*omega - mu).^2 ./ max(omega.*(1-omega),eps);
[~,idx] = max(sigma2); t = (idx-1)/(nbins-1);
end
end
% ---------- FFT 初值（估计 P 与 beta） ----------
function [P_est, beta_est] = guess_period_phase_fft(x, y, P0, guard)
% x: mm, y: mm；guard = [lo hi] 针对 P/P0 的约束
x = x(:); y = y(:);
% 先线性去趋势
p = polyfit(x,y,1); y_d = y - polyval(p,x);
% 频谱（均匀化采样）
Nx = 1024;
xu = linspace(min(x), max(x), Nx)';
yu = interp1(x, y_d, xu, 'linear', 'extrap');
Y = fft(yu - mean(yu));
f = (0:Nx-1)'/ (xu(end)-xu(1)); % 周期频率 (1/mm)
% 只取前半谱
half = 2:floor(Nx/2);
[~,pk] = max(abs(Y(half)));
f0 = f(half(pk));
P_est = 1/max(f0, 1e-6);
P_est = min(max(P_est, guard(1)*P0), guard(2)*P0);
% 相位：用最小二乘在已知 P 下解 beta 和 C、R
omega = 2*pi/P_est;
A = [sin(omega*x), cos(omega*x), ones(size(x))];
coef = A\y;
beta_est = atan2(coef(2), coef(1)); % 因 y = a*sin + b*cos
end
% ---------- 一阶谐波
拟合（P 可微调） ----------
function [Rfit,Pfit,betafit,Cfit] = fit_sine_firstharm(x,y,p0,lb,ub)
% 目标：y ≈ C + R*sin(2πx/P + β)
% 用简单的坐标下降：交替优化 {P} 与 {R,β,C}
x = x(:); y = y(:);
Rfit=p0(1); Pfit=p0(2); betafit=p0(3); Cfit=p0(4);
for it=1:15
% 1) 固定 P，线性最小二乘解 {R,β,C}
w = 2*pi/Pfit;
A = [sin(w*x), cos(w*x), ones(size(x))];
c = A\y; a=c(1); b=c(2); Cfit=c(3);
Rfit = hypot(a,b);
betafit = atan2(b, a); % y=a*sin+ b*cos = R*sin(w x + beta) with beta=atan2(b,a)
% 2) 微调 P（线性化一维搜索）
Pcand = linspace(max(lb(2),0.95*Pfit), min(ub(2),1.05*Pfit), 21);
err = zeros(size(Pcand));
for i=1:numel(Pcand)
wi = 2*pi/Pcand(i);
Ai = [sin(wi*x), cos(wi*x), ones(size(x))];
ci = Ai\y;
yi = Ai*ci;
err(i) = mean((y-yi).^2);
end
[~,m] = min(err); Pfit = Pcand(m);
end
end
% ---------- 等间距重采样 ----------
function [xr,yr] = resample_profile(x,y,N)
xr = linspace(min(x), max(x), N);
yr = interp1(x,y,xr,'linear','extrap');
end
% ---------- 曲率自适应重采样 ----------
function [xr,yr] = adaptive_resample_curvature(x,y,N)
x = x(:).'; y = y(:).';
dx = gradient(x); dy = gradient(y);
ddx = gradient(dx); ddy = gradient(dy);
kappa = abs(ddy)./( (dx.^2 + dy.^2).^(3/2) + eps ); % 近似曲率
w = kappa + 1e-3; w = w/sum(w);
s = [0, cumsum(w)]; s = s/s(end);
sQ = linspace(0,1,N);
xr = interp1(linspace(0,1,numel(x)), x, sQ,'linear','extrap');
yr = interp1(linspace(0,1,numel(y)), y, sQ,'linear','extrap');
end
% ---------- 面积加权重采样（|Δy| 大处加密） ----------
function [xr,yr] = area_weighted_resample(x,y,N)
x = x(:).'; y = y(:).';
w = [0 abs(diff(y))]; w = w + mean(w); w = w/sum(w);
s = [0, cumsum(w)]; s = s/s(end);
sQ = linspace(0,1,N);
xr = interp1(linspace(0,1,numel(x)), x, sQ,'linear','extrap');
yr = interp1(linspace(0,1,numel(y)), y, sQ,'linear','extrap');
end
% ---------- Z2 指标（按文中式(3)含 L 的常见写法） ----------
function Z2 = z2_metric(x,y,L)
% x,y 为等距/非等距
离散点（mm），L 取钻孔周长
dx = diff(x); dy = diff(y);
Z2 = sum( (dy.^2) ./ max(dx,1e-9) ) / max(L,1e-9);
end
”针对2025 年中国研究生数学建模竞赛 C 题文件内容对上述设计的代码方案进行优化提升，最后给出正确完整的代码
本文介绍了LOWESS平滑算法的实现原理与步骤。通过遍历散点集，使用局部加权最小二乘法拟合数据，进而实现线性和非线性数据的平滑处理。该文还详细说明了MATLAB中lowess函数的具体实现流程。
最低0.47元/天 解锁文章
扫一扫
2885
被折叠的 条评论
为什么被折叠?
到【灌水乐园】发言
