voronoi diagram (Tyson polygon) application – Good Manners

Welcome to pay attention to more exciting things
Follow me to learn common algorithms and data structures, solve multiple problems, and combat dimensionality reduction.

voronoi diagram solutionClick to go

Question link: https://vjudge.net/problem/URAL-1504

The main idea of the title

There is a table, which is round in shape.

There are many cakes on the table, each cake has its own size.

People can sit anywhere on the edge of the table. After sitting down, they can take the piece of cake closest to them. Find the position where you can get the largest cake as possible.

Thinking &Analysis

First find the voronoi area where each cake is located.

Check whether there is any intersection between the area and the circle.

Traverse the cakes with intersection points and select the intersection point that can get the largest cake.

Algorithm details

1) Determine whether the area intersects the circle
You can find all the vertices of the area and determine whether any vertices are on or outside the circle.

2) Find the intersection of a straight line and a circle

You can use projection to find point P first, and use the right triangle to know the length of PB’
From point P along the AB direction, you can get the B’ coordinate. Similarly, there is an intersection point in the opposite direction.

3) Special circumstances
When there are only 2 cakes or all cakes are on a line, triangulation cannot be performed.

At this time, select the largest cake and find the intersection point of the vertical line and the circle.

Algorithm process

step 1 Perform delaunay triangulation on boundary points
Step 2: Traverse each boundary point, collect adjacent edges, and sort them counterclockwise.
Solve for the mid-perpendicular of each adjacent side, and find intersections with the boundary and other perpendiculars respectively (only adjacent mid-perpendiculars have intersections) to determine whether there is an intersection with the circle.
Step 3 Find the intersection point with the circle and determine whether the intersection point is within the voronoi area.
step 4 Select the location where you can get the biggest cake

Code

#include<stdio.h>
#include<cmath>
#include <algorithm>
#include <vector>
#include <list>
#include <cstring>
#include <utility>


using namespace std;
const double EPS = 1e-8;

const int N = 1e6 + 10;
const int M = 1e6 + 10;

int cmp(double d) {<!-- -->
if (abs(d) < EPS)return 0;
if (d > 0)return 1;
return -1;
}

class Point {<!-- -->
public:
double x, y;
int id;

Point() {<!-- -->}
Point(double a, double b) :x(a), y(b) {<!-- -->}
Point(const Point & amp; p) :x(p.x), y(p.y), id(p.id) {<!-- -->}

void in() {<!-- -->
scanf("%lf %lf", & amp;x, & amp;y);
}
void out() {<!-- -->
printf("%f %f\
", x, y);
}

double dis() {<!-- -->
return sqrt(x * x + y * y);
}

double dis2() {<!-- -->
return x * x + y * y;
}

Point operator -() const {<!-- -->
return Point(-x, -y);
}

Point operator -(const Point & amp; p) const {<!-- -->
return Point(x - p.x, y - p.y);
}

Point operator + (const Point & amp; p) const {<!-- -->
return Point(x + p.x, y + p.y);
}
Point operator *(double d)const {<!-- -->
return Point(x * d, y * d);
}

Point operator /(double d)const {<!-- -->
return Point(x / d, y / d);
}


void operator -=(Point & amp; p) {<!-- -->
x -= p.x;
y -= p.y;
}

void operator + =(Point & amp; p) {<!-- -->
x + = p.x;
y + = p.y;
}
void operator *=(double d) {<!-- -->
x *= d;
y*=d;
}

void operator /=(double d) {<!-- -->
this ->operator*= (1 / d);
}

bool operator<(const Point & amp; a) const {<!-- -->
return x < a.x || (abs(x - a.x) < EPS & amp; & y < a.y);
}

bool operator==(const Point & amp; a) const {<!-- -->
return abs(x - a.x) < EPS & amp; & amp; abs(y - a.y) < EPS;
}
};

//Vector operations

double cross(const Point & amp; a, const Point & amp; b) {<!-- -->
return a.x * b.y - a.y * b.x;
}

double dot(const Point & amp; a, const Point & amp; b) {<!-- -->
return a.x * b.x + a.y * b.y;
}


