src/intersect.c

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