[코드모음]-기하,

/* 기하 Basic structure */

struct Point{

        double x, y;

 

        Point(){}

        Point(double _x, double _y){ x = _x; y = _y; }

 

        Point operator + (const Point& rhs) const{ return Point(x+rhs.x, y+rhs.y); }

        Point operator - (const Point& rhs) const{ return Point(x - rhs.x, y - rhs.y); }

        bool operator < (const Point& rhs){   return (x != rhs.x) ? (x < rhs.x) : (y < rhs.y); }

        Point operator * (double a) { return Point(a * x, a * y); }

 

        double len(){ return sqrt((double)x*x + y*y); }

        Point normalize(){

               double len = sqrt((double)x*x + y*y);

               return Point(x/len, y/len);

        }

        Point projection(Point& b){

               Point n = b.normalize();

               return n * dot(n,*this);

        }

};

double dot(Point& a, Point& b){ return a.x * b.x + a.y * b.y; }

double cross(Point &a, Point& b){ return a.x * b.y - a.y * b.x; }

 

/* 라인 교차여부 */

nt cross(P a, P b) {      return a.x * b.y - a.y*b.x; }

int ccw(P a, P b, P c) {return cross(b - a, c - a);} 

bool isInter(P a, P b, P c, P d) {

         int ab = ccw(a, b, c) * ccw(a, b, d);

         int cd = ccw(c, d, a) * ccw(c, d, b);

         if (ab == 0 && cd == 0) {

                  if (b < a) swap(a, b);

                  if (d < c) swap(c, d);

                  // [a..b]인 경우만 true. 즉, !( d<- a~b -> c )

                  return !(d<a || b<c);

         }

         return ab <= 0 && cd <= 0;

}

 

/* 내외부 판별 */

bool isInside(P q, vector<P>& p) {

         int cnt = 0;

         int x, x1, x2, y, y1, y2;

 

         x = q.x; y = q.y;

         for (int i = 0; i < p.size(); ++i) {

                  int j = (i + 1) % p.size();

                  x1 = p[i].x; y1 = p[i].y;

                  x2 = p[j].x; y2 = p[j].y;

                 

                  //핵심. [low..high) 임. 낮은쪽 점과 그 위쪽 라인을 지날 때만. 평행의 경우는 해당 안됨

                  if ((y1 > y) != (y2 > y)) {

                           //두 점을 이루는 직선상에서, y점에 대한 x값

                           double atX = (x2 - x1) * (y - y1) / (y2 - y1) + x1;

                           if (x < atX) ++cnt;

                  }

         }

         return (cnt % 2) > 0; //홀수이면 inside

}

 

/* convex hull */

ll ccw(Point p1, Point p2, Point p3) {

        return (p2 - p1).cross(p3 - p1);

}

 

double dist(Point a, Point b) {

        return (b.x - a.x) * (b.x - a.x) + (b.y - a.y)*(b.y - a.y);

}

 

bool comp(Point a, Point b) {

        ll cr = ccw(mini, a, b);

        if (cr > 0 || (cr == 0 && dist(mini, a) < dist(mini, b))) return true;

        return false;

}

 

int graham(){

        //1. find most below y. if same y, choose min x -> O(n)

        mini = a[1];

        for (int i = 2; i <= n; i++) {

               if (mini.y > a[i].y || mini.y == a[i].y && mini.x > a[i].x) mini = a[i];

        }

 

        //2. sort       

        sort(a + 1, a + n + 1, comp);

 

        //3. a[n]->a[0]에서부터 시작해서,

        int st_top = 0;

        st[0] = a[n]; //주의!!!

        st[++st_top] = a[1];

 

        for (int i = 2; i <= n; ++i) { //n-1이 아닌 n까지라는 것에 유의

               while (ccw(st[st_top - 1], st[st_top], a[i]) <= 0) st_top--; //주의!!

               st[++st_top] = a[i];

        }

 

        return st_top;

}

 

/* 회전하는 캘리퍼스 - 컨벡스헐의 가장 긴 거리 구하기*/

- 볼록다각형에 버니아캘리퍼스를 대고, 돌려가면서 가장 거리가 긴 두 점을 찾는 원리 이용

1. x축을 기준으로, 다각형의 가장 왼쪽에 있는 점과, 가장 오른쪽에 있는 점을 찾는다. 그 점의 인덱스: left, right

2. 다각형의 각 선분에 대한 단위벡터를 구해 둔다. 

        vector<V> u(N+1);

        for (int i = 1; i <= N; ++i) u[i] = (p[(i + 1) % N] - p[i]).normalize();

3. a=left에서 b=right까지 단위벡터 선분과 캘리퍼스와의 각을 구하면서 작은 각도쪽으로 회전.

   즉, cosA= dot(diaA, u[a])와 cosB를 통해 A와 B의 각도를 비교하고, 작은쪽 선분으로 dia 치환

4. 위의 루틴을 돌면서, p[a]~p[b] 사잇거리가 가장 긴 것 찾아낸다.

double diameter() {

         //1. 가장 왼쪽, 오른쪽 점 찾는다. idx:left, right

         int leftMax = p[1].x, rightMax = p[1].x;

         int left = 1, right = 1;

         for (int i = 2; i <= N; ++i) {

                  if (p[i].x < leftMax) {

                           leftMax = p[i].x; left = i;

                  }

                  if (p[i].x > rightMax) {

                           rightMax = p[i].x; right = i;

                  }

         }

 

         //2. 도형의 선분을 단위벡터로 구해둔다. u[i]: i ~ (i+1) 사이 단위벡터

         vector<V> u(N+1);

         for (int i = 1; i <= N; ++i)        u[i] = (p[(i + 1) % N] - p[i]).normalize();

        

 

         //3. a=left에서 b=right까지 단위벡터 선분과 calipers와의 각을 구하면서 작은 각도쪽으로 회전.

         //   즉, cosA= dot(diaA, u[a])와 cosB를 통해 A와 B의 각도를 비교하고, 작은쪽 선분으로 dia 치환

         double maxLen = (p[right] - p[left]).length();

         V diaA = V(0, 1);

         V diaB = V(0, -1);

         int a = left, b = right;

         while (a!=right || b!=left) {

                  double cosA = dot(diaA, u[a]);

                  double cosB = dot(diaB, u[b]);

 

                  if (cosA > cosB) { // ==> A < B

                           diaA = u[a];

                           diaB = V(diaA.x * -1, -diaA.y * -1);

                           a = (a + 1) % N;

                  }

                  else {

                           diaB = u[b];

                           diaA = V(diaB.x * -1, diaB.y * -1);

                           b = (b + 1) % N;

                  }

                  double len = (p[b] - p[a]).length();

                  if (len > maxLen) maxLen = len;

         }

 

         return maxLen;

}

 

-끝-