class Point3D {<!-- -->
public:
double x, y, z;

Point3D() {<!-- -->}
Point3D(double a, double b, double c) :x(a), y(b), z(c) {<!-- -->}
Point3D(const Point3D & amp; p) :x(p.x), y(p.y), z(p.z) {<!-- -->}

double dis() {<!-- -->
return sqrt(x * x + y * y + z * z);
}

double dis2() {<!-- -->
return x * x + y * y + z * z;
}

Point3D operator -(const Point3D & amp; p) const {<!-- -->
return Point3D(x - p.x, y - p.y, z - p.z);
}

Point3D operator + (const Point3D & amp; p) const {<!-- -->
return Point3D(x + p.x, y + p.y, z + p.z);
}
Point3D operator *(double d)const {<!-- -->
return Point3D(x * d, y * d, z * d);
}

Point3D operator /(double d)const {<!-- -->
return Point3D(x / d, y / d, z / d);
}


void operator -=(Point3D & amp; p) {<!-- -->
x -= p.x;
y -= p.y;
z -= p.z;
}

void operator + =(Point3D & amp; p) {<!-- -->
x + = p.x;
y + = p.y;
z + = p.z;
}
void operator *=(double d) {<!-- -->
x *= d;
y*=d;
z *= d;
}

void operator /=(double d) {<!-- -->
this ->operator*= (1 / d);
}
};

//Vector operations
Point3D cross(const Point3D & amp; a, const Point3D & amp; b) {<!-- -->
return Point3D(a.y * b.z - a.z * b.y, -a.x * b.z + a.z * b.x,
a.x * b.y - a.y * b.x);
}

double dot(const Point3D & amp; a, const Point3D & amp; b) {<!-- -->
return a.x * b.x + a.y * b.y + a.z * b.z;
}


class Line {<!-- -->
public:
Point front, tail;
Line() {<!-- -->}
Line(Point a, Point b) :front(a), tail(b) {<!-- -->}
};

/*
0 disjoint
1 intersect
0 Parallel/Coincident
*/
int cross(const Line & amp; a, const Line & amp; b) {<!-- -->
Point dir1 = a.front - a.tail;
Point dir2 = b.front - b.tail;
if (cmp(cross(dir1, dir2)) == 0) {<!-- -->
return 0;
}

if (cmp(cross(a.front - b.tail, dir2)) * cmp(cross(a.tail - b.tail, dir2)) >= 0)return 0;
if (cmp(cross(b.front - a.tail, dir1)) * cmp(cross(b.tail - a.tail, dir1)) >= 0)return 0;
return 1;
}


int inCircle(Point p0, Point p1, Point p2, Point p3) {<!-- -->
Point d1 = p1 - p0;
Point d2 = p2 - p0;
if (cross(d1, d2) < 0)return inCircle(p0, p2, p1, p3); // Ensure that the normal direction of the plane is upward

// Build mapping point
Point3D lift0(p0.x, p0.y, p0.dis2());
Point3D lift1(p1.x, p1.y, p1.dis2());
Point3D lift2(p2.x, p2.y, p2.dis2());
Point3D lift3(p3.x, p3.y, p3.dis2());

Point3D z1(lift1 - lift0), z2(lift2 - lift0);
Point3D normal = cross(z1, z2); // Calculate plane normal direction
double project = dot(normal, lift3 - lift0); // Calculate the distance from the point to the plane

return cmp(project);
}



class EdgeDelaunay {<!-- -->
public:
int id;
std::list<EdgeDelaunay>::iterator c;
EdgeDelaunay(int id = 0) {<!-- --> this->id = id; }
};

