第17节 实例-倾斜摄影数据挖填方

致谢

感谢网友提出这个问题。大家有问题也可以在评论区提出，有问必有答。

问题描述

在第12节（这是第17节，往前翻翻）挖填方的基础上，仍有网友提问怎么做倾斜摄影的挖填方，仿照大疆的一个什么软件的结果。如图：

下图是本程序的对比（结果大差不差）：

各方确认各点的精确值并没有，因此上面的点我只能大差不差的点，而上图还一个最关键信息是：基准面：最低点

也就是以点出的点的最低点为基准面来计算挖填方。同时也有网友提出来我要在这个图形上挖个沟，沟的模型有，地表的模型有，求挖沟的挖方与填方。这与本节中所说的方法都是一样的。以供大家参考。

资源下载（这一节是有模型的）

本文集包括本节所有资源包括模型代码都在此下载，按节的序号有文件或文件夹：

【击此打开网盘资源链接】

本节程序功能

1. 程序加载倾斜摄影，调整到要测试挖填方的位置后，按住左shift键，鼠标在地图上单击，则会出现小球。超过三个点则会出现小球围成的面。
2. 点完之后，点击s键，就会进行挖填方的计算。即使用原始点围面的面来计算，也使用以最低点为基准面围成的面来计算。

关键思想

1. 比之第12节，我们求交发生了变化，首先与我们点的面来求交，其次再与地形求交。比较其z值的大小，看谁在谁之上。若面在地形之上，则是填方，若地形在面之上则是挖方。
2. 若依最低点为基准面，则只需要遍历所有的顶点的z值，找到最小的，在上一步的基础上将与点的面的求交点的z值设置为最低点的z值。再按照上一步的思路求一下就OK。

以下为全部代码，模型在网盘中下载

#include <osgViewer/Viewer>
#include <osgDB/ReadFile>
#include <osgGA/EventHandler>
#include <osg/ShapeDrawable>
#include <osg/Shape>

osg::Node* _base = nullptr;//地形
osg::Group* _root = new osg::Group;//场景根结点


osg::Geode* _mask = new osg::Geode;//鼠标点出来的面
osg::Geometry* _maskGeom = new osg::Geometry;
osg::DrawArrays* _da = new osg::DrawArrays;

osg::Group* _mark = new osg::Group;//小球

osg::Node* createSphere(osg::Vec3 center)
{
    osg::Geode* gnode = new osg::Geode;
    osg::ShapeDrawable* sd = new osg::ShapeDrawable(new osg::Sphere(center, 0.2));
    gnode->addDrawable(sd);
    gnode->getOrCreateStateSet()->setMode(GL_LIGHTING, osg::StateAttribute::OFF);
    return gnode;
}

osg::Geode* CreateBox(osg::Vec3 from, osg::Vec3 to, float stepSize, osg::Vec4 color)
{
    osg::Geode* gnode = new osg::Geode;
    osg::ShapeDrawable* sd = new osg::ShapeDrawable(new osg::Box((from + to) / 2, stepSize, stepSize, std::abs(to.z() - from.z())));
    sd->setColor(color);
    gnode->addDrawable(sd);
    return gnode;
}

class MyEventHandler : public osgGA::GUIEventHandler
{
public:
    MyEventHandler()
    {
        vertexArray = new osg::Vec3Array;
        _maskGeom->setVertexArray(vertexArray);

        osg::Vec4Array* color = new osg::Vec4Array;
        color->push_back(osg::Vec4(1.0, 0.0, 0.0, 0.4));
        _maskGeom->setColorArray(color, osg::Array::BIND_OVERALL);
    }

