home *** CD-ROM | disk | FTP | other *** search
/ Usenet 1994 January / usenetsourcesnewsgroupsinfomagicjanuary1994.iso / sources / x / volume19 / xfig / part07 / u_geom.c < prev   
Encoding:
C/C++ Source or Header  |  1993-05-27  |  8.9 KB  |  396 lines
  1. /*
  2.  * FIG : Facility for Interactive Generation of figures
  3.  * Copyright (c) 1985 by Supoj Sutanthavibul
  4.  *
  5.  * "Permission to use, copy, modify, distribute, and sell this software and its
  6.  * documentation for any purpose is hereby granted without fee, provided that
  7.  * the above copyright notice appear in all copies and that both the copyright
  8.  * notice and this permission notice appear in supporting documentation. 
  9.  * No representations are made about the suitability of this software for 
  10.  * any purpose.  It is provided "as is" without express or implied warranty."
  11.  */
  12.  
  13. /*
  14.  * Routines dealing with geometry under the following headings:
  15.  *    COMPUTE NORMAL, CLOSE TO VECTOR, COMPUTE ARC CENTER,
  16.  *    COMPUTE ANGLE, COMPUTE DIRECTION, LATEX LINE ROUTINES.
  17.  */
  18.  
  19. #include "fig.h"
  20. #include "resources.h"
  21. #include "object.h"
  22.  
  23. /*************************** COMPUTE NORMAL **********************
  24.  
  25. Input arguments :
  26.     (x1,y1)(x2,y2) : the vector
  27.     direction : direction of the normal vector to (x1,y1)(x2,y2)
  28. Output arguments :
  29.     (*x,*y)(x2,y2) : a normal vector.
  30. Return value : none
  31.  
  32. ******************************************************************/
  33.  
  34. compute_normal(x1, y1, x2, y2, direction, x, y)
  35.     float        x1, y1;
  36.     int            x2, y2, direction, *x, *y;
  37. {
  38.     if (direction) {        /* counter clockwise  */
  39.     *x = round(x2 - (y2 - y1));
  40.     *y = round(y2 - (x1 - x2));
  41.     } else {
  42.     *x = round(x2 + (y2 - y1));
  43.     *y = round(y2 + (x1 - x2));
  44.     }
  45. }
  46.  
  47. /******************** CLOSE TO VECTOR **************************
  48.  
  49. Input arguments:
  50.     (x1,y1)(x2,y2) : the vector
  51.     (xp,yp) : the point
  52.     d : tolerance (max. allowable distance from the point to the vector)
  53.     dd : d * d
  54. Output arguments:
  55.     (*px,*py) : a point on the vector which is not far from (xp,yp)
  56.         by more than d. Normally the vector (*px,*py)(xp,yp)
  57.         is normal to vector (x1,y1)(x2,y2) except when (xp,yp)
  58.         is within d from (x1,y1) or (x2,y2), in which cases,
  59.         (*px,*py) = (x1,y1) or (x2,y2) respectively.
  60. Return value :
  61.     0 : No point on the vector is within d from (xp, yp)
  62.     1 : (*px, *py) is such a point.
  63.  
  64. ******************************************************************/
  65.  
  66. close_to_vector(x1, y1, x2, y2, xp, yp, d, dd, px, py)
  67.     int            x1, y1, x2, y2, xp, yp, d;
  68.     float        dd;
  69.     int           *px, *py;
  70. {
  71.     int            xmin, ymin, xmax, ymax;
  72.     float        x, y, slope, D2, dx, dy;
  73.  
  74.     if (abs(xp - x1) <= d && abs(yp - y1) <= d) {
  75.     *px = x1;
  76.     *py = y1;
  77.     return (1);
  78.     }
  79.     if (abs(xp - x2) <= d && abs(yp - y2) <= d) {
  80.     *px = x2;
  81.     *py = y2;
  82.     return (1);
  83.     }
  84.     if (x1 < x2) {
  85.     xmin = x1 - d;
  86.     xmax = x2 + d;
  87.     } else {
  88.     xmin = x2 - d;
  89.     xmax = x1 + d;
  90.     }
  91.     if (xp < xmin || xmax < xp)
  92.     return (0);
  93.  
  94.     if (y1 < y2) {
  95.     ymin = y1 - d;
  96.     ymax = y2 + d;
  97.     } else {
  98.     ymin = y2 - d;
  99.     ymax = y1 + d;
  100.     }
  101.     if (yp < ymin || ymax < yp)
  102.     return (0);
  103.  
  104.     if (x2 == x1) {
  105.     x = x1;
  106.     y = yp;
  107.     } else if (y1 == y2) {
  108.     x = xp;
  109.     y = y1;
  110.     } else {
  111.     slope = ((float) (x2 - x1)) / ((float) (y2 - y1));
  112.     y = (slope * (xp - x1 + slope * y1) + yp) / (1 + slope * slope);
  113.     x = ((float) x1) + slope * (y - y1);
  114.     }
  115.     dx = ((float) xp) - x;
  116.     dy = ((float) yp) - y;
  117.     D2 = dx * dx + dy * dy;
  118.     if (D2 < dd) {
  119.     *px = (int) (x + .5);
  120.     *py = (int) (y + .5);
  121.     return (1);
  122.     }
  123.     return (0);
  124. }
  125.  
  126. /********************* COMPUTE ARC CENTER ******************
  127.  
  128. Input arguments :
  129.     p1, p2, p3 : 3 points on the arc
  130. Output arguments :
  131.     (*x,*y) : Center of the arc
  132. Return value :
  133.     0 : if p1, p2, p3 are co-linear.
  134.     1 : if they are not.
  135.  
  136. *************************************************************/
  137.  
  138. int
  139. compute_arccenter(p1, p2, p3, x, y)
  140.     F_pos        p1, p2, p3;
  141.     float       *x, *y;
  142. {
  143.     float        s12, s13, len1, len2, len3, dx12, dy12, dx13, dy13;
  144.  
  145.     if (p1.x == p3.x && p1.y == p3.y)
  146.     return 0;
  147.  
  148.     dx12 = p1.x - p2.x;
  149.     dy12 = p1.y - p2.y;
  150.     dx13 = p1.x - p3.x;
  151.     dy13 = p1.y - p3.y;
  152.  
  153.     s12 = asin((double) (dy12 / sqrt((double) (dx12 * dx12 + dy12 * dy12))));
  154.     s13 = asin((double) (dy13 / sqrt((double) (dx13 * dx13 + dy13 * dy13))));
  155.     if (fabs(s12 - s13) < .01)
  156.     return 0;
  157.  
  158.     len1 = p1.x * p1.x + p1.y * p1.y;
  159.     len2 = p2.x * p2.x + p2.y * p2.y;
  160.     len3 = p3.x * p3.x + p3.y * p3.y;
  161.     *y = (dx12 * (len3 - len1) - dx13 * (len2 - len1)) /
  162.     (2 * (dx13 * dy12 - dx12 * dy13));
  163.     if (p1.x != p3.x)
  164.     *x = (len3 + 2 * (*y) * dy13 - len1) / (2 * (-dx13));
  165.     else
  166.     *x = (len2 + 2 * (*y) * dy12 - len1) / (2 * (-dx12));
  167.     return 1;
  168. }
  169.  
  170. /********************* COMPUTE ANGLE ************************
  171.  
  172. Input arguments :
  173.     (dx,dy) : the vector (0,0)(dx,dy)
  174. Output arguments : none
  175. Return value : the angle of the vector in the range [0, 2PI)
  176.  
  177. *************************************************************/
  178.  
  179. float
  180. compute_angle(dx, dy)        /* compute the angle between 0 to 2PI  */
  181.     float        dx, dy;
  182. {
  183.     float        alpha;
  184.  
  185.     if (dx == 0) {
  186.     if (dy > 0)
  187.         alpha = M_PI_2;
  188.     else
  189.         alpha = 3 * M_PI_2;
  190.     } else if (dy == 0) {
  191.     if (dx > 0)
  192.         alpha = 0;
  193.     else
  194.         alpha = M_PI;
  195.     } else {
  196.     alpha = atan((double) (dy / dx));    /* range = -PI/2 to PI/2 */
  197.     if (dx < 0)
  198.         alpha += M_PI;
  199.     else if (dy < 0)
  200.         alpha += M_2PI;
  201.     }
  202.     return (alpha);
  203. }
  204.  
  205.  
  206. /********************* COMPUTE DIRECTION ********************
  207.  
  208. Input arguments :
  209.     p1, p2, p3 : 3 points of an arc with p1 the first and p3 the last.
  210. Output arguments : none
  211. Return value :
  212.     0 : if the arc passes p1, p2 and p3 (in that order) in
  213.         clockwise direction
  214.     1 : if direction is counterclockwise
  215.  
  216. *************************************************************/
  217.  
  218. int
  219. compute_direction(p1, p2, p3)
  220.     F_pos        p1, p2, p3;
  221. {
  222.     float        diff, dx, dy, alpha, theta;
  223.  
  224.     dx = p2.x - p1.x;
  225.     dy = p1.y - p2.y;        /* because origin of the screen is on the
  226.                  * upper left corner */
  227.  
  228.     alpha = compute_angle(dx, dy);
  229.  
  230.     dx = p3.x - p2.x;
  231.     dy = p2.y - p3.y;
  232.     theta = compute_angle(dx, dy);
  233.  
  234.     diff = theta - alpha;
  235.     if ((0 < diff && diff < M_PI) || diff < -M_PI) {
  236.     return (1);        /* counterclockwise */
  237.     }
  238.     return (0);            /* clockwise */
  239. }
  240.  
  241. /*********************** LATEX LINE ROUTINES ***************************/
  242.  
  243. int
  244. pgcd(a, b)
  245.     int            a, b;
  246.  
  247. /*
  248.  * compute greatest common divisor, assuming 0 < a <= b
  249.  */
  250. {
  251.     b = b % a;
  252.     return (b) ? gcd(b, a) : a;
  253. }
  254.  
  255. int
  256. gcd(a, b)
  257.     int            a, b;
  258.  
  259. /*
  260.  * compute greatest common divisor
  261.  */
  262. {
  263.     if (a < 0)
  264.     a = -a;
  265.     if (b < 0)
  266.     b = -b;
  267.     return (a <= b) ? pgcd(a, b) : pgcd(b, a);
  268. }
  269.  
  270.  
  271. int
  272. lcm(a, b)
  273.     int            a, b;
  274.  
  275. /*
  276.  * compute least common multiple
  277.  */
  278. {
  279.     return abs(a * b) / gcd(a, b);
  280. }
  281.  
  282.  
  283. double        rad2deg = 57.295779513082320877;
  284.  
  285. struct angle_table {
  286.     int            x, y;
  287.     double        angle;
  288. };
  289.  
  290. struct angle_table line_angles[25] =
  291. {{0, 1, 90.0},
  292. {1, 0, 0.0},
  293. {1, 1, 45.0},
  294. {1, 2, 63.434948822922010648},
  295. {1, 3, 71.565051177077989351},
  296. {1, 4, 75.963756532073521417},
  297. {1, 5, 78.690067525979786913},
  298. {1, 6, 80.537677791974382609},
  299. {2, 1, 26.565051177077989351},
  300. {2, 3, 56.309932474020213086},
  301. {2, 5, 68.198590513648188229},
  302. {3, 1, 18.434948822922010648},
  303. {3, 2, 33.690067525979786913},
  304. {3, 4, 53.130102354155978703},
  305. {3, 5, 59.036243467926478582},
  306. {4, 1, 14.036243467926478588},
  307. {4, 3, 36.869897645844021297},
  308. {4, 5, 51.340191745909909396},
  309. {5, 1, 11.309932474020213086},
  310. {5, 2, 21.801409486351811770},
  311. {5, 3, 30.963756532073521417},
  312. {5, 4, 38.659808254090090604},
  313. {5, 6, 50.194428907734805993},
  314. {6, 1, 9.4623222080256173906},
  315. {6, 5, 39.805571092265194006}
  316. };
  317.  
  318. struct angle_table arrow_angles[13] =
  319. {{0, 1, 90.0},
  320. {1, 0, 0.0},
  321. {1, 1, 45.0},
  322. {1, 2, 63.434948822922010648},
  323. {1, 3, 71.565051177077989351},
  324. {1, 4, 75.963756532073521417},
  325. {2, 1, 26.565051177077989351},
  326. {2, 3, 56.309932474020213086},
  327. {3, 1, 18.434948822922010648},
  328. {3, 2, 33.690067525979786913},
  329. {3, 4, 53.130102354155978703},
  330. {4, 1, 14.036243467926478588},
  331. {4, 3, 36.869897645844021297},
  332. };
  333.  
  334. get_slope(dx, dy, sxp, syp, arrow)
  335.     int            dx, dy, *sxp, *syp, arrow;
  336. {
  337.     double        angle;
  338.     int            i, s, max;
  339.     double        d, d1;
  340.     struct angle_table *st;
  341.  
  342.     if (dx == 0) {
  343.     *sxp = 0;
  344.     *syp = signof(dy);
  345.     return;
  346.     }
  347.     angle = atan((double) abs(dy) / (double) abs(dx)) * rad2deg;
  348.     if (arrow) {
  349.     st = arrow_angles;
  350.     max = 13;
  351.     } else {
  352.     st = line_angles;
  353.     max = 25;
  354.     }
  355.     s = 0;
  356.     d = 9.9e9;
  357.     for (i = 0; i < max; i++) {
  358.     d1 = fabs(angle - st[i].angle);
  359.     if (d1 < d) {
  360.         s = i;
  361.         d = d1;
  362.     }
  363.     }
  364.     *sxp = st[s].x;
  365.     if (dx < 0)
  366.     *sxp = -*sxp;
  367.     *syp = st[s].y;
  368.     if (dy < 0)
  369.     *syp = -*syp;
  370. }
  371.  
  372. latex_endpoint(x1, y1, x2, y2, xout, yout, arrow, magnet)
  373.     int            x1, y1, x2, y2;
  374.     int           *xout, *yout;
  375.     int            arrow, magnet;
  376. {
  377.     int            dx, dy, sx, sy, ds, dsx, dsy;
  378.  
  379.     dx = x2 - x1;
  380.     dy = y2 - y1;
  381.     get_slope(dx, dy, &sx, &sy, arrow);
  382.     if (abs(sx) >= abs(sy)) {
  383.     ds = lcm(sx, magnet * gcd(sx, magnet));
  384.     dsx = (2 * abs(dx) / ds + 1) / 2;
  385.     dsx = (dx >= 0) ? dsx * ds : -dsx * ds;
  386.     *xout = x1 + dsx;
  387.     *yout = y1 + dsx * sy / sx;
  388.     } else {
  389.     ds = lcm(sy, magnet * gcd(sy, magnet));
  390.     dsy = (2 * abs(dy) / ds + 1) / 2;
  391.     dsy = (dy >= 0) ? dsy * ds : -dsy * ds;
  392.     *yout = y1 + dsy;
  393.     *xout = x1 + dsy * sx / sy;
  394.     }
  395. }
  396.