lib/geometry.c

   1 /* Dia -- an diagram creation/manipulation program
   2  * Copyright (C) 1998 Alexander Larsson
   3  *
   4  * Matrix code from GIMP:app/transform_tool.c
   5  * Copyright (C) 1995 Spencer Kimball and Peter Mattis
   6  *
   7  * This program is free software; you can redistribute it and/or modify
   8  * it under the terms of the GNU General Public License as published by
   9  * the Free Software Foundation; either version 2 of the License, or
  10  * (at your option) any later version.
  11  *
  12  * This program is distributed in the hope that it will be useful,
  13  * but WITHOUT ANY WARRANTY; without even the implied warranty of
  14  * MERCHANTABILITY or FITNESS FOR A PARTICULAR PURPOSE.  See the
  15  * GNU General Public License for more details.
  16  *
  17  * You should have received a copy of the GNU General Public License
  18  * along with this program; if not, write to the Free Software
  19  * Foundation, Inc., 59 Temple Place - Suite 330, Boston, MA 02111-1307, USA.
  20  */
  21
  22 /* include glib.h first, so we don't pick up its inline functions as well */
  23 #include <config.h>
  24
  25 #include <glib.h>
  26 /* include normal versions of the inlined functions here ... */
  27 #undef G_INLINE_FUNC
  28 #define G_INLINE_FUNC extern
  29 #define G_CAN_INLINE 1
  30 #include "geometry.h"
  31
  32 void
  33 rectangle_union(Rectangle *r1, const Rectangle *r2)
  34 {
  35   r1->top = MIN( r1->top, r2->top );
  36   r1->bottom = MAX( r1->bottom, r2->bottom );
  37   r1->left = MIN( r1->left, r2->left );
  38   r1->right = MAX( r1->right, r2->right );
  39 }
  40
  41 void
  42 int_rectangle_union(IntRectangle *r1, const IntRectangle *r2)
  43 {
  44   r1->top = MIN( r1->top, r2->top );
  45   r1->bottom = MAX( r1->bottom, r2->bottom );
  46   r1->left = MIN( r1->left, r2->left );
  47   r1->right = MAX( r1->right, r2->right );
  48 }
  49
  50 void
  51 rectangle_intersection(Rectangle *r1, const Rectangle *r2)
  52 {
  53   r1->top = MAX( r1->top, r2->top );
  54   r1->bottom = MIN( r1->bottom, r2->bottom );
  55   r1->left = MAX( r1->left, r2->left );
  56   r1->right = MIN( r1->right, r2->right );
  57
  58   /* Is rectangle empty? */
  59   if ( (r1->top >= r1->bottom) ||
  60        (r1->left >= r1->right) ) {
  61     r1->top = r1->bottom = r1->left = r1->right = 0.0;
  62   }
  63 }
  64
  65 int
  66 rectangle_intersects(const Rectangle *r1, const Rectangle *r2)
  67 {
  68   if ( (r1->right < r2->left) ||
  69        (r1->left > r2->right) ||
  70        (r1->top > r2->bottom) ||
  71        (r1->bottom < r2->top) )
  72     return FALSE;
  73
  74   return TRUE;
  75 }
  76
  77 int
  78 point_in_rectangle(const Rectangle* r, const Point *p)
  79 {
  80   if ( (p->x < r->left) ||
  81        (p->x > r->right) ||
  82        (p->y > r->bottom) ||
  83        (p->y < r->top) )
  84     return FALSE;
  85
  86   return TRUE;
  87 }
  88
  89 int
  90 rectangle_in_rectangle(const Rectangle* outer, const Rectangle *inner)
  91 {
  92   if ( (inner->left < outer->left) ||
  93        (inner->right > outer->right) ||
  94        (inner->top < outer->top) ||
  95        (inner->bottom > outer->bottom))
  96     return FALSE;
  97
  98   return TRUE;
  99 }
 100
 101 void
 102 rectangle_add_point(Rectangle *r, const Point *p)
 103 {
 104   if (p->x < r->left)
 105     r->left = p->x;
 106   else if (p->x > r->right)
 107     r->right = p->x;
 108
 109   if (p->y < r->top)
 110     r->top = p->y;
 111   else if (p->y > r->bottom)
 112     r->bottom = p->y;
 113 }
 114
 115
 116 /*
 117  * This function estimates the distance from a point to a rectangle.
 118  * If the point is in the rectangle, 0.0 is returned. Otherwise the
 119  * distance in a manhattan metric from the point to the nearest point
 120  * on the rectangle is returned.
 121  */
 122 real
 123 distance_rectangle_point(const Rectangle *rect, const Point *point)
 124 {
 125   real dx = 0.0;
 126   real dy = 0.0;
 127
 128   if (point->x < rect->left ) {
 129     dx = rect->left - point->x;
 130   } else if (point->x > rect->right ) {
 131     dx = point->x - rect->right;
 132   }
 133
 134   if (point->y < rect->top ) {
 135     dy = rect->top - point->y;
 136   } else if (point->y > rect->bottom ) {
 137     dy = point->y - rect->bottom;
 138   }
 139
 140   return dx + dy;
 141 }
 142
 143 /*
 144  * This function estimates the distance from a point to a line segment
 145  * specified by two endpoints.
 146  * If the point is on the line segment, 0.0 is returned. Otherwise the
 147  * distance in the R^2 metric from the point to the nearest point
 148  * on the line segment is returned. Does one sqrt per call.
 149  * Philosophical bug: line_width is ignored iff point is beyond
 150  * end of line segment.
 151  */
 152 real
 153 distance_line_point(const Point *line_start, const Point *line_end,
 154                     real line_width, const Point *point)
 155 {
 156   Point v1, v2;
 157   real v1_lensq;
 158   real projlen;
 159   real perp_dist;
 160
 161   v1 = *line_end;
 162   point_sub(&v1,line_start);
 163
 164   v2 = *point;
 165   point_sub(&v2, line_start);
 166
 167   v1_lensq = point_dot(&v1,&v1);
 168
 169   if (v1_lensq<0.000001) {
 170     return sqrt(point_dot(&v2,&v2));
 171   }
 172
 173   projlen = point_dot(&v1,&v2) / v1_lensq;
 174
 175   if (projlen<0.0) {
 176     return sqrt(point_dot(&v2,&v2));
 177   }
 178
 179   if (projlen>1.0) {
 180     Point v3 = *point;
 181     point_sub(&v3,line_end);
 182     return sqrt(point_dot(&v3,&v3));
 183   }
 184
 185   point_scale(&v1, projlen);
 186   point_sub(&v1, &v2);
 187
 188   perp_dist = sqrt(point_dot(&v1,&v1));
 189
 190   perp_dist -= line_width / 2.0;
 191   if (perp_dist < 0.0) {
 192     perp_dist = 0.0;
 193   }
 194
 195   return perp_dist;
 196 }
 197
 198 /* returns 1 iff the line crosses the ray from (-Inf, rayend.y) to rayend */
 199 static guint
 200 line_crosses_ray(const Point *line_start,
 201                  const Point *line_end, const Point *rayend)
 202 {
 203   coord xpos;
 204
 205   /* swap end points if necessary */
 206   if (line_start->y > line_end->y) {
 207     const Point *tmp;
 208
 209     tmp = line_start;
 210     line_start = line_end;
 211     line_end = tmp;
 212   }
 213   /* if y coords of line do not include rayend.y */
 214   if (line_start->y > rayend->y || line_end->y < rayend->y)
 215     return 0;
 216   /* Avoid division by zero */
 217   if (line_end->y - line_start->y < 0.00000000001) {
 218     return line_end->y - rayend->y < 0.00000000001;
 219   }
 220   xpos = line_start->x + (rayend->y - line_start->y) *
 221     (line_end->x - line_start->x) / (line_end->y - line_start->y);
 222   return xpos <= rayend->x;
 223 }
 224
 225 real
 226 distance_polygon_point(const Point *poly, guint npoints, real line_width,
 227                        const Point *point)
 228 {
 229   guint i, last = npoints - 1;
 230   real line_dist = G_MAXFLOAT;
 231   guint crossings = 0;
 232
 233   /* calculate ray crossings and line distances */
 234   for (i = 0; i < npoints; i++) {
 235     real dist;
 236
 237     crossings += line_crosses_ray(&poly[last], &poly[i], point);
 238     dist = distance_line_point(&poly[last], &poly[i], line_width, point);
 239     line_dist = MIN(line_dist, dist);
 240     last = i;
 241   }
 242   /* If there is an odd number of ray crossings, we are inside the polygon.
 243    * Otherwise, return the minium distance from a line segment */
 244   if (crossings % 2 == 1)
 245     return 0.0;
 246   else
 247     return line_dist;
 248 }
 249
 250 /* number of segments to use in bezier curve approximation */
 251 #define NBEZ_SEGS 10
 252
 253 /* if cross is not NULL, it will be incremented for each ray crossing */
 254 static real
 255 bez_point_distance_and_ray_crosses(const Point *b1,
 256                                    const Point *b2, const Point *b3,
 257                                    const Point *b4,
 258                                    real line_width, const Point *point,
 259                                    guint *cross)
 260 {
 261   static gboolean calculated_coeff = FALSE;
 262   static real coeff[NBEZ_SEGS+1][4];
 263   int i;
 264   real line_dist = G_MAXFLOAT;
 265   Point prev, pt;
 266
 267   if (!calculated_coeff) {
 268     for (i = 0; i <= NBEZ_SEGS; i++) {
 269       real t1 = ((real)i)/NBEZ_SEGS, t2 = t1*t1, t3 = t1*t2;
 270       real it1 = 1-t1, it2 = it1*it1, it3 = it1*it2;
 271
 272       coeff[i][0] = it3;
 273       coeff[i][1] = 3 * t1 * it2;
 274       coeff[i][2] = 3 * t2 * it1;
 275       coeff[i][3] = t3;
 276     }
 277   }
 278   calculated_coeff = TRUE;
 279
 280   prev.x = coeff[0][0] * b1->x + coeff[0][1] * b2->x +
 281            coeff[0][2] * b3->x + coeff[0][3] * b4->x;
 282   prev.y = coeff[0][0] * b1->y + coeff[0][1] * b2->y +
 283            coeff[0][2] * b3->y + coeff[0][3] * b4->y;
 284   for (i = 1; i <= NBEZ_SEGS; i++) {
 285     real dist;
 286
 287     pt.x = coeff[i][0] * b1->x + coeff[i][1] * b2->x +
 288            coeff[i][2] * b3->x + coeff[i][3] * b4->x;
 289     pt.y = coeff[i][0] * b1->y + coeff[i][1] * b2->y +
 290            coeff[i][2] * b3->y + coeff[i][3] * b4->y;
 291
 292     dist = distance_line_point(&prev, &pt, line_width, point);
 293     line_dist = MIN(line_dist, dist);
 294     if (cross)
 295       *cross += line_crosses_ray(&prev, &pt, point);
 296
 297     prev = pt;
 298   }
 299   return line_dist;
 300 }
 301
 302 real
 303 distance_bez_seg_point(const Point *b1, const Point *b2,
 304                        const Point *b3, const Point *b4,
 305                        real line_width, const Point *point)
 306 {
 307   return bez_point_distance_and_ray_crosses(b1, b2, b3, b4,
 308                                             line_width, point, NULL);
 309 }
 310
 311 real
 312 distance_bez_line_point(const BezPoint *b, guint npoints,
 313                         real line_width, const Point *point)
 314 {
 315   Point last;
 316   guint i;
 317   real line_dist = G_MAXFLOAT;
 318
 319   g_return_val_if_fail(b[0].type == BEZ_MOVE_TO, -1);
 320
 321   last = b[0].p1;
 322
 323   for (i = 1; i < npoints; i++) {
 324     real dist;
 325
 326     switch (b[i].type) {
 327     case BEZ_MOVE_TO:
 328       /*g_assert_not_reached();*/
 329       g_warning("BEZ_MOVE_TO found half way through a bezier line");
 330       break;
 331     case BEZ_LINE_TO:
 332       dist = distance_line_point(&last, &b[i].p1, line_width, point);
 333       line_dist = MIN(line_dist, dist);
 334       last = b[i].p1;
 335       break;
 336     case BEZ_CURVE_TO:
 337       dist = bez_point_distance_and_ray_crosses(&last, &b[i].p1, &b[i].p2,
 338                                                 &b[i].p3, line_width, point,
 339                                                 NULL);
 340       line_dist = MIN(line_dist, dist);
 341       last = b[i].p3;
 342       break;
 343     }
 344   }
 345   return line_dist;
 346 }
 347
 348 real
 349 distance_bez_shape_point(const BezPoint *b, guint npoints,
 350                          real line_width, const Point *point)
 351 {
 352   Point last;
 353   guint i;
 354   real line_dist = G_MAXFLOAT;
 355   guint crossings = 0;
 356
 357   g_return_val_if_fail(b[0].type == BEZ_MOVE_TO, -1);
 358
 359   last = b[0].p1;
 360
 361   for (i = 1; i < npoints; i++) {
 362     real dist;
 363
 364     switch (b[i].type) {
 365     case BEZ_MOVE_TO:
 366       /*g_assert_not_reached();*/
 367       g_warning("BEZ_MOVE_TO found half way through a bezier shape");
 368       break;
 369     case BEZ_LINE_TO:
 370       dist = distance_line_point(&last, &b[i].p1, line_width, point);
 371       crossings += line_crosses_ray(&last, &b[i].p1, point);
 372       line_dist = MIN(line_dist, dist);
 373       last = b[i].p1;
 374       break;
 375     case BEZ_CURVE_TO:
 376       dist = bez_point_distance_and_ray_crosses(&last, &b[i].p1, &b[i].p2,
 377                                                 &b[i].p3, line_width, point,
 378                                                 &crossings);
 379       line_dist = MIN(line_dist, dist);
 380       last = b[i].p3;
 381       break;
 382     }
 383   }
 384   /* If there is an odd number of ray crossings, we are inside the polygon.
 385    * Otherwise, return the minium distance from a line segment */
 386   if (crossings % 2 == 1)
 387     return 0.0;
 388   else
 389     return line_dist;
 390 }
 391
 392 real
 393 distance_ellipse_point(const Point *centre, real width, real height,
 394                        real line_width, const Point *point)
 395 {
 396   real w2 = width*width, h2 = height*height;
 397   real scale, rad, dist;
 398   Point pt;
 399
 400   /* find the radius of the ellipse in the appropriate direction */
 401
 402   /* find the point of intersection between line (x=cx+(px-cx)t; y=cy+(py-cy)t)
 403    * and ellipse ((x-cx)^2)/(w/2)^2 + ((y-cy)^2)/(h/2)^2 = 1 */
 404   /* radius along ray is sqrt((px-cx)^2 * t^2 + (py-cy)^2 * t^2) */
 405
 406   pt = *point;
 407   point_sub(&pt, centre);
 408   pt.x *= pt.x;
 409   pt.y *= pt.y;  /* pt = (point - centre).^2 */
 410
 411   scale = w2 * h2 / (4*h2*pt.x + 4*w2*pt.y);
 412   rad = sqrt((pt.x + pt.y)*scale) + line_width/2;
 413
 414   dist = sqrt(pt.x + pt.y);
 415
 416   if (dist <= rad)
 417     return 0.0;
 418   else
 419     return dist - rad;
 420 }
 421
 422
 423 void
 424 transform_point (Matrix m, Point *src, Point *dest)
 425 {
 426   real xx, yy, ww;
 427
 428   xx = m[0][0] * src->x + m[0][1] * src->y + m[0][2];
 429   yy = m[1][0] * src->x + m[1][1] * src->y + m[1][2];
 430   ww = m[2][0] * src->x + m[2][1] * src->y + m[2][2];
 431
 432   if (!ww)
 433     ww = 1.0;
 434
 435   dest->x = xx / ww;
 436   dest->y = yy / ww;
 437 }
 438
 439 void
 440 mult_matrix (Matrix m1, Matrix m2)
 441 {
 442   Matrix result;
 443   int i, j, k;
 444
 445   for (i = 0; i < 3; i++)
 446     for (j = 0; j < 3; j++)
 447       {
 448         result [i][j] = 0.0;
 449         for (k = 0; k < 3; k++)
 450           result [i][j] += m1 [i][k] * m2[k][j];
 451       }
 452
 453   /*  copy the result into matrix 2  */
 454   for (i = 0; i < 3; i++)
 455     for (j = 0; j < 3; j++)
 456       m2 [i][j] = result [i][j];
 457 }
 458
 459 void
 460 identity_matrix (Matrix m)
 461 {
 462   int i, j;
 463
 464   for (i = 0; i < 3; i++)
 465     for (j = 0; j < 3; j++)
 466       m[i][j] = (i == j) ? 1 : 0;
 467
 468 }
 469
 470 void
 471 translate_matrix (Matrix m, real x, real y)
 472 {
 473   Matrix trans;
 474
 475   identity_matrix (trans);
 476   trans[0][2] = x;
 477   trans[1][2] = y;
 478   mult_matrix (trans, m);
 479 }
 480
 481 void
 482 scale_matrix (Matrix m, real x, real y)
 483 {
 484   Matrix scale;
 485
 486   identity_matrix (scale);
 487   scale[0][0] = x;
 488   scale[1][1] = y;
 489   mult_matrix (scale, m);
 490 }
 491
 492 void
 493 rotate_matrix (Matrix m, real theta)
 494 {
 495   Matrix rotate;
 496   real cos_theta, sin_theta;
 497
 498   cos_theta = cos (theta);
 499   sin_theta = sin (theta);
 500
 501   identity_matrix (rotate);
 502   rotate[0][0] = cos_theta;
 503   rotate[0][1] = -sin_theta;
 504   rotate[1][0] = sin_theta;
 505   rotate[1][1] = cos_theta;
 506   mult_matrix (rotate, m);
 507 }
 508
 509 void
 510 xshear_matrix (Matrix m, real shear)
 511 {
 512   Matrix shear_m;
 513
 514   identity_matrix (shear_m);
 515   shear_m[0][1] = shear;
 516   mult_matrix (shear_m, m);
 517 }
 518
 519 void
 520 yshear_matrix (Matrix m, real shear)
 521 {
 522   Matrix shear_m;
 523
 524   identity_matrix (shear_m);
 525   shear_m[1][0] = shear;
 526   mult_matrix (shear_m, m);
 527 }
 528
 529 /*  find the determinate for a 3x3 matrix  */
 530 static real
 531 determinate (Matrix m)
 532 {
 533   int i;
 534   double det = 0;
 535
 536   for (i = 0; i < 3; i ++)
 537     {
 538       det += m[0][i] * m[1][(i+1)%3] * m[2][(i+2)%3];
 539       det -= m[2][i] * m[1][(i+1)%3] * m[0][(i+2)%3];
 540     }
 541
 542   return det;
 543 }
 544
 545
 546
 547
 548
 549