计算几何模板

来源:互联网 发布:淘宝俄区key 编辑:程序博客网 时间:2024/06/07 13:15


多圆面积交

typedef long long LL;  typedef unsigned long long ULL;  typedef vector <int> VI;  const int INF = 0x3f3f3f3f;  const double eps = 1e-10;  const int MOD = 100000007;  const int MAXN = 1000010;  const double PI = acos(-1.0);  #define sqr(x) ((x)*(x))  const int N = 1010;  double area[N];  int n;    int dcmp(double x)  {      if (x < -eps) return -1;      else return x > eps;  }    struct cp  {      double x, y, r, angle;      int d;      cp() {}      cp(double xx, double yy, double ang = 0, int t = 0)      {          x = xx;          y = yy;          angle = ang;          d = t;      }      void get()      {          scanf("%lf%lf%lf", &x, &y, &r);          d = 1;      }  } cir[N], tp[N * 2];    double dis(cp a, cp b)  {      return sqrt(sqr(a.x - b.x) + sqr(a.y - b.y));  }    double cross(cp p0, cp p1, cp p2)  {      return (p1.x - p0.x) * (p2.y - p0.y) - (p1.y - p0.y) * (p2.x - p0.x);  }    int CirCrossCir(cp p1, double r1, cp p2, double r2, cp &cp1, cp &cp2)  {      double mx = p2.x - p1.x, sx = p2.x + p1.x, mx2 = mx * mx;      double my = p2.y - p1.y, sy = p2.y + p1.y, my2 = my * my;      double sq = mx2 + my2, d = -(sq - sqr(r1 - r2)) * (sq - sqr(r1 + r2));      if (d + eps < 0) return 0;      if (d < eps) d = 0;      else d = sqrt(d);      double x = mx * ((r1 + r2) * (r1 - r2) + mx * sx) + sx * my2;      double y = my * ((r1 + r2) * (r1 - r2) + my * sy) + sy * mx2;      double dx = mx * d, dy = my * d;      sq *= 2;      cp1.x = (x - dy) / sq;      cp1.y = (y + dx) / sq;      cp2.x = (x + dy) / sq;      cp2.y = (y - dx) / sq;      if (d > eps) return 2;      else return 1;  }    bool circmp(const cp& u, const cp& v)  {      return dcmp(u.r - v.r) < 0;  }    bool cmp(const cp& u, const cp& v)  {      if (dcmp(u.angle - v.angle)) return u.angle < v.angle;      return u.d > v.d;  }    double calc(cp cir, cp cp1, cp cp2)  {      double ans = (cp2.angle - cp1.angle) * sqr(cir.r)                   - cross(cir, cp1, cp2) + cross(cp(0, 0), cp1, cp2);      return ans / 2;  }    void CirUnion(cp cir[], int n)  {      cp cp1, cp2;      sort(cir, cir + n, circmp);      for (int i = 0; i < n; ++i)          for (int j = i + 1; j < n; ++j)              if (dcmp(dis(cir[i], cir[j]) + cir[i].r - cir[j].r) <= 0)                  cir[i].d++;      for (int i = 0; i < n; ++i)      {          int tn = 0, cnt = 0;          for (int j = 0; j < n; ++j)          {              if (i == j) continue;              if (CirCrossCir(cir[i], cir[i].r, cir[j], cir[j].r,                              cp2, cp1) < 2) continue;              cp1.angle = atan2(cp1.y - cir[i].y, cp1.x - cir[i].x);              cp2.angle = atan2(cp2.y - cir[i].y, cp2.x - cir[i].x);              cp1.d = 1;              tp[tn++] = cp1;              cp2.d = -1;              tp[tn++] = cp2;              if (dcmp(cp1.angle - cp2.angle) > 0) cnt++;          }          tp[tn++] = cp(cir[i].x - cir[i].r, cir[i].y, PI, -cnt);          tp[tn++] = cp(cir[i].x - cir[i].r, cir[i].y, -PI, cnt);          sort(tp, tp + tn, cmp);          int p, s = cir[i].d + tp[0].d;          for (int j = 1; j < tn; ++j)          {              p = s;              s += tp[j].d;              area[p] += calc(cir[i], tp[j - 1], tp[j]);          }      }  }    void solve()  {      scanf("%d", &n);      for (int i = 0; i < n; ++i)          cir[i].get();      memset(area, 0, sizeof(area));      CirUnion(cir, n);      //去掉重复计算的      for (int i = 1; i <= n; ++i)      {          area[i] -= area[i + 1];      }      //area[i]为重叠了i次的面积      //tot 为总面积      double tot = 0;      for(int i=1; i<=n; i++) tot += area[i];      printf("%f\n", tot);  }    int main()  {      //freopen("input.txt", "r", stdin);      return 0;  }  

多圆面积并