class Delaunay {<!-- -->
public:
std::list<EdgeDelaunay> head[N]; // graph
Point p[N];
int n = 0;

void init(int psize, Point ps[]) {<!-- -->
this->n = psize;
memcpy(this->p, ps, sizeof(Point) * n);
std::sort(this->p, this->p + n);
divide(0, n - 1);
}

void addEdge(int u, int v) {<!-- -->
head[u].push_front(EdgeDelaunay(v));
head[v].push_front(EdgeDelaunay(u));
head[u].begin()->c = head[v].begin();
head[v].begin()->c = head[u].begin();
}

void divide(int l, int r) {<!-- -->
if (r - l <= 1) {<!-- --> // #point <= 2
for (int i = l; i <= r; i + + )
for (int j = i + 1; j <= r; j + + ) addEdge(i, j);
return;
}
int mid = (l + r) / 2;
divide(l, mid);
divide(mid + 1, r);

std::list<EdgeDelaunay>::iterator it;
int nowl = l, nowr = r;

for (int update = 1; update;) {<!-- -->
// Find the lowest line position on the left
update = 0;
Point ptL = p[nowl], ptR = p[nowr];
for (it = head[nowl].begin(); it != head[nowl].end(); it + + ) {<!-- -->
Point t = p[it->id];
double v = cross(ptL - ptR, t - ptR);
if (cmp(v) > 0 || (cmp(v) == 0 & amp; & amp; (t - ptR).dis() < (ptL - ptR).dis())) {<!-- -->
nowl = it->id, update = 1;
break;
}
}
if (update) continue;
// Find the lowest line position on the right
for (it = head[nowr].begin(); it != head[nowr].end(); it + + ) {<!-- -->
Point t = p[it->id];
double v = cross(ptR - ptL, t - ptL);
if (cmp(v) < 0 || (cmp(v) == 0 & amp; & amp; (t - ptL).dis() < (ptL - ptR).dis())) {<!-- -->
nowr = it->id, update = 1;
break;
}
}
}

addEdge(nowl, nowr); //Add baseline

for (; true;) {<!-- -->
Point ptL = p[nowl], ptR = p[nowr];
int ch = -1, side = 0;
for (it = head[nowl].begin(); it != head[nowl].end(); it + + ) {<!-- -->
if (cmp(cross(ptR - ptL, p[it->id] - ptL)) <= 0)continue; // Determine whether the angle is less than 180
if (ch == -1 || inCircle(ptL, ptR, p[ch], p[it->id]) < 0) {<!-- -->
ch = it->id, side = -1;
}
}
for (it = head[nowr].begin(); it != head[nowr].end(); it + + ) {<!-- -->
if (cmp(cross(p[it->id] - ptR, ptL - ptR)) <= 0) continue; // Determine whether the angle is less than 180
if (ch == -1 || inCircle(ptL, ptR, p[ch], p[it->id]) < 0) {<!-- -->
ch = it->id, side = 1;
}
}
if (ch == -1) break; // All lines have been added
if (side == -1) {<!-- -->
for (it = head[nowl].begin(); it != head[nowl].end();) {<!-- -->
// Determine whether they intersect. Edges are not considered to intersect.
if (cross(Line(ptL, p[it->id]), Line(ptR, p[ch]))) {<!-- -->
head[it->id].erase(it->c);
head[nowl].erase(it + + );
}
else {<!-- -->
it++;
}
}
nowl = ch;
addEdge(nowl, nowr);
}
else {<!-- -->
for (it = head[nowr].begin(); it != head[nowr].end();) {<!-- -->
// Determine whether they intersect. Edges are not considered to intersect.
if (cross(Line(ptR, p[it->id]), Line(ptL, p[ch]))) {<!-- -->
head[it->id].erase(it->c);
head[nowr].erase(it + + );
}
else {<!-- -->
it++;
}
}
nowr = ch;
addEdge(nowl, nowr);
}
}
}

std::vector<std::pair<int, int> > getEdge() {<!-- -->
std::vector<std::pair<int, int> > ret;
ret.reserve(n);
std::list<EdgeDelaunay>::iterator it;
for (int i = 0; i < n; i + + ) {<!-- -->
for (it = head[i].begin(); it != head[i].end(); it + + ) {<!-- -->
ret.push_back(std::make_pair(p[i].id, p[it->id].id));
}
}
return ret;
}
};

