C++
Siv3D
Siv3DDay 23

GEOSによる図形演算

Siv3D Advent Calendar 2017 23 日目の記事です。

導入

GEOS - Geometry Engine, Open Source とは、ジオメトリ計算を行うためのライブラリです。
主に GIS(地理情報システム)に用いられているそうです。

ダウンロードは https://trac.osgeo.org/geos から行えます。

Siv3D 用にビルドするのが少し面倒なので、VS2015 でビルド済みのライブラリファイルを GitHub に上げました(GEOS のヘッダファイルは入っていません)。
https://github.com/agehama/geos_build_VS2015

0. 図形の作成

geos_example0.cpp
#include "GeosSiv3DHelper.hpp"

void Main()
{
    while (System::Update())
    {
        gg::Polygon* poly = ToGeosPoly(Rect(100, 50).setCenter(Window::Center()).asPolygon());

        ToSivPoly(*poly).draw();

        delete poly;
    }
}

ex0.png

1. 図形の接触判定

図形の接触判定はほぼ Siv3D と同じですが(形状の種類では円に対応している Siv3D の方が強い)、GEOS には DE-9IM と呼ばれる9種類の接触関係を一度に計算できるという機能があります(Intersects、Contains などそれぞれの関係のみを指定して調べることもできます)。

geos_example1.cpp
#include "GeosSiv3DHelper.hpp"

void Main()
{
    Font font(30);

    gg::Polygon* star = ToGeosPoly(Geometry2D::CreateNStar(5, 100, 50).movedBy(Window::Size()*0.5));

    while (System::Update())
    {
        gg::Polygon* plus = ToGeosPoly(Geometry2D::CreatePlus(30.0, 10.0).movedBy(Mouse::Pos()));

        gg::IntersectionMatrix* relation = star->relate(*plus);

        font(CharacterSet::Widen(relation->toString())).draw();

        ToSivPoly(*star).draw(Color(255, 0, 0));
        ToSivPoly(*plus).draw(Color(0, 255, 0));

        delete plus;
        delete relation;
    }

    delete star;
}

example1.gif

2. ブーリアン演算

論理和(Union)、論理積(intersection)、差分(difference)、排他的論理和(symDifference)の四つの演算機能があります。

geos_example2.cpp
#include "GeosSiv3DHelper.hpp"

void Main()
{
    gg::Polygon* star = ToGeosPoly(Geometry2D::CreateNStar(5, 100, 50).movedBy(Window::Size()*0.5));

    double angle = 0.0;
    while (System::Update())
    {
        angle += 0.01;
        gg::Polygon* plus = ToGeosPoly(Geometry2D::CreatePlus(30.0, 10.0).rotated(angle).movedBy(Window::Size()*0.5));

        gg::Geometry* resultGPoly = plus->symDifference(star);
        const MultiPolygon resultSPoly = ToSivPoly(*resultGPoly);

        resultSPoly.draw();
        resultSPoly.drawWireframe(1.0, Palette::Black);

        ToSivPoly(*plus).drawFrame();
        ToSivPoly(*star).drawFrame();

        delete plus;
        delete resultGPoly;
    }

    delete star;
}

example2.gif

3. バッファ演算

geos_example3.cpp
#include "GeosSiv3DHelper.hpp"

void Main()
{
    const double pos = 150.0;

    gg::Polygon* star = ToGeosPoly(Geometry2D::CreateNStar(5, 100, 50).movedBy(pos, Window::Height()*0.5));
    gg::Polygon* plus = ToGeosPoly(Geometry2D::CreatePlus(100.0, 50.0).movedBy(Window::Width() - pos, Window::Height()*0.5));

    while (System::Update())
    {
        gg::Geometry* resultGStar = star->buffer(15.0, 10);
        gg::Geometry* resultGPlus = plus->buffer(-15.0, 10);

        const MultiPolygon resultSStar = ToSivPoly(*resultGStar);
        const MultiPolygon resultSPlus = ToSivPoly(*resultGPlus);

        resultSStar.draw();
        resultSStar.drawWireframe(1.0, Palette::Gray);

        resultSPlus.draw();
        resultSPlus.drawWireframe(1.0, Palette::Gray);

        ToSivPoly(*star).drawFrame(1.0, Palette::Cyan);
        ToSivPoly(*plus).drawFrame(1.0, Palette::Cyan);

        delete resultGStar;
        delete resultGPlus;
    }

    delete star;
    delete plus;
}

ex3.png

4. 折れ線からポリゴンの生成

geos_example4.cpp
#include "GeosSiv3DHelper.hpp"