#include<iostream>#include<cstdio>#include<cmath>#include<algorithm>#include<cstring>#define ld double#define eps 1e-13using namespace std;int n,top,st,ed;ld xl[1001],xr[1001];ld ans;bool del[1001];struct data{ld x,y,r;}t[1001],sk[1001];struct line{ld l,r;}p[1001];ld dis(data a,data b){return sqrt((a.x-b.x)*(a.x-b.x)+(a.y-b.y)*(a.y-b.y));}bool cmp1(data a,data b){return a.r<b.r;}bool cmp2(data a,data b){return a.x-a.r<b.x-b.r;}bool cmp3(line a,line b){return a.l<b.l;}void ini(){    scanf("%d",&n);    for(int i=1;i<=n;i++)       {scanf("%lf%lf%lf",&t[i].x,&t[i].y,&t[i].r);}   sort(t+1,t+n+1,cmp1);   for(int i=1;i<=n;i++)      for(int j=i+1;j<=n;j++)         if(dis(t[i],t[j])<=t[j].r-t[i].r)               {del[i]=1;break;}   for(int i=1;i<=n;i++)if(!del[i])sk[++top]=t[i];n=top;   sort(sk+1,sk+n+1,cmp2);}ld getf(ld x){    int sz=0,i,j;ld r,len=0,dis;    for(i=st;i<=ed;i++)   {       if(x<=xl[i]||x>=xr[i])continue;       dis=sqrt(sk[i].r-(x-sk[i].x)*(x-sk[i].x));       p[++sz].l=sk[i].y-dis;p[sz].r=sk[i].y+dis;   }    sort(p+1,p+sz+1,cmp3);    for(i=1;i<=sz;i++)    {        r=p[i].r;        for(j=i+1;j<=sz;j++)        {            if(p[j].l>r)break;            if(r<p[j].r)r=p[j].r;        }        len+=r-p[i].l;i=j-1;    }    return len;}ld cal(ld l,ld fl,ld fmid,ld fr){return (fl+fmid*4+fr)*l/6;}ld simpson(ld l,ld mid,ld r,ld fl,ld fmid,ld fr,ld s){    ld m1=(l+mid)/2,m2=(r+mid)/2;    ld f1=getf(m1),f2=getf(m2);    ld g1=cal(mid-l,fl,f1,fmid),g2=cal(r-mid,fmid,f2,fr);    if(fabs(g1+g2-s)<eps)return g1+g2;   return simpson(l,m1,mid,fl,f1,fmid,g1)+simpson(mid,m2,r,fmid,f2,fr,g2);}void work(){    int i,j;ld l,r,mid,fl,fr,fmid;    for(i=1;i<=n;i++){xl[i]=sk[i].x-sk[i].r;xr[i]=sk[i].x+sk[i].r;sk[i].r*=sk[i].r;}    for(i=1;i<=n;i++)    {        l=xl[i];r=xr[i];        for(j=i+1;j<=n;j++)        {            if(xl[j]>r)break;         if(xr[j]>r)r=xr[j];        }        st=i;ed=j-1;i=j-1;      mid=(l+r)/2;      fl=getf(l);fr=getf(r);fmid=getf(mid);      ans+=simpson(l,mid,r,fl,fmid,fr,cal(r-l,fl,fmid,fr));    }}int main(){    ini();    work();    printf("%.3lf",ans);    return 0;}

给你一堆点,找锐角三角形。

TWO POINTER 思想。

统计出所有锐角和直=钝角的数目。

做法是这样的:对每个点对所有点极角排序,然后TWO POINTER计算每一个锐角(一个边上有好几个点也会被统计好几次),直角钝角。然后ans=(锐角个数-直角钝角个数*2)/3;因为每一个角度可能也只可能出现在一个三角形中。

但是一上来没弄懂,所以队友想了一个思路也能过,先统计直角和钝角的个数。然后C(n,3)统计所有情况,不能构成三角形的是三点共线的情况,所以n^2logn统计三点共线(TWO POINTER)。然后再减去直角钝角个数即可。两种方法均可过。

#include <cstdio>#include <cstring>#include <vector>#include <cmath>#include <stack>#include <cstdlib>#include <queue>#include <map>#include <iostream>#include <algorithm>using namespace std;typedef long long LL;struct Point{    int x,y;    Point(int _x = 0,int _y = 0)    {        x=_x;        y=_y;    }    Point operator -(const Point &a) const    {        return Point(x-a.x,y-a.y);    }}p[2010];int n;LL xmul(const Point &a,const Point &b){    return a.x * 1ll * b.y - a.y * 1ll * b.x;}LL dot(const Point &a,const Point &b){    return a.x * 1ll * b.x + a.y * 1ll * b.y;}bool cmp(const Point &a,const Point &b){    if (a.y * 1ll * b.y <= 0)    {        if (a.y > 0 || b.y > 0) return a.y < b.y;        if (a.y == 0 && b.y == 0) return a.x < b.x;    }    return xmul(a,b) > 0;}vector<Point>vec;int main(){    while (scanf ("%d",&n)!=EOF)    {        LL zhidun=0,rui=0,ans=0;        for (int i=0;i<n;++i)        {            scanf ("%d%d",&p[i].x,&p[i].y);        }        for (int cor=0;cor<n;++cor)        {            vec.clear();            for (int i=0;i<n;++i)            {                if (i!=cor)                vec.emplace_back(p[i]-p[cor]);            }            sort(vec.begin(),vec.end(),cmp);            vec.insert(vec.end(),vec.begin(),vec.end());            int j=0,k=0,r=0;            for (int i=0;i<n-1;++i)            {                while (j<i+n-1&&xmul(vec[i],vec[j])==0&&dot(vec[i],vec[j])>0) ++j;                k=max(k,j);                while (k<i+n-1&&xmul(vec[i],vec[k])>0&&dot(vec[i],vec[k])>0) ++k;                r=max(r,k);                while (r<i+n-1&&xmul(vec[i],vec[r])>0) ++r;                rui+=k-j;                zhidun+=r-k;            }        }        ans=(rui-2*zhidun)/3;        printf ("%I64d\n",ans);    }    return 0;}





原创粉丝点击