CodeArchive CircleBoundaryTraversal

Circle Boundary Traversal

Don't mess up with comments.

#include <iostream>
#include <cmath>
#include <complex>
using namespace std;
typedef long double ld;
typedef complex < ld > point;
const ld eps = 1e-9;
const ld PI = 2*acos((ld)0);
int vis[128], done[128][128], ind[128][128], sz[128], cnt, n;
point c[128];
ld out[128][128], in[128][128], r[128], x, y;

void dfs(int now){
  if(vis[now]) return;
  vis[now] = 1;
  for(int i = 0; i < n; i++)
    if(!vis[i] && abs(c[now] - c[i]) < r[now] + r[i]) dfs(i);
}

int main(){

  int t; cin >> t;
  while(t-- && cin >> n){
    for(int i = 0; i < n && cin >> x >> y >> r[i]; i++) c[i] = point(x, y);
    memset(vis, 0, sizeof(vis)), memset(sz, 0, sizeof(sz)), memset(done, 0, sizeof(done));

    for(int i = 0; i < n; i++)
      for(int j = 0; j < n; j++){
        ld d = abs(c[i] - c[j]);
        if(i != j && abs(r[i] - r[j]) < d && d < r[i] + r[j]){
          point v = c[j] - c[i];
          ld cosa = (r[i]*r[i] + d*d - r[j]*r[j])/(2*r[i]*d);
          ld sina = sqrt(1 - cosa*cosa);
          v *= point(cosa, -sina); // v is direction to intersection 
          v /= abs(v), v *= r[i];
          point inter = v + c[i]; // now inter is the intersection point;


          bool good = true;
          for(int k = 0; k < n; k++) good &= (abs(inter - c[k]) + eps > r[k]);

          if(good){  // only if not inside a circle
            out[i][sz[i]] = atan2(v.imag(), v.real());
            if(out[i][sz[i]] < 0) out[i][sz[i]] += 2*PI;
            ind[i][sz[i]] = j;

            v = c[i] - c[j];
            cosa = (r[j]*r[j] + d*d - r[i]*r[i])/(2*r[j]*d);
            sina = sqrt(1 - cosa*cosa);

            v *= point(cosa, sina); // the other turn
            in[i][sz[i]++] = atan2(v.imag(), v.real());
            if(in[i][sz[i] - 1] < 0) in[i][sz[i] - 1] += 2*PI;
          } 
        }
      }

    cnt = 1;
    for(int i = 0; i < n; i++){
      if(sz[i] && !vis[i]) cnt--, dfs(i);
      for(int j = 0; j < sz[i]; j++)
        if(!done[i][j]){
          cnt++;
          ld angle = in[i][j] - eps;
          int now = ind[i][j];
          while(true){
            ld best = 1e99;
            int p = -1;
            for(int k = 0; k < sz[now]; k++){
              ld dif = out[now][k] - angle; if(dif < 0) dif += PI*2;
              if(dif < best) best = dif, p = k;
            }
            if(done[now][p]) break;
            done[now][p] = 1;
            angle = in[now][p] - eps;
            now = ind[now][p];
          }
        }
    }
    cout << cnt - 1 << endl;
  }

  return 0;
}

-- Main.MatthewChan - 29 Mar 2006