1. 程式人生 > >P1742 最小圓覆蓋(計算幾何)

P1742 最小圓覆蓋(計算幾何)

line int ase 化簡 當前 std def 圓心 同時

體驗過\(O(n^3)\)\(10^5\)嗎?快來體驗一波當\(wys\)的快感吧\(QAQ\)

前置芝士1:二元一次方程組求解


\[\begin{cases}a1 * x + b1*y=c1\\a2 * x + b2*y=c2\end{cases}\]
(其中\(a1,a2,b1,b2,c1,c2\)為已知量)

\(②\)式得:
\[x=\frac{c2-b2*y}{a2}\]

帶入\(①\)式並化簡得:
\[y=\frac{c1-\frac{a1*c2}{a2}}{b1-\frac{a1*b2}{a2}}\]

分子分母同時乘以\(a2\)得:
\[y=\frac{a2*c1-a1*c2}{a2*b1-a1*b2}\]

同理可得(把\(a,b\)互換即可):
\[x=\frac{b2*c1-b1*c2}{b2*a1-b1*a2}\]

前置芝士2:三點定圓

給出三個點,求出圓心&半徑

\[\begin{cases}x1^2-2x1*x0+x0^2+y1^2-2y1*y0+y0^2=r^2\\x2^2-2x2*x0+x0^2+y2^2-2y2*y0+y0^2=r^2\\x3^2-2x3*x0+x0^2+y3^2-2y3*y0+y0^2=r^2\end{cases}\]

\(②-①\)\(③-①\),並化簡得:

\[\begin{cases}2*(x2-x1)x+2*(y2-y1)y=x2^2-x1^2+y2^2-y1^2\\2*(x3-x1)x+2*(y3-y1)y=x3^2-x1^2+y3^2-y1^2\end{cases}\]

我們將三點定圓的柿子對應二元一次方程組中,可知:

\[a1=x2-x1,\quad a2=x3-x1\]
\[b1=y2-y1,\quad b2=y3-y1\]
\[c1=\frac{x2^2-x1^2+y2^2-y1^2}{2},\quad c2=\frac{x3^2-x1^2+y3^2-y1^2}{2}\]

然後就可以根據三個點求出圓心和半徑了

正文

跟據前置芝士,我們知道對於任意三個不共線的點,我們可以求出三點定的圓,所以一個明顯的想法就是枚舉三個點

我們先枚舉第一個點,有兩種情況

①:當前點在當前外面,即\(dis(\)圓心,該點\()>r\)那麽我們不管這個點

②:不是情況①的情況,那麽我們就需要重新構造這個圓來包含所有的點了

怎麽構造呢?我們重新枚舉兩外兩個已經遍歷過的點,組成三個點。同理,若重新構造的圓包括了三個點,那麽就不管,若有任意一個在圓外,那麽我們根據前置芝士重新確定圓心和半徑即可

PS:本題出題人過於duliu,故意構造數據卡掉了上述解法,所以我們需要一個神奇的東西:隨\((da)\)\((luan)\)\((shu)\)\((ju)\)法,來防止掉精度

代碼如下:

#include<bits/stdc++.h>
using namespace std;
#define il inline
#define re register
#define D double
il int read()
{
    re int x = 0, f = 1; re char c = getchar();
    while(c < '0' || c > '9') { if(c == '-') f = -1; c = getchar();}
    while(c >= '0' && c <= '9') x = x * 10 + c - 48, c = getchar();
    return x * f;
}
#define rep(i, s, t) for(re int i = s; i <= t; ++ i)
#define eps 1e-12
#define maxn 100005
#define ff(x) (x) * (x)
int n, m;
D r;
struct node
{
    D x, y;
}o, e[maxn];
il D dis(node a, node b){return sqrt(ff(a.x - b.x) + ff(a.y - b.y));}
il void get(node a, node b, node c)
{
    D a1 = b.x - a.x, a2 = c.x - a.x, b1 = b.y - a.y, b2 = c.y - a.y;
    D c1 = (ff(b.x) - ff(a.x) + ff(b.y) - ff(a.y));
    D c2 = (ff(c.x) - ff(a.x) + ff(c.y) - ff(a.y));
    o = (node){(b2 * c1 - b1 * c2) / (b2 * a1 * 2 - b1 * a2 * 2), 
               (a2 * c1 - a1 * c2) / (a2 * b1 * 2 - a1 * b2 * 2)};
    r = dis(a, o);
}
il void work()
{
    o = e[1], r = 0;
    rep(i, 2, n)
    {
        if(dis(o, e[i]) > r + eps)
        {
            o = e[i], r = 0;
            rep(j, 1, i - 1)
            {
                if(dis(o, e[j]) > r + eps)
                {
                    o.x = (e[i].x + e[j].x) / 2, o.y = (e[i].y + e[j].y) / 2;
                    r = dis(o, e[j]);
                    rep(k, 1, j - 1) if(dis(o, e[k]) > r + eps) get(e[i], e[j], e[k]);
                }
            }
        }
    //  printf("%.10lf\n%.10lf %.10lf\n", r, o.x, o.y);
    }
}
int main()
    n = read();
    rep(i, 1, n) scanf("%lf%lf", &e[i].x, &e[i].y);
    random_shuffle(e + 1, e + n + 1);
    work();
    printf("%.10lf\n%.10lf %.10lf", r, o.x, o.y);
    return 0;
}

P1742 最小圓覆蓋(計算幾何)