Разбить выпуклый многоугольник на треугольники

#include <cassert>
#include <cmath>
#include <vector>
#include <algorithm>
#include <limits>
#include <iostream>
 
struct Point 
{ 
  int x, y; 
};
 
double distance(const Point &p1, const Point &p2)
{
  Point delta = { p2.x - p1.x, p2.y - p1.y };
  return std::sqrt((double) delta.x * delta.x + (double) delta.y * delta.y);
}
 
struct PartialSolution 
{
  unsigned k = -1;
  double weight;
};
 
typedef std::vector<std::vector<PartialSolution>> PartialSolutions;
 
double MWT(const Point poly[], unsigned i, unsigned j, PartialSolutions &ps)
{
  unsigned n = ps.size();
  assert(n > 0 && i < j && j < n);
 
  if (j == i + 1)
    // Исходное ребро - возвращаем вес 0
    return 0;
 
  PartialSolution &solution = ps[i][j];
  if (solution.k != -1)
    // Такую подзадачу мы уже решали - возвращаем готовое решение
    return solution.weight;
 
  // Такую подзадачу мы еще не решали - решаем ее полностью
 
  // Вес диагонали
  double weight_ij = i > 0 || j < n - 1 ? distance(poly[i], poly[j]) : 0;
  // Не забываем, что диагональ [0, n - 1] - это исходное ребро и ее вес равен 0
 
  double min_weight = std::numeric_limits<double>::max();
  unsigned min_k;
 
  for (unsigned k = i + 1; k < j; ++k)
  {
    double weight = MWT(poly, i, k, ps) + MWT(poly, k, j, ps);
    if (weight < min_weight)
    {
      min_weight = weight;
      min_k = k;
    }
  }
 
  // Решение нашей подзадачи
  double weight = weight_ij + min_weight;
  
  // Запоминаем решение на будущее
  solution.k = min_k;
  solution.weight = weight;
 
  return weight;
}
 
double MWT(const Point poly[], PartialSolutions &ps)
{
  unsigned n = ps.size();
  assert(n > 0);
  return MWT(poly, 0, n - 1, ps);
}
 
void MWT_diagonals(const Point poly[], unsigned i, unsigned j, const PartialSolutions &ps)
{
  unsigned n = ps.size();
  assert(n > 0 && i < j && j < n);
 
  if (j == i + 1)
    return;
    
  if (i > 0 || j < n - 1)
    std::cout << "(" << i << " - " << j << ") - " << distance(poly[i], poly[j]) << std::endl;
 
  const PartialSolution &solution = ps[i][j];
  assert(solution.k != -1);
 
  MWT_diagonals(poly, i, solution.k, ps);
  MWT_diagonals(poly, solution.k, j, ps);
}
 
void MWT_diagonals(const Point poly[], const PartialSolutions &ps)
{
  unsigned n = ps.size();
  assert(n > 0);
  MWT_diagonals(poly, 0, n - 1, ps);
}
 
int main()
{
  Point poly[] = { { 0, 0 }, { 5, -50 }, {20, -50}, { 200, 0 }, { 5, 50 } };
  unsigned n = sizeof poly / sizeof *poly;
 
  PartialSolutions ps(n);
  for (auto &s : ps) s.resize(n);
 
  std::cout << MWT(poly, ps) << std::endl;
 
  MWT_diagonals(poly, ps);
}