切片法计算点云体积时,首先对点云按照每层切片,但是对于地物而言,有些地方存在缺失导致切片时某些切片层可能计算不准,为此可以从下到上依次将点云到每一层的点统计出来并基于此作为切片点云数据。然后对每层切片计算面积,再通过多边体计算公式实现点云体积计算。
式中为顶处棱锥体积,
为沿Z轴方向最高分割点处的切片投影面面积,
,
为当前投影面面积与上层平面投影面面积。因为最高处只有一个点,所以最高处据最高处最近的切片层计算公式得变换为棱锥计算公式。
double compute_volume(const pcl::PointCloud<pcl::PointXYZ>::ConstPtr& cloud)
{
pcl::PointXYZ minPt, maxPt;
pcl::getMinMax3D(*cloud, minPt, maxPt);
double delta = maxPt.z - minPt.z;
double h{ 0.1 };
int levels = delta / h;
double residual_delta = delta / h - levels;
vector<double> areacount(levels);
for (int i{ 0 }; i < levels; i++)
{
pcl::PointCloud<pcl::PointXYZ>::Ptr compute_cloud(new pcl::PointCloud<pcl::PointXYZ>);
for (int j{ 0 }; j < cloud->size(); j++)
{
if (cloud->points[j].z > i * h + minPt.z)
{
compute_cloud->push_back(cloud->at(j));
}
}
double comarea = compute_area(compute_cloud);
areacount[i] = comarea;
}
double volume{ 0 };
for (int i = 0; i < areacount.size(); i++)
{
if (i + 1 == areacount.size())
break;
double volume_ = ((areacount[i] + areacount[i + 1]) * h) / 2;
volume += volume_;
}
volume += (1 / 3) * h * (areacount[areacount.size() - 1]);
return volume;
}