読者です 読者をやめる 読者になる 読者になる

土下座しながら探索中

主に競技プログラミング

AOJ 1301 : Malfatti Circles

AOJ 幾何

問題リンク:Malfatti Circles | Aizu Online Judge

問題概要:

解法:
3つの円の半径はそれぞれ以下のサイトにのっている式で求めることが可能

http://forumgeom.fau.edu/FG2003volume3/FG200308.pdf

2つの円の半径を2分探索をして求める方法もあるみたい

2辺に接している半径rの円の中心点を求める方法は次のとおり

コード:

//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;
}