polygon.cc

   1
   2 // Written and (C) by Francois Fleuret
   3 // Contact <francois.fleuret@idiap.ch> for comments & bug reports
   4
   5 #include <cmath>
   6 #include "polygon.h"
   7
   8 static const scalar_t dl = 20.0;
   9 static const scalar_t repulsion_constant = 0.2;
  10 static const scalar_t dissipation = 0.5;
  11
  12 Polygon::Polygon(scalar_t mass,
  13                  scalar_t red, scalar_t green, scalar_t blue,
  14                  scalar_t *x, scalar_t *y,
  15                  int nv) : _mass(mass),
  16                            _moment_of_inertia(0), _radius(0),
  17                            _relative_x(new scalar_t[nv]), _relative_y(new scalar_t[nv]),
  18                            _total_nb_dots(0),
  19                            _nb_dots(new int[nv]),
  20                            _effecting_edge(0),
  21                            _length(new scalar_t[nv]),
  22                            _triangles(new Triangle[nv-2]),
  23                            _initialized(false), _nailed(false),
  24                            _nb_vertices(nv),
  25                            _x(new scalar_t[nv]), _y(new scalar_t[nv]),
  26                            _red(red), _green(green), _blue(blue) {
  27   _last_center_x = 0;
  28   _last_center_y = 0;
  29   _last_theta = 0;
  30
  31   if(x) for(int i = 0; i < nv; i++) _relative_x[i] = x[i];
  32   if(y) for(int i = 0; i < nv; i++) _relative_y[i] = y[i];
  33 }
  34
  35 Polygon::~Polygon() {
  36   delete[] _relative_x;
  37   delete[] _relative_y;
  38   delete[] _x;
  39   delete[] _y;
  40   delete[] _nb_dots;
  41   delete[] _length;
  42   delete[] _triangles;
  43   delete[] _effecting_edge;
  44 }
  45
  46 Polygon *Polygon::clone() {
  47   return new Polygon(_mass, _red, _green, _blue, _relative_x, _relative_y, _nb_vertices);
  48 }
  49
  50 void Polygon::color_xfig(XFigTracer *tracer) {
  51   tracer->add_color(int(255 * _red), int(255 * _green), int(255 * _blue));
  52 }
  53
  54 void Polygon::print_xfig(XFigTracer *tracer) {
  55   tracer->draw_polygon(int(255 * _red), int(255 * _green), int(255 * _blue),
  56                        _nb_vertices, _x, _y);
  57 }
  58
  59 void Polygon::draw(SimpleWindow *window) {
  60   window->color(_red, _green, _blue);
  61   int x[_nb_vertices], y[_nb_vertices];
  62   for(int n = 0; n < _nb_vertices; n++) {
  63     x[n] = int(_x[n]);
  64     y[n] = int(_y[n]);
  65   }
  66   window->fill_polygon(_nb_vertices, x, y);
  67 }
  68
  69 void Polygon::draw_contours(SimpleWindow *window) {
  70   int x[_nb_vertices], y[_nb_vertices];
  71   for(int n = 0; n < _nb_vertices; n++) {
  72     x[n] = int(_x[n]);
  73     y[n] = int(_y[n]);
  74   }
  75   // window->color(0.0, 0.0, 0.0);
  76   window->color(1.0, 1.0, 1.0);
  77   for(int n = 0; n < _nb_vertices; n++) {
  78     window->draw_line(x[n], y[n], x[(n+1)%_nb_vertices], y[(n+1)%_nb_vertices]);
  79   }
  80 }
  81
  82 #ifdef CAIRO_SUPPORT
  83 void Polygon::draw(cairo_t* context_resource) {
  84   cairo_set_line_width(context_resource, 1.0);
  85   cairo_set_source_rgb (context_resource, _red, _green, _blue);
  86   cairo_move_to(context_resource, _x[0], _y[0]);
  87   for(int n = 0; n < _nb_vertices; n++) {
  88     cairo_line_to(context_resource, _x[n], _y[n]);
  89   }
  90   cairo_close_path(context_resource);
  91   cairo_stroke_preserve(context_resource);
  92   cairo_fill(context_resource);
  93 }
  94
  95 void Polygon::draw_contours(cairo_t* context_resource) {
  96   cairo_set_line_width(context_resource, 1.0);
  97   // cairo_set_source_rgb (context_resource, 0.0, 0.0, 0.0);
  98   cairo_set_source_rgb (context_resource, 1.0, 1.0, 1.0);
  99   cairo_move_to(context_resource, _x[0], _y[0]);
 100   for(int n = 0; n < _nb_vertices; n++) {
 101     cairo_line_to(context_resource, _x[n], _y[n]);
 102   }
 103   cairo_close_path(context_resource);
 104   cairo_stroke(context_resource);
 105 }
 106 #endif
 107
 108 void Polygon::set_vertex(int k, scalar_t x, scalar_t y) {
 109   _relative_x[k] = x;
 110   _relative_y[k] = y;
 111 }
 112
 113 void Polygon::set_position(scalar_t center_x, scalar_t center_y, scalar_t theta) {
 114   _center_x = center_x;
 115   _center_y = center_y;
 116   _theta = theta;
 117 }
 118
 119 void Polygon::set_speed(scalar_t dcenter_x, scalar_t dcenter_y, scalar_t dtheta) {
 120   _dcenter_x = dcenter_x;
 121   _dcenter_y = dcenter_y;
 122   _dtheta = dtheta;
 123 }
 124
 125 bool Polygon::contain(scalar_t x, scalar_t y) {
 126   for(int t = 0; t < _nb_vertices-2; t++) {
 127     scalar_t xa = _x[_triangles[t].a], ya = _y[_triangles[t].a];
 128     scalar_t xb = _x[_triangles[t].b], yb = _y[_triangles[t].b];
 129     scalar_t xc = _x[_triangles[t].c], yc = _y[_triangles[t].c];
 130     if(prod_vect(x - xa, y - ya, xb - xa, yb - ya) <= 0 &&
 131        prod_vect(x - xb, y - yb, xc - xb, yc - yb) <= 0 &&
 132        prod_vect(x - xc, y - yc, xa - xc, ya - yc) <= 0) return true;
 133   }
 134   return false;
 135 }
 136
 137 void Polygon::triangularize(int &nt, int nb, int *index) {
 138   if(nb == 3) {
 139
 140     if(nt >= _nb_vertices-2) {
 141       cerr << "Error type #1 in triangularization." << endl;
 142       exit(1);
 143     }
 144
 145     _triangles[nt].a = index[0];
 146     _triangles[nt].b = index[1];
 147     _triangles[nt].c = index[2];
 148     nt++;
 149
 150   } else {
 151     int best_m = -1, best_n = -1;
 152     scalar_t best_split = -1, det, s = -1, t = -1;
 153
 154     for(int n = 0; n < nb; n++) for(int m = 0; m < n; m++) if(n > m+1 && m+nb > n+1) {
 155       bool no_intersection = true;
 156       for(int k = 0; no_intersection & (k < nb); k++)
 157         if(k != n && k != m && (k+1)%nb != n && (k+1)%nb != m) {
 158           intersection(_relative_x[index[n]], _relative_y[index[n]],
 159                        _relative_x[index[m]], _relative_y[index[m]],
 160                        _relative_x[index[k]], _relative_y[index[k]],
 161                        _relative_x[index[(k+1)%nb]], _relative_y[index[(k+1)%nb]], det, s, t);
 162           no_intersection = det == 0 || s < 0 || s > 1 || t < 0 || t > 1;
 163         }
 164
 165       if(no_intersection) {
 166         scalar_t a1 = 0, a2 = 0;
 167         for(int k = 0; k < nb; k++) if(k >= m && k < n)
 168           a1 += prod_vect(_relative_x[index[k]] - _relative_x[index[m]],
 169                           _relative_y[index[k]] - _relative_y[index[m]],
 170                           _relative_x[index[k+1]] - _relative_x[index[m]],
 171                           _relative_y[index[k+1]] - _relative_y[index[m]]);
 172         else
 173           a2 += prod_vect(_relative_x[index[k]] - _relative_x[index[m]],
 174                           _relative_y[index[k]] - _relative_y[index[m]],
 175                           _relative_x[index[(k+1)%nb]] - _relative_x[index[m]],
 176                           _relative_y[index[(k+1)%nb]] - _relative_y[index[m]]);
 177
 178         if((a1 * a2 > 0 && best_split < 0) || (abs(a1 - a2) < best_split)) {
 179           best_n = n; best_m = m;
 180           best_split = abs(a1 - a2);
 181         }
 182       }
 183     }
 184
 185     if(best_n >= 0 && best_m >= 0) {
 186       int index_neg[nb], index_pos[nb];
 187       int neg = 0, pos = 0;
 188       for(int k = 0; k < nb; k++) {
 189         if(k >= best_m && k <= best_n) index_pos[pos++] = index[k];
 190         if(k <= best_m || k >= best_n) index_neg[neg++] = index[k];
 191       }
 192       if(pos < 3 || neg < 3) {
 193         cerr << "Error type #2 in triangularization." << endl;
 194         exit(1);
 195       }
 196       triangularize(nt, pos, index_pos);
 197       triangularize(nt, neg, index_neg);
 198     } else {
 199       cerr << "Error type #3 in triangularization." << endl;
 200       exit(1);
 201     }
 202   }
 203 }
 204
 205 void Polygon::initialize(int nb_polygons) {
 206   scalar_t a;
 207
 208   _nb_polygons = nb_polygons;
 209
 210   a = _relative_x[_nb_vertices - 1] * _relative_y[0] - _relative_x[0] * _relative_y[_nb_vertices - 1];
 211   for(int n = 0; n < _nb_vertices - 1; n++)
 212     a += _relative_x[n] * _relative_y[n+1] - _relative_x[n+1] * _relative_y[n];
 213   a *= 0.5;
 214
 215   // Reorder the vertices
 216
 217   if(a < 0) {
 218     a = -a;
 219     scalar_t x, y;
 220     for(int n = 0; n < _nb_vertices / 2; n++) {
 221       x = _relative_x[n];
 222       y = _relative_y[n];
 223       _relative_x[n] = _relative_x[_nb_vertices - 1 - n];
 224       _relative_y[n] = _relative_y[_nb_vertices - 1 - n];
 225       _relative_x[_nb_vertices - 1 - n] = x;
 226       _relative_y[_nb_vertices - 1 - n] = y;
 227     }
 228   }
 229
 230   // Compute the center of mass and moment of inertia
 231
 232   scalar_t sx, sy, w;
 233   w = 0;
 234   sx = 0;
 235   sy = 0;
 236   for(int n = 0; n < _nb_vertices; n++) {
 237     int np = (n+1)%_nb_vertices;
 238     w =_relative_x[n] * _relative_y[np] - _relative_x[np] * _relative_y[n];
 239     sx += (_relative_x[n] + _relative_x[np]) * w;
 240     sy += (_relative_y[n] + _relative_y[np]) * w;
 241   }
 242   sx /= 6 * a;
 243   sy /= 6 * a;
 244
 245   _radius = 0;
 246   for(int n = 0; n < _nb_vertices; n++) {
 247     _relative_x[n] -= sx;
 248     _relative_y[n] -= sy;
 249     scalar_t r = sqrt(sq(_relative_x[n]) + sq(_relative_y[n]));
 250     if(r > _radius) _radius = r;
 251   }
 252
 253   scalar_t num = 0, den = 0;
 254   for(int n = 0; n < _nb_vertices; n++) {
 255     int np = (n+1)%_nb_vertices;
 256     den += abs(prod_vect(_relative_x[np], _relative_y[np], _relative_x[n], _relative_y[n]));
 257     num += abs(prod_vect(_relative_x[np], _relative_y[np], _relative_x[n], _relative_y[n])) *
 258       (_relative_x[np] * _relative_x[np] + _relative_y[np] * _relative_y[np] +
 259        _relative_x[np] * _relative_x[n] + _relative_y[np] * _relative_y[n] +
 260        _relative_x[n] * _relative_x[n] + _relative_y[n] * _relative_y[n]);
 261   }
 262
 263   _moment_of_inertia = num / (6 * den);
 264
 265   scalar_t vx = cos(_theta), vy = sin(_theta);
 266
 267   for(int n = 0; n < _nb_vertices; n++) {
 268     _x[n] = _center_x + _relative_x[n] * vx + _relative_y[n] * vy;
 269     _y[n] = _center_y - _relative_x[n] * vy + _relative_y[n] * vx;
 270   }
 271
 272   scalar_t length;
 273   _total_nb_dots = 0;
 274   for(int n = 0; n < _nb_vertices; n++) {
 275     length = sqrt(sq(_relative_x[n] - _relative_x[(n+1)%_nb_vertices]) +
 276                   sq(_relative_y[n] - _relative_y[(n+1)%_nb_vertices]));
 277     _length[n] = length;
 278     _nb_dots[n] = int(length / dl + 1);
 279     _total_nb_dots += _nb_dots[n];
 280   }
 281
 282   delete[] _effecting_edge;
 283   _effecting_edge = new int[_nb_polygons * _total_nb_dots];
 284   for(int p = 0; p < _nb_polygons * _total_nb_dots; p++) _effecting_edge[p] = -1;
 285
 286   int nt = 0;
 287   int index[_nb_vertices];
 288   for(int n = 0; n < _nb_vertices; n++) index[n] = n;
 289   triangularize(nt, _nb_vertices, index);
 290
 291   _initialized = true;
 292 }
 293
 294 bool Polygon::update(scalar_t dt) {
 295   if(!_nailed) {
 296     _center_x += _dcenter_x * dt;
 297     _center_y += _dcenter_y * dt;
 298     _theta += _dtheta * dt;
 299   }
 300
 301   scalar_t d = exp(log(dissipation) * dt);
 302   _dcenter_x *= d;
 303   _dcenter_y *= d;
 304   _dtheta *= d;
 305
 306   scalar_t vx = cos(_theta), vy = sin(_theta);
 307
 308   for(int n = 0; n < _nb_vertices; n++) {
 309     _x[n] = _center_x + _relative_x[n] * vx + _relative_y[n] * vy;
 310     _y[n] = _center_y - _relative_x[n] * vy + _relative_y[n] * vx;
 311   }
 312
 313   if(abs(_center_x - _last_center_x) +
 314      abs(_center_y - _last_center_y) +
 315      abs(_theta - _last_theta) * _radius > 0.1) {
 316     _last_center_x = _center_x;
 317     _last_center_y = _center_y;
 318     _last_theta = _theta;
 319     return true;
 320   } else return false;
 321 }
 322
 323 scalar_t Polygon::relative_x(scalar_t ax, scalar_t ay) {
 324   return (ax - _center_x) * cos(_theta) - (ay - _center_y) * sin(_theta);
 325 }
 326
 327 scalar_t Polygon::relative_y(scalar_t ax, scalar_t ay) {
 328   return (ax - _center_x) * sin(_theta) + (ay - _center_y) * cos(_theta);
 329 }
 330
 331 scalar_t Polygon::absolute_x(scalar_t rx, scalar_t ry) {
 332   return _center_x + rx * cos(_theta) + ry * sin(_theta);
 333 }
 334
 335 scalar_t Polygon::absolute_y(scalar_t rx, scalar_t ry) {
 336   return _center_y - rx * sin(_theta) + ry * cos(_theta);
 337 }
 338
 339 void Polygon::apply_force(scalar_t dt, scalar_t x, scalar_t y, scalar_t fx, scalar_t fy) {
 340   _dcenter_x += fx / _mass * dt;
 341   _dcenter_y += fy / _mass * dt;
 342   _dtheta -= prod_vect(x - _center_x, y - _center_y, fx, fy) / (_mass * _moment_of_inertia) * dt;
 343 }
 344
 345 void Polygon::apply_border_forces(scalar_t dt, scalar_t xmax, scalar_t ymax) {
 346   for(int v = 0; v < _nb_vertices; v++) {
 347     int vp = (v+1)%_nb_vertices;
 348     for(int d = 0; d < _nb_dots[v]; d++) {
 349       scalar_t s = scalar_t(d * dl)/_length[v];
 350       scalar_t x = _x[v] * (1 - s) + _x[vp] * s;
 351       scalar_t y = _y[v] * (1 - s) + _y[vp] * s;
 352       scalar_t vx = 0, vy = 0;
 353       if(x < 0) vx = x;
 354       else if(x > xmax) vx = x - xmax;
 355       if(y < 0) vy = y;
 356       else if(y > ymax) vy = y - ymax;
 357       apply_force(dt, x, y, - dl * vx * repulsion_constant, - dl * vy * repulsion_constant);
 358     }
 359   }
 360 }
 361
 362 void Polygon::apply_collision_forces(scalar_t dt, int n_polygon, Polygon *p) {
 363   scalar_t closest_x[_total_nb_dots], closest_y[_total_nb_dots];
 364   bool inside[_total_nb_dots];
 365   scalar_t distance[_total_nb_dots];
 366   int _new_effecting_edge[_total_nb_dots];
 367
 368   int first_dot = 0;
 369
 370   for(int v = 0; v < _nb_vertices; v++) {
 371     int vp = (v+1)%_nb_vertices;
 372     scalar_t x = _x[v], y = _y[v], xp = _x[vp], yp = _y[vp];
 373
 374     for(int d = 0; d < _nb_dots[v]; d++) {
 375       inside[d] = false;
 376       distance[d] = FLT_MAX;
 377     }
 378
 379     // First, we tag the dots located inside the polygon p
 380
 381     for(int t = 0; t < p->_nb_vertices-2; t++) {
 382       scalar_t min = 0, max = 1;
 383       scalar_t xa = p->_x[p->_triangles[t].a], ya = p->_y[p->_triangles[t].a];
 384       scalar_t xb = p->_x[p->_triangles[t].b], yb = p->_y[p->_triangles[t].b];
 385       scalar_t xc = p->_x[p->_triangles[t].c], yc = p->_y[p->_triangles[t].c];
 386       scalar_t den, num;
 387
 388       const scalar_t eps = 1e-6;
 389
 390       den = prod_vect(xp - x, yp - y, xb - xa, yb - ya);
 391       num = prod_vect(xa - x, ya - y, xb - xa, yb - ya);
 392       if(den > eps) {
 393         if(num / den < max) max = num / den;
 394       } else if(den < -eps) {
 395         if(num / den > min) min = num / den;
 396       } else {
 397         if(num < 0) { min = 1; max = 0; }
 398       }
 399
 400       den = prod_vect(xp - x, yp - y, xc - xb, yc - yb);
 401       num = prod_vect(xb - x, yb - y, xc - xb, yc - yb);
 402       if(den > eps) {
 403         if(num / den < max) max = num / den;
 404       } else if(den < -eps) {
 405         if(num / den > min) min = num / den;
 406       } else {
 407         if(num < 0) { min = 1; max = 0; }
 408       }
 409
 410       den = prod_vect(xp - x, yp - y, xa - xc, ya - yc);
 411       num = prod_vect(xc - x, yc - y, xa - xc, ya - yc);
 412       if(den > eps) {
 413         if(num / den < max) max = num / den;
 414       } else if(den < -eps) {
 415         if(num / den > min) min = num / den;
 416       } else {
 417         if(num < 0) { min = 1; max = 0; }
 418       }
 419
 420       for(int d = 0; d < _nb_dots[v]; d++) {
 421         scalar_t s = scalar_t(d * dl)/_length[v];
 422         if(s >= min && s <= max) inside[d] = true;
 423       }
 424     }
 425
 426     // Then, we compute for each dot what is the closest point on
 427     // the boundary of p
 428
 429     for(int m = 0; m < p->_nb_vertices; m++) {
 430       int mp = (m+1)%p->_nb_vertices;
 431       scalar_t xa = p->_x[m], ya = p->_y[m];
 432       scalar_t xb = p->_x[mp], yb = p->_y[mp];
 433       scalar_t gamma0 = ((x - xa) * (xb - xa) + (y - ya) * (yb - ya)) / sq(p->_length[m]);
 434       scalar_t gamma1 = ((xp - x) * (xb - xa) + (yp - y) * (yb - ya)) / sq(p->_length[m]);
 435       scalar_t delta0 = (prod_vect(xb - xa, yb - ya, x - xa, y - ya)) / p->_length[m];
 436       scalar_t delta1 = (prod_vect(xb - xa, yb - ya, xp - x, yp - y)) / p->_length[m];
 437
 438       for(int d = 0; d < _nb_dots[v]; d++) if(inside[d]) {
 439         int r = _effecting_edge[(first_dot + d) * _nb_polygons + n_polygon];
 440
 441         // If there is already a spring, we look only at the
 442         // vertices next to the current one
 443
 444         if(r < 0 || m == r || m == (r+1)%p->_nb_vertices || (m+1)%p->_nb_vertices == r) {
 445
 446           scalar_t s = scalar_t(d * dl)/_length[v];
 447           scalar_t delta = abs(s * delta1 + delta0);
 448           if(delta < distance[d]) {
 449             scalar_t gamma = s * gamma1 + gamma0;
 450             if(gamma < 0) {
 451               scalar_t l = sqrt(sq(x * (1 - s) + xp * s - xa) + sq(y * (1 - s) + yp * s - ya));
 452               if(l < distance[d]) {
 453                 distance[d] = l;
 454                 closest_x[d] = xa;
 455                 closest_y[d] = ya;
 456                 _new_effecting_edge[first_dot + d] = m;
 457               }
 458             } else if(gamma > 1) {
 459               scalar_t l = sqrt(sq(x * (1 - s) + xp * s - xb) + sq(y * (1 - s) + yp * s - yb));
 460               if(l < distance[d]) {
 461                 distance[d] = l;
 462                 closest_x[d] = xb;
 463                 closest_y[d] = yb;
 464                 _new_effecting_edge[first_dot + d] = m;
 465               }
 466             } else {
 467               distance[d] = delta;
 468               closest_x[d] = xa * (1 - gamma) + xb * gamma;
 469               closest_y[d] = ya * (1 - gamma) + yb * gamma;
 470               _new_effecting_edge[first_dot + d] = m;
 471             }
 472           }
 473         }
 474       } else _new_effecting_edge[first_dot + d] = -1;
 475     }
 476
 477     // Apply forces
 478
 479     for(int d = 0; d < _nb_dots[v]; d++) if(inside[d]) {
 480       scalar_t s = scalar_t(d * dl)/_length[v];
 481       scalar_t x = _x[v] * (1 - s) + _x[vp] * s;
 482       scalar_t y = _y[v] * (1 - s) + _y[vp] * s;
 483       scalar_t vx = x - closest_x[d];
 484       scalar_t vy = y - closest_y[d];
 485
 486       p->apply_force(dt,
 487                      closest_x[d], closest_y[d],
 488                      dl * vx  * repulsion_constant, dl * vy * repulsion_constant);
 489
 490       apply_force(dt,
 491                   x, y,
 492                   - dl * vx * repulsion_constant, - dl * vy * repulsion_constant);
 493     }
 494
 495     first_dot += _nb_dots[v];
 496   }
 497
 498   for(int d = 0; d < _total_nb_dots; d++)
 499     _effecting_edge[d * _nb_polygons + n_polygon] = _new_effecting_edge[d];
 500
 501 }
 502
 503 bool Polygon::collide(Polygon *p) {
 504   for(int n = 0; n < _nb_vertices; n++) {
 505     int np = (n+1)%_nb_vertices;
 506     for(int m = 0; m < p->_nb_vertices; m++) {
 507       int mp = (m+1)%p->_nb_vertices;
 508       scalar_t det, s = -1, t = -1;
 509       intersection(_x[n], _y[n], _x[np], _y[np],
 510                    p->_x[m], p->_y[m], p->_x[mp], p->_y[mp], det, s, t);
 511       if(det != 0 && s>= 0 && s <= 1&& t >= 0 && t <= 1) return true;
 512     }
 513   }
 514
 515   for(int n = 0; n < _nb_vertices; n++) if(p->contain(_x[n], _y[n])) return true;
 516   for(int n = 0; n < p->_nb_vertices; n++) if(contain(p->_x[n], p->_y[n])) return true;
 517
 518   return false;
 519 }