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-1520
The main idea of the title
There is a city whose shape is circular.
There are many chemical plants in the city.
Now I want to place explosives with the same explosion damage radius in each chemical plant. The explosion can affect the entire city. Ask what is the minimum radius.
Thinking &Analysis
First find the edges of the voronoi area where each chemical plant is located.
Determine whether the intersection point of the sides is within the circle.
The intersection point between the edge of the voronoi region and the circle is used to determine whether it is within the edge area of the voronoi region.
The intersection of the diameter and the circle where the chemical plant is located may also be the farthest point.
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
There is only 1 chemical plant, take the furthest distance from the edge of the circle.
Triangulation is not possible when there are only 2 chemical plants or when all plants are on one line.
At this time, sort the coordinates, and take the mid-perpendicular line as the dividing line between each two closest points. It is also necessary to consider whether the intersection of the diameter and the circle is the farthest point.
- The distance from a point p in a circle to the circle is monotonic on the same side of the diameter
Special case: when p is the center of the circle, the distance from p to the circle is the radius.
Normal situation:
In the above figure, p is a point within the great circle, points A and B are the two ends of the diameter where p is located, Ap
Algorithm process
1 point, 2 points, special treatment for collinearity
generally:
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 Find the intersection of the diameter and the circle and determine whether it is within the voronoi area
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;n, & amp;r);
int a, b;
for (int i = 0; i < n; + + i) {<!-- -->
//scanf("%lf%lf", & amp;oiPs[i].x, & amp;oiPs[i].y);
scanf("%d%d", & amp;a, & amp;b);
oiPs[i].x = a;
oiPs[i].y = b;
oiPs[i].id = i;
}
if (n == 1) {<!-- -->
printf("%.10f\\
", r + oiPs[0].dis());
return;
}
double maxRadius = 0;
// Determine whether it is collinear
if (oneLine(n)) {<!-- -->
// Arrange a line and separate each two points with a vertical bisector.
sort(oiPs, oiPs + n);
for (int i = 0; i < n; + + i) {<!-- -->
vector<pair<Point, Point>> midLines;
if (i) {<!-- -->
Point mid = (oiPs[i] + oiPs[i-1]) / 2;
Point dir = oiPs[i] - oiPs[i - 1];
dir = {<!-- --> -dir.y, dir.x }; // Rotate 90 degrees
dir /= dir.dis();
midLines.push_back({<!-- --> mid, dir });
}
if (i + 1 < n) {<!-- -->
Point mid = (oiPs[i] + oiPs[i + 1]) / 2;
Point dir = oiPs[i] - oiPs[i + 1];
dir = {<!-- --> -dir.y, dir.x }; // Rotate 90 degrees
dir /= dir.dis();
midLines.push_back({<!-- --> mid, dir });
}
int k = 0;
// Find the intersection point of diameter and circle
if (oiPs[i].dis() > EPS) {<!-- -->
tmpPs[k + + ] = oiPs[i] / oiPs[i].dis() * r;
tmpPs[k] = -tmpPs[k - 1];
k++;
}
for (auto & amp; midLine : midLines) {<!-- -->
// Find the intersection point with the circle
auto cirs=insectCircle(midLine.first, midLine.second, r);
for (auto & amp; c : cirs) {<!-- -->
maxRadius = max(maxRadius, (c - oiPs[i]).dis());
}
//Exclude diameter intersections that are not within the region
for (int j = 0; j < k; + + j) {<!-- -->
// Determine whether the point is in the same half plane as Lowpoint
if (cmp(cross(oiPs[i] - midLine.first, midLine.second) * cross(midLine.second, tmpPs[j] - midLine.first)) > 0) {<!-- -->
swap(tmpPs[k - 1], tmpPs[j]);
k--;
j--;
continue;
}
}
}
// Get the maximum value of the diameter intersection
for (int j = 0; j < k; + + j) {<!-- -->
maxRadius = max(maxRadius, (tmpPs[j] - oiPs[i]).dis());
}
}
printf("%.10f\\
", maxRadius);
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);
}
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
int k = 0;
// Find the intersection point between voronoi boundaries
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
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) {<!-- -->
// Find the intersection point
auto pr = intersection(mid, dir, mid2, dir2);
Point ablePoint = mid2 + dir2 * pr.second;
if (cmp(ablePoint.dis() - r) <= 0) maxRadius = max(maxRadius, (ablePoint - lowPoint).dis());
}
// Find the intersection point of diameter and circle
if (lowPoint.dis() > EPS) {<!-- -->
tmpPs[k + + ] = lowPoint / lowPoint.dis() * r;
tmpPs[k] = -tmpPs[k - 1];
k++;
}
dir /= dir.dis();
// Find the intersection point with the circle
auto cirs = insectCircle(mid, dir, r);
for (auto & amp; c : cirs) {<!-- -->
tmpPs[k + + ] = c;
}
}
//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) {<!-- -->
maxRadius = max(maxRadius, (tmpPs[j] - lowPoint).dis());
}
}
printf("%.10f\\
", maxRadius);
//printf("%d\\
", maxSize);
//printf("%.10f\\
", seat.dis());
}
int main() {<!-- -->
solve();
return 0;
}
/*
5 10
0 0
1 0
0 1
0 -1
-1 0
5 10
-1 0
0 0
-1 1
-1 -1
-2 0
5 10
0 0
10 0
0 10
0 -10
-10 0
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
2 5
2 0
-2 0
2 5
twenty two
-twenty two
6 10
1 0
-1 0
2 0
-2 0
3 0
-3 0
6 10
1 -1
-1 -1
twenty one
-twenty one
3-1
-3 -1
6 10
-4 -1
-1 -1
-5 -1
-twenty one
-6 -1
-3 -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.
