简介MATLAB空间雕刻三维重建项目面向计算机视觉与图像处理学习者围绕序列图像恢复物体三维几何这一核心问题提供从视图配准、遮挡分析到体素模型构建与后处理的完整实现参考。压缩包共99个文件包含jpg与ppm序列图像、m脚本、mat数据、txt说明及pdf资料约44.68MB便于按图像、代码和文档三类内容检索学习。已有787人学习下载说明其在相关课程实验与自学场景中具有一定参考价值。读者可结合恐龙模型的序列图像理解空间雕刻的可见性判定与体素雕刻流程并通过脚本和中间数据复现实验、观察参数变化对重建结果的影响积累三维建模与MATLAB编程的排错经验。1. 从一个空盒子开始SpaceCarving 为什么靠“删”而不是“建”完成三维重建拿到一叠标定好的序列图像想把里面的物体还原成三维模型多数人第一反应是找特征点、做匹配、三角化再拼点云。SpaceCarving 走的是完全相反的路先在空间里摆一个足够大的空盒子切成成千上万个体素然后逐张图去问一句——这个体素投到这张图上落不落在物体轮廓里只要有一张图回答“不”这个体素就被永久删掉。全程不做特征匹配也不估深度最后活下来的那坨体素就是模型。这套思路用 MATLAB 实现核心循环不到一百行却能稳定跑出带凹槽的实体模型。真正值得搞懂的不是代码而是它凭什么敢“删”——删错一次就没有后悔药它靠什么保证不删到真东西。2. 视觉外壳剪掉的是“不可能”留下的是“不能排除”2.1 轮廓锥的交集视觉外壳到底“包”住了什么把一张图的物体轮廓反投影回三维空间得到的不是一个面而是一个以相机光心为顶点、以轮廓曲线为准线的锥体——轮廓内部任意一个像素都对应一条从光心射出去的射线这些射线的并集就是锥体。物理上能确定的事情只有一件物体上任何一个真实存在的点它在每张图里的投影必然落在该图的轮廓内部因此必然落在这个锥体里。于是把所有相机的锥体求交集得到的区域一定包含真实物体。这个交集叫视觉外壳visual hullSpaceCarving 雕刻出来的就是它。注意这个包含关系是单向的视觉外壳 ⊇ 真实物体多出来的部分是“所有相机都没能排除掉”的空间不是错误而是信息不足的必然结果。整个算法的正确性就建立在这一条单向包含上——它只删除“至少被一张图证明是背景”的体素绝不凭空增加体素。这也解释了为什么它一旦删错就无法挽回。集合交集没有反操作一个体素被判死之后后面的任何一张图都无法把它救回来。常见的实现里这个“不可逆”被当成特性而不是缺陷只要轮廓是保守的宁可画大一点删除就永远不会伤到真实表面。反过来说如果某张图的轮廓分割把物体切掉了一块那块区域会被整条锥体贯穿式地削掉模型上就会出现一个平滑的切口这种错误在结果里非常显眼也极难修补。2.2 凹面去哪了视觉外壳的理论天花板视觉外壳的紧致程度完全取决于视角分布。对于凸的、或者接近广义圆柱的物体几十个视角就能贴得非常准对于凹陷情况完全不同。设想一个碗从正上方看轮廓是一个圆从侧面看轮廓是一个 U 形但碗的内壁在所有视角下都被碗沿挡住没有任何一张图的轮廓能“凹”进碗的内部。结果就是这个碗会被填成一个实心的半球壳块。这不是实现问题是原理问题。判断一个体素该不该保留只看它在图像上的投影是否落在轮廓内凹面内部的点投影恰好落在轮廓内于是被保留。要想把凹面挖出来必须引入轮廓之外的信息——像素颜色。这就是后面要讲的 photo-consistency 精修也是 SpaceCarving 这个名字里“carving”之外常被忽略的另一半。工程上有个实用推论如果你的物体是凸的或者只有浅凹陷纯轮廓法足够如果物体有明显的镂空、内腔、把手纯轮廓法必然失败别在分辨率和阈值上反复折腾那是玄学不是调参能解决的。2.3 体素化把连续的交集变成布尔数组连续空间里的交集没法直接算必须离散化。做法是在一个包围盒里均匀撒网格点每个点代表一个小立方体体素用体素是否“存活”来表示它是否属于视觉外壳。这里有三个必须想清楚的细节。第一是判定用体素中心还是体素整体。用中心点判定最简单一个体素投成一个像素查一次轮廓表就行。但这样会有半个体素量级的抖动靠近轮廓边界的体素中心点稍微偏一点判定就翻转。更保守的做法是把体素投成一个小四边形四个角都落在轮廓内部才保留代价是结果整体胖一圈。实际项目里常用折中方案——用轮廓的距离变换给出一个内缩余量体素中心到背景的距离大于某个半径才算数这样既保守又只查一次表。第二是分辨率与精度的关系。包围盒边长 L、每轴 N 个体素时体素边长是 L/N而体素对角线是 L·√3/N。模型表面的定位误差量级就在这个对角线上所以把 N 从 128 提到 256精度只提升一倍内存却涨八倍。这个账要在动手之前算清楚。第三是数据结构。三维逻辑数组logical(N,N,N)是最省事的选择占 N³ 字节128³ 才 2MB256³ 是 16MB都能接受。但如果为了加速把每个体素的齐次坐标都存成 double内存立刻爆掉——那是 4×N³ 个双精度数128³ 就是 64MB256³ 是 512MB再乘上分块拷贝的临时变量很容易把机器拖死。2.4 剪枝顺序交换律带来的确定性与效率差所有相机的判定是求交集而交集满足交换律和结合律所以处理顺序不影响最终结果。这一点在调试时非常有用如果你怀疑某张图有问题把它单独拿掉重跑结果的变化是可解释的不会因为顺序变化产生蝴蝶效应。但顺序极大影响速度。最朴素的做法是对每个相机遍历全部体素复杂度是 O(C·N³)如果先用手工指定的一两个“正面”视角粗剪一遍把大部分空体素干掉后续相机只需要处理存活体素实际计算量能掉一个数量级。更系统的做法是八叉树把包围盒当成一个根节点判定它八个角或整个块是否全在某个锥体外——全在外就整块删掉全在内部就整块保留待细分跨边界才继续分。% 由粗到细的八叉树式 carving 骨架carveBlock / refineBlock 为局部函数 occ true(N0, N0, N0); % 粗层N0 一般取 8 或 16 for L 1:log2(N / N0) occ carveBlock(occ, cams, bbox); % 在当前层对所有存活块做一次全量判定 occ refineBlock(occ); % 每个存活块细分成 2×2×2 个子块 endcarveBlock内部对每个块只做一次投影用块的中心和半径构造一个保守的包围圆去查轮廓——圆完全在轮廓外就整块删完全在内就整块留跨边界的标记为“待细分”。refineBlock把标记块拆成八个子块未标记的块整体复制。相比全分辨率遍历这个结构在物体只占包围盒一小部分时优势巨大因为粗层就能砍掉 90% 以上的空间。3. 投影链路一个体素是怎么落到像素上的3.1 投影矩阵 P把四个变换压成一次矩阵乘从世界坐标到像素坐标要经过三步世界坐标乘外参变成相机坐标相机坐标除以深度做透视除法变成归一化坐标再乘内参变成像素坐标。这三步可以合成一个 3×4 的矩阵 PXc R * Xw t % 外参旋转 平移 u fx * Xc(1)/Xc(3) cx v fy * Xc(2)/Xc(3) cy写成齐次形式就是[u; v; w] P * [Xw; 1]其中P K * [R t]K [fx 0 cx; 0 fy cy; 0 0 1]。用 P 的好处不只是少写几行矩阵乘之后得到的第三个分量 w 恰好等于相机坐标下的深度 Zc用它一次就能筛掉“在相机背后”的点。这一步漏掉投影会翻车得莫名其妙——背后的点被透视除法除成负数u、v 跑到图像另一侧恰好落在轮廓里于是本该删掉的体素被保留模型朝着相机背后长出一大片鬼影。参数含义需要记牢fx、fy 是以像素为单位的焦距cx、cy 是主点一般接近图像中心但不完全等于R 是 3×3 旋转t 是 3×1 平移。MATLAB 的图像坐标系 y 轴向下标定工具给出的 K 一般已经按这个约定组织好自己手搓的时候最容易在这里反号。3.2 一次算完所有体素向量化投影的写法MATLAB 的循环慢是出了名的把 N³ 个体素逐个 for 循环投影128³ 就是两百万次循环跑一晚上都算不完。正确写法是把所有体素坐标堆成一个 4×M 的矩阵一次矩阵乘搞定所有投影。% V: 4×M 齐次坐标M N^3每一列是一个体素 [U, Vv, W] ndgrid(linspace(bb(1), bb(2), N), ... linspace(bb(3), bb(4), N), ... linspace(bb(5), bb(6), N)); Pts [U(:); Vv(:); W(:); ones(1, numel(U))]; % 4×Msingle 可省一半内存 for k 1:numel(cams) pr cams(k).P * Pts; % 3×M z pr(3, :); ok z 1e-6; % 只保留相机前方的体素 u pr(1, ok) ./ z(ok); v pr(2, ok) ./ z(ok); in u 1 u W_img v 1 v H_img; % 落在图像范围内 idx sub2ind([H_img, W_img], round(v(in)), round(u(in))); hit false(1, size(Pts, 2)); hit(ok) false; % 先全判为未命中 tmp false(1, nnz(ok)); tmp(in) cams(k).sil(idx); % 查该相机的轮廓表 hit(ok) tmp; alive alive hit; % 一张图说背景永久删除 end关键点逐个说。pr P * Pts是一次 3×4 乘 4×MBLAS 会把它优化到接近内存带宽的极限。ok z 1e-6是前面强调的深度筛选阈值取 1e-6 而不是 0是因为浮点误差可能让本该为正的深度变成 -1e-17除出来是一个巨大的负数round之后索引越界直接报错。in筛掉投影跑到图像外的体素——这些体素无法判断按保守原则应当保留所以用一个hit全 false 的向量接住只把真正落在图像内的位置填上查表结果。最后alive hit完成交集注意是“与”不是“或”任意一张图判定为背景就删除。cams(k).sil是预先二值化好的轮廓图逻辑矩阵用sub2ind把 (v,u) 转成线性索引。这里round到最近像素等价于用体素中心判定如果想保守一点把round换成floor并在查表时对 2×2 邻域取“与”就能实现前面说的四角判定。3.3 轮廓查表距离变换比硬阈值抗噪直接用二值轮廓查表有个隐患轮廓边界上的像素归属本身就不确定。前景分割的边缘往往有一两个像素的抖动体素刚好投在边界上时判定结果完全由这一两个像素决定最终模型表面会出现一层随机的锯齿和孤立小碎块。常见的补救办法是不查二值图改查距离变换。bwdist(~sil)给出每个前景像素到最近背景像素的距离单位是像素。判定条件从“落在轮廓内”改成“到背景的距离大于 r”r 取 1 到 2 个像素。这在几何上等价于把轮廓整体向内收缩 r 个像素也就是把一个可能偏大的轮廓修正回真实边界附近。它的好处是把硬判决变成了带余量的判决分割边缘的随机抖动被平滑掉模型表面明显干净。代价是结果会瘦一圈。r 取大了会把细结构比如手指、细杆整段削掉因为这些部位在图像上的宽度可能只有三五个像素。经验做法是 r 取 1然后靠提高体素分辨率找补精度如果发现细结构断裂先去看是不是轮廓本身就没分好而不是急着调 r。3.4 遮挡问题为什么轮廓法不需要 z-buffer光度法必须知道每个体素在每张图里是不是可见的被挡住的像素颜色不能拿来比较否则一个被前景遮住的背景点会被误判成“颜色不一致”而删除。求可见性就要渲染深度图也就是 z-buffer成本不低。轮廓法完全绕开了这个问题。它不比较颜色只比较“投影像素是不是背景”。一个体素即使被别的东西挡住它在图像上的投影位置仍然落在物体轮廓内部判定依然是“保留”。遮挡只会让判定变得更宽松该删的没删掉绝不会误删这正好符合算法的保守性原则。这是 SpaceCarving 相比体素着色voxel coloring最大的工程优势不需要可见性计算不需要按深度排序所有相机的判定可以完全并行、任意顺序。代价就是凹面恢复不了。想要凹面就得回到可见性 光度一致性那条路上而那一步的复杂度是另一个量级。4. MATLAB 工程实现体素数组、向量化与等值面收尾4.1 先算内存账N³ 到底吃多少动手之前把这张表算一遍能省掉大量“跑到一半内存不足”的调试时间。每轴分辨率 N体素总数 N³齐次坐标矩阵4×M double逻辑体素数组single 坐标矩阵6426.2 万8.4 MB0.26 MB4.2 MB128209.7 万67 MB2.1 MB33.5 MB2561677.7 万537 MB16.8 MB268 MB5121.34 亿4.3 GB134 MB2.1 GB逻辑体素数组本身很小真正吃内存的是那个 4×M 的坐标矩阵以及矩阵乘法过程中产生的临时变量P*Pts还会再生成一个 3×M。所以瓶颈不在分辨率而在“一次性把所有体素都摊平”。256³ 以上就必须分块把 M 切成几万到几十万一组循环处理每块算完立刻把结果写回逻辑数组。分块的大小按可用内存调一般让单块的坐标矩阵控制在 100MB 以内比较安全。另一个省内存的招是用single。投影精度损失在小数点后第七位左右对像素级判定完全够用内存直接减半。注意cams(k).P也要转成 single否则 mixed 运算会自动升回 double白忙一场。4.2 分块循环的写法与索引对齐分块最容易出错的地方是索引对不上坐标矩阵按块切了逻辑数组还是整块的写完一块忘了把块内索引映射回全局索引结果所有块都写到同一个位置。chunk 2e5; % 每块 20 万个体素 M N^3; alive true(1, M, like, single(1)); % 全存活标记线性索引 for s 1:chunk:M e min(s chunk - 1, M); idx s:e; % 本块的全局线性索引 Pts [Xg(idx); Yg(idx); Zg(idx); ones(1, numel(idx), single)]; keep true(1, numel(idx), like, single(1)); for k 1:numel(cams) keep keep projectAndTest(cams(k), Pts, W_img, H_img); end alive(idx) keep; % 用全局索引写回绝不能写成 alive(1:e-s1) end alive reshape(alive, N, N, N);Xg、Yg、Zg是用ndgrid提前算好的三个 N³ 向量每个 N³ 个 single256³ 时各 67MB可以接受比每次重新linspacendgrid快得多也避免了三维数组展开顺序的混乱。alive(idx) keep这一行是整个分块逻辑的正确性核心idx必须是全局索引块内偏移只能用在Pts和keep上。踩过这个坑的人都知道写错之后结果不会报错只是模型少了一大半而且缺的位置很有规律——沿着某条轴整齐地缺看到这种规律就该怀疑索引。projectAndTest就是 3.2 节那段投影 查表的封装返回本块内每个体素是否通过该相机的测试。4.3 从体素到场isosurface 与顶点回映射逻辑体素数组是 0/1 的三维场直接扔给isosurface会在 0.5 处提取等值面得到三角网格。但有两个细节会让结果很难看。第一逻辑场的等值面是方块状的直接提取出来是台阶状的锯齿表面。常规做法是先smooth3做一次三维高斯平滑把 0/1 变成连续场再提等值面表面立刻圆润。平滑核不能太大否则细结构会被抹掉smooth3(vol, gaussian, 3)是个常用起点。第二isosurface返回的顶点坐标是数组下标1 到 N 的浮点数必须手动映射回世界坐标否则模型尺寸和位置全错。vol smooth3(single(alive), gaussian, 3); % 0/1 场 → 连续场 fv isosurface(vol, 0.5); % 提取 0.5 等值面 % 下标 → 世界坐标注意 linspace 的端点对应下标 1 和 N fv.vertices(:,1) bb(1) (fv.vertices(:,1) - 1) * (bb(2) - bb(1)) / (N - 1); fv.vertices(:,2) bb(3) (fv.vertices(:,2) - 1) * (bb(4) - bb(3)) / (N - 1); fv.vertices(:,3) bb(5) (fv.vertices(:,3) - 1) * (bb(6) - bb(5)) / (N - 1); patch(Faces, fv.faces, Vertices, fv.vertices, ... FaceColor, [0.82 0.84 0.88], EdgeColor, none, ... FaceLighting, gouraud);映射公式里的(N-1)和-1是配套的ndgrid生成的第 1 个采样点对应包围盒下界第 N 个对应上界中间等距。写成bb(1) (v-1)*(range)/(N-1)就对了写成bb(1) v*range/N会整体偏移半个到一个体素模型看起来“差一点点”非常难查。4.4 可选精修用光度一致性把凹面挖出来纯轮廓法的结果就是视觉外壳凹面填实。要往下走一步就得在存活的体素上再跑一遍光度一致性检验对每个体素把它投到所有能看到它的相机上取这些像素的颜色如果颜色差异超过阈值就删掉。凹面内部的空间从不同视角看到的颜色分别来自凹陷的两侧壁差异很大于是被正确地挖出来。这一步的难点全在可见性。需要为每个相机维护一张深度图可以用体素投影的 z 值做 splatting 得到近似判定体素在某个相机里是否可见不可见的相机要跳过。实现成本比轮廓法高好几倍收益是凹面细节。一个折中的渐进策略是先用轮廓法在粗分辨率上拿到视觉外壳再只在存活体素附近做一次细分 光度检验把计算量集中砸在有希望的区域。5. 排查清单投影对不上、结果膨胀与内存爆掉的五种典型翻车5.1 模型整个消失或者缩成一个空壳现象跑完发现alive几乎全 false或者只剩薄薄一层视觉上什么都没有。原因几乎总是投影落到了图像外导致所有体素都被判成“背景”。细分下去有三种包围盒的世界坐标单位和相机标定的单位不一致一个用毫米一个用米差一千倍外参的坐标系定义反了有些标定输出的是相机到世界的变换需要求逆再用像素索引基准差一MATLAB 从 1 开始某些工具输出从 0 开始round之后整体偏移一个像素在图像边缘就会越界。解决不要一上来就跑全量。挑三五个已知世界坐标的点比如标定板上某个角点用 P 投影回去和它在图像上的实际像素位置比误差应该在 1 个像素以内。这一步对不上后面全是白搭。这个检查我一般写成一个小函数每次改标定就顺手跑一遍。5.2 结果膨胀成一个圆滚滚的“肥皂”现象模型轮廓大致对但所有凹陷都被填平边缘异常肥厚看起来像裹了一层保鲜膜。原因是轮廓太宽松。常见于背景不均匀、阈值用全局固定值、或者分割后做了imfill填洞。轮廓一旦比真实物体大一圈反投影的锥体就大一圈交集自然胖。另外如果查表用的是体素中心判定而没有内缩余量边界体素会系统性地多保留半个。解决先肉眼比对轮廓图和原图看边缘是不是明显在外面用 Otsu 或自适应阈值替代固定阈值对轮廓做 1 到 2 个像素的腐蚀等价于 3.3 节的距离变换余量。注意腐蚀不能过头否则细结构会断。5.3 表面有整齐的切口或莫名空洞现象模型上出现一块平滑的斜面切口或者某个区域被挖穿。原因是某一张图的轮廓分割失败——物体在那张图里被切掉了一部分对应的锥体把真实表面削掉了。这种情况在结果里非常典型切口是平面状的因为锥体的边界是平面而且沿着某张图的视线方向。另一种可能是相机外参标定错了那张图的锥体整个偏了。解决逐相机做 leave-one-out——去掉某张图重跑如果模型只在去掉某张图时变好问题就在这张图上。便宜一点的检查是把每张图的轮廓叠加回原图看一眼人工扫一遍几十张图比调参快得多。5.4 表面全是台阶和毛刺网格还特别大现象isosurface出来的面片数量爆炸渲染卡顿表面有明显的方块状阶梯。原因是直接对 0/1 逻辑场提等值面或者分辨率太低。逻辑场的等值面天然是方块拼的isosurface会为每个方块边界生成面片数量是连续场的两三倍。解决先smooth3再提面或者干脆把alive转成single后做一次轻微的高斯模糊。如果模型本身分辨率就不够比如物体只占 30 个体素宽平滑只能让锯齿变圆细节回不来该提分辨率还是得提。另外可以用reducepatch做一次面片简化把对精度没贡献的三角形合并掉。5.5 跑到一半内存不足或者跑了一整晚没动静现象Out of memory或者进度条几个小时不动。原因是全分辨率一次性摊平或者用了逐体素的 for 循环。前者是 4×N³ 的矩阵太大后者是 MATLAB 解释器的循环开销。解决分块4.2 节并且用single把最内层的循环换成矩阵运算开启parfor对相机维度并行相机之间完全独立这是最容易并行的方向。如果这些都做了还是慢说明分辨率定得太高先降到 128 看到结果确认流程正确再往上加。6. 取舍轮廓法的天花板在哪里什么时候该换一条路视觉外壳的上界是数学事实不是实现缺陷。这意味着在纯轮廓法这条路上投入再多工程优化也只能逼近这个上界永远突破不了。判断要不要换方法看的不是模型好不好看而是你的物体有多少“轮廓看不见的凹”。一个粗略的判断标准如果从任何方向看过去物体的轮廓都能反映出它的凹陷那视觉外壳就很接近真值如果存在“从所有视角看都被外沿包住”的内腔那部分必然丢失。前者比如浅浮雕、石块、多面体零件后者比如带把手的杯子、镂空的框架、有内壁的容器。方法需要的信息凹面能力计算成本典型失效场景纯轮廓 carving仅二值轮廓 标定无只有轮廓可见的凹低N³×C内腔、把手、镂空轮廓 光度精修轮廓 颜色 可见性有取决于纹理中高需深度图无纹理表面、镜面反射稠密匹配 MVS特征 颜色 极线搜索强高弱纹理、重复纹理实际项目里我倾向于分两阶段第一阶段用轮廓法在 128³ 上快速拿到视觉外壳这一步几乎是确定性的跑出来什么样心里有数第二阶段只在存活体素附近做一次局部细分和光度检验把计算预算花在真正需要的地方。这样即使光度那一步效果不好退回纯轮廓也还是个能用的模型不会全盘翻车。验证也有讲究。手里没有真值模型的时候最实用的检验是重投影一致性把重建出来的网格重新投影到每张图上和原图的轮廓做交并比。IoU 高不代表模型对视觉外壳天然重投影一致但 IoU 低一定说明某处出了问题。另一个办法是留一验证——留几张图不参与重建用它们做检验如果模型在留出的视角下轮廓明显对不上说明前面的锥体求交有系统偏差。最后一个血泪经验这套流程里最容易出错的从来不是 carving 本身而是它前面的标定和分割。carving 的代码写死了也就一百行跑一次几秒钟而轮廓分错一个像素你可能要花两天才能定位到是哪张图。从那以后我每次跑 carving 之前都强制先把所有轮廓叠回原图逐张过一遍再挑五个已知点做投影验证两步都过了才允许开跑。希望帮到你。本文还有配套的精品资源点击获取