AOJ 1301 : Malfatti Circles
問題リンク:Malfatti Circles | Aizu Online Judge
問題概要:
略
解法:
3つの円の半径はそれぞれ以下のサイトにのっている式で求めることが可能
http://forumgeom.fau.edu/FG2003volume3/FG200308.pdf
2つの円の半径を2分探索をして求める方法もあるみたい
コード:
//rad は角度をラジアンで持たせること Point rotate(Point a,double rad) { return Point(cos(rad)*a.x - sin(rad)*a.y,sin(rad)*a.x + cos(rad)*a.y); } // 度をラジアンに変換 double toRad(double agl) { return agl*M_PI/180.0; } double getArea(vector<Point>& vec) { double sum = 0; for(int i=0;i<vec.size();i++) sum += cross(vec[i],vec[(i+1)%vec.size()]); return fabs(sum)/2.0; } double getDist(Point a,Point b) { return sqrt(norm(a-b)); } int ccw(Point p0,Point p1,Point p2) { Point a = p1-p0; Point b = p2-p0; if(cross(a,b) > EPS)return COUNTER_CLOCKWISE; if(cross(a,b) < -EPS)return CLOCKWISE; if(dot(a,b) < -EPS)return ONLINE_BACK; if(norm(a) < norm(b))return ONLINE_FRONT; return ON_SEGMENT; } struct P { }; bool isIntersect(Point p1,Point p2,Point p3,Point p4) { return (ccw(p1,p2,p3) * ccw(p1,p2,p4) <= 0 && ccw(p3,p4,p1) * ccw(p3,p4,p2) <= 0 ); } bool isIntersect(Segment s1,Segment s2) { return isIntersect(s1.p1,s1.p2,s2.p1,s2.p2); } Point getCrossPoint(Point a1,Point a2,Point b1,Point b2) { assert(isIntersect(a1,a2,b1,b2));//交差していない これを外してもある程度(またはしっかり?)はうごいてくれる Point base = b2 - b1; double d1 = abs(cross(base,a1-b1)); double d2 = abs(cross(base,a2-b1)); double t = d1/(d1+d2); return a1+(a2-a1)*t; } // Line Line 上と同じ assertなし Point getCrossPointLines( Line s1, Line s2){ Point a = s1.p2 - s1.p1; Point base = s2.p2 -s2.p1; return s1.p1 + a * (cross(base, s2.p1 - s1.p1)/cross(base, a)); } //三角形ABCの内接円半径を求める double getInradius(Point A,Point B,Point C) { vector<Point> vec(3); vec[0] = A, vec[1] = B, vec[2] = C; double S = getArea(vec); double a = getDist(B,C); double b = getDist(A,C); double c = getDist(A,B); return 2.0*S/(a+b+c); } //segをrだけ平行にスライド Segment slideSeg(Segment seg,double r) { Point p = (seg.p2-seg.p1)*r/abs(seg.p2-seg.p1); p = rotate(p,toRad(90.0)); seg = Segment(seg.p1+p,seg.p2+p); if(fabs(seg.p1.x) < EPS)seg.p1.x = 0; if(fabs(seg.p1.y) < EPS)seg.p1.y = 0; if(fabs(seg.p2.x) < EPS)seg.p2.x = 0; if(fabs(seg.p2.y) < EPS)seg.p2.y = 0; return seg; } //三角形ABCの内心を求める Point getIncenter(Point A,Point B,Point C) { double inradius = getInradius(A,B,C); Segment seg1 = Segment(A,B); Segment seg2 = Segment(B,C); seg1 = slideSeg(seg1,inradius); seg2 = slideSeg(seg2,inradius); return getCrossPointLines(seg1,seg2); } vector<Point> getMalfattieCircles(Point A,Point B,Point C) { Point I = getIncenter(A,B,C); double r = getInradius(A,B,C); double a = getDist(B,C); double b = getDist(C,A); double c = getDist(A,B); double s = (a+b+c)/2.0; double IA = getDist(I,A); double IB = getDist(I,B); double IC = getDist(I,C); /* r1,r2,r3はそれぞれ三角形内の円の半径 チェック済み */ double r1 = ( r/(2.0*(s-a)) )*( s - r - (IB + IC - IA)); double r2 = ( r/(2.0*(s-b)) )*( s - r - (IC + IA - IB)); double r3 = ( r/(2.0*(s-c)) )*( s - r - (IA + IB - IC)); cout << setiosflags(ios::fixed) << setprecision(6) << r1 << " " << r2 << " " << r3 << endl; /* 以下は未チェック(ちょっとだけチェック済み) 三角形内の3つの円の中心点 */ vector<Point> vec(3); Segment seg1 = Segment(C,A); Segment seg2 = Segment(A,B); seg1 = slideSeg(seg1,r1); seg2 = slideSeg(seg2,r1); vec[0] = getCrossPointLines(seg1,seg2); //cout << vec[0].x << "," << vec[0].y << endl; seg1 = Segment(A,B); seg2 = Segment(B,C); seg1 = slideSeg(seg1,r2); seg2 = slideSeg(seg2,r2); vec[1] = getCrossPointLines(seg1,seg2); //cout << vec[1].x << "," << vec[1].y << endl; seg1 = Segment(B,C); seg2 = Segment(C,A); seg1 = slideSeg(seg1,r3); seg2 = slideSeg(seg2,r3); vec[2] = getCrossPointLines(seg1,seg2); //cout << vec[2].x << "," << vec[2].y << endl; return vec; } int main() { while(true) { Point A,B,C; cin >> A.x >> A.y >> B.x >> B.y >> C.x >> C.y; if(A.x+A.y+B.x+B.y+C.x+C.y == 0)break; getMalfattieCircles(A,B,C); } return 0; }