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   Coord left, right;
  70   int covered;
  71   Coord area;
  72 }
  73 SegmentTreeNode;
  74
  75 typedef struct
  76 {
  77   SegmentTreeNode *nodes;
  78   int size;
  79 }
  80 SegmentTree;
  81
  82 typedef struct
  83 {
  84   Coord *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   Coord last;
  97   int i, n;
  98   /* create sorted list of Y coordinates */
  99   yCoords.size = 2 * boxlist->BoxN;
 100   yCoords.p = (Coord *)calloc (yCoords.size, sizeof (*yCoords.p));
 101   for (i = 0; i < boxlist->BoxN; i++)
 102     {
 103       yCoords.p[2 * i] = boxlist->Box[i].Y1;
 104       yCoords.p[2 * i + 1] = boxlist->Box[i].Y2;
 105     }
 106   qsort (yCoords.p, yCoords.size, sizeof (*yCoords.p), comparepos);
 107   /* count uniq y coords */
 108   last = 0;
 109   for (n = 0, i = 0; i < yCoords.size; i++)
 110     if (i == 0 || yCoords.p[i] != last)
 111       yCoords.p[n++] = last = yCoords.p[i];
 112   yCoords.size = n;
 113   return yCoords;
 114 }
 115
 116 /* ---------------------------------------------------------------------------
 117  * Create an empty segment tree from the given sorted list of uniq y coords.
 118  */
 119 static SegmentTree
 120 createSegmentTree (Coord * yCoords, int N)
 121 {
 122   SegmentTree st;
 123   int i;
 124   /* size is twice the nearest larger power of 2 */
 125   st.size = 2 * nextpwrof2 (N);
 126   st.nodes = (SegmentTreeNode *)calloc (st.size, sizeof (*st.nodes));
 127   /* initialize the rightmost leaf node */
 128   st.nodes[st.size - 1].left = (N > 0) ? yCoords[--N] : 10;
 129   st.nodes[st.size - 1].right = st.nodes[st.size - 1].left + 1;
 130   /* initialize the rest of the leaf nodes */
 131   for (i = st.size - 2; i >= st.size / 2; i--)
 132     {
 133       st.nodes[i].right = st.nodes[i + 1].left;
 134       st.nodes[i].left = (N > 0) ? yCoords[--N] : st.nodes[i].right - 1;
 135     }
 136   /* initialize the internal nodes */
 137   for (; i > 0; i--)
 138     {                           /* node 0 is not used */
 139       st.nodes[i].right = st.nodes[2 * i + 1].right;
 140       st.nodes[i].left = st.nodes[2 * i].left;
 141     }
 142   /* done! */
 143   return st;
 144 }
 145
 146 void
 147 insertSegment (SegmentTree * st, int n, Coord Y1, Coord Y2)
 148 {
 149   Coord discriminant;
 150   if (st->nodes[n].left >= Y1 && st->nodes[n].right <= Y2)
 151     {
 152       st->nodes[n].covered++;
 153     }
 154   else
 155     {
 156       assert (n < st->size / 2);
 157       discriminant = st->nodes[n * 2 + 1 /*right */ ].left;
 158       if (Y1 < discriminant)
 159         insertSegment (st, n * 2, Y1, Y2);
 160       if (discriminant < Y2)
 161         insertSegment (st, n * 2 + 1, Y1, Y2);
 162     }
 163   /* fixup area */
 164   st->nodes[n].area = (st->nodes[n].covered > 0) ?
 165     (st->nodes[n].right - st->nodes[n].left) :
 166     (n >= st->size / 2) ? 0 :
 167     st->nodes[n * 2].area + st->nodes[n * 2 + 1].area;
 168 }
 169
 170 void
 171 deleteSegment (SegmentTree * st, int n, Coord Y1, Coord Y2)
 172 {
 173   Coord discriminant;
 174   if (st->nodes[n].left >= Y1 && st->nodes[n].right <= Y2)
 175     {
 176       assert (st->nodes[n].covered);
 177       --st->nodes[n].covered;
 178     }
 179   else
 180     {
 181       assert (n < st->size / 2);
 182       discriminant = st->nodes[n * 2 + 1 /*right */ ].left;
 183       if (Y1 < discriminant)
 184         deleteSegment (st, n * 2, Y1, Y2);
 185       if (discriminant < Y2)
 186         deleteSegment (st, n * 2 + 1, Y1, Y2);
 187     }
 188   /* fixup area */
 189   st->nodes[n].area = (st->nodes[n].covered > 0) ?
 190     (st->nodes[n].right - st->nodes[n].left) :
 191     (n >= st->size / 2) ? 0 :
 192     st->nodes[n * 2].area + st->nodes[n * 2 + 1].area;
 193 }
 194
 195 /* ---------------------------------------------------------------------------
 196  * Compute the area of the intersection of the given rectangles; that is,
 197  * the area covered by more than one rectangle (counted twice if the area is
 198  * covered by three rectangles, three times if covered by four rectangles,
 199  * etc.).
 200  * Runs in O(N ln N) time.
 201  */
 202 double
 203 ComputeIntersectionArea (BoxListTypePtr boxlist)
 204 {
 205   Cardinal i;
 206   double area = 0.0;
 207   /* first get the aggregate area. */
 208   for (i = 0; i < boxlist->BoxN; i++)
 209     area += (double) (boxlist->Box[i].X2 - boxlist->Box[i].X1) *
 210       (double) (boxlist->Box[i].Y2 - boxlist->Box[i].Y1);
 211   /* intersection area is aggregate - union. */
 212   return area * 0.0001 - ComputeUnionArea (boxlist);
 213 }
 214
 215 /* ---------------------------------------------------------------------------
 216  * Compute the area of the union of the given rectangles.
 217  * O(N ln N) time.
 218  */
 219 double
 220 ComputeUnionArea (BoxListTypePtr boxlist)
 221 {
 222   BoxTypePtr *rectLeft, *rectRight;
 223   Cardinal i, j;
 224   LocationList yCoords;
 225   SegmentTree segtree;
 226   Coord lastX;
 227   double area = 0.0;
 228
 229   if (boxlist->BoxN == 0)
 230     return 0.0;
 231   /* create sorted list of Y coordinates */
 232   yCoords = createSortedYList (boxlist);
 233   /* now create empty segment tree */
 234   segtree = createSegmentTree (yCoords.p, yCoords.size);
 235   free (yCoords.p);
 236   /* create sorted list of left and right X coordinates of rectangles */
 237   rectLeft = (BoxTypePtr *)calloc (boxlist->BoxN, sizeof (*rectLeft));
 238   rectRight = (BoxTypePtr *)calloc (boxlist->BoxN, sizeof (*rectRight));
 239   for (i = 0; i < boxlist->BoxN; i++)
 240     {
 241       assert (boxlist->Box[i].X1 <= boxlist->Box[i].X2);
 242       assert (boxlist->Box[i].Y1 <= boxlist->Box[i].Y2);
 243       rectLeft[i] = rectRight[i] = &boxlist->Box[i];
 244     }
 245   qsort (rectLeft, boxlist->BoxN, sizeof (*rectLeft), compareleft);
 246   qsort (rectRight, boxlist->BoxN, sizeof (*rectRight), compareright);
 247   /* sweep through x segments from left to right */
 248   i = j = 0;
 249   lastX = rectLeft[0]->X1;
 250   while (j < boxlist->BoxN)
 251     {
 252       assert (i <= boxlist->BoxN);
 253       /* i will step through rectLeft, j will through rectRight */
 254       if (i == boxlist->BoxN || rectRight[j]->X2 < rectLeft[i]->X1)
 255         {
 256           /* right edge of rectangle */
 257           BoxTypePtr b = rectRight[j++];
 258           /* check lastX */
 259           if (b->X2 != lastX)
 260             {
 261               assert (lastX < b->X2);
 262               area += (double) (b->X2 - lastX) * segtree.nodes[1].area;
 263               lastX = b->X2;
 264             }
 265           /* remove a segment from the segment tree. */
 266           deleteSegment (&segtree, 1, b->Y1, b->Y2);
 267         }
 268       else
 269         {
 270           /* left edge of rectangle */
 271           BoxTypePtr b = rectLeft[i++];
 272           /* check lastX */
 273           if (b->X1 != lastX)
 274             {
 275               assert (lastX < b->X1);
 276               area += (double) (b->X1 - lastX) * segtree.nodes[1].area;
 277               lastX = b->X1;
 278             }
 279           /* add a segment from the segment tree. */
 280           insertSegment (&segtree, 1, b->Y1, b->Y2);
 281         }
 282     }
 283   free (rectLeft);
 284   free (rectRight);
 285   free (segtree.nodes);
 286   return area * 0.0001;
 287 }
 288 static int
 289 compareleft (const void *ptr1, const void *ptr2)
 290 {
 291   BoxTypePtr *b1 = (BoxTypePtr *) ptr1, *b2 = (BoxTypePtr *) ptr2;
 292   return (*b1)->X1 - (*b2)->X1;
 293 }
 294 static int
 295 compareright (const void *ptr1, const void *ptr2)
 296 {
 297   BoxTypePtr *b1 = (BoxTypePtr *) ptr1, *b2 = (BoxTypePtr *) ptr2;
 298   return (*b1)->X2 - (*b2)->X2;
 299 }
 300 static int
 301 comparepos (const void *ptr1, const void *ptr2)
 302 {
 303   return *((Coord *) ptr1) - *((Coord *) ptr2);
 304 }
 305 static int
 306 nextpwrof2 (int n)
 307 {
 308   int r = 1;
 309   while (n != 0)
 310     {
 311       n /= 2;
 312       r *= 2;
 313     }
 314   return r;
 315 }