void Main()
{
    gg::CoordinateArraySequence points;
    std::vector<LineString> lines;

    while (System::Update())
    {
        if (Input::MouseL.clicked)
        {
            points.add(gg::Coordinate(Mouse::Pos().x, Mouse::Pos().y));
        }

        if (Input::MouseL.pressed && !points.empty() && 30.0 < points.back().distance(gg::Coordinate(Mouse::Pos().x, Mouse::Pos().y)))
        {
            points.add(gg::Coordinate(Mouse::Pos().x, Mouse::Pos().y));
        }

        if (Input::MouseL.released)
        {
            gg::GeometryFactory::unique_ptr factory = gg::GeometryFactory::create();

            gg::PrecisionModel model(gg::PrecisionModel::Type::FLOATING);

            gob::BufferParameters param(10, gob::BufferParameters::CAP_ROUND, gob::BufferParameters::JOIN_BEVEL, 10.0);
            //gob::BufferParameters param(10, gob::BufferParameters::CAP_FLAT, gob::BufferParameters::JOIN_BEVEL, 10.0);

            gob::OffsetCurveBuilder builder(&model, param);

            std::vector<gg::CoordinateSequence*> result;

            builder.getLineCurve(&points, 15.0, result);
            //builder.getSingleSidedLineCurve(&points, 10.0, result, true, false);

            for (gg::CoordinateSequence* p : result)
            {
                std::vector<Vec2> ps;
                for (size_t i = 0; i < p->size(); ++i)
                {
                    ps.emplace_back(p->getX(i), p->getY(i));
                }

                lines.emplace_back(ps);
            }

            points.clear();
        }

        if (!points.empty())
        {
            for (size_t i = 0; i+1 < points.size(); ++i)
            {
                const double x0 = points.getX(i);
                const double y0 = points.getY(i);
                const double x1 = points.getX(i + 1);
                const double y1 = points.getY(i + 1);
                Line(x0, y0, x1, y1).draw();
            }
        }

        for (const auto& l : lines)
        {
            l.draw();
        }
    }
}

example4.gif

5. 図形同士の距離、最近接点

geos_example5.cpp
#include "GeosSiv3DHelper.hpp"

void Main()
{
    gg::CoordinateArraySequence points;
    std::vector<LineString> lines;

    gg::Polygon* star = ToGeosPoly(Geometry2D::CreateNStar(5, 100, 50).movedBy(Window::Size()*0.5));

    while (System::Update())
    {
        std::vector<Vec2> points;

        gg::Polygon* plus = ToGeosPoly(Geometry2D::CreatePlus(40.0, 20.0).movedBy(Mouse::Pos()));
        god::DistanceOp distanceOp(*star, *plus);

        gg::CoordinateSequence* ps = distanceOp.closestPoints();

        for (size_t i = 0; i < ps->size(); ++i)
        {
            const double x0 = ps->getX(i);
            const double y0 = ps->getY(i);
            points.emplace_back(x0, y0);
        }

        LineString(points).draw(1.0, Palette::Red);

        ToSivPoly(*star).drawFrame();
        ToSivPoly(*plus).drawFrame();

        delete plus;
    }

    delete star;
}

example5.gif

6. 凸包、ボロノイ分割

geos_example6.cpp
#include "GeosSiv3DHelper.hpp"

void Main()
{
    double angle = 0.0;

    while (System::Update())
    {
        angle += 0.02;
        gg::Polygon* star = ToGeosPoly(Geometry2D::CreatePlus(100.0, 20.0).rotated(angle).movedBy(Window::Size()*0.5));

        gt::VoronoiDiagramBuilder builder;
        builder.setSites(*star->getCoordinates());

        gg::GeometryFactory::unique_ptr factory = gg::GeometryFactory::create();
        std::auto_ptr<gg::Geometry> edges = builder.getDiagramEdges(*factory);

        gg::CoordinateSequence* cs = edges->getCoordinates();

        ToSivPoly(*star).draw();

        for (size_t i = 0; i + 1 < cs->size(); ++i)
        {
            const double x0 = cs->getX(i);
            const double y0 = cs->getY(i);
            const double x1 = cs->getX(i + 1);
            const double y1 = cs->getY(i + 1);
            Line(x0, y0, x1, y1).draw(1.0, Palette::Cyan);
        }

        ToSivPoly(*edges->getEnvelope()).drawFrame(1.0, Palette::Cyan);

        delete star;
    }
}

example6.gif

終わりに

一昨日触り始めたばかりですが便利かもしれません。