【动态规划、计算几何】最优三角剖分

整理的算法模板合集： ACM模板

最优三角剖分

问题描述：

给一个有n个顶点的凸多边形，有很多方法进行三角剖分(polygon triangulation) 。给每个三角形规定一个权函数w(i,j,k)w(i,j,k)w(i,j,k)（比如三角形的周长或者三顶点的权和或者三角形的面积等等），求让所有三角形权和最大的方案。

　这个问题的关键在于状态的定义，因为如果允许随意切割，显然任意“半成品” 多边形的各个顶点可以是原多边形中随意选取的，很难简洁的定义成状态。但我们又可以发现，对于同一种切割方法，我们可以有多种切割顺序，但切割方法就已经决定了这个方法产生的结果，我们只要计算其中一种切割顺序就好了，所以不如把决策顺序规范化。

定义 d[i,j] 为从顶点 i 开始到顶点 j所构成的多边形的最大三角剖分权和（编号是逆时针分配的），则边i→j 在此多边形内一定恰好属于一个三角形 i→j→k ，所以多边形就被分成了三部分：三角形ikj 及其左右两部分的小多边形。这个切割方法保证了两个子多边形它们的顶点编号是连续的。所以我们只需要每次将代表这个多边形的两个顶点（始、末）作为分割后的其中一个三角形的一条边，然后枚举第三个顶点 k 的情况，再依次分割下去直到分割完成。

状态转移方程：d[i,j]=max{d[i,k]+d[k,j]+w(i,j,k)}d[i,j] = max\{d[i,k]+d[k,j]+w(i,j,k)\}d[i,j]=max{d[i,k]+d[k,j]+w(i,j,k)}
时间复杂度：O(n3)O(n^3)O(n3)
最终的答案：d[1][n−1]d[1][n - 1]d[1][n−1]

下面给出例题：

UVA1331 最大面积最小的三角剖分 Minimax Triangulation

本题要求的是最大面积最小的三角剖分，所以转移方程有所变化。

需要注意的是，本题并没有指明这个多边形的凹凸性，所以对于凹多边形的情况，需要加上一个判断，即判断Δijk是否是合法的三角形：如果Δijk内部含有多边形的一个顶点，那么就是不合法的三角形。

注意本题还会卡精度，所以需要判断精度误差

要注意我们这里使用的判断三角形是否合法使用的是O(n)O(n)O(n)的做法（枚举点 i(i≠a,i≠b,i≠c)i(i\neq a, i\neq b, i\neq c)i(i=a,i=b,i=c)，若△a,b,i+S△a,c,i+S△b,c,i=S△a,b,c\triangle_{a,b,i}+S\triangle_{a,c,i}+S\triangle_{b,c,i} = S\triangle_{a,b,c}△a,b,i+S△a,c,i+S△b,c,i=S△a,b,c，则点 i 位于 △a,b,c\triangle_{a,b,c}△a,b,c内部，则三角形 a,b,c 是不合法的。），使得总体的时间复杂度达到了O(n4)O(n^4)O(n4)。

