首先將橢圓的圓心移動到原點處，同時矩形也進行相應的移動。

在次之後，使用仿射變換，選定任意一個軸變換使其變成圓，那麼相應的矩形也同等比例的變換。

變換之後，來一個板子把他們整合起來，結束。

  1 int dcmp(double x) {
  2     if (fabs(x) < eps)return 0;
  3     return x > 0 ? 1 : -1;
  4 }
  5 struct Point {
  6     double x, y;
  7     Point(double _x = 0, double _y = 0) {
  8         x = _x;y = _y;
 
  9     }
 10 };
 11 Point operator + (const Point& a, const Point& b) {
 12     return Point(a.x + b.x, a.y + b.y);
 13 }
 14 Point operator - (const Point& a, const Point& b) {
 15     return Point(a.x - b.x, a.y - b.y);
 16 }
 17 Point operator * (const Point& a, const double& p) {
 
 18     return Point(a.x * p, a.y * p);
 19 }
 20 Point operator / (const Point& a, const double& p) {
 21     return Point(a.x / p, a.y / p);
 22 }
 23 bool operator < (const Point& a, const Point& b) {
 24     return a.x < b.x || (dcmp(a.x - b.x) == 0 && a.y < b.y);
 25
 
 }
 26 bool operator == (const Point& a, const Point& b) {
 27     return dcmp(a.x - b.x) == 0 && dcmp(a.y - b.y) == 0;
 28 }
 29 double Dot(Point  a, Point b) {
 30     return a.x * b.x + a.y * b.y;
 31 }
 32 double Length(Point a) {
 33     return sqrt(Dot(a, a));
 34 }
 35 double Angle(Point a, Point b) {
 36     return acos(Dot(a, b) / Length(a) / Length(b));
 37 }
 38 double angle(Point a) {
 39     return atan2(a.y, a.x);
 40 }
 41 double Cross(Point a, Point b) {
 42     return a.x * b.y - a.y * b.x;
 43 }
 44 Point vecunit(Point a) {
 45     return a / Length(a);
 46 }
 47 Point Normal(Point a) {
 48     return Point(-a.y, a.x) / Length(a);
 49 }
 50 Point Rotate(Point a, double rad) {
 51     return Point(a.x * cos(rad) - a.y * sin(rad), a.x * sin(rad) + a.y * cos(rad));
 52 }
 53 double Area2(Point a, Point b, Point c) {
 54     return Length(Cross(b - a, c - a));
 55 }
 56 bool OnSegment(Point p, Point a1, Point a2) {
 57     return dcmp(Cross(a1 - p, a2 - p)) == 0 && dcmp(Dot(a1 - p, a2 - p)) <= 0;
 58 }
 59 struct Line {
 60     Point p, v;
 61     double ang;
 62     Line() {};
 63     Line(Point p, Point v) :p(p), v(v) {
 64         ang = atan2(v.y, v.x);
 65     }
 66     bool operator < (const Line& L) const {
 67         return ang < L.ang;
 68     }
 69     Point point(double d) {
 70         return p + (v * d);
 71     }
 72 };
 73 bool OnLeft(const Line& L, const Point& p) {
 74     return Cross(L.v, p - L.p) >= 0;
 75 }
 76 Point GetLineIntersection(Point p, Point v, Point q, Point w) {
 77     Point u = p - q;
 78     double t = Cross(w, u) / Cross(v, w);
 79     return p + v * t;
 80 }
 81 Point GetLineIntersection(Line a, Line b) {
 82     return GetLineIntersection(a.p, a.v, b.p, b.v);
 83 }
 84 double PolyArea(vector<Point> p) {
 85     int n = p.size();
 86     double ans = 0;
 87     for (int i = 1;i < n - 1;i++)
 88         ans += Cross(p[i] - p[0], p[i + 1] - p[0]);
 89     return fabs(ans) / 2;
 90 }
 91 struct Circle
 92 {
 93     Point c;
 94     double r;
 95     Circle() {}
 96     Circle(Point c, double r) :c(c), r(r) {}
 97     Point point(double a) //根據圓心角求點座標      
 98     {
 99         return Point(c.x + cos(a) * r, c.y + sin(a) * r);
100     }
101 };
102 
103 bool InCircle(Point x, Circle c) {
104     return dcmp(c.r - Length(c.c - x)) >= 0;
105 }
106 bool OnCircle(Point x, Circle c) {
107     return dcmp(c.r - Length(c.c - x)) == 0;
108 }
109 int getSegCircleIntersection(Line L, Circle C, Point* sol) {
110     Point nor = Normal(L.v);
111     Line p1 = Line(C.c, nor);
112     Point ip = GetLineIntersection(p1, L);
113     double dis = Length(ip - C.c);
114     if (dcmp(dis - C.r) > 0)return 0;
115     Point dxy = vecunit(L.v) * sqrt(C.r * C.r - dis * dis);
116     int ret = 0;
117     sol[ret] = ip + dxy;
118     if (OnSegment(sol[ret], L.p, L.point(1)))ret++;
119     sol[ret] = ip - dxy;
120     if (OnSegment(sol[ret], L.p, L.point(1)))ret++;
121     return ret;
122 }
123 double SegCircleArea(Circle C, Point a, Point b) {
124     double a1 = angle(a - C.c);
125     double a2 = angle(b - C.c);
126     double da = fabs(a1 - a2);
127     if (da > pi)da = pi * 2 - da;
128     return dcmp(Cross(b - C.c, a - C.c)) * da * C.r * C.r / 2.0;
129 }
130 double PolyCircleArea(Circle C, Point* p, int n) {
131     double ret = 0;
132     Point sol[2];
133     p[n] = p[0];
134     for (int i = 0;i < n;i++) {
135         double t1, t2;
136         int cnt = getSegCircleIntersection(Line(p[i], p[i + 1] - p[i]), C, sol);  //判斷線段與圓有幾個交點，  
137     //  cout<<"cnt="<<cnt<<" "<<p[i].x<<" "<<p[i].y<<" "<<p[i+1].x<<" "<<p[i+1].y<<endl;  
138     //  cout<<"C: "<<C.c.x<<" "<<C.c.y<<" "<<C.r<<endl;  
139     //  cout<<"sol ";for(int j=0;j<cnt;j++)cout<<sol[j].x<<" "<<sol[j].y<<" ";cout<<endl; 
140         //cout << cnt << endl;
141         if (cnt == 0) {  //0個交點，判斷線段在多邊形內部還是外部。  
142             if (!InCircle(p[i], C) || !InCircle(p[i + 1], C))ret += SegCircleArea(C, p[i], p[i + 1]); //外部直接計算圓弧面積   
143             else ret += Cross(p[i + 1] - C.c, p[i] - C.c) / 2;  //內部計算三角形面積。
144         }
145         if (cnt == 1) {
146             if (InCircle(p[i], C) && (!InCircle(p[i + 1], C) || OnCircle(p[i + 1], C)))ret += Cross(sol[0] - C.c, p[i] - C.c) / 2, ret += SegCircleArea(C, sol[0], p[i + 1]);//,cout<<"jj-1"<<endl;  
147             else ret += SegCircleArea(C, p[i], sol[0]), ret += Cross(p[i + 1] - C.c, sol[0] - C.c) / 2;//,cout<<"jj-2"<<endl;  
148         }
149         if (cnt == 2) {  //兩個交點  
150             if ((p[i] < p[i + 1]) ^ (sol[0] < sol[1]))swap(sol[0], sol[1]);
151             ret += SegCircleArea(C, p[i], sol[0]);
152             ret += Cross(sol[1] - C.c, sol[0] - C.c) / 2;
153             ret += SegCircleArea(C, sol[1], p[i + 1]);
154         }
155         // cout<<ret<<endl;      
156     }
157     return fabs(ret);
158 }