src/intersect.c

   1 /*
   2  *                            COPYRIGHT
   3  *
   4  *  PCB, interactive printed circuit board design
   5  *  Copyright (C) 1994,1995,1996 Thomas Nau
   6  *  Copyright (C) 1998,1999,2000,2001 harry eaton
   7  *
   8  *  This program is free software; you can redistribute it and/or modify
   9  *  it under the terms of the GNU General Public License as published by
  10  *  the Free Software Foundation; either version 2 of the License, or
  11  *  (at your option) any later version.
  12  *
  13  *  This program is distributed in the hope that it will be useful,
  14  *  but WITHOUT ANY WARRANTY; without even the implied warranty of
  15  *  MERCHANTABILITY or FITNESS FOR A PARTICULAR PURPOSE.  See the
  16  *  GNU General Public License for more details.
  17  *
  18  *  You should have received a copy of the GNU General Public License
  19  *  along with this program; if not, write to the Free Software
  20  *  Foundation, Inc., 675 Mass Ave, Cambridge, MA 02139, USA.
  21  *
  22  *  Contact addresses for paper mail and Email:
  23  *  harry eaton, 6697 Buttonhole Ct, Columbia, MD 21044 USA
  24  *  haceaton@aplcomm.jhuapl.edu
  25  *
  26  */
  27
  28 /* this file, intersect.c, was written and is
  29  * Copyright (c) 2001 C. Scott Ananian
  30  */
  31
  32 /* rectangle intersection/union routines.
  33  */
  34
  35 #ifdef HAVE_CONFIG_H
  36 #include "config.h"
  37 #endif
  38
  39 #include "global.h"
  40
  41 #include <assert.h>
  42
  43 #include "data.h"
  44 #include "intersect.h"
  45 #include "mymem.h"
  46
  47 #ifdef HAVE_LIBDMALLOC
  48 #include <dmalloc.h>
  49 #endif
  50
  51 /* ---------------------------------------------------------------------------
  52  * some local prototypes
  53  */
  54 static int compareleft (const void *ptr1, const void *ptr2);
  55 static int compareright (const void *ptr1, const void *ptr2);
  56 static int comparepos (const void *ptr1, const void *ptr2);
  57 static int nextpwrof2 (int i);
  58
  59 /* ---------------------------------------------------------------------------
  60  * some local types
  61  */
  62 typedef struct
  63 {
  64   Coord left, right;
  65   int covered;
  66   Coord area;
  67 }
  68 SegmentTreeNode;
  69
  70 typedef struct
  71 {
  72   SegmentTreeNode *nodes;
  73   int size;
  74 }
  75 SegmentTree;
  76
  77 typedef struct
  78 {
  79   Coord *p;
  80   int size;
  81 }
  82 LocationList;
  83
  84 /* ---------------------------------------------------------------------------
  85  * Create a sorted list of unique y coords from a BoxList.
  86  */
  87 static LocationList
  88 createSortedYList (BoxListType *boxlist)
  89 {
  90   LocationList yCoords;
  91   Coord last;
  92   int i, n;
  93   /* create sorted list of Y coordinates */
  94   yCoords.size = 2 * boxlist->BoxN;
  95   yCoords.p = (Coord *)calloc (yCoords.size, sizeof (*yCoords.p));
  96   for (i = 0; i < boxlist->BoxN; i++)
  97     {
  98       yCoords.p[2 * i] = boxlist->Box[i].Y1;
  99       yCoords.p[2 * i + 1] = boxlist->Box[i].Y2;
 100     }
 101   qsort (yCoords.p, yCoords.size, sizeof (*yCoords.p), comparepos);
 102   /* count uniq y coords */
 103   last = 0;
 104   for (n = 0, i = 0; i < yCoords.size; i++)
 105     if (i == 0 || yCoords.p[i] != last)
 106       yCoords.p[n++] = last = yCoords.p[i];
 107   yCoords.size = n;
 108   return yCoords;
 109 }
 110
 111 /* ---------------------------------------------------------------------------
 112  * Create an empty segment tree from the given sorted list of uniq y coords.
 113  */
 114 static SegmentTree
 115 createSegmentTree (Coord * yCoords, int N)
 116 {
 117   SegmentTree st;
 118   int i;
 119   /* size is twice the nearest larger power of 2 */
 120   st.size = 2 * nextpwrof2 (N);
 121   st.nodes = (SegmentTreeNode *)calloc (st.size, sizeof (*st.nodes));
 122   /* initialize the rightmost leaf node */
 123   st.nodes[st.size - 1].left = (N > 0) ? yCoords[--N] : 10;
 124   st.nodes[st.size - 1].right = st.nodes[st.size - 1].left + 1;
 125   /* initialize the rest of the leaf nodes */
 126   for (i = st.size - 2; i >= st.size / 2; i--)
 127     {
 128       st.nodes[i].right = st.nodes[i + 1].left;
 129       st.nodes[i].left = (N > 0) ? yCoords[--N] : st.nodes[i].right - 1;
 130     }
 131   /* initialize the internal nodes */
 132   for (; i > 0; i--)
 133     {                           /* node 0 is not used */
 134       st.nodes[i].right = st.nodes[2 * i + 1].right;
 135       st.nodes[i].left = st.nodes[2 * i].left;
 136     }
 137   /* done! */
 138   return st;
 139 }
 140
 141 void
 142 insertSegment (SegmentTree * st, int n, Coord Y1, Coord Y2)
 143 {
 144   Coord discriminant;
 145   if (st->nodes[n].left >= Y1 && st->nodes[n].right <= Y2)
 146     {
 147       st->nodes[n].covered++;
 148     }
 149   else
 150     {
 151       assert (n < st->size / 2);
 152       discriminant = st->nodes[n * 2 + 1 /*right */ ].left;
 153       if (Y1 < discriminant)
 154         insertSegment (st, n * 2, Y1, Y2);
 155       if (discriminant < Y2)
 156         insertSegment (st, n * 2 + 1, Y1, Y2);
 157     }
 158   /* fixup area */
 159   st->nodes[n].area = (st->nodes[n].covered > 0) ?
 160     (st->nodes[n].right - st->nodes[n].left) :
 161     (n >= st->size / 2) ? 0 :
 162     st->nodes[n * 2].area + st->nodes[n * 2 + 1].area;
 163 }
 164
 165 void
 166 deleteSegment (SegmentTree * st, int n, Coord Y1, Coord Y2)
 167 {
 168   Coord discriminant;
 169   if (st->nodes[n].left >= Y1 && st->nodes[n].right <= Y2)
 170     {
 171       assert (st->nodes[n].covered);
 172       --st->nodes[n].covered;
 173     }
 174   else
 175     {
 176       assert (n < st->size / 2);
 177       discriminant = st->nodes[n * 2 + 1 /*right */ ].left;
 178       if (Y1 < discriminant)
 179         deleteSegment (st, n * 2, Y1, Y2);
 180       if (discriminant < Y2)
 181         deleteSegment (st, n * 2 + 1, Y1, Y2);
 182     }
 183   /* fixup area */
 184   st->nodes[n].area = (st->nodes[n].covered > 0) ?
 185     (st->nodes[n].right - st->nodes[n].left) :
 186     (n >= st->size / 2) ? 0 :
 187     st->nodes[n * 2].area + st->nodes[n * 2 + 1].area;
 188 }
 189
 190 /* ---------------------------------------------------------------------------
 191  * Compute the area of the intersection of the given rectangles; that is,
 192  * the area covered by more than one rectangle (counted twice if the area is
 193  * covered by three rectangles, three times if covered by four rectangles,
 194  * etc.).
 195  * Runs in O(N ln N) time.
 196  */
 197 double
 198 ComputeIntersectionArea (BoxListType *boxlist)
 199 {
 200   Cardinal i;
 201   double area = 0.0;
 202   /* first get the aggregate area. */
 203   for (i = 0; i < boxlist->BoxN; i++)
 204     area += (double) (boxlist->Box[i].X2 - boxlist->Box[i].X1) *
 205       (double) (boxlist->Box[i].Y2 - boxlist->Box[i].Y1);
 206   /* intersection area is aggregate - union. */
 207   return area * 0.0001 - ComputeUnionArea (boxlist);
 208 }
 209
 210 /* ---------------------------------------------------------------------------
 211  * Compute the area of the union of the given rectangles.
 212  * O(N ln N) time.
 213  */
 214 double
 215 ComputeUnionArea (BoxListType *boxlist)
 216 {
 217   BoxType **rectLeft, **rectRight;
 218   Cardinal i, j;
 219   LocationList yCoords;
 220   SegmentTree segtree;
 221   Coord lastX;
 222   double area = 0.0;
 223
 224   if (boxlist->BoxN == 0)
 225     return 0.0;
 226   /* create sorted list of Y coordinates */
 227   yCoords = createSortedYList (boxlist);
 228   /* now create empty segment tree */
 229   segtree = createSegmentTree (yCoords.p, yCoords.size);
 230   free (yCoords.p);
 231   /* create sorted list of left and right X coordinates of rectangles */
 232   rectLeft = (BoxType **)calloc (boxlist->BoxN, sizeof (*rectLeft));
 233   rectRight = (BoxType **)calloc (boxlist->BoxN, sizeof (*rectRight));
 234   for (i = 0; i < boxlist->BoxN; i++)
 235     {
 236       assert (boxlist->Box[i].X1 <= boxlist->Box[i].X2);
 237       assert (boxlist->Box[i].Y1 <= boxlist->Box[i].Y2);
 238       rectLeft[i] = rectRight[i] = &boxlist->Box[i];
 239     }
 240   qsort (rectLeft, boxlist->BoxN, sizeof (*rectLeft), compareleft);
 241   qsort (rectRight, boxlist->BoxN, sizeof (*rectRight), compareright);
 242   /* sweep through x segments from left to right */
 243   i = j = 0;
 244   lastX = rectLeft[0]->X1;
 245   while (j < boxlist->BoxN)
 246     {
 247       assert (i <= boxlist->BoxN);
 248       /* i will step through rectLeft, j will through rectRight */
 249       if (i == boxlist->BoxN || rectRight[j]->X2 < rectLeft[i]->X1)
 250         {
 251           /* right edge of rectangle */
 252           BoxType *b = rectRight[j++];
 253           /* check lastX */
 254           if (b->X2 != lastX)
 255             {
 256               assert (lastX < b->X2);
 257               area += (double) (b->X2 - lastX) * segtree.nodes[1].area;
 258               lastX = b->X2;
 259             }
 260           /* remove a segment from the segment tree. */
 261           deleteSegment (&segtree, 1, b->Y1, b->Y2);
 262         }
 263       else
 264         {
 265           /* left edge of rectangle */
 266           BoxType *b = rectLeft[i++];
 267           /* check lastX */
 268           if (b->X1 != lastX)
 269             {
 270               assert (lastX < b->X1);
 271               area += (double) (b->X1 - lastX) * segtree.nodes[1].area;
 272               lastX = b->X1;
 273             }
 274           /* add a segment from the segment tree. */
 275           insertSegment (&segtree, 1, b->Y1, b->Y2);
 276         }
 277     }
 278   free (rectLeft);
 279   free (rectRight);
 280   free (segtree.nodes);
 281   return area * 0.0001;
 282 }
 283 static int
 284 compareleft (const void *ptr1, const void *ptr2)
 285 {
 286   BoxType **b1 = (BoxType **) ptr1;
 287   BoxType **b2 = (BoxType **) ptr2;
 288   return (*b1)->X1 - (*b2)->X1;
 289 }
 290 static int
 291 compareright (const void *ptr1, const void *ptr2)
 292 {
 293   BoxType **b1 = (BoxType **) ptr1;
 294   BoxType **b2 = (BoxType **) ptr2;
 295   return (*b1)->X2 - (*b2)->X2;
 296 }
 297 static int
 298 comparepos (const void *ptr1, const void *ptr2)
 299 {
 300   return *((Coord *) ptr1) - *((Coord *) ptr2);
 301 }
 302 static int
 303 nextpwrof2 (int n)
 304 {
 305   int r = 1;
 306   while (n != 0)
 307     {
 308       n /= 2;
 309       r *= 2;
 310     }
 311   return r;
 312 }