/*
Point p to p + r represents line segment 1
Point q to q + s represents line segment 2
For a point on line segment 1, use p' = p + t*r (0<=t<=1)
For a point on line segment 2, use q' = q + u*s (0<=u<=1)
Let the two equations be equal to find the intersection point p + t*r = q + u*s
Cross multiply s on both sides
(p + t*r)Xs = (q + u*s)Xs
pXs + t*rXs = qXs
t = (q-p)Xs/(rXs)
In the same way,
u = (p-q)Xr/(sXr) -> u = (q-p)Xr/(rXs)

There are 4 situations as follows:
1. Collinear, sXr==0 & amp; & amp; (q-p)Xr==0, calculate the proportion t0 of the projection of (q-p) on r on the length of r,
Calculate the proportion t1 of the projection of (q + s-p) on r on the length of r, and check whether [t0, t1] intersects with the range [0, 1].
If t0>t1, compare [t1, t0] to see if it intersects with the range [0, 1].
t0 = (q-p)*r/(r*r)
t1 = (q + s-p)*r/(r*r) = t0 + s · r / (r · r)
2. Parallel sXr==0 & amp; & amp; (q-p)Xr!=0
3. 0<=u<=1 & amp; & amp; 0<=t<=1 has an intersection
4. Other u and t are not in the range of 0 to 0, and there is no intersection point.
*/
pair<double, double> intersection(const Point & amp; q, const Point & amp; s, const Point & amp; p, const Point & amp; r) {<!-- -->
// Calculate (q-p)Xr
auto qpr = cross(q - p, r);
auto qps = cross(q - p, s);

auto rXs = cross(r, s);
if (cmp(rXs) == 0)return {<!-- --> -1, -1 }; // Parallel or collinear
//Solve for t, u
// t = (q-p)Xs/(rXs)
auto t = qps / rXs;

// u = (q-p)Xr/(rXs)
auto u = qpr / rXs;

return {<!-- --> u, t };
}


Point oiPs[N];
Delaunay de;
Point lowPoint;
int ind[M];
Point tmpPs[N]; //Storage the intersection point with the boundary
int cakeSize[N];

vector<Point> insectCircle(const Point & amp; A, const Point & amp; dir, double r) {<!-- -->
vector<Point> ans;
Point P = A + dir * dot(dir, -A);
double op = abs(cross(dir, A));
if (cmp(op - r) > 0)return ans;

double Pb = sqrt(r * r - op * op);
if (cmp(Pb) == 0) ans.push_back(P);
else {<!-- -->
ans.push_back(P + dir * Pb);
ans.push_back(P - dir * Pb);
}

return ans;
}

// Sort by polar coordinates
bool sortcmp(int i, int j) {<!-- -->
Point pi = oiPs[i] - lowPoint;
Point pj = oiPs[j] - lowPoint;

// On different sides of the upper and lower halves, the upper half takes priority
if (cmp(pi.y * pj.y) < 0) return pi.y > pj.y;
pi /= pi.dis();
pj /= pj.dis();
// There is an optimization for 1, 0, x
if (cmp(pi.x - 1) == 0 || 0 == cmp(pj.x - 1)) return pi.x > pj.x;

double d = cmp(cross(pi, pj)); // Determine whether the same side rotates counterclockwise
return d > 0;
}

bool oneLine(int n) {<!-- -->
if (n < 3)return true;
for (int i = 2; i < n; + + i) {<!-- -->
if (cmp(cross(oiPs[1] - oiPs[0], oiPs[i] - oiPs[0])) != 0) return false;
}

return true;
}

