src/region_utils.c

   1 /*
   2  * Grace - Graphics for Exploratory Data Analysis
   3  *
   4  * Home page: http://plasma-gate.weizmann.ac.il/Grace/
   5  *
   6  * Copyright (c) 1991-1995 Paul J Turner, Portland, OR
   7  * Copyright (c) 1996-2003 Grace Development Team
   8  *
   9  * Maintained by Evgeny Stambulchik
  10  *
  11  *
  12  *                           All Rights Reserved
  13  *
  14  *    This program is free software; you can redistribute it and/or modify
  15  *    it under the terms of the GNU General Public License as published by
  16  *    the Free Software Foundation; either version 2 of the License, or
  17  *    (at your option) any later version.
  18  *
  19  *    This program is distributed in the hope that it will be useful,
  20  *    but WITHOUT ANY WARRANTY; without even the implied warranty of
  21  *    MERCHANTABILITY or FITNESS FOR A PARTICULAR PURPOSE.  See the
  22  *    GNU General Public License for more details.
  23  *
  24  *    You should have received a copy of the GNU General Public License
  25  *    along with this program; if not, write to the Free Software
  26  *    Foundation, Inc., 675 Mass Ave, Cambridge, MA 02139, USA.
  27  */
  28
  29 /*
  30  *
  31  * routines to allocate, manipulate, and return
  32  * information about regions.
  33  *
  34  */
  35
  36 #include <config.h>
  37
  38 #include <stdio.h>
  39 #include <stdlib.h>
  40
  41 #include "core_utils.h"
  42 #include "protos.h"
  43
  44 /*
  45  * routines to determine if a point lies in a polygon
  46 */
  47 static int intersect_to_left(const WPoint *wp,
  48     const WPoint *wp1, const WPoint *wp2)
  49 {
  50     double xtmp, m, b;
  51
  52     /* ignore horizontal lines */
  53     if (wp1->y == wp2->y) {
  54         return 0;
  55     }
  56     /* not contained vertically */
  57     if (((wp->y < wp1->y) && (wp->y < wp2->y)) ||
  58         ((wp->y > wp1->y) && (wp->y > wp2->y))) {
  59         return 0;
  60     }
  61     /* none of the above, compute the intersection */
  62     if ((xtmp = wp2->x - wp1->x) != 0.0) {
  63         m = (wp2->y - wp1->y) / xtmp;
  64         b = wp1->y - m * wp1->x;
  65         xtmp = (wp->y - b) / m;
  66     } else {
  67         xtmp = wp1->x;
  68     }
  69     if (xtmp <= wp->x) {
  70         /* check for intersections at a vertex */
  71         /* if this is the max ordinate then accept */
  72         if (wp->y == wp1->y) {
  73             if (wp1->y > wp2->y) {
  74                 return 1;
  75             } else {
  76                 return 0;
  77             }
  78         }
  79         /* check for intersections at a vertex */
  80         if (wp->y == wp2->y) {
  81             if (wp2->y > wp1->y) {
  82                 return 1;
  83             } else {
  84                 return 0;
  85             }
  86         }
  87         /* no vertices intersected */
  88         return 1;
  89     }
  90     return 0;
  91 }
  92
  93 /*
  94  * determine if (x,y) is in the polygon xlist[], ylist[]
  95  */
  96 static int inbound(const WPoint *wp, const WPoint *wps, int n)
  97 {
  98     int i, l = 0;
  99
 100     for (i = 0; i < n; i++) {
 101         l += intersect_to_left(wp, &wps[i], &wps[(i + 1) % n]);
 102     }
 103     return l % 2;
 104 }
 105
 106 static int isleft(const WPoint *wp, const WPoint *wp1, const WPoint *wp2)
 107 {
 108     if ((wp2->x - wp1->x)*(wp->y - wp2->y) -
 109         (wp2->y - wp1->y)*(wp->x - wp2->x) >= 0) {
 110         return TRUE;
 111     } else {
 112         return FALSE;
 113     }
 114 }
 115
 116 static int inband(const WPoint *wp, const WPoint *wp1, const WPoint *wp2)
 117 {
 118     double s, d2;
 119
 120     s  = (wp->x - wp1->x)*(wp2->x - wp1->x) +
 121          (wp->y - wp1->y)*(wp2->y - wp1->y);
 122     d2 = (wp2->x - wp1->x)*(wp2->x - wp1->x) +
 123          (wp2->y - wp1->y)*(wp2->y - wp1->y);
 124     if (s < 0 || s > d2) {
 125         return FALSE;
 126     } else {
 127         return TRUE;
 128     }
 129 }
 130
 131 int inregion(const Quark *q, const WPoint *wp)
 132 {
 133     region *r = region_get_data(q);
 134
 135     if (!r) {
 136         return FALSE;
 137     }
 138
 139     switch (r->type) {
 140     case REGION_POLYGON:
 141         if (r->n > 2) {
 142             return inbound(wp, r->wps, r->n);
 143         } else
 144         if (r->n == 2) {
 145             return isleft(wp, &r->wps[0], &r->wps[1]);
 146         } else {
 147             return FALSE;
 148         }
 149         break;
 150     case REGION_BAND:
 151         if (r->n == 2) {
 152             return inband(wp, &r->wps[0], &r->wps[1]);
 153         } else {
 154             return FALSE;
 155         }
 156         break;
 157     }
 158
 159     return FALSE;
 160 }