1.算法仿真效果
matlab2022a仿真结果如下：
2.算法涉及理论知识概要
双目立体视觉(Binocular StereoVision)是机器视觉的一种重要形式，它是基于视差原理并利用成像设备从不同的位置获取被测物体的两幅图像，通过计算图像对应点间的位置偏差，来获取物体三维几何信息的方法。融合两只眼睛获得的图像并观察它们之间的差别，使我们可以获得明显的深度感，建立特征间的对应关系，将同一空间物理点在不同图像中的映像点对应起来，这个差别，我们称作视差(Disparity)。双目图像就是通过左右两个摄像头拍摄的图片，存在一定视差。深度就是指实际物体到摄像头的距离。
立体视觉意味着人工智能可以通过一对相机来感知图像的深度以及物体的距离。大多数三维相机模型都是基于立体视觉理论和技术的。两台摄像机之间设置一定的距离，这样它们就可以从不同的角度“看”物体。评估两个图像之间的对应关系，人工智能确定到目标的距离，分析，并建立目标的3D结构。
为了度量体积，需要进行下面的步骤：
1、对一组连续的帧进行累积，这将增加错误恢复的弹性，用于对3D场景进行平均或细化。
2、只在场景中选择产品定义的点。这是通过使用颜色分割、模板匹配或神经网络语义分割来实现的。最快的方法是颜色分割。该方法的缺点是将设置绑定到特定的产品上，如果背景颜色和对象不是很清楚，可能会产生不好的结果。如果GPU优化是可能的，那么使用U形卷积神经网络如U-net和高级U-net，或全卷积神经网络会有较高的性能和分割精度。
3、对一个由产品定义的3D点组成的场景进行聚类。每个簇是一个对象。
4、为每个簇形成一个凸多边形，消除3d场景对象的边缘缺陷
5、使用线性插值恢复缺失的3D点
6、通过基于场景的几何聚类区域积分来计算单个对象的体积。
7、最后，计算所有物体的总体积。
3.MATLAB核心程序
%立体参数的可视化
figure;
showExtrinsics(stereoParams);
for ij = 1:10
..............................................................
figure;
subplot(3,2,1)
imshow(stereoAnaglyph(frameLeftRect, frameRightRect));
title('Rectified Frames');
frameLeftGray = rgb2gray(frameLeftRect);
frameRightGray = rgb2gray(frameRightRect);
%% 三维重建
points3D = reconstructScene(disparityMap, stereoParams);
.......................................................................
subplot(3,2,6)
mesh(X, Y, Z);
%% 体积计算
Heights = sum(sum(abs(Z)));
S = abs((y0(2) - y0(1)) * (x0(2) - x0(1)));
V(ij) = Heights * S ; % 体积
clear frameLeft frameRight
clear frameLeftRect frameRightRect
clear frameLeftGray frameRightGray disparityMap points3D ptCloud
clear ptCloudA Temp i j ptCloudB ptCloudC
clear model inlierIndices outlierIndices ptCloudPlane ptCloudD ptCloudE
clear x y z x0 y0 X Y Z Heights S
ij
end
figure;
plot(V,'-bs',...
'LineWidth',1,...
'MarkerSize',6,...
'MarkerEdgeColor','k',...
'MarkerFaceColor',[0.0,0.9,0.0]);
hold on
plot(Vreal,'-bo',...
'LineWidth',1,...
'MarkerSize',6,...
'MarkerEdgeColor','k',...
'MarkerFaceColor',[0.9,0.0,0.0]);
xlabel('物体编号');
ylabel('体积mm3');
legend('算法计算值','真实值');
......................................................
modelNormal = model.Normal;
X_normal = [1 0 0];
Y_normal = [0 1 0];
% Z_normal = [0 0 1];
alpha = acos(dot(modelNormal, X_normal) / (norm(modelNormal) * norm(X_normal)));
beta = acos(dot(modelNormal, Y_normal) / (norm(modelNormal) * norm(Y_normal)));
theta = 0;
alpha = -(alpha - pi / 2);
beta = (beta - pi / 2);
% 旋转矩阵
R = [cos(beta) * cos(theta) cos(beta) * sin(theta) -sin(beta);
-cos(alpha) * sin(theta) + sin(alpha) * sin(beta) * cos(theta) cos(alpha) * cos(theta) + sin(alpha) * sin(beta) * sin(theta) sin(alpha) * cos(beta);
sin(alpha) * sin(theta) + cos(alpha) * sin(beta) * cos(theta) -sin(alpha) * cos(theta) + cos(alpha) * sin(beta) * sin(theta) cos(alpha) * cos(beta)];
Temp(:, 1) = pcIn.Location(:, 1);
Temp(:, 2) = pcIn.Location(:, 2);
Temp(:, 3) = pcIn.Location(:, 3);
% 旋转
Temp = Temp * R;
% Temp(:, 3) = Temp(:, 3) - min(Temp(:, 3));
% 平移
a = model.Parameters(1);
b = model.Parameters(2);
c = model.Parameters(3);
d = model.Parameters(4);
X = -150 : 0.1 : 150;
Y = -150 : 0.1 : 150;
panelZ = -(a * X + b * Y + d) / c;
Temp(:, 3) = Temp(:, 3) - median(panelZ);
pcOut = pointCloud(Temp);