Advertisement
istinishat

GeometryTemplate

Aug 23rd, 2017
264
0
Never
Not a member of Pastebin yet? Sign Up, it unlocks many cool features!
  1. #include <iostream>
  2. #include <vector>
  3. #include <cmath>
  4. #include <cassert>
  5.  
  6. using namespace std;
  7.  
  8. double INF = 1e100;
  9. double EPS = 1e-12;
  10.  
  11. struct PT {
  12.   double x, y;
  13.   PT() {}
  14.   PT(double x, double y) : x(x), y(y) {}
  15.   PT(const PT &p) : x(p.x), y(p.y)    {}
  16.   PT operator + (const PT &p)  const { return PT(x+p.x, y+p.y); }
  17.   PT operator - (const PT &p)  const { return PT(x-p.x, y-p.y); }
  18.   PT operator * (double c)     const { return PT(x*c,   y*c  ); }
  19.   PT operator / (double c)     const { return PT(x/c,   y/c  ); }
  20. };
  21.  
  22. double dot(PT p, PT q)     { return p.x*q.x+p.y*q.y; }
  23. double dist2(PT p, PT q)   { return dot(p-q,p-q); }
  24. double cross(PT p, PT q)   { return p.x*q.y-p.y*q.x; }
  25. ostream &operator<<(ostream &os, const PT &p) {
  26.   os << "(" << p.x << "," << p.y << ")";
  27. }
  28.  
  29. istream &operator>>(istream &is,PT &p) {
  30.   is>>p.x>>p.y;
  31. }
  32.  
  33. // rotate a point CCW or CW around the origin
  34. PT RotateCCW90(PT p)   { return PT(-p.y,p.x); }
  35. PT RotateCW90(PT p)    { return PT(p.y,-p.x); }
  36. PT RotateCCW(PT p, double t) {
  37.   return PT(p.x*cos(t)-p.y*sin(t), p.x*sin(t)+p.y*cos(t));
  38. }
  39.  
  40. // project point c onto line through a and b
  41. // assuming a != b
  42. PT ProjectPointLine(PT a, PT b, PT c) {
  43.   return a + (b-a)*dot(c-a, b-a)/dot(b-a, b-a);
  44. }
  45.  
  46. // project point c onto line segment through a and b
  47. PT ProjectPointSegment(PT a, PT b, PT c) {
  48.   double r = dot(b-a,b-a);
  49.   if (fabs(r) < EPS) return a;
  50.   r = dot(c-a, b-a)/r;
  51.   if (r < 0) return a;
  52.   if (r > 1) return b;
  53.   return a + (b-a)*r;
  54. }
  55.  
  56. // compute distance from c to segment between a and b
  57. double DistancePointSegment(PT a, PT b, PT c) {
  58.   return sqrt(dist2(c, ProjectPointSegment(a, b, c)));
  59. }
  60.  
  61. // compute distance between point (x,y,z) and plane ax+by+cz=d
  62. double DistancePointPlane(double x, double y, double z,
  63.                           double a, double b, double c, double d)
  64. {
  65.   return fabs(a*x+b*y+c*z-d)/sqrt(a*a+b*b+c*c);
  66. }
  67.  
  68. // determine if lines from a to b and c to d are parallel or collinear
  69. bool LinesParallel(PT a, PT b, PT c, PT d) {
  70.   return fabs(cross(b-a, c-d)) < EPS;
  71. }
  72.  
  73. bool LinesCollinear(PT a, PT b, PT c, PT d) {
  74.   return LinesParallel(a, b, c, d)
  75.       && fabs(cross(a-b, a-c)) < EPS
  76.       && fabs(cross(c-d, c-a)) < EPS;
  77. }
  78.  
  79. // determine if line segment from a to b intersects with
  80. // line segment from c to d
  81. bool SegmentsIntersect(PT a, PT b, PT c, PT d) {
  82.   if (LinesCollinear(a, b, c, d)) {
  83.     if (dist2(a, c) < EPS || dist2(a, d) < EPS ||
  84.       dist2(b, c) < EPS || dist2(b, d) < EPS) return true;
  85.     if (dot(c-a, c-b) > 0 && dot(d-a, d-b) > 0 && dot(c-b, d-b) > 0)
  86.       return false;
  87.     return true;
  88.   }
  89.   if (cross(d-a, b-a) * cross(c-a, b-a) > 0) return false;
  90.   if (cross(a-c, d-c) * cross(b-c, d-c) > 0) return false;
  91.   return true;
  92. }
  93.  
  94. // compute intersection of line passing through a and b
  95. // with line passing through c and d, assuming that unique
  96. // intersection exists; for segment intersection, check if
  97. // segments intersect first
  98. PT ComputeLineIntersection(PT a, PT b, PT c, PT d) {
  99.   b=b-a; d=c-d; c=c-a;
  100.   assert(dot(b, b) > EPS && dot(d, d) > EPS);
  101.   return a + b*cross(c, d)/cross(b, d);
  102. }
  103.  
  104. // compute center of circle given three points
  105. PT ComputeCircleCenter(PT a, PT b, PT c) {
  106.   b=(a+b)/2;
  107.   c=(a+c)/2;
  108.   return ComputeLineIntersection(b, b+RotateCW90(a-b), c, c+RotateCW90(a-c));
  109. }
  110.  
  111. // determine if point is in a possibly non-convex polygon (by William
  112. // Randolph Franklin); returns 1 for strictly interior points, 0 for
  113. // strictly exterior points, and 0 or 1 for the remaining points.
  114. // Note that it is possible to convert this into an *exact* test using
  115. // integer arithmetic by taking care of the division appropriately
  116. // (making sure to deal with signs properly) and then by writing exact
  117. // tests for checking point on polygon boundary
  118. bool PointInPolygon(const vector<PT> &p, PT q) {
  119.   bool c = 0;
  120.   for (int i = 0; i < p.size(); i++){
  121.     int j = (i+1)%p.size();
  122.     if ((p[i].y <= q.y && q.y < p[j].y ||
  123.       p[j].y <= q.y && q.y < p[i].y) &&
  124.       q.x < p[i].x + (p[j].x - p[i].x) * (q.y - p[i].y) / (p[j].y - p[i].y))
  125.       c = !c;
  126.   }
  127.   return c;
  128. }
  129.  
  130. // determine if point is on the boundary of a polygon
  131. bool PointOnPolygon(const vector<PT> &p, PT q) {
  132.   for (int i = 0; i < p.size(); i++)
  133.     if (dist2(ProjectPointSegment(p[i], p[(i+1)%p.size()], q), q) < EPS)
  134.       return true;
  135.     return false;
  136. }
  137.  
  138. // compute intersection of line through points a and b with
  139. // circle centered at c with radius r > 0
  140. vector<PT> CircleLineIntersection(PT a, PT b, PT c, double r) {
  141.   vector<PT> ret;
  142.   b = b-a;
  143.   a = a-c;
  144.   double A = dot(b, b);
  145.   double B = dot(a, b);
  146.   double C = dot(a, a) - r*r;
  147.   double D = B*B - A*C;
  148.   if (D < -EPS) return ret;
  149.   ret.push_back(c+a+b*(-B+sqrt(D+EPS))/A);
  150.   if (D > EPS)
  151.     ret.push_back(c+a+b*(-B-sqrt(D))/A);
  152.   return ret;
  153. }
  154.  
  155. // compute intersection of circle centered at a with radius r
  156. // with circle centered at b with radius R
  157. vector<PT> CircleCircleIntersection(PT a, PT b, double r, double R) {
  158.   vector<PT> ret;
  159.   double d = sqrt(dist2(a, b));
  160.   if (d > r+R || d+min(r, R) < max(r, R)) return ret;
  161.   double x = (d*d-R*R+r*r)/(2*d);
  162.   double y = sqrt(r*r-x*x);
  163.   PT v = (b-a)/d;
  164.   ret.push_back(a+v*x + RotateCCW90(v)*y);
  165.   if (y > 0)
  166.     ret.push_back(a+v*x - RotateCCW90(v)*y);
  167.   return ret;
  168. }
  169.  
  170. // This code computes the area or centroid of a (possibly nonconvex)
  171. // polygon, assuming that the coordinates are listed in a clockwise or
  172. // counterclockwise fashion.  Note that the centroid is often known as
  173. // the "center of gravity" or "center of mass".
  174. double ComputeSignedArea(const vector<PT> &p) {
  175.   double area = 0;
  176.   for(int i = 0; i < p.size(); i++) {
  177.     int j = (i+1) % p.size();
  178.     area += p[i].x*p[j].y - p[j].x*p[i].y;
  179.   }
  180.   return area / 2.0;
  181. }
  182.  
  183. double ComputeArea(const vector<PT> &p) {
  184.   return fabs(ComputeSignedArea(p));
  185. }
  186.  
  187. PT ComputeCentroid(const vector<PT> &p) {
  188.   PT c(0,0);
  189.   double scale = 6.0 * ComputeSignedArea(p);
  190.   for (int i = 0; i < p.size(); i++){
  191.     int j = (i+1) % p.size();
  192.     c = c + (p[i]+p[j])*(p[i].x*p[j].y - p[j].x*p[i].y);
  193.   }
  194.   return c / scale;
  195. }
  196.  
  197. // tests whether or not a given polygon (in CW or CCW order) is simple
  198. bool IsSimple(const vector<PT> &p) {
  199.   for (int i = 0; i < p.size(); i++) {
  200.     for (int k = i+1; k < p.size(); k++) {
  201.       int j = (i+1) % p.size();
  202.       int l = (k+1) % p.size();
  203.       if (i == l || j == k) continue;
  204.       if (SegmentsIntersect(p[i], p[j], p[k], p[l]))
  205.         return false;
  206.     }
  207.   }
  208.   return true;
  209. }
  210.  
  211.  
  212. double areaOfTriangle(double ax,double ay,double bx,double by,double cx,double cy)
  213. {
  214.     double t=(ax*(by-cy)+bx*(cy-ay)+cx*(ay-by))/2;
  215.     if(t<0)t*=-1.00;
  216.     return t;
  217. }
  218.  
  219. pair<double,double> ratio_point(double x1,double y1,double x2,double y2,double m,double n)
  220. {
  221.     //AB=2BC
  222.     //A=x1,y1  C=x2,y2
  223.     //m=2,n=1 , find B
  224.     double x=(n*x1+m*x2)/(m+n);
  225.     double y=(n*y1+m*y2)/(m+n);
  226.     return mp(x,y);
  227. }
  228.  
  229. int main() {
  230.  
  231.  
  232. }
Advertisement
Add Comment
Please, Sign In to add comment
Advertisement