前言
这道题目第一个测资比较简单,但第二个测资我很坚持求公式解但还是失败了,但至少有得到角度,还好最后发现官方提供一些讯息(比赛结束的解析),但发现其实题目就有说算凸包了,虽然还有更快的作法但没有去实作它,最后也多学到了一个新的演算法。
原题目
题目
一驾神秘的太空船出现在这个天空中,天空是三维空间的平面,在y = -3 km与xz平面平行。外星人的船是一个边长1公里的实心立方体,中心(0km,0km,0km),八个角(+/- 0.5km,+ / - 0.5km,+ / - 0.5km) 。这艘船投影了一个影子,阴影是立方体在平面上的正交投影。(太阳是沿着y轴在天空上方无限远的点。)
只要外星人满足军方的要求,军方就愿意容忍这艘船:阴影面积覆盖只可接近Akm^2。我们了解船只不能改变尺寸,传的中心不能移动但是传可以任意旋转到位。
请找到一个方式,使阴影面积最接近外星人可旋转角度。使用三个点表达你的旋转。
只要外星人满足他们的官僚要求,军方就愿意容忍这艘船:阴影必须覆盖可接受地接近A km 2的飞机区域(参见输出部分以获得精确的定义)。他们聘请了一位几何语言学专家来向外国人传达这种需求。在你的通讯中,到目前为止,你已经了解到船只不能改变尺寸,船的中心不能移动,但船只能够任意旋转到位。
请找到一种方式,使阴影的面积接近外星人可以旋转角度。(x.y.z的座标)。
输入
输入第一行T为测试笔数,A为所需阴影面积km^2单位,小数点六位。
至少会有一种方法可以旋转太空船。
输出
每个测试,首先输出一行Case #x:,x是测试编号(从1开始)。然后,输出另外三行,如上所述,三个提供的面中心之一的x,y和z坐标。可使用小数(例如,0.000123456)或科学记号(例如,1.23456e-4)。
当且仅当满足以下所有条件时,您的答案才是是正确的:
1.从每个点到原点的距离(以km为单位)必须介于0.5 - 10^-6和0.5 + 10^6之间。
2.连接原点到每个点的段之间的角度(以弧度表示)必须在π / 2 - 10^6和π / 2 + 10^-6之间。
3.阴影区域(以km^2为单位),通过将所有8个顶点投影到y = -3平面上并找到那些投影点的凸包区域计算得出,必须介于A - 10^-6和 A + 10^-6之间。(后面不太了解意思大概是说它答案检测是用向量加法之类)
请注意,您可能需在小数点后输出超过6位数才能通过。如果有多个可接受的答案,您可以输出其中任何一个。
範围
1 ≤ T ≤ 100.
时间限制: 30秒.
记忆体限制: 1GB.
测试1
1.000000 ≤ A ≤ 1.414213
测试2
1.000000 ≤ A ≤ 1.732050
範例
输入
2
1.000000
1.414213
输出
Case #1:
0.5 0 0
0 0.5 0
0 0 0.5
Case #2:
0.3535533905932738 0.3535533905932738 0
-0.3535533905932738 0.3535533905932738 0
0 0 0.5
解析
这题主要就是要我们求它x.y.z旋转矩阵,图片来源3D小画家。
所需公式
三角函数:
cos(角度(弧度)) = 邻边 / 对边
角度弧度 = acos(邻边 / 对边)
旋转矩阵:
这部分因为题目是投影到xz平面因此我的想法是把z和y旋转矩阵做交换(应该可以用原本的不用这样做)。
依序上到下左到右都是x.y.z位置。
绕着x轴 = [1, 0, 0]
[0, cos(度), sin(度)]
[0, -sin(度), cos(度)]
绕着y轴 = [cos(度), 0, -sin(度)]
[0, 1, 0]
[sin(度), 0, cos(度)]
绕着z轴 = [cos(度), sin(度), 0]
[-sin(度), cos(度), 0]
[0, 0, 1]
旋转矩阵
凸包位置搜寻:
利用向量叉积来取得外围的点,下面演算法会详细说明。
礼物包装演算法
凸多边形面积计算:
主要是计算两点x.y的差,就可以取的三角形的底和高,即可计算面积。
面积推算过程
1.测试一,我们可以如下图来得知最大面积为sqrt(2)。
图一未旋转度,黑色斜线的右下角往左画过去(垂直)为45度角,长度为1km。图二旋转25度,黑色斜线的右下角往左画过去(垂直)为(45-25)度角,长度为1.3289...km。图三旋转45度,黑色斜线的右下角画过去(垂直)为(45 - 45)度,长度为sqrt(2)。第一种算法:
我们乘上对边就可得到邻边。邻边 = cos(45 - 角度) * sqrt(2)
会发现阴影面积的长*宽 = 邻边 * 1 = A面积。
因此我们可以求反三角函数求出角度求出 角度 = 45度 -acos(邻边 / 对边)。
第二种算法(可跳过):
虽然这算法多此一举就是了,但当时我忘记用45度扣掉所以用另一个角去算,临时想到的解决方法但不推荐。
首先图四灰色线=A,黄色线=B,则我们可以知道A+B=面积,和毕氏定理a^2 + b^2 = sqrt(2)^2,用二元一次方程式解,求得:
a = abs(((-2 * 面积) + sqrt(pow((-2 * 面积), 2) - 8 * (面积 * 面积 - 1))) / 4),算出来的为灰色线,在使用反三角函数得 角度 = asin(A / sqrt(2))
再将角度带入旋转公式即可(这部分答案用原始的旋转公式也可)
2.测试二,测试二部分我们计算每个角的位置在投影在XZ平面上,使用礼物包装法点连起来(六边形),计算面积,使用二分法取得精準度到最少小数点6位。
这部分我延续上面的图3的位置在往上转N度会等于最大面积sqrt(3),如图五,因为最大面积是一个正六边形,用面积推算回去得到的边长是1 / sqrt(6)即是图五黑色底线,这时候得知 角度 = acos(邻边 / 斜边),得到最大角度为35.2644度,因此就可以用二分法从0度到35.2644取的我们的面积来比较。
完整程式码
#include <iostream>#include <vector>#include "math.h"#include <algorithm>using namespace std;typedef long double LONG_D;typedef vector<vector<LONG_D>> Matrix;enum MatrixType{DEFAULT,RotatX,RotatY,RotatZ,};struct Point{Point() : _x(0.0), _y(0.0) {}Point(const LONG_D x, const LONG_D y) : _x(x), _y(y) {}LONG_D _x;LONG_D _y;};const LONG_D PI = atan(1) * 4;//F, G, H, E, B, C, D, Aconst Matrix _cubicPoint = { {0.5, 0.5, 0.5}, { 0.5, 0.5, -0.5 }, { -0.5, 0.5, -0.5 }, {-0.5, 0.5, 0.5 },{ 0.5, -0.5, 0.5 }, { 0.5, -0.5, -0.5 }, { -0.5, -0.5, -0.5 }, { -0.5, -0.5, 0.5, } };inline const LONG_D Ang2Rad(const LONG_D& angle){return angle / 180.0 * PI;}inline const LONG_D Rad2Ang(const LONG_D& radian){return radian / PI * 180.0;}/* 向量叉积向量。,大于180度叉积 > 0、大于180度叉积 < 0、等于 180° 叉积 = 0|a|向量 |b| 向量|a| x |b| = {{|i|, |k|, |j|},{x1, y1, 0},{x2, y2, 0} }解行列式 |j| * (x1 * y2 - x2 * y2) + 0 * (...) + 0 * (....)推算出来公式: |j| * (x1 * y2 - x2 * y2) 内积大于0度(小于等于则是超过180度,刚好可拿来使用)*/inline const bool Cross(const Point& p1, const Point& vertex, const Point& p2){return ((p2._x - vertex._x) * (p1._y - vertex._y) - (p2._y - vertex._y) * (p1._x - vertex._x)) > 0.0;}Matrix GetMatrix(const LONG_D& radian, const MatrixType type){Matrix matrix;switch (type){case RotatX:matrix = { { 1, 0, 0 },{ 0, cos(radian), sin(radian) },{ 0, -sin(radian), cos(radian) } };break;case RotatY:matrix = { { cos(radian), 0, -sin(radian) },{ 0, 1, 0 },{ sin(radian), 0, cos(radian) } };break;case RotatZ:matrix = { { cos(radian), sin(radian), 0 },{ -sin(radian), cos(radian), 0 },{ 0, 0, 1 } };break;case DEFAULT:matrix = { { 0.5, 0, 0 },{ 0, 0.5, 0 },{ 0, 0, 0.5 } };break;}return matrix;}Matrix MulMatrix(const Matrix& matrix1, const Matrix& matrix2){Matrix mul(matrix1.size(), vector<LONG_D>(matrix2[0].size(), 0));for (size_t row = 0; row < matrix1.size(); row++){for (size_t col = 0; col < matrix2[0].size(); col++){for (size_t index = 0; index < matrix2.size(); index++){mul[row][col] += (matrix1[row][index] * matrix2[index][col]);}}}return mul;}vector<LONG_D> MulPoint(const Matrix& matrix1, const vector<LONG_D>& point){vector<LONG_D> mul(matrix1.size(), 0);for (size_t row = 0; row < matrix1.size(); row++){for (size_t col = 0; col < point.size(); col++){mul[row] += (matrix1[row][col] * point[col]);}}return mul;}vector<Point> GetHullPoint(const vector<Point>& rotatePoint){vector<Point> hullPoint;int start = 0;for (size_t index = 0; index < rotatePoint.size(); index++){if (rotatePoint[index]._x < rotatePoint[start]._x|| (rotatePoint[index]._x == rotatePoint[start]._x && rotatePoint[index]._y < rotatePoint[start]._y)){start = index;}}int findEnd = start;while (true){hullPoint.push_back(rotatePoint[findEnd]);int find = start;for (size_t index = 0; index < rotatePoint.size(); index++){// 判断凸型顶点if (find == findEnd || Cross(rotatePoint[index], hullPoint[hullPoint.size() - 1], rotatePoint[find])){find = index;}}findEnd = find;if (start == find){break;}}return hullPoint;}LONG_D GetHullArea(const vector<Point>& point){LONG_D area = 0;area += abs(point[point.size() - 1]._x * point[0]._y - point[point.size() - 1]._y * point[0]._x);for (size_t index = 1; index < point.size(); index++){area += abs(point[index - 1]._x * point[index]._y - point[index - 1]._y * point[index]._x);}area *= 0.5;return area;}LONG_D GetArea(const Matrix& rotateZ, const LONG_D& radian){const Matrix rotateX = GetMatrix(radian, RotatX);vector<Point> rotatePoint(_cubicPoint.size());// 原始座标, 对Z轴心旋转45度, 对X轴心旋转radian弧度, 投影在X.Z, 取得2D座标for (size_t index = 0; index < _cubicPoint.size(); index++){vector<LONG_D> point = MulPoint(rotateX, rotateZ[index]);rotatePoint[index] = Point(point[0], point[2]);}vector<Point> hullPoint = GetHullPoint(rotatePoint);LONG_D area = GetHullArea(hullPoint);return area;}Matrix Compute(const LONG_D& area){if (area <= sqrt(2)){//const LONG_D aLen = abs(((-2 * area) + sqrt(pow((-2 * area), 2) - 8 * (area * area - 1))) / 4); // a^2.....//const LONG_D radian = asin(aLen);LONG_D radian = Ang2Rad(45.0) - abs(acos(area / sqrt(2)));const Matrix rotateZ = GetMatrix(radian, RotatZ);return MulMatrix(rotateZ, GetMatrix(0, DEFAULT));}const Matrix rotateZ = GetMatrix(Ang2Rad(45.0), RotatZ);Matrix rotateZPoint(8);for (size_t index = 0; index < _cubicPoint.size(); index++){rotateZPoint[index] = MulPoint(rotateZ, _cubicPoint[index]);}LONG_D minRadian = Ang2Rad(0);LONG_D maxRadian = Ang2Rad(35.2644);LONG_D nowRadian = maxRadian;LONG_D nowArea = GetArea(rotateZPoint, maxRadian);while (abs(area - nowArea) > 0.00000000001){const LONG_D middleRadian = (minRadian + maxRadian) / 2.0;const LONG_D middleArea = GetArea(rotateZPoint, middleRadian);if (area >= middleArea){minRadian = middleRadian;}else{maxRadian = middleRadian;}nowRadian = middleRadian;nowArea = middleArea;}Matrix rotateX = GetMatrix(nowRadian, RotatX);Matrix matrix (3);for (int index = 0; index < 3; index++){vector<LONG_D> center(3, 0);center[index] = 0.5;matrix[index] = MulPoint(rotateX, MulPoint(rotateZ, center));}return matrix;}void Print(Matrix& matrix, const int index){cout.precision(16);cout << "Case #" << index + 1 << ":" << endl;for (size_t row = 0; row < matrix.size(); row++){for (size_t col = 0; col < matrix.size(); col++){cout << matrix[row][col] << " ";}cout << endl;}}void Input(){int testCount = 0;while (cin >> testCount){LONG_D area = 0L;for (int index = 0; index < testCount; index++){cin >> area;Matrix result = Compute(area);Print(result, index);}}}int main(){Input();return 0;}
程式码解析
前置作业
1.宣告两个别名LONG_D和Matrix,方便观看。
2.宣告旋转类型MatrixType。
3.建立一个结构Point用来记录2D座标。
4.计算PI的弧度(因C++都是用弧度运算)。
5.宣告和指派8个顶点位置给_cubicPoint如下图
注:题目说八个角为+-0.5km。
6.Ang2Rad角度转弧度。
7.Rad2Ang弧度转角度。
8.GetMatrix部分为得取旋转矩阵或原始x.y.z座标矩阵。
注:原始为何是0.5? 题目说原点是(0,0,0)是顶点中心,也就是顶点到顶点 / 2 = 原始座标。
我放大到200如下方程式码的x.y.z的座标,并画出它在2D中的显示,如下图。
3D转2D公式:
X = 原点X - cos(45度) * x(X水平投影) + y大小
Y = 原点Y + sin(45度) * x(X在垂直向上的投影) - z大小
因为Z和X如果垂直那Y就会往45度方向所以才会这样计算,可以自己推算看看。
----------C#--3D转2D可视化START---------
//x = origin - cos(45) * x + y; //原点X - cos(45度) * x(X水平投影) + y大小//y = origin + sin(45) * x - z; //原点Y + sin(45度) * x(X在垂直向上的投影) - z大小public Point2D GetPoint2D(Point3D origin){ int x = (int)(origin.X - Math.Cos(45) * _x + _y); int y = (int)(origin.Y + Math.Sin(45) * _x - _z); return new Point2D(x, y);}Point3D origin = new Point3D(300, 300, 300);Point3D xEnd = new Point3D(200, 0, 0);Point3D yEnd = new Point3D(0, 200, 0);Point3D zEnd = new Point3D(0, 0, 200);Point2D xLine = xEnd.GetPoint2D(origin);e.Graphics.DrawLine(pen, origin.X, origin.Y, xLine.X, xLine.Y);Point2D yLine = yEnd.GetPoint2D(origin);e.Graphics.DrawLine(pen, origin.X, origin.Y, yLine.X, yLine.Y);Point2D zLine = zEnd.GetPoint2D(origin);e.Graphics.DrawLine(pen, origin.X, origin.Y, zLine.X, zLine.Y);
----------C#--3D转2D可视化END---------
9.MulMatrix为矩阵相乘,左矩阵的水平列每一个对上右矩阵垂直行两两相乘再加总。
如下範例:
[1, 2, 3] [1, 3, 4]
A = [4, 5, 6] X [2, 2, 2]
[4, 5, 6] [3, 1, 5]
A[0][0] = (1 * 1) + (2 * 2) + (3 * 3)
A[0][1] = (1 * 3) + (2 * 2) + (3 * 1)
A[0][2] = (1 * 4) + (2 * 2) + (3 * 5)
A[1][0] = (4 * 1) + (5 * 2) + (6 * 3)
.....略。
10.MulPoint是计算顶点在矩阵旋转后的座标如上一样,只是右矩阵只有一列水平列。
typedef long double LONG_D;typedef vector<vector<LONG_D>> Matrix;enum MatrixType{DEFAULT,RotatX,RotatY,RotatZ,};struct Point{Point() : _x(0.0), _y(0.0) {}Point(const LONG_D x, const LONG_D y) : _x(x), _y(y) {}LONG_D _x;LONG_D _y;};const LONG_D PI = atan(1) * 4;//F, G, H, E, B, C, D, A正方体顶点位置const Matrix _cubicPoint = { {0.5, 0.5, 0.5}, { 0.5, 0.5, -0.5 }, { -0.5, 0.5, -0.5 }, {-0.5, 0.5, 0.5 },{ 0.5, -0.5, 0.5 }, { 0.5, -0.5, -0.5 }, { -0.5, -0.5, -0.5 }, { -0.5, -0.5, 0.5, } };// 角度转弧度inline const LONG_D Ang2Rad(const LONG_D& angle){return angle / 180.0 * PI;}// 弧度转角度inline const LONG_D Rad2Ang(const LONG_D& radian){return radian / PI * 180.0;}// 得取旋转矩阵或原点位置Matrix GetMatrix(const LONG_D& radian, const MatrixType type){Matrix matrix;switch (type){case RotatX:matrix = { { 1, 0, 0 },{ 0, cos(radian), sin(radian) },{ 0, -sin(radian), cos(radian) } };break;case RotatY:matrix = { { cos(radian), 0, -sin(radian) },{ 0, 1, 0 },{ sin(radian), 0, cos(radian) } };break;case RotatZ:matrix = { { cos(radian), sin(radian), 0 },{ -sin(radian), cos(radian), 0 },{ 0, 0, 1 } };break;case DEFAULT:matrix = { { 0.5, 0, 0 },{ 0, 0.5, 0 },{ 0, 0, 0.5 } };break;}return matrix;}// 矩阵相乘Matrix MulMatrix(const Matrix& matrix1, const Matrix& matrix2){Matrix mul(matrix1.size(), vector<LONG_D>(matrix2[0].size(), 0));for (size_t row = 0; row < matrix1.size(); row++){for (size_t col = 0; col < matrix2[0].size(); col++){for (size_t index = 0; index < matrix2.size(); index++){mul[row][col] += (matrix1[row][index] * matrix2[index][col]);}}}return mul;}// 顶点乘上旋转矩阵vector<LONG_D> MulPoint(const Matrix& matrix1, const vector<LONG_D>& point){vector<LONG_D> mul(matrix1.size(), 0);for (size_t row = 0; row < matrix1.size(); row++){for (size_t col = 0; col < point.size(); col++){mul[row] += (matrix1[row][col] * point[col]);}}return mul;}
输入
1.testCount测试资料T笔数。
2.area为面积大小。
3.Compute函数开始计算旋转矩阵。
4.Print输出回传结果。
void Input(){int testCount = 0;while (cin >> testCount){LONG_D area = 0;for (int index = 0; index < testCount; index++){cin >> area;Matrix result = Compute(area);Print(result, index);}}}
计算测资1
1.如果面积小于等于sqrt(2)(测资一)可以用我们分析的公式直接解决。
2.先把45度转弧度再带入公式得到角度。
3.我们选择绕着Z轴转,得取Z轴的样板rotateZ。
4.MulMatrix将Z轴样板乘上原始的座标即是我们所要的答案(这部分其实应该座标对旋转矩阵乘)。
Matrix Compute(const LONG_D& area){if (area <= sqrt(2)){//const LONG_D aLen = abs(((-2 * area) + sqrt(pow((-2 * area), 2) - 8 * (area * area - 1))) / 4); // a^2.....//const LONG_D radian = asin(aLen);LONG_D radian = Ang2Rad(45.0) - abs(acos(area / sqrt(2)));const Matrix rotateZ = GetMatrix(radian, RotatZ);return MulMatrix(rotateZ, GetMatrix(0, DEFAULT));}}
计算测资2
若大于sqrt(2)我们就要旋转两次,这部分可以优化(官方的分析有影片)。
1.首先我先直接取得我们刚刚的旋转Z轴模板=第一次旋转。
2.Z轴模板乘上rotateZPoint 8个顶点,得取目前8个顶点的三维座标。
3.设定目前最大最小的弧度与目前弧度(minRadian.maxRadian.nowRadian)和面积(nowArea)。
4.使用二分搜寻,取得两弧度中间,GetArea计算面积(若要理解下方还有说明函数)。
若A面积>中间面积则範围缩小为中间值~最大面积,反之,A面积<中间面积则範围缩小为最小面积~中间值,再将目前面积设定为中间面积,直到精準度到达6位。
5.使用我们计算好的弧度得取旋转X轴的模板,这里使用迴圈跑三次将x.y.z原本的三个点先乘上Z再乘上X(不得交换先旋转哪一个就必须先乘上哪个矩阵)。
Matrix Compute(const LONG_D& area){const Matrix rotateZ = GetMatrix(Ang2Rad(45.0), RotatZ);Matrix rotateZPoint(8);for (size_t index = 0; index < _cubicPoint.size(); index++){rotateZPoint[index] = MulPoint(rotateZ, _cubicPoint[index]);}LONG_D minRadian = Ang2Rad(0);LONG_D maxRadian = Ang2Rad(35.2644);LONG_D nowRadian = maxRadian;LONG_D nowArea = GetArea(rotateZPoint, maxRadian);while (abs(area - nowArea) > 0.00000000001){const LONG_D middleRadian = (minRadian + maxRadian) / 2.0;const LONG_D middleArea = GetArea(rotateZPoint, middleRadian);if (area >= middleArea){minRadian = middleRadian;}else{maxRadian = middleRadian;}nowRadian = middleRadian;nowArea = middleArea;}Matrix rotateX = GetMatrix(nowRadian, RotatX);Matrix matrix (3);for (int index = 0; index < 3; index++){vector<LONG_D> center(3, 0);center[index] = 0.5;matrix[index] = MulPoint(rotateX, MulPoint(rotateZ, center));}return matrix;}
得取旋转顶点座标
1.我们用传入的已旋转45度Z的八个顶点座标乘上目前我们预测radian的X旋转样板,因为我们投影在XZ平面所以我们只需要取的X和Z的做边转换到2D平面。
2.GetHullPoint取得投影在2D的凸多边型顶点座标(若要理解下方还有说明函数)。
3.GetHullArea取得凸多边形面积(若要理解下方还有说明函数)。
LONG_D GetArea(const Matrix& rotateZ, const LONG_D& radian){const Matrix rotateX = GetMatrix(radian, RotatX);vector<Point> rotatePoint(_cubicPoint.size());// 原始座标, 对Z轴心旋转45度, 对X轴心旋转radian弧度, 投影在X.Z, 取得2D座标for (size_t index = 0; index < _cubicPoint.size(); index++){vector<LONG_D> point = MulPoint(rotateX, rotateZ[index]);rotatePoint[index] = Point(point[0], point[2]);}vector<Point> hullPoint = GetHullPoint(rotatePoint);LONG_D area = GetHullArea(hullPoint);return area;}
得取凸多边形座标
1.GetHullPoint参照一个容器2D座标过去即可计算。
2.先宣告一个容器纪录顶点,和宣告并指派start用来记录第一个最左边的。
3.跑回圈搜寻目前最左边的点,就是搜寻最小x和y座标,并在指派给start。
4.在宣告一个findEnd来记录每次寻找凸点最左边的索引位置,先暂时指派最左边的索引值。
5.开始跑回圈直到我们跑回第一个最左边的点。
6.先把上一个我们找到的findEnd点加入,在宣告一个我们目前找到左边的索引位置find,暂时指派第一个最左边位置 给他。
7.跑回圈寻找最左边的点。
8.先判断find == findEnd因为上一个已经是顶点所以将当下索引设定为我们目前找到的在继续比较。
例如上一轮找到5当索引值到5时我们比较叉积时候他又在左边find就变成5,然而5已经在上一轮被设定过,我们必须跳过他先将6设定为我们找到的。
9.在判断叉积Cross传入两个2D座标和一个vertex顶点(上一个凸点)方可计算叉积(注解有打公式),下方我用我的方式解释。
假如vertex(1,2)是我们找到的上一个最左边的点。
* * * * * *
* * x * * *
* * * x * *
* 1 * * * *
* * * * * *
* * * * * *
公式 ((p2._x - vertex._x) * (p1._y - vertex._y)
- (p2._y - vertex._y) * (p1._x - vertex._x))
1.vertex = (1,2)
p1 = (3,3)
p2 = (2,4)
(2 - 1) * (3 - 2) - (4 - 2) * (3 - 2) = -1
2.vertex = (1,2)
p1 = (2,4)
p2 = (3,3)
(3 - 1) * (4 -2) - (3 - 2) * (2 - 2) = 2 * 2 = 4
第一假如(3,3)是我们现在找到的点,(2,4)是上个被找到的点,算出结果为负数。
第二假如(2,4)是我们现在找到的点,(3,3)是上个被找到的点,算出结果为正数。
这是运用叉积运算来计算出向量的差,因为180度的向量不会是负的。
若还需要更详细正确的原理可上网搜寻。
10.若已经搜寻完我们将findEnd指派find给他,若find等于我们第一个最左边的代表凸点都已经搜寻到,即可结束。
vector<Point> GetHullPoint(const vector<Point>& rotatePoint){vector<Point> hullPoint;int start = 0;for (size_t index = 0; index < rotatePoint.size(); index++){if (rotatePoint[index]._x < rotatePoint[start]._x|| (rotatePoint[index]._x == rotatePoint[start]._x && rotatePoint[index]._y < rotatePoint[start]._y)){start = index;}}int findEnd = start;while (true){hullPoint.push_back(rotatePoint[findEnd]);int find = start;for (size_t index = 0; index < rotatePoint.size(); index++){// 判断凸型顶点if (find == findEnd || Cross(rotatePoint[index], hullPoint[hullPoint.size() - 1], rotatePoint[find])){find = index;}}findEnd = find;if (start == find){break;}}return hullPoint;}/* 向量叉积向量。,大于180度叉积 > 0、大于180度叉积 < 0、等于 180° 叉积 = 0|a|向量 |b| 向量|a| x |b| = {{|i|, |k|, |j|},{x1, y1, 0},{x2, y2, 0} }解行列式 |j| * (x1 * y2 - x2 * y2) + 0 * (...) + 0 * (....)推算出来公式: |j| * (x1 * y2 - x2 * y2) 内积大于0度(小于等于则是超过180度,刚好可拿来使用)*/inline const bool Cross(const Point& p1, const Point& vertex, const Point& p2){return ((p2._x - vertex._x) * (p1._y - vertex._y) - (p2._y - vertex._y) * (p1._x - vertex._x)) > 0.0;}
计算凸多边型
1.GetHullArea传入凸点座标。
2.我们去计算点到点的x.y在算出三角形面积,上面有推倒文章公式。
例如,我们有A座标B座标C座标D座标。
开始计算A->B两点的x和y所组成的三角形,在计算B->C两点,在计算C->D两点,最后在计算D->A两点。
你可以想像它计算方式就是切割出很多个长方形在计算三角形面积。
若想更了解可以详看过程,就可得知为何要这样做。
面积推算过程
LONG_D GetHullArea(const vector<Point>& point){LONG_D area = 0;area += abs(point[point.size() - 1]._x * point[0]._y - point[point.size() - 1]._y * point[0]._x);for (size_t index = 1; index < point.size(); index++){area += abs(point[index - 1]._x * point[index]._y - point[index - 1]._y * point[index]._x);}area *= 0.5;return area;}
结语
这题目最快解法大概是计算向量或是运用正方体最高点的距离差去计算,但这些都牵扯到满多数学,对我平常很少在用几何的人可能要花更多时间去了解,但最开心的是多学到一种演算法(礼物包装演算法)。
这篇文章实在有够长,仔细看其实每个部份就是一小块演算法组合起来,打那么长
有笔误或理解错的希望大家包涵。