GEOSによる図形演算

[Siv3D Advent Calendar 2017](https://qiita.com/advent-calendar/2017/siv3d) 23 日目の記事です。

# 導入
[GEOS - Geometry Engine, Open Source](https://live.osgeo.org/ja/overview/geos_overview.html) とは、ジオメトリ計算を行うためのライブラリです。
主に GIS（地理情報システム）に用いられているそうです。

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

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


# 0. 図形の作成

```cpp: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](https://qiita-image-store.s3.amazonaws.com/0/100367/9f52743c-323e-2aad-dfe0-4c39059046a7.png)


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

```cpp: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](https://qiita-image-store.s3.amazonaws.com/0/100367/7c1fc797-673e-0e1d-17ab-ad6a41098e30.gif)


# 2. ブーリアン演算
論理和(Union)、論理積(intersection)、差分(difference)、排他的論理和(symDifference)の四つの演算機能があります。

```cpp: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](https://qiita-image-store.s3.amazonaws.com/0/100367/6a1d0f48-93dc-aa1d-ac48-0180148f9f50.gif)


# 3. バッファ演算

```cpp: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](https://qiita-image-store.s3.amazonaws.com/0/100367/5e165068-3100-815e-8b85-9821080a7438.png)


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

```cpp: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](https://qiita-image-store.s3.amazonaws.com/0/100367/84640607-2b20-6ae5-afe3-db1744cf6ce3.gif)


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

```cpp: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](https://qiita-image-store.s3.amazonaws.com/0/100367/5526afa8-c677-2ef1-8c1b-67d31f93cf96.gif)



# 6. 凸包、ボロノイ分割

```cpp: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](https://qiita-image-store.s3.amazonaws.com/0/100367/81b91e4d-95a1-adaa-c343-05a0e95a1710.gif)


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