void solve() {<!-- -->
int n, r;

scanf("%d%d", & amp;r, & amp;n);
for (int i = 0; i < n; + + i) {<!-- -->
scanf("%lf%lf", & amp;oiPs[i].x, & amp;oiPs[i].y);
scanf("%d", cakeSize + i);
oiPs[i].id = i;
}
// Determine whether it is collinear
if (oneLine(n)) {<!-- -->
//Select the largest
int indbig = 0;
for (int i = 1; i < n; + + i) {<!-- -->
if (cakeSize[i] > cakeSize[indbig])indbig = i;
}

Point dir = oiPs[n - 1] - oiPs[0];
dir = {<!-- --> -dir.y, dir.x };
dir /= dir.dis();

auto insectPs = insectCircle(oiPs[indbig], dir, r);
Point seat = insectPs[0];
printf("%.10f %.10f\
", seat.x, seat.y);
return;
}

oiPs[n] = oiPs[0];

de.init(n, oiPs);
auto oiedges = de.getEdge();
vector<vector<int>> link(n, vector<int>());
for (auto oie : oiedges) {<!-- -->
link[oie.first].push_back(oie.second);
}

int maxSize = 0;
Point seat;

for (int i = 0; i < n; + + i) {<!-- -->
//Traverse each boundary point, collect adjacent edges, and sort them counterclockwise.
int len = 0;
for (auto to : link[i]) {<!-- -->
ind[len + + ] = to;
}

lowPoint = oiPs[i];
sort(ind, ind + len, sortcmp);
ind[len] = ind[0]; // Add loop optimization
// Find the intersection point between voronoi boundaries
bool isInsect = false; // Mark whether the vonoroi cell intersects with the circle
for (int i = 0; i < len & amp; & amp; !isInsect; + + i) {<!-- -->
Point mid = (lowPoint + oiPs[ind[i]]) / 2;
Point dir = oiPs[ind[i]] - lowPoint;
dir = {<!-- --> -dir.y, dir.x }; // Rotate 90 degrees

Point mid2 = (lowPoint + oiPs[ind[i + 1]]) / 2;
Point dir2 = oiPs[ind[i + 1]] - lowPoint;
dir2 = {<!-- --> -dir2.y, dir2.x }; // Rotate 90 degrees

// Determine whether it is the boundary of the capital (the angle cannot be greater than 180)
if (cmp(cross(dir, dir2)) <= 0) {<!-- -->
isInsect = true;
break;
}

// Find the intersection point
auto pr = intersection(mid, dir, mid2, dir2);

Point ablePoint = mid2 + dir2 * pr.second;
if (cmp(ablePoint.dis() - r) >= 0)isInsect = true;
}



int k = 0;
// Find the intersection point with the circle
for (int i = 0; i < len & amp; & amp; isInsect; + + i) {<!-- -->
Point mid = (lowPoint + oiPs[ind[i]]) / 2;
Point dir = oiPs[ind[i]] - lowPoint;
dir = {<!-- --> -dir.y, dir.x }; // Rotate 90 degrees

dir /= dir.dis();

auto insectPs = insectCircle(mid, dir, r);
for (auto & amp; c : insectPs) {<!-- -->
tmpPs[k + + ] = c;
//c.out();
}
}

//Exclude invalid intersections
for (int i = 0; i < len; + + i) {<!-- -->
Point mid = (lowPoint + oiPs[ind[i]]) / 2;
Point dir = oiPs[ind[i]] - lowPoint;
dir = {<!-- --> -dir.y, dir.x }; // Rotate 90 degrees
for (int j = 0; j < k; + + j) {<!-- -->
// Determine whether the point is in the same half plane as Lowpoint
if (cmp(cross(lowPoint - mid, dir) * cross(dir, tmpPs[j] - mid)) > 0) {<!-- -->
swap(tmpPs[k - 1], tmpPs[j]);
k--;
j--;
continue;
}
}
}
for (int j = 0; j < k; + + j) {<!-- -->
if (cakeSize[i] > maxSize) {<!-- -->
maxSize = cakeSize[i];
seat = tmpPs[j];
//seat.out();
}
}
}

printf("%.10f %.10f\
", seat.x, seat.y);
//printf("%d\
", maxSize);
//printf("%.10f\
", seat.dis());
}



int main() {<!-- -->
solve();
return 0;

}

/*
10 5
0 0 100
1 0 1
0 1 2
0 -1 3
-1 0 4


10 3
1 -1 100
2 2 200
-2.5 -2.56 1

10 5
0.2 0 100
1 0 1
2 0 2
3 0 1
4 0 4

10 5
9.89 0 100
1 0 1
2 0 2
3 0 1
4 0 4

10 2
10 0 100
1 0 1


*/

I am a coder, and I hope that through my sharing, it will be easier for everyone to learn computer knowledge. Creation is not easy, please help click on the link to the official account.