#include<cstdio>
#include<cmath>
#include<algorithm>
#include<iostream>
#include<cstring>using namespace std;
typedef long long ll;
const int N = 507, M = 5000007, INF = 0x3f3f3f3f;
const double eps = 1e-6;struct Point{double x, y;
}a[N];double f[N][N];int n;double Area(double a, double b, double c)
{double p = (a + b + c) / 2.0;return sqrt(p * (p - a) * (p - b) * (p - c));
}double get_dist(Point a, Point b)
{return sqrt((a.x - b.x) * (a.x - b.x) + (a.y - b.y) * (a.y - b.y));
}bool check(int A, int B, int C)//判断凹多边形中切出来的是不是合法的三角形
{for(int i = 1; i <= n; ++ i){if(i != A && i != B && i != C){double s = Area(get_dist(a[A], a[B]), get_dist(a[B], a[C]), get_dist(a[A], a[C]));double p = Area(get_dist(a[A], a[B]), get_dist(a[B], a[i]), get_dist(a[A], a[i])) + Area(get_dist(a[A], a[i]), get_dist(a[i], a[C]), get_dist(a[A], a[C])) + Area(get_dist(a[i], a[B]), get_dist(a[B], a[C]), get_dist(a[i], a[C]));if(fabs(s - p) < eps)return false;}}return true;
}int main()
{int t;scanf("%d", &t);while(t -- ){scanf("%d", &n);for(int i = 1; i <= n; ++ i){scanf("%lf%lf", &a[i].x, &a[i].y);}for(int i = 1; i <= n - 2; ++ i){f[i][i + 2] = Area(get_dist(a[i], a[i + 1]), get_dist(a[i + 1], a[i + 2]), get_dist(a[i + 2], a[i]));}for(int len = 3; len < n; ++ len){//至少长度为3 至少三个点for(int i = 1; i <= n - len; ++ i){int j = i + len;f[i][j] = INF;for(int k = i + 1; k < j; ++ k){if(check(i, j, k)){//题目要求的是最大三角形最小，所以找到最大三角形转移f[i][j] = min(f[i][j], max(max(f[i][k], f[k][j]), Area(get_dist(a[i], a[k]), get_dist(a[k], a[j]), get_dist(a[j], a[i]))));}}}}printf("%.1f\n", f[1][n]);}
}

另一种方法是根据三角形三个点的两个向量的叉积是否小于0，如果小于0说明不合法，有点玄学，我代码也没调出来，就放一个那位想出来这个思路的大佬的代码吧。如果可以O(1)O(1)O(1)判断的话整体的时间复杂度就降到了O(n3)O(n^3)O(n3)，虽然数据只有 50…
大佬博客的链接：https://www.luogu.com.cn/blog/user17952/solution-uva1331

#include <iostream>
#include <cstdio>
#include <algorithm>using namespace std;const int Maxn = 0x3f3f3f3f;
const int N = 55;
int nxt[N], f[N][N], n, m, Ans;struct point
{int x, y;point() {}point(int X, int Y):x(X), y(Y) {}friend inline point operator - (const point &a, const point &b) {return point(b.x - a.x, b.y - a.y); }friend inline int operator * (const point &a, const point &b){return a.x * b.y - b.x * a.y;}
}a[N];inline int Max(int x, int y) {return x > y ? x : y;}
inline int Min(int x, int y) {return x < y ? x : y;}
inline void CkMin(int &x, int y) {if (x > y) x = y;}inline int dp(int l, int r)
{if (f[l][r] != -1) return f[l][r];if (r == nxt[l]) return f[l][r] = 0;int res = Maxn;for (int i = nxt[l]; i != r; i = nxt[i]){int tmp = (a[i] - a[r]) * (a[i] - a[l]);if (tmp > 0) CkMin(res, Max(Max(dp(l, i), tmp), dp(i, r)));}return f[l][r] = res;
}int main()
{while (scanf("%d", &m) != EOF) {while (m--) { scanf("%d", &n);for (int i = 1; i <= n; ++i) scanf("%d%d", &a[i].x, &a[i].y);for (int i = 1; i < n; ++i)nxt[i] = i + 1; nxt[n] = 1;int sum = 0; for (int i = 1; i < n; ++i)sum += a[i] * a[i + 1]; sum += a[n] * a[1];if (sum < 0)for (int i = 1, im = n >> 1; i <= im; ++i) swap(a[i], a[n - i + 1]);for (int i = 1; i <= n; ++i)for (int j = 1; j <= n; ++j)f[i][j] = -1;for (int i = 1; i <= n; ++i)for (int j = 1; j <= n; ++j)    if (i != j) dp(i, j);Ans = Maxn;for (int i = 1; i <= n; ++i)for (int j = 1; j <= n; ++j)if (i != j) CkMin(Ans, Max(f[i][j], f[j][i]));printf("%.1lf\n", (double)Ans / 2.0);}}return 0;
}