版权声明:转载请标注源地址。 https://blog.csdn.net/rkexy/article/details/82145812
由于我的极差记忆力,我打算把这个破玩意先记下来。因为以后会有改动(Delaunay三角网生成算法),我不想把一个好的东西改坏了。。。
好吧……
凸包生成算法,:
1.先在指定的宽(width
)高(height
)范围内生成一堆随机点;
1.1. 生成N个不重复的正整数,使用洗牌算法让生成的数字不重复;
1.2. 将每个数字分解成坐标。可以设想一个二维数组,每个数字依次填进数组内。那么,对于数字A来说,它能够生成的坐标则为:
x = A % width;
y = (A% width== 0) ? A / width : A / width+ 1);
2. 对所有随机点进行排序;
2.1.找到 y 值最小的点
pointYmin
;
2.2. 计算每个离散点
pt
、
pointYmin
、和X轴组成角的弧度
radian
(或者角度
angle
也可以)。其中,点
pointYmin
为角心;
2.3. 依照上面计算出的弧度进行排序。
3. 开始找凸包
3.1.(我懒得写了……) https://www.cnblogs.com/aiguona/p/7232243.html
3.2. 大体上就是,首先离散点集之中入栈两个点,然后需要计算每三个点
ptA
(栈内点,次顶点),
ptO
(栈内点,栈顶点),
ptB
(栈外点,当前点)之间的逆时针夹角,如果夹角是小于180(π)的,那么则
ptB
入栈,否则
ptO
出栈。直到所有离散点都计算完成。
以下是代码……C++11
File: stdinc.h
#include <cmath>
#include <random>
#include <stdarg.h>
#include <string>
#include <exception>
#include <vector>
#include <stack>
#include <map>
#include <set>
#pragma once
#ifdef _DEBUG
#define ErrorThrow(errMsg) { LogThrow::ErrorOccur(errMsg, __FILE__, __LINE__);}
#include <cassert>
#else
#define ErrorThrow(errMsg) {}
#endif // _DEBUG
#include <iostream>
static const double PI = 3.14159265358980/*...80... -> ...7932384626...*/;
static const double precision = 1000.0;
static const double doubleESP = 1e-008;
namespace LogThrow
{
#include <windows.h>
static void ErrorOccur(const std::string& errMsg, const std::string& fileName, int line)
{
std::string msg("[ERROR] " + errMsg + '\n' + "File: " + fileName + '\n' + "Line: " + std::to_string(line));
std::logic_error(msg.c_str());
MessageBox(NULL, msg.c_str(), NULL, MB_ICONINFORMATION | MB_OK);
ExitProcess(0);
}
};
namespace MyMathTools
{
inline static bool LessThan(double a, double b)
{
if (a - b < doubleESP)
{
return true;
}
return false;
}
inline static bool IsEqual(double a, double b)
{
if (std::abs(a - b) <= doubleESP)
{
return true;
}
return false;
}
};
File: Point2D.cpp
#include "stdinc.h"
#pragma once
//Point & Vector
template<typename T>
class Point2D
{
private:
public:
typedef std::vector<const Point2D<int>> Points;
T x, y;
public:
Point2D(T nx, T ny) :x(nx), y(ny)
{
}
Point2D() :Point2D{ 0, 0 }
{
}
Point2D(const Point2D& np) :Point2D{ np.x, np.y }
{
}
void Init(T x, T y)
{
this->x = x;
this->y = y;
}
void Transposition()
{
T total = x + y;
x = total - x;
y = total - y;
}
void operator= (const Point2D& np/*new point*/)
{
Init(np.x, np.y);
}
bool operator== (const Point2D& pointOther)
{
if ((pointOther.x == this->x) &&
(pointOther.y == this->y))
{
return true;
}
return false;
}
friend std::ostream& operator<<(std::ostream& out, const Point2D& pt)
{
out << "(" << pt.x << "," << pt.y << ")";
return out;
}
public:
//[min, max)
/*const */Point2D& RandomPoint(int maxX, int maxY)
{
if (maxX <= 0 || maxY <= 0)
{
std::logic_error("The less then 0!");
}
std::random_device rd;
std::mt19937 mt(rd());
this->x = mt() % maxX;
this->y = mt() % maxY;
//I can modify the left value.
return *this;
}
template<typename T1, typename T2>
static double Distance(const Point2D<T1>& ptA, const Point2D<T2>& ptB)
{
double cutX = ptB.x - ptA.x;
double cutY = ptB.y - ptA.y;
return std::pow((cutX * cutX + cutY * cutY), .5);
}
template<typename T1>
double Distance(const Point2D<T1>& ptOther)
{
return Point2D::Distance<T, T1>(*this, ptOther);
}
//Without repetition: Shuffle algorithm.
static void RandomPoints(int count, int maxX, int maxY, std::vector<const Point2D<int>>& __out points)
{
//auto LcmFunc = [maxX, maxY](int numA, int numB)
//{
// int max = (maxX > maxY) ? maxX : maxY;
// int min = maxX + maxY - max;
// long lcm = -1;
// for (int i = 1; i <= min; ++i)
// {
// lcm = static_cast<long>(max * i);
// if (lcm % min == 0)
// {
// break;
// }
// }
// return lcm;
//};
//int lcm = static_cast<char32_t>(LcmFunc(maxX, maxY));
//if (lcm == -1)
//{
// return false;
//}
if (count >= maxX * maxY * 0.5)
{
ErrorThrow("Error Count.");
}
points.clear();
auto GetX = [maxX](int num)
{
return num % maxX;
};
auto GetY = [maxX](int num)
{
return (num % maxX == 0) ? num / maxX : num / maxX + 1;
};
static std::random_device rd;
static std::mt19937 mt(rd());
const long total = maxX * maxY;
std::vector<long> nums;
nums.resize(total);
points.resize(count);
for (long i = 0; i < total; ++i)
{
nums[i] = i;
}
for (int i = 0; i < count; ++i)
{
//It is array nums index.
long randomNum = mt() % (total - i) + i;
std::swap(nums[i], nums[randomNum]);
//Swap over, the num[i] is new random number.
//Get the point.
points[i] = { GetX(nums[i]), GetY(nums[i]) };
}
}
};
template<typename T>
class Vector2D final : public Point2D<T>
{
public:
Vector2D() : Point2D()
{
}
Vector2D(T nx, T ny) : Point2D(nx, ny)
{
}
Vector2D(const Vector2D& np) :Point2D(np)
{
}
double Length() const
{
return std::pow((x * x + y * y), .5);
}
//Vector AB
void Vector(const Point2D& ptA, const Point2D& ptB)
{
Init(ptB.x - ptA.x, ptB.y - ptA.y);
}
double Radian(const Vector2D& otherVec, bool isQuadrant = true) const
{
return Vector2D::Radian(*this, otherVec, isQuadrant);
}
//anti-clockwise angle.
//If the return value less than 0, the Radian is less than π.
static double Radian(const Vector2D& vec1, const Vector2D& vec2, bool isQuadrant = true)
{
if (isQuadrant)
{
//CHECK IT:
double radian = std::atan2(vec2.y, vec2.x) - std::atan2(vec1.y, vec1.x); //return [-2π, 2π]
if (radian < 0)
{
radian += 2 * PI;
}
return radian;
}
return std::acos(Vector2D::Dot(vec1, vec2) / (vec1.Length() * vec2.Length())); //return [0, π]
}
double Angle(const Vector2D& otherVec, bool isQuadrant = true) const
{
return Vector2D::Angle(*this, otherVec, isQuadrant);
}
//The same as function Radian.
static double Angle(const Vector2D& vec1, const Vector2D& vec2, bool isQuadrant = true)
{
return Vector2D::Radian(vec1, vec2, isQuadrant) / PI * 180;
}
T Dot(const Vector2D& otherVec) const
{
return Vector2D::Dot(*this, otherVec);
}
static T Dot(const Vector2D& vec1, const Vector2D& vec2)
{
return vec1.x * vec2.x + vec1.y * vec2.y;
}
T Cross(const Vector2D& otherVec) const
{
return Vector2D::Cross(*this, otherVec);
}
//2D Vector have no cross. But we can set the Z is 0.
//So , return a value is (0, 0, Ax*By - Ay*Bx);
static T Cross(const Vector2D& vec1, const Vector2D& vec2)
{
return vec1.x*vec2.y - vec1.y*vec2.x;
}
};
File: ConvexHull.h
//#include "Line.cpp"
#include "Point2D.cpp"
#include "stdinc.h"
#pragma once
class ConvexHull final
{
private:
typedef std::vector<const Point2D<int>> Points;
typedef Points::const_iterator PointItor;
Points m_points;
Points m_CHpoints;
public:
void operator= (const ConvexHull&) = delete;
ConvexHull();
ConvexHull(Points& __out points);
~ConvexHull();
public:
void TestPoints();
void AddPoint(const Point2D<int>& __in newPoint);
void AddRandomPoints(int count, int maxX, int maxY);
void GetConvexHullPoints(Points& __out points);
void GetAllPoints(Points& __out points);
void Generate();
private:
void SortPoints();
};
File: ConvexHull.cpp
#include "ConvexHull.h"
#include <algorithm>
ConvexHull::ConvexHull()
{
}
ConvexHull::ConvexHull(Points& __out points):
m_points(points)
{
}
ConvexHull::~ConvexHull()
{
m_points.clear();
m_CHpoints.clear();
}
void ConvexHull::AddPoint(const Point2D<int>& __in newPoint)
{
m_points.push_back(newPoint);
}
void ConvexHull::GetConvexHullPoints(Points& __out points)
{
points = m_CHpoints;
}
void ConvexHull::GetAllPoints(Points& __out points)
{
points = m_points;
}
//Without repetition: Shuffle algorithm.
void ConvexHull::AddRandomPoints(int count, int maxX, int maxY)
{
std::vector<const Point2D<int>> pointsArray;
Point2D<int>::RandomPoints(count, maxX, maxY, pointsArray);
m_points.clear();
m_points.resize(count);
for (unsigned int i = 0; i < count; ++i)
{
const Point2D<int>& eachpt = pointsArray.at(i);
m_points[i] = eachpt;
}
}
void ConvexHull::SortPoints()
{
if (m_points.size() == 0)
{
ErrorThrow("Empty Points.");
return;
}
//Selecting a point what Y min.
Point2D<int> pointYmin = *m_points.begin();
for each (const Point2D<int>& each in m_points)
{
const Point2D<int>& eachpt = each;
if (eachpt.y < pointYmin.y)
{
pointYmin = eachpt;
}
else if (eachpt.y == pointYmin.y)
{
if (eachpt.x < pointYmin.x)
{
pointYmin = eachpt;
}
}
}
auto SortRuleFunc = [&pointYmin](const Point2D<int>& ptA, const Point2D<int>& ptB)
{
static const Vector2D<int> baseVecX = { 1, 0 };
Vector2D<int> vecA;
vecA.Vector(pointYmin, ptA);
Vector2D<int> vecB;
vecB.Vector(pointYmin, ptB);
double radianA = Vector2D<int>::Radian(baseVecX, vecA);
double radianB = Vector2D<int>::Radian(baseVecX, vecB);
//If collinear occurs, the nearest is before.
if (std::abs(radianA - radianB) <= doubleESP)
{
return vecA.Length() < vecB.Length();
}
//return (radianA - radianB) < doubleESP; //ERROR: ATTENTION: Check the Strict Weak Ordering...
return radianA < radianB; //Ascending order.
};
std::sort(m_points.begin(), m_points.end(), SortRuleFunc);
}
void ConvexHull::Generate()
{
const int pointsCount = m_points.size();
if (pointsCount < 3)
{
ErrorThrow("Points count too less.");
return;
}
SortPoints();
std::stack<const Point2D<int>> stackCHpt;
stackCHpt.push(m_points[0]);
stackCHpt.push(m_points[1]);
int ptIndex = 2; //Total is m_points.size().
while (ptIndex < pointsCount)
{
const Point2D<int>& ptO = stackCHpt.top();
stackCHpt.pop();
const Point2D<int>& ptA = stackCHpt.top();
stackCHpt.push(ptO);
const Point2D<int>& ptB = m_points[ptIndex]; //Current point.
Vector2D<int> vecA, vecB;
vecA.Vector(ptO, ptA);
vecB.Vector(ptO, ptB);
//Vector B turn to Vector A.
double angle = Vector2D<int>::Angle(vecB, vecA);
if ((angle - 0 >= 0) && (angle - 180 <= doubleESP))
{
stackCHpt.push(ptB);
ptIndex++;
}
else
{
stackCHpt.pop(); //pop point O.
}
if (stackCHpt.size() < 2)
{
//If the sort rule is not ok...
ErrorThrow("Error stackCHpt size.");
return;
}
}
/*Over Generation.*/
//
int lengthCHPoint = stackCHpt.size();
m_CHpoints.clear();
for (int i = 0; i < lengthCHPoint; ++i)
{
const Point2D<int>& pt = stackCHpt.top();
m_CHpoints.push_back(pt);
stackCHpt.pop();
}
}
void ConvexHull::TestPoints()
{
AddPoint({ 0, 0 });
AddPoint({ 5, 1 });
AddPoint({ 9, 2 });
AddPoint({ 1, 1 });
AddPoint({ 13, 6 });
AddPoint({ 12, 7 });
}
File: Line.cpp (暂时没用到)
#include "Point2D.cpp"
#include <array>
#include <iostream>
#pragma once
template<typename T>
class Line
{
public:
struct Equation
{
public:
T A, B, C;
public:
//Ax + By + C = 0
Equation()
{
}
Equation(const Point2D<T>& ptA, const Point2D<T>& ptB)
{
(*this)(ptA, ptB);
}
public:
void operator() (const Point2D<T>& ptA, const Point2D<T>& ptB)
{
A = ptB.y - ptA.y;
B = ptA.x - ptB.x;
C = ptB.x * ptA.y - ptA.x * ptB.y;
}
double GetY(double X)const
{
if (B == 0)
{
return INT_MAX;
}
return -1.0 * (C + A * X) / B;
}
double GetX(double Y) const
{
if (A == 0)
{
return INT_MAX;
}
return -1.0 * (C + B * Y) / A;
}
void GetSlope(Vector2D<double>& __out slope) const
{
if (B == 0) //X === -1.0 * C / A;
{
slope = { 0, -1.0 * C / A };
return;
}
slope = {1, -1.0 * A / B};
}
};
private:
bool CheckLine(const Point2D<T>& ptA, const Point2D<T>& ptB) const
{
if (std::abs(ptA.x - ptB.x) <= doubleESP &&
std::abs(ptA.y - ptB.y) <= doubleESP)
{
return false;
}
return true;
}
public:
Line() : Line({0, 0}, {1, 0})
{
}
Line(const Point2D<T>& a, const Point2D<T>& b) : ptA(a), ptB(b)
{
if (false == CheckLine(a, b))
{
ErrorThrow("They are the same point.");
return;
}
equation(a, b);
}
Line(const Line& nl) :Line{ nl.ptA, nl.ptB }
{
}
Line(T xA, T yA, T xB, T yB) :Line({ xA, yA }, { xB, yB })
{
}
void Init(const Point2D<T>& a, const Point2D<T>& b)
{
if (!CheckLine(a, b))
{
ErrorThrow("They are the same point.");
return;
}
ptA = a;
ptB = b;
equation(a, b);
}
void Init(T xA, T yA, T xB, T yB)
{
if (!CheckLine({ xA, yA }, { xB, yB }))
{
ErrorThrow("They are the same point.");
return;
}
ptA = { xA, yA };
ptB = { xB, yB };
equation(ptA, ptB);
}
void TransVector(Vector2D<T>& __out vec, bool isReverse/*B -> A*/ = false) const
{
if (isReverse) /*B -> A*/
vec.Init(ptA.x - ptB.x, ptA.y - ptB.y);
else /*A -> B*/
vec.Init(ptB.x - ptA.x, ptB.y - ptA.y);
}
void operator= (const Line& nl)
{
Init(nl.ptA, nl.ptB);
}
void Midperpendicular(Line<double>& __out resLine) const
{
Point2D<double> midPoint = { 0.5 * (GetPointA().x + GetPointB().x), 0.5 * (GetPointA().y + GetPointB().y) };
Vector2D<double> slope;
equation.GetSlope(slope);
Point2D<double> ptA, ptB;
if (0 == slope.x)
{
//Perpendicular to the X-axis.
ptA = { midPoint.x - 1.0, midPoint.y };
}
else
{
double k2 = -1.0 / (slope.y / slope.x);
double b2 = midPoint.y - midPoint.x * k2;
ptA = { midPoint.x + 1.0, k2 * (midPoint.x + 1.0) + b2 }; //Y = kX + b;
}
ptB = { midPoint.x, midPoint.y };
resLine.Init(ptA, ptB);
}
static void Intersection(const Line& __in lineA, const Line& __in lineB, Point2D<double>& __out pt)
{
if (Line::IsCollinear(4, lineA.ptA, lineA.ptB, lineB.ptA, lineB.ptB) == true)
{
pt = { INT_MAX, INT_MAX };
return;
}
T aA, aB, aC;
T bA, bB, bC;
aA = lineA.equation.A;
aB = lineA.equation.B;
aC = lineA.equation.C;
bA = lineB.equation.A;
bB = lineB.equation.B;
bC = lineB.equation.C;
double k = static_cast<double>(aA * bB - aB * bA);
if (MyMathTools::IsEqual(k, 0.0))
{
pt = { INT_MAX, INT_MAX };
return;
}
pt = { (bB* -1 * aC - (-1) * bC * aB) / k, (aA * -1 * bC - (-1) * aC * bA) / k };
}
void Intersection(const Line& __in lineOther, Point2D<double>& __out pt) const
{
Line::Intersection(*this, lineOther, pt);
}
static bool IsCollinear(int ptCount, const Point2D<T>& pt1, const Point2D<T>& pt2, const Point2D<T>& pt3, ...)
{
Vector2D<T> vecAB, vecBC;
vecAB.Vector(pt1, pt2);
vecBC.Vector(pt2, pt3);
if (!MyMathTools::IsEqual(vecAB.x * 1.0 / vecBC.x, vecAB.y * 1.0 / vecBC.y))
{
return false;
}
va_list aps;
va_start(aps, pt3);
const Point2D<T>* pNextArg = &pt3;
bool res = true;
while (ptCount - 3 > 0)
{
const Point2D<T> ptX = va_arg(aps, const Point2D<T>);
Vector2D<T> vecAx;
vecAx.Vector(pt1, ptX);
if (!MyMathTools::IsEqual(vecAB.x * 1.0 / vecAx.x, vecAB.y * 1.0 / vecAx.y))
{
res = false;
break;
}
--ptCount;
}
va_end(aps);
return res;
}
const Point2D<T>& GetPointA() const
{
return ptA;
}
const Point2D<T>& GetPointB() const
{
return ptB;
}
private:
Point2D<T> ptA, ptB;
public:
Equation equation;
};
File: Triangle.h (暂时没用到)
#include "Point2D.cpp"
#include "Line.cpp"
#include "Circle.cpp"
#pragma once
class Triangle final
{
public:
Triangle();
Triangle(const Point2D<int>& ptA, const Point2D<int>& ptB, const Point2D<int>& ptC);
~Triangle();
public:
bool IsPointInTriangle(const Point2D<int>& pt) const;
void Circumcircle(Circle<double, double>& __out circle) const;
double Area() const;
private:
std::array<Point2D<int>, 3> m_vertices /*m_vertex*/;
std::array<Line<int>, 3> m_lines;
};
File: Triangle.cpp (暂时没用到)
#include "Triangle.h"
Triangle::Triangle() :Triangle({ 0, 0 }, { 0, 3 }, { 4, 0 })
{
}
Triangle::Triangle(const Point2D<int>& ptA, const Point2D<int>& ptB, const Point2D<int>& ptC)
{
if (Line<int>::IsCollinear(3, ptA, ptB, ptC))
{
ErrorThrow("It is Collinear.");
return;
}
m_vertices[0] = ptA;
m_vertices[1] = ptB;
m_vertices[2] = ptC;
m_lines[0].Init(ptA, ptB);
m_lines[1].Init(ptB, ptC);
m_lines[2].Init(ptC, ptA);
}
Triangle::~Triangle()
{
}
void Triangle::Circumcircle(Circle<double, double>& __out circle) const
{
Line<int> lAB, lBC;
Line<double> lAB_mp, lBC_mp; /*_mp: Midperpendicular.*/
lAB = m_lines[0];
lBC = m_lines[1];
lAB.Midperpendicular(lAB_mp);
lBC.Midperpendicular(lBC_mp);
Point2D<double> intersection;
Line<double>::Intersection(lAB_mp, lBC_mp, intersection);
if (MyMathTools::IsEqual(intersection.x, INT_MAX) &&
MyMathTools::IsEqual(intersection.y, INT_MAX))
{
ErrorThrow("Error Triangle: It's not a triangle perhaps.");
return;
}
circle.center = intersection;
circle.radius = intersection.Distance<int>(m_vertices[0]);
}
double Triangle::Area() const
{
Vector2D<int> vecBA, vecBC;
m_lines[0].TransVector(vecBA, true);
m_lines[1].TransVector(vecBC);
return 0.5 * Vector2D<int>::Cross(vecBA, vecBC);
}
bool Triangle::IsPointInTriangle(const Point2D<int>& pt) const
{
if (Line<int>::IsCollinear(3, pt, m_vertices[0], m_vertices[1]) ||
Line<int>::IsCollinear(3, pt, m_vertices[1], m_vertices[2]) ||
Line<int>::IsCollinear(3, pt, m_vertices[0], m_vertices[2])
)
{
return true;
}
Triangle trianglePAB(pt, m_vertices[0], m_vertices[1]);
Triangle trianglePBC(pt, m_vertices[1], m_vertices[2]);
Triangle trianglePAC(pt, m_vertices[0], m_vertices[2]);
const double thisArea = Area();
const double totalArea = trianglePAB.Area() + trianglePBC.Area() + trianglePAC.Area();
if (MyMathTools::IsEqual(thisArea, totalArea))
{
return true;
}
return false;
}
File: Circle.cpp (暂时没用到)
#include "Point2D.cpp"
#pragma once
template<typename TC, typename TR>
class Circle final
{
public:
Circle() :Circle({0,0}, 1)
{
}
Circle(const Point2D<TC>& centerPt, TR radius) :
center(centerPt),
radius(radius)
{
}
double Area()
{
return PI * radius * radius;
}
double Perimeter()
{
return 2 * PI * radius;
}
template<typename TP>
bool IsPointInCircle(const Point2D<TP>& pt)
{
return !MyMathTools::LessThan(radius, Point2D<TP>::Distance(center, pt));
}
public:
TR radius;
Point2D<TC> center;
};
File: Win32Project1.cpp (win32默认程序,测试用)
case WM_PAINT:
hdc = BeginPaint(hWnd, &ps);
{
ConvexHull ch;
const int count =215;
ch.AddRandomPoints(count, 1000, 600);
//ch.TestPoints();
ch.Generate();
std::vector<const Point2D<int>> points;
ch.GetAllPoints(points);
const int pointWidth = 2;
int i = 0;
for each (const Point2D<int>& eachPt in points)
{
const Point2D<int>& pt = eachPt;
Rectangle(hdc, pt.x - pointWidth, pt.y - pointWidth, pt.x + pointWidth, pt.y + pointWidth);
i++;
}
HBRUSH hBrush = CreateSolidBrush(RGB(255, 0, 0));
HBRUSH hOldBrush = (HBRUSH)SelectObject(hdc, hBrush);
std::vector<const Point2D<int>> chpoints;
ch.GetConvexHullPoints(chpoints);
const int pointWidthCH = 3;
for each (const Point2D<int>& eachPt in chpoints)
{
const Point2D<int>& pt = eachPt;
Rectangle(hdc, pt.x - pointWidthCH, pt.y - pointWidthCH, pt.x + pointWidthCH, pt.y + pointWidthCH);
}
for (size_t i = 0; i < chpoints.size() - 1; ++i)
{
const Point2D<int>& eachPt = chpoints[i];
MoveToEx(hdc, eachPt.x, eachPt.y, NULL);
LineTo(hdc, chpoints[i + 1].x, chpoints[i + 1].y);
}
MoveToEx(hdc, (chpoints.end() - 1)->x, (chpoints.end() - 1)->y, NULL);
LineTo(hdc, chpoints[0].x, chpoints[0].y);
SelectObject(hdc, hOldBrush);
DeleteObject(hBrush);
}
//{
// Triangle trg;
// Circle<double, double> cc;
// trg.Circumcircle(cc);
//}
EndPaint(hWnd, &ps);
break;