XL2 &#wbDTypeV2ObjIDDDirPL1zBEcAlgoEcMEcNEcBSizeEcIndexEcDistCSumAlgoPartNumsPartETags 211a14125370b7d28e56d66cd2e2b99dPartSizes"PartASizes"PartIdxSize"MTime#wbDMetaSysx-minio-internal-inline-datatruex-rustfs-internal-inline-datatrueMetaUsrcontent-typeapplication/octet-streametag 211a14125370b7d28e56d66cd2e2b99dv\null" s c=ݩ4JB,żk]clc; clear; close all; %% ========================================================= % 二维两步基线校正 % 第1步：竖向(Y方向) 小窗口 block-min + pchip % 第2步：横向(X方向) 大窗口 block-min + pchip %% ========================================================= %% 1. 载入数据 load('.\data\pc1.mat'); assert(exist('pc1','var')==1, 'pc1.mat 中没有变量 pc1'); assert(isa(pc1, 'pointCloud'), 'pc1 必须是 pointCloud 对象'); P = double(pc1.Location); assert(size(P,2)==3, 'pc1.Location 必须是 N x 3'); Xw = P(:,1); Yw = P(:,2); Zw = P(:,3); % valid = isfinite(Xw) & isfinite(Yw) & isfinite(Zw); valid = isfinite(Xw) & isfinite(Yw) & isfinite(Zw) ... & ~(Xw == 0 & Yw == 0 & Zw == 0); Xw = Xw(valid); Yw = Yw(valid); Zw = Zw(valid); fprintf('有效点数: %d\n', numel(Xw)); %% 2. 参数设置 scanStep = 0.1; % Y方向扫描步进 xGridStep = 0.2; % X方向网格步长 % ===== 第1遍：竖向 ===== blockWidthY = 1; % 小一些 minPtsPerBlockY = 2; % ===== 第2遍：横向 ===== blockWidthX = 5; % 大一些 minPtsPerBlockX = 3; minPtsPerColumn = 5; minPtsPerRow = 10; showSampleProfiles = true; %% 3. 构造规则二维高度图 Zmap(y,x) Yq = round(Yw / scanStep) * scanStep; yGrid = unique(Yq(:)); % 列向量 ny = numel(yGrid); xMin = min(Xw); xMax = max(Xw); xGrid = (xMin:xGridStep:xMax); % 行向量 nx = numel(xGrid); fprintf('构造规则网格: ny=%d, nx=%d\n', ny, nx); % ---- 把原始点映射到二维网格 ---- xBin = discretize(Xw, [xGrid, xGrid(end)+xGridStep]); % yGrid 是严格按 0.1 间隔来的，直接用索引公式更稳 yBin = round((Yq - yGrid(1)) / scanStep) + 1; validBin = ~isnan(xBin) & (yBin >= 1) & (yBin <= ny); xBin2 = xBin(validBin); yBin2 = yBin(validBin); zUse = Zw(validBin); linIdx = sub2ind([ny, nx], yBin2, xBin2); [uniqIdx, ~, ic] = unique(linIdx); zMean = accumarray(ic, zUse, [], @mean, NaN); Zmap = nan(ny, nx); Zmap(uniqIdx) = zMean; % 补少量空洞 Zmap = fillmissing(Zmap, 'nearest', 2); Zmap = fillmissing(Zmap, 'nearest', 1); %% 4. 第1遍：先对每条扫描线（按行）做 block-min + pchip fprintf('开始第1遍：逐条扫描线基线...\n'); B1 = nan(ny, nx); % 每一行就是一条扫描线，所以这里按行循环 tic; for i = 1:ny zrow = Zmap(i, :); xrow = xGrid; validRow = isfinite(zrow); if nnz(validRow) < minPtsPerRow continue; end B1(i, :) = baseline_block_min_pchip_1d( ... xrow(validRow)', zrow(validRow)', xGrid', ... blockWidthX, minPtsPerBlockX)'; end toc; B1 = fillmissing(B1, 'nearest', 2); B1 = fillmissing(B1, 'nearest', 1); R1 = Zmap - B1; %% 5. 第2遍：再沿跨扫描线方向（按列）在 R1 上做 block-min + pchip fprintf('开始第2遍：跨扫描线方向基线...\n'); B2 = nan(ny, nx); % 这里才是按列循环 for j = 1:nx zcol = R1(:, j); ycol = yGrid; validCol = isfinite(zcol); if nnz(validCol) < minPtsPerColumn continue; end B2(:, j) = baseline_block_min_pchip_1d( ... ycol(validCol), zcol(validCol), ycol, ... blockWidthY, minPtsPerBlockY); end B2 = fillmissing(B2, 'nearest', 2); B2 = fillmissing(B2, 'nearest', 1); %% 6. 合成总基线与最终残差 baselineTotal = B1 + B2; detailMap = Zmap - baselineTotal; %% 7. 把二维总基线插值回原始点 fprintf('开始把基线插值回原始点云...\n'); [XG, YG] = meshgrid(xGrid, yGrid); baseline_each_pt = interp2(XG, YG, baselineTotal, Xw, Yq, 'linear', NaN); validInterp = isfinite(baseline_each_pt); detail_each_pt = nan(size(Zw)); detail_each_pt(validInterp) = Zw(validInterp) - baseline_each_pt(validInterp); X_final = Xw(validInterp); Y_final = Yw(validInterp); detail_final = detail_each_pt(validInterp); xyz_detail_final = [X_final, Y_final, detail_final]; dMin = min(detail_final); dMax = max(detail_final); D8 = uint8(255 * (detail_final - dMin) / (dMax - dMin + eps)); C_detail = repmat(D8, 1, 3); pc_detail_final = pointCloud(xyz_detail_final, 'Color', C_detail); %% 8. 可视化 ind = find(xyz_detail_final(:,3)>-0.1); k = xyz_detail_final; k = k(ind,:); ind = find(k(:,3)<1.5); k = k(ind,:); figure,pcshow(k); figure('Color','w','Name','原始 Zmap'); imagesc(xGrid, yGrid, Zmap); axis image; set(gca, 'YDir', 'normal'); colormap turbo; colorbar; xlabel('X'); ylabel('Y'); title('原始规则高度图 Zmap'); figure('Color','w','Name','第1遍竖向基线 Bv'); imagesc(xGrid, yGrid, Bv); axis image; set(gca, 'YDir', 'normal'); colormap turbo; colorbar; xlabel('X'); ylabel('Y'); title('第1遍：竖向基线 Bv'); figure('Color','w','Name','第一次残差 R1'); imagesc(xGrid, yGrid, R1); axis image; set(gca, 'YDir', 'normal'); colormap turbo; colorbar; xlabel('X'); ylabel('Y'); title('第一次残差 R1 = Zmap - Bv'); figure('Color','w','Name','第2遍横向基线 Bh'); imagesc(xGrid, yGrid, Bh); axis image; set(gca, 'YDir', 'normal'); colormap turbo; colorbar; xlabel('X'); ylabel('Y'); title('第2遍：横向基线 Bh'); figure('Color','w','Name','最终 detailMap'); imagesc(xGrid, yGrid, detailMap); axis image; set(gca, 'YDir', 'normal'); colormap turbo; colorbar; xlabel('X'); ylabel('Y'); title('最终残差图 detailMap'); figure('Color','w','Name','最终 detail 点云'); showN = min(500000, size(xyz_detail_final,1)); rng(1); idxShow = randperm(size(xyz_detail_final,1), showN); pcshow(pointCloud(xyz_detail_final(idxShow,:), 'Color', C_detail(idxShow,:))); xlabel('X'); ylabel('Y'); zlabel('detail'); title('所有处理后点的点云 pc_detail_final'); if showSampleProfiles rowIdx = round(ny/2); colIdx = round(nx/2); figure('Color','w','Name','样例行剖面'); plot(xGrid, Zmap(rowIdx,:), 'k-', 'LineWidth', 1.2); hold on; grid on; plot(xGrid, baselineTotal(rowIdx,:), 'r-', 'LineWidth', 2); plot(xGrid, detailMap(rowIdx,:), 'b-', 'LineWidth', 1); legend({'原始Z','总基线','最终残差'}); xlabel('X'); ylabel('Z / detail'); title(sprintf('样例行剖面, rowIdx=%d', rowIdx)); figure('Color','w','Name','样例列剖面'); plot(yGrid, Zmap(:,colIdx), 'k-', 'LineWidth', 1.2); hold on; grid on; plot(yGrid, Bv(:,colIdx), 'r-', 'LineWidth', 2); plot(yGrid, R1(:,colIdx), 'b-', 'LineWidth', 1); legend({'原始Z','竖向基线','第一次残差'}); xlabel('Y'); ylabel('Z / detail'); title(sprintf('样例列剖面, colIdx=%d', colIdx)); end %% 9. 保存结果 result = struct(); result.scanStep = scanStep; result.xGridStep = xGridStep; result.blockWidthY = blockWidthY; result.minPtsPerBlockY = minPtsPerBlockY; result.blockWidthX = blockWidthX; result.minPtsPerBlockX = minPtsPerBlockX; result.xGrid = xGrid; result.yGrid = yGrid; result.Zmap = Zmap; result.Bv = Bv; result.R1 = R1; result.Bh = Bh; result.baselineTotal = baselineTotal; result.detailMap = detailMap; result.xyz_detail_final = xyz_detail_final; result.pc_detail_final = pc_detail_final; save('two_pass_block_min_pchip_result.mat', 'result', '-v7.3'); fprintf('结果已保存到 two_pass_block_min_pchip_result.mat\n'); fprintf('处理完成。\n'); %% ========================================================= % 局部函数：一维 block-min + pchip %% ========================================================= function baseline = baseline_block_min_pchip_1d(s, z, sQuery, blockWidth, minPtsPerBlock) s = s(:); z = z(:); sQuery = sQuery(:); valid = isfinite(s) & isfinite(z); s = s(valid); z = z(valid); baseline = nan(size(sQuery)); if numel(s) < 2 return; end [s, idx] = sort(s); z = z(idx); sMin = min(s); sMax = max(s); edges = sMin:blockWidth:sMax; if numel(edges) < 2 || edges(end) < sMax edges = [edges, sMax]; end numBlocks = numel(edges) - 1; sc = []; zb = []; for k = 1:numBlocks sL = edges(k); sR = edges(k+1); if k < numBlocks idxk = (s >= sL) & (s < sR); else idxk = (s >= sL) & (s <= sR); end if nnz(idxk) < minPtsPerBlock continue; end sc(end+1,1) = 0.5 * (sL + sR); %#ok % zb(end+1,1) = min(z(idxk)); %#ok zb(end+1,1) = prctile(z(idxk), 10); end if numel(sc) < 2 return; end inRange = (sQuery >= min(sc)) & (sQuery <= max(sc)); baseline(inRange) = interp1(sc, zb, sQuery(inRange), 'pchip'); baseline(sQuery < min(sc)) = zb(1); baseline(sQuery > max(sc)) = zb(end); end