1.算法仿真效果
matlab2022a仿真结果如下：
2.算法涉及理论知识概要
visual hull用于简单、快速地获取三维模型。目前,获取三维模型的方式有:
利用传统几何构造技术直接构造模型.利用三维扫描设备对真实物体进行扫描,进而重建出模型.利用从各个视角拍摄的真实物体的多幅图像重建模型.由图像重建三维模型技术又可分为两类:一类是通过多幅深度图像重建模型,另一类是通过多幅照片生成物体的可见外壳visual hull.
Laurentini 最早提出了可见外壳(visual hulls )的概念,即可见外壳就是由空间物体的所有已知侧影轮廓线决定的该物体的凸包，是一种由形状到轮廓的3D重建技术来构建几何对象。具体方法是，假定场景中的前景物体可以和背景相分离。基于这样的假设，原始图像可以经过阈值化前/背景二值图像，得到轮廓图像（silhouette image）。前景掩膜，所谓的轮廓，是相应3D前景物体的2D投影。伴随着相机外参，轮廓定义了反向投影的通用椎体，该锥体包含这个前景物体。
当利用透视投影的方式从多个视角观察某一空间物体时,在每一个视角上都会得到一条该物体的侧影轮廓线。这条侧影轮廓线和对应的透视投影中心将共同决定一个三维空间中一般形状的锥体，轮廓椎体（silhouette cone）。而由所有已知侧影轮廓线和对应的透视投影中心决定的锥体的intersection相交点，最终将决定一个包含该物体的凸包，即物体的可见外壳。在多数情况下,空间物体的可见外壳是该物体的合理逼近.
Visual hull（可视外壳）是由Laurentini引入的有形状到轮廓的3D重建技术构建的几何实体。这种技术假定场景中的前景物体可以和背景相分离。基于这样的假设，原始图像可以经过阈值化处理成前/背景二值图像，我们将这种图像称作轮廓图像（silhouette image）。所谓的轮廓，是相应3D前景物体的2D投影。伴随着相机观看参数，轮廓定义了后向投影的涵盖实际物体的通用锥体。这个锥体被称作轮廓锥体（silhouette cone）。两个锥体相交的地方被称作可视外壳，它是实际3D物体的几何边界。
所以说，可视外壳就是由空间物体的所有已知侧影轮廓线决定的该物体的凸包。当利用透视投影的方式从多个视角观察某一空间物体时，在每一个视角上都会得到一条该物体的侧影轮廓线。这条侧影轮廓线和对应的透视投影中心将共同决定一个三维空间中一般形状的锥体。显然，该物体必将落于这个锥体之内;而由所有已知侧影轮廓线和对应的透视投影中心决定的锥体的交最终将决定一个包含该物体的凸包，这个凸包就是物体的可视外壳。在多数情况下，空间物体的可视外壳是该物体的合理逼近。
Martin和Aggarwal 最早提出了用多幅照片生成可视外壳的方法。他们利用真实物体在照片上的侧影轮廓线和相关的相机参数实现了物体的三维重建。Chien等人采用八叉树结构来表示物体的可视外壳。这种八叉树结构是通过预先在平行投影照片上生成表示物体的四叉树的基础上建立的。Potmesil和Szeliski同样建立了可视外壳的八叉树表示，但输入改为从任意视角拍摄的透视投影照片。Garcia 等人也提出了一种改进的八叉树表示方法。
3.MATLAB核心程序 `object_points3D = [voxels(:,1)'; voxels(:,2)'; voxels(:,3)'; ones(1, length(voxels))];
voxels(:, 4) = zeros(size(voxels(:, 4)));
voxel_inds = {};
img_vals = {};
dmin = depth_range(1);
dmax = depth_range(2);
img_size = size(silhouettes);
for i = 1:size(M,3)
r = M(1:3, 1:3, i);
t = M(1:3, 4, i);
if t == [0;0;0]
continue;
end
center = [-(r')* t; 1];
cam_center = repmat(center, 1, size(object_points3D,2));
KM = [K]*M(1:3, 1:4, i);
points2D = KM*object_points3D;
points2D = points2D./[points2D(3,:); points2D(3,:); points2D(3,:)];
cur_silhouette = silhouettes(:,:,i);
pixThresh = mean(cur_silhouette(:)) + std(cur_silhouette(:));
cur_mask = cur_silhouette > pixThresh;
saliencyFlag = zeros(1,size(points2D,2),'logical');
for ii = 1:size(points2D,2)
if points2D(2,ii) > 1 && points2D(2,ii) < img_size(1) && ...
points2D(1,ii) > 1 && points2D(1,ii) < img_size(2)
saliencyFlag(ii) = cur_mask(round(points2D(2,ii)),round(points2D(1,ii)));
end
end
saliencyFlag = ones(1,size(points2D,2),'logical');
d = sqrt(sum((cam_center - object_points3D).^2, 1));
pts_ind = find(... %d > dmin & d < dmax ...
points2D(2,:) > 1 & points2D(2,:) < img_size(1) ...
& points2D(1,:) > 1 & points2D(1,:) < img_size(2) ...
& saliencyFlag);
pi = points2D(1,pts_ind);
pj = points2D(2,pts_ind);
[img_val, ind, object_points_cam] = GetSilhouetVals([pi;pj], cur_silhouette);
ind = pts_ind(ind);
voxel_inds{i} = ind;
img_vals{i} = img_val;
if(display_projected_voxels)
figure(fid),
imagesc(cur_silhouette);title(i);hold on
plot(object_points_cam(1,:), object_points_cam(2,:), '.g');hold off
end
end`