    virtual bool handle(const osgGA::GUIEventAdapter& ea, osgGA::GUIActionAdapter& aa)
    {
        if (ea.getEventType() == ea.PUSH)
        {
            //是左键且左shift按下了
            if ((ea.getButton() == ea.LEFT_MOUSE_BUTTON) && (ea.getModKeyMask() == ea.MODKEY_LEFT_SHIFT))
            {
                osgUtil::LineSegmentIntersector::Intersections intersections;
                osgViewer::View* view = dynamic_cast<osgViewer::View*>(&aa);
                if (view->computeIntersections(ea, intersections))
                {
                    //有交点，先画个球
                    osg::Vec3 hitPoint = intersections.begin()->getWorldIntersectPoint();
                    _mark->addChild(createSphere(hitPoint));

                    vertexArray->push_back(hitPoint);

                    if (vertexArray->size() >= 3)
                    {
                        _da->set(GL_POLYGON, 0, vertexArray->size());
                        _maskGeom->dirtyDisplayList();//redraw
                    }
                    
                }

                
            }
        }

        if (ea.getEventType() == ea.KEYDOWN)
        {
            if((ea.getKey() == 'x')|| (ea.getKey() == 'X'))
            {
                _da->set(GL_POLYGON, 0, 0);
                _maskGeom->dirtyDisplayList();//redraw

                _mark->removeChildren(0, _mark->getNumChildren());

                vertexArray->clear();
            }

            if ((ea.getKey() == 's') || (ea.getKey() == 'S'))
            {
                //开始计算挖方量与填方量
                //与点围成的面先求是否有交点，有交点再求与地形交点，两个交点之间的距离就是需要挖/填的方
                //其中与面的距离在上，则是填方，与地形交点距离在上，则是挖方，使用z值大小来判断距离
                //求要求的场景的boudingbox
                _maskGeom->dirtyBound();//因为顶点改变，必须重新计算包围盒
                osg::BoundingSphere bs = _maskGeom->getBound();
                //求出xyz的最小值和最大值，然后根据间隔来求
                osg::Vec3 center = bs.center();
                float r = bs.radius();

                //以点的原点为基准面
                //挖方量
                float volumesW = 0.0;
                //填方量
                float volumesT = 0.0;

                //以最低点为基准面
                //挖方量
                float volumesWL = 0.0;
                //填方量
                float volumesTL = 0.0;
                //查找基准面的z值的最小值
                float zmin = FLT_MAX;
                for (int i = 0; i < vertexArray->size(); i++)
                {
                    if (zmin > vertexArray->at(i).z())
                    {
                        zmin = vertexArray->at(i).z();
                    }
                }

                //采样一个方向上采50份
                float stepSize = r / 50;
                for (float fromx = center.x() - r; fromx <= center.x() + r; fromx += stepSize)
                {
                    for (float fromy = center.y() - r; fromy <= center.y() + r; fromy += stepSize)
                    {
                        //乘1000是为了确保与地形相交，与面相交不需要乘1000
                        osg::Vec3 from = osg::Vec3(fromx, fromy, center.z() - r*1000);
                        osg::Vec3 to = osg::Vec3(fromx, fromy, center.z() + r*1000);
                        osgUtil::LineSegmentIntersector::Intersections intersections;
                        osg::ref_ptr<osgUtil::LineSegmentIntersector> ls = new osgUtil::LineSegmentIntersector(from, to);
                        osg::ref_ptr<osgUtil::IntersectionVisitor> iv = new osgUtil::IntersectionVisitor(ls);
                        _mask->accept(*iv.get());

                        if (ls->containsIntersections()) //与点出来的面有交点
                        {
                            intersections = ls->getIntersections();

                            //与地形必有交点
                            osgUtil::LineSegmentIntersector::Intersections intersections0;
                            osg::ref_ptr<osgUtil::LineSegmentIntersector> ls0 = new osgUtil::LineSegmentIntersector(from, to);
                            osg::ref_ptr<osgUtil::IntersectionVisitor> iv0 = new osgUtil::IntersectionVisitor(ls0);
                            _base->accept(*iv0.get());

                            if (ls0->containsIntersections()) //与地形必有交点
                            {
                                intersections0 = ls0->getIntersections();


                                osg::Vec3 maskInter = intersections.begin()->getWorldIntersectPoint();
                                osg::Vec3 terrainInter = intersections0.begin()->getWorldIntersectPoint();

                                if (maskInter.z() > terrainInter.z())
                                {
                                    volumesT += (stepSize * stepSize * std::abs(maskInter.z() - terrainInter.z()));
                                    _mark->addChild(CreateBox(maskInter, terrainInter, stepSize, osg::Vec4(0.0, 1.0, 0.0, 1.0)));
                                }
                                else
                                {
                                    volumesW += (stepSize * stepSize * std::abs(maskInter.z() - terrainInter.z()));
                                    _mark->addChild(CreateBox(maskInter, terrainInter, stepSize, osg::Vec4(1.0, 0.0, 0.0, 1.0)));
                                }

                                maskInter.z() = zmin;
                                if (maskInter.z() > terrainInter.z())
                                {
                                    volumesTL += (stepSize * stepSize * std::abs(maskInter.z() - terrainInter.z()));
                                }
                                else
                                {
                                    volumesWL += (stepSize * stepSize * std::abs(maskInter.z() - terrainInter.z()));
                                }

                            }                            
                        }

                    }
                }

                std::cout << "采用原点为基准面：" << std::endl;
                std::cout << "挖方量：" << volumesW << std::endl;
                std::cout << "填方量：" << volumesT << std::endl;
                std::cout << "采用最低点为基准面：" << std::endl;
                std::cout << "挖方量：" << volumesWL << std::endl;
                std::cout << "填方量：" << volumesTL << std::endl;
                
            }
        }

        return false;
    }

    osg::Vec3Array* vertexArray;
};


int main()
{
    osgViewer::Viewer viewer;

    _base = osgDB::readNodeFile("model/Block/Block.osgb");
    _base->getOrCreateStateSet()->setMode(GL_LIGHTING, osg::StateAttribute::OFF);
    _root->addChild(_base);

    _mask->addDrawable(_maskGeom);
    _maskGeom->addPrimitiveSet(_da);
    _mask->getOrCreateStateSet()->setMode(GL_BLEND, osg::StateAttribute::ON);
    _mask->getOrCreateStateSet()->setRenderingHint(osg::StateSet::TRANSPARENT_BIN);

    _root->addChild(_mark);
    _root->addChild(_mask);

    viewer.setSceneData(_root);
    viewer.addEventHandler(new MyEventHandler());
    return viewer.run();
}