partition.c

Go to the documentation of this file.

/*************************************************************************
 * Copyright (c) 2011 AT&T Intellectual Property
 * All rights reserved. This program and the accompanying materials
 * are made available under the terms of the Eclipse Public License v2.0
 * which accompanies this distribution, and is available at
 * https://www.eclipse.org/org/documents/epl-2.0/EPL-2.0.html
 *
 * Contributors: Details at https://graphviz.org
 *************************************************************************/
 
#include "config.h"
#include <common/boxes.h>
#include <common/geomprocs.h>
#include <math.h>
#include <ortho/partition.h>
#include <ortho/trap.h>
#include <stdbool.h>
#include <stdint.h>
#include <stdio.h>
#include <stdlib.h>
#include <util/alloc.h>
#include <util/bitarray.h>
#include <util/gv_math.h>
#include <util/list.h>
#include <util/prisize_t.h>
#include <util/random.h>
 
#ifndef DEBUG
#define DEBUG 0
#endif
 
#define NPOINTS 4 /* only rectangles */
 
#define TR_FROM_UP 1 /* for traverse-direction */
#define TR_FROM_DN 2
 
#define SP_SIMPLE_LRUP 1 /* for splitting trapezoids */
#define SP_SIMPLE_LRDN 2
#define SP_2UP_2DN 3
#define SP_2UP_LEFT 4
#define SP_2UP_RIGHT 5
#define SP_2DN_LEFT 6
#define SP_2DN_RIGHT 7
#define SP_NOSPLIT -1
 
#define DOT(v0, v1) ((v0).x * (v1).x + (v0).y * (v1).y)
#define CROSS_SINE(v0, v1) ((v0).x * (v1).y - (v1).x * (v0).y)
#define LENGTH(v0) hypot((v0).x, (v0).y)
 
typedef struct {
  size_t vnum;
  size_t next;
  size_t prev;
} monchain_t;
 
typedef struct {
  size_t next; 
  size_t pos;  
} vertex_t;
 
typedef struct {
  pointf pt;
  LIST(vertex_t) vertex;
} vertexchain_t;
 
static size_t chain_idx;
static size_t mon_idx;
/* contains position of any vertex in */
/* the monotone chain for the polygon */
static size_t *mon;
 
/* return a new mon structure from the table */
#define newmon() (++mon_idx)
/* return a new chain element from the table */
#define new_chain_element() (++chain_idx)
 
static void convert(boxf bb, int flip, int ccw, pointf *pts) {
  pts[0] = bb.LL;
  pts[2] = bb.UR;
  if (ccw) {
    pts[1].x = bb.UR.x;
    pts[1].y = bb.LL.y;
    pts[3].x = bb.LL.x;
    pts[3].y = bb.UR.y;
  } else {
    pts[1].x = bb.LL.x;
    pts[1].y = bb.UR.y;
    pts[3].x = bb.UR.x;
    pts[3].y = bb.LL.y;
  }
  if (flip) {
    int i;
    for (i = 0; i < NPOINTS; i++) {
      pts[i] = perp(pts[i]);
    }
  }
}
 
static size_t store(segment_t *seg, size_t first, pointf *pts) {
  const size_t last = first + NPOINTS - 1;
 
  for (size_t i = first, j = 0; i <= last; i++, j++) {
    if (i == first) {
      seg[i].next = first + 1;
      seg[i].prev = last;
    } else if (i == last) {
      seg[i].next = first;
      seg[i].prev = last - 1;
    } else {
      seg[i].next = i + 1;
      seg[i].prev = i - 1;
    }
    seg[i].is_inserted = false;
    seg[seg[i].prev].v1 = seg[i].v0 = pts[j];
  }
  return (last + 1);
}
 
static void genSegments(cell *cells, size_t ncells, boxf bb, segment_t *seg,
                        int flip) {
  size_t i = 1;
  pointf pts[4];
 
  convert(bb, flip, 1, pts);
  i = store(seg, i, pts);
  for (size_t j = 0; j < ncells; j++) {
    convert(cells[j].bb, flip, 0, pts);
    i = store(seg, i, pts);
  }
}
 
/* Generate a random permutation of the segments 1..n */
static size_t *generateRandomOrdering(size_t n) {
  size_t *const permute = gv_permutation(n);
 
  // make indexing 1-based
  for (size_t i = 0; i < n; i++) {
    ++permute[i];
  }
 
  return permute;
}
 
/* Function returns true if the trapezoid lies inside the polygon */
static bool inside_polygon(trap_t *t, segment_t *seg) {
  size_t rseg = t->rseg;
 
  if (!t->is_valid)
    return false;
 
  if (t->lseg == 0 || t->rseg == 0)
    return false;
 
  if ((!is_valid_trap(t->u0) && !is_valid_trap(t->u1)) ||
      (!is_valid_trap(t->d0) && !is_valid_trap(t->d1))) // triangle
    return greater_than(seg[rseg].v1, seg[rseg].v0);
 
  return false;
}
 
static double get_angle(pointf vp0, pointf vpnext, pointf vp1) {
  const pointf v0 = sub_pointf(vpnext, vp0);
  const pointf v1 = sub_pointf(vp1, vp0);
 
  if (CROSS_SINE(v0, v1) >= 0) /* sine is positive */
    return DOT(v0, v1) / LENGTH(v0) / LENGTH(v1);
  return -1.0 * DOT(v0, v1) / LENGTH(v0) / LENGTH(v1) - 2;
}
 
static void get_vertex_positions(const vertexchain_t *vert, size_t v0,
                                 size_t v1, size_t *ip, size_t *iq) {
  double angle, temp;
 
  const vertexchain_t *const vp0 = &vert[v0];
  const vertexchain_t *const vp1 = &vert[v1];
 
  /* p is identified as follows. Scan from (v0, v1) rightwards till */
  /* you hit the first segment starting from v0. That chain is the */
  /* chain of our interest */
 
  angle = -4.0;
  *ip = 0;
  for (size_t i = 0; i < LIST_SIZE(&vp0->vertex); i++) {
    const size_t next = LIST_GET(&vp0->vertex, i).next;
    if (next == 0)
      continue;
    if ((temp = get_angle(vp0->pt, vert[next].pt, vp1->pt)) > angle) {
      angle = temp;
      *ip = i;
    }
  }
 
  /* Do similar actions for q */
 
  angle = -4.0;
  *iq = 0;
  for (size_t i = 0; i < LIST_SIZE(&vp1->vertex); i++) {
    const size_t next = LIST_GET(&vp1->vertex, i).next;
    if (next == 0)
      continue;
    if ((temp = get_angle(vp1->pt, vert[next].pt, vp0->pt)) > angle) {
      angle = temp;
      *iq = i;
    }
  }
}
 
/* v0 and v1 are specified in anti-clockwise order with respect to
 * the current monotone polygon mcur. Split the current polygon into
 * two polygons using the diagonal (v0, v1)
 *
 * @param vert Chain initial information
 * @param chain Monotone polygon chain to operate on
 */
static size_t make_new_monotone_poly(vertexchain_t *vert, monchain_t *chain,
                                     size_t mcur, size_t v0, size_t v1) {
  size_t ip, iq;
  const size_t mnew = newmon();
  vertexchain_t *vp0, *vp1;
 
  vp0 = &vert[v0];
  vp1 = &vert[v1];
 
  get_vertex_positions(vert, v0, v1, &ip, &iq);
 
  // if a vertex has no chain yet, synthesize a zeroed element
  if (ip == 0 && LIST_IS_EMPTY(&vp0->vertex)) {
    LIST_APPEND(&vp0->vertex, (vertex_t){0});
  }
  if (iq == 0 && LIST_IS_EMPTY(&vp1->vertex)) {
    LIST_APPEND(&vp1->vertex, (vertex_t){0});
  }
 
  const size_t p = LIST_GET(&vp0->vertex, ip).pos;
  const size_t q = LIST_GET(&vp1->vertex, iq).pos;
 
  /* At this stage, we have got the positions of v0 and v1 in the */
  /* desired chain. Now modify the linked lists */
 
  const size_t i = new_chain_element(); // for the new list
  const size_t j = new_chain_element();
 
  chain[i].vnum = v0;
  chain[j].vnum = v1;
 
  chain[i].next = chain[p].next;
  chain[chain[p].next].prev = i;
  chain[i].prev = j;
  chain[j].next = i;
  chain[j].prev = chain[q].prev;
  chain[chain[q].prev].next = j;
 
  chain[p].next = q;
  chain[q].prev = p;
 
  LIST_AT(&vp0->vertex, ip)->next = v1;
 
  LIST_APPEND(&vp0->vertex,
              (vertex_t){.next = chain[chain[i].next].vnum, .pos = i});
  LIST_APPEND(&vp1->vertex, (vertex_t){.next = v0, .pos = j});
 
#if DEBUG > 0
  fprintf(stderr,
          "make_poly: mcur = %" PRISIZE_T ", (v0, v1) = (%" PRISIZE_T
          ", %" PRISIZE_T ")\n",
          mcur, v0, v1);
  fprintf(stderr, "next posns = (p, q) = (%" PRISIZE_T ", %" PRISIZE_T ")\n", p,
          q);
#endif
 
  mon[mcur] = p;
  mon[mnew] = i;
  return mnew;
}
 
static void traverse_polygon(vertexchain_t *vert, bitarray_t *visited,
                             boxes_t *decomp, segment_t *seg, traps_t *tr,
                             size_t mcur, size_t trnum, size_t from, int flip,
                             int dir, monchain_t *chain) {
  size_t mnew;
 
  if (!is_valid_trap(trnum) || bitarray_get(*visited, trnum))
    return;
 
  trap_t *t = LIST_AT(tr, trnum);
 
  bitarray_set(visited, trnum, true);
 
  if (t->hi.y > t->lo.y + C_EPS &&
      fp_equal(seg[t->lseg].v0.x, seg[t->lseg].v1.x) &&
      fp_equal(seg[t->rseg].v0.x, seg[t->rseg].v1.x)) {
    boxf newbox = {0};
    if (flip) {
      newbox.LL.x = t->lo.y;
      newbox.LL.y = -seg[t->rseg].v0.x;
      newbox.UR.x = t->hi.y;
      newbox.UR.y = -seg[t->lseg].v0.x;
    } else {
      newbox.LL.x = seg[t->lseg].v0.x;
      newbox.LL.y = t->lo.y;
      newbox.UR.x = seg[t->rseg].v0.x;
      newbox.UR.y = t->hi.y;
    }
    LIST_APPEND(decomp, newbox);
  }
 
  /* We have much more information available here. */
  /* rseg: goes upwards   */
  /* lseg: goes downwards */
 
  /* Initially assume that dir = TR_FROM_DN (from the left) */
  /* Switch v0 and v1 if necessary afterwards */
 
  /* special cases for triangles with cusps at the opposite ends. */
  /* take care of this first */
  if (!is_valid_trap(t->u0) && !is_valid_trap(t->u1)) {
    if (is_valid_trap(t->d0) &&
        is_valid_trap(t->d1)) { // downward opening triangle
      const size_t v0 = LIST_GET(tr, t->d1).lseg;
      const size_t v1 = t->lseg;
      if (from == t->d1) {
        mnew = make_new_monotone_poly(vert, chain, mcur, v1, v0);
        traverse_polygon(vert, visited, decomp, seg, tr, mcur, t->d1, trnum,
                         flip, TR_FROM_UP, chain);
        traverse_polygon(vert, visited, decomp, seg, tr, mnew, t->d0, trnum,
                         flip, TR_FROM_UP, chain);
      } else {
        mnew = make_new_monotone_poly(vert, chain, mcur, v0, v1);
        traverse_polygon(vert, visited, decomp, seg, tr, mcur, t->d0, trnum,
                         flip, TR_FROM_UP, chain);
        traverse_polygon(vert, visited, decomp, seg, tr, mnew, t->d1, trnum,
                         flip, TR_FROM_UP, chain);
      }
    } else {
      /* Just traverse all neighbours */
      traverse_polygon(vert, visited, decomp, seg, tr, mcur, t->u0, trnum, flip,
                       TR_FROM_DN, chain);
      traverse_polygon(vert, visited, decomp, seg, tr, mcur, t->u1, trnum, flip,
                       TR_FROM_DN, chain);
      traverse_polygon(vert, visited, decomp, seg, tr, mcur, t->d0, trnum, flip,
                       TR_FROM_UP, chain);
      traverse_polygon(vert, visited, decomp, seg, tr, mcur, t->d1, trnum, flip,
                       TR_FROM_UP, chain);
    }
  }
 
  else if (!is_valid_trap(t->d0) && !is_valid_trap(t->d1)) {
    if (is_valid_trap(t->u0) &&
        is_valid_trap(t->u1)) { // upward opening triangle
      const size_t v0 = t->rseg;
      const size_t v1 = LIST_GET(tr, t->u0).rseg;
      if (from == t->u1) {
        mnew = make_new_monotone_poly(vert, chain, mcur, v1, v0);
        traverse_polygon(vert, visited, decomp, seg, tr, mcur, t->u1, trnum,
                         flip, TR_FROM_DN, chain);
        traverse_polygon(vert, visited, decomp, seg, tr, mnew, t->u0, trnum,
                         flip, TR_FROM_DN, chain);
      } else {
        mnew = make_new_monotone_poly(vert, chain, mcur, v0, v1);
        traverse_polygon(vert, visited, decomp, seg, tr, mcur, t->u0, trnum,
                         flip, TR_FROM_DN, chain);
        traverse_polygon(vert, visited, decomp, seg, tr, mnew, t->u1, trnum,
                         flip, TR_FROM_DN, chain);
      }
    } else {
      /* Just traverse all neighbours */
      traverse_polygon(vert, visited, decomp, seg, tr, mcur, t->u0, trnum, flip,
                       TR_FROM_DN, chain);
      traverse_polygon(vert, visited, decomp, seg, tr, mcur, t->u1, trnum, flip,
                       TR_FROM_DN, chain);
      traverse_polygon(vert, visited, decomp, seg, tr, mcur, t->d0, trnum, flip,
                       TR_FROM_UP, chain);
      traverse_polygon(vert, visited, decomp, seg, tr, mcur, t->d1, trnum, flip,
                       TR_FROM_UP, chain);
    }
  }
 
  else if (is_valid_trap(t->u0) && is_valid_trap(t->u1)) {
    if (is_valid_trap(t->d0) &&
        is_valid_trap(t->d1)) { // downward + upward cusps
      const size_t v0 = LIST_GET(tr, t->d1).lseg;
      const size_t v1 = LIST_GET(tr, t->u0).rseg;
      if ((dir == TR_FROM_DN && t->d1 == from) ||
          (dir == TR_FROM_UP && t->u1 == from)) {
        mnew = make_new_monotone_poly(vert, chain, mcur, v1, v0);
        traverse_polygon(vert, visited, decomp, seg, tr, mcur, t->u1, trnum,
                         flip, TR_FROM_DN, chain);
        traverse_polygon(vert, visited, decomp, seg, tr, mcur, t->d1, trnum,
                         flip, TR_FROM_UP, chain);
        traverse_polygon(vert, visited, decomp, seg, tr, mnew, t->u0, trnum,
                         flip, TR_FROM_DN, chain);
        traverse_polygon(vert, visited, decomp, seg, tr, mnew, t->d0, trnum,
                         flip, TR_FROM_UP, chain);
      } else {
        mnew = make_new_monotone_poly(vert, chain, mcur, v0, v1);
        traverse_polygon(vert, visited, decomp, seg, tr, mcur, t->u0, trnum,
                         flip, TR_FROM_DN, chain);
        traverse_polygon(vert, visited, decomp, seg, tr, mcur, t->d0, trnum,
                         flip, TR_FROM_UP, chain);
        traverse_polygon(vert, visited, decomp, seg, tr, mnew, t->u1, trnum,
                         flip, TR_FROM_DN, chain);
        traverse_polygon(vert, visited, decomp, seg, tr, mnew, t->d1, trnum,
                         flip, TR_FROM_UP, chain);
      }
    } else /* only downward cusp */
    {
      if (equal_to(t->lo, seg[t->lseg].v1)) {
        const size_t v0 = LIST_GET(tr, t->u0).rseg;
        const size_t v1 = seg[t->lseg].next;
 
        if (dir == TR_FROM_UP && t->u0 == from) {
          mnew = make_new_monotone_poly(vert, chain, mcur, v1, v0);
          traverse_polygon(vert, visited, decomp, seg, tr, mcur, t->u0, trnum,
                           flip, TR_FROM_DN, chain);
          traverse_polygon(vert, visited, decomp, seg, tr, mnew, t->d0, trnum,
                           flip, TR_FROM_UP, chain);
          traverse_polygon(vert, visited, decomp, seg, tr, mnew, t->u1, trnum,
                           flip, TR_FROM_DN, chain);
          traverse_polygon(vert, visited, decomp, seg, tr, mnew, t->d1, trnum,
                           flip, TR_FROM_UP, chain);
        } else {
          mnew = make_new_monotone_poly(vert, chain, mcur, v0, v1);
          traverse_polygon(vert, visited, decomp, seg, tr, mcur, t->u1, trnum,
                           flip, TR_FROM_DN, chain);
          traverse_polygon(vert, visited, decomp, seg, tr, mcur, t->d0, trnum,
                           flip, TR_FROM_UP, chain);
          traverse_polygon(vert, visited, decomp, seg, tr, mcur, t->d1, trnum,
                           flip, TR_FROM_UP, chain);
          traverse_polygon(vert, visited, decomp, seg, tr, mnew, t->u0, trnum,
                           flip, TR_FROM_DN, chain);
        }
      } else {
        const size_t v0 = t->rseg;
        const size_t v1 = LIST_GET(tr, t->u0).rseg;
        if (dir == TR_FROM_UP && t->u1 == from) {
          mnew = make_new_monotone_poly(vert, chain, mcur, v1, v0);
          traverse_polygon(vert, visited, decomp, seg, tr, mcur, t->u1, trnum,
                           flip, TR_FROM_DN, chain);
          traverse_polygon(vert, visited, decomp, seg, tr, mnew, t->d1, trnum,
                           flip, TR_FROM_UP, chain);
          traverse_polygon(vert, visited, decomp, seg, tr, mnew, t->d0, trnum,
                           flip, TR_FROM_UP, chain);
          traverse_polygon(vert, visited, decomp, seg, tr, mnew, t->u0, trnum,
                           flip, TR_FROM_DN, chain);
        } else {
          mnew = make_new_monotone_poly(vert, chain, mcur, v0, v1);
          traverse_polygon(vert, visited, decomp, seg, tr, mcur, t->u0, trnum,
                           flip, TR_FROM_DN, chain);
          traverse_polygon(vert, visited, decomp, seg, tr, mcur, t->d0, trnum,
                           flip, TR_FROM_UP, chain);
          traverse_polygon(vert, visited, decomp, seg, tr, mcur, t->d1, trnum,
                           flip, TR_FROM_UP, chain);
          traverse_polygon(vert, visited, decomp, seg, tr, mnew, t->u1, trnum,
                           flip, TR_FROM_DN, chain);
        }
      }
    }
  } else if (is_valid_trap(t->u0) || is_valid_trap(t->u1)) { // no downward cusp
    if (is_valid_trap(t->d0) && is_valid_trap(t->d1)) {      // only upward cusp
      if (equal_to(t->hi, seg[t->lseg].v0)) {
        const size_t v0 = LIST_GET(tr, t->d1).lseg;
        const size_t v1 = t->lseg;
        if (!(dir == TR_FROM_DN && t->d0 == from)) {
          mnew = make_new_monotone_poly(vert, chain, mcur, v1, v0);
          traverse_polygon(vert, visited, decomp, seg, tr, mcur, t->u1, trnum,
                           flip, TR_FROM_DN, chain);
          traverse_polygon(vert, visited, decomp, seg, tr, mcur, t->d1, trnum,
                           flip, TR_FROM_UP, chain);
          traverse_polygon(vert, visited, decomp, seg, tr, mcur, t->u0, trnum,
                           flip, TR_FROM_DN, chain);
          traverse_polygon(vert, visited, decomp, seg, tr, mnew, t->d0, trnum,
                           flip, TR_FROM_UP, chain);
        } else {
          mnew = make_new_monotone_poly(vert, chain, mcur, v0, v1);
          traverse_polygon(vert, visited, decomp, seg, tr, mcur, t->d0, trnum,
                           flip, TR_FROM_UP, chain);
          traverse_polygon(vert, visited, decomp, seg, tr, mnew, t->u0, trnum,
                           flip, TR_FROM_DN, chain);
          traverse_polygon(vert, visited, decomp, seg, tr, mnew, t->u1, trnum,
                           flip, TR_FROM_DN, chain);
          traverse_polygon(vert, visited, decomp, seg, tr, mnew, t->d1, trnum,
                           flip, TR_FROM_UP, chain);
        }
      } else {
        const size_t v0 = LIST_GET(tr, t->d1).lseg;
        const size_t v1 = seg[t->rseg].next;
 
        if (dir == TR_FROM_DN && t->d1 == from) {
          mnew = make_new_monotone_poly(vert, chain, mcur, v1, v0);
          traverse_polygon(vert, visited, decomp, seg, tr, mcur, t->d1, trnum,
                           flip, TR_FROM_UP, chain);
          traverse_polygon(vert, visited, decomp, seg, tr, mnew, t->u1, trnum,
                           flip, TR_FROM_DN, chain);
          traverse_polygon(vert, visited, decomp, seg, tr, mnew, t->u0, trnum,
                           flip, TR_FROM_DN, chain);
          traverse_polygon(vert, visited, decomp, seg, tr, mnew, t->d0, trnum,
                           flip, TR_FROM_UP, chain);
        } else {
          mnew = make_new_monotone_poly(vert, chain, mcur, v0, v1);
          traverse_polygon(vert, visited, decomp, seg, tr, mcur, t->u0, trnum,
                           flip, TR_FROM_DN, chain);
          traverse_polygon(vert, visited, decomp, seg, tr, mcur, t->d0, trnum,
                           flip, TR_FROM_UP, chain);
          traverse_polygon(vert, visited, decomp, seg, tr, mcur, t->u1, trnum,
                           flip, TR_FROM_DN, chain);
          traverse_polygon(vert, visited, decomp, seg, tr, mnew, t->d1, trnum,
                           flip, TR_FROM_UP, chain);
        }
      }
    } else /* no cusp */
    {
      if (equal_to(t->hi, seg[t->lseg].v0) &&
          equal_to(t->lo, seg[t->rseg].v0)) {
        const size_t v0 = t->rseg;
        const size_t v1 = t->lseg;
        if (dir == TR_FROM_UP) {
          mnew = make_new_monotone_poly(vert, chain, mcur, v1, v0);
          traverse_polygon(vert, visited, decomp, seg, tr, mcur, t->u0, trnum,
                           flip, TR_FROM_DN, chain);
          traverse_polygon(vert, visited, decomp, seg, tr, mcur, t->u1, trnum,
                           flip, TR_FROM_DN, chain);
          traverse_polygon(vert, visited, decomp, seg, tr, mnew, t->d1, trnum,
                           flip, TR_FROM_UP, chain);
          traverse_polygon(vert, visited, decomp, seg, tr, mnew, t->d0, trnum,
                           flip, TR_FROM_UP, chain);
        } else {
          mnew = make_new_monotone_poly(vert, chain, mcur, v0, v1);
          traverse_polygon(vert, visited, decomp, seg, tr, mcur, t->d1, trnum,
                           flip, TR_FROM_UP, chain);
          traverse_polygon(vert, visited, decomp, seg, tr, mcur, t->d0, trnum,
                           flip, TR_FROM_UP, chain);
          traverse_polygon(vert, visited, decomp, seg, tr, mnew, t->u0, trnum,
                           flip, TR_FROM_DN, chain);
          traverse_polygon(vert, visited, decomp, seg, tr, mnew, t->u1, trnum,
                           flip, TR_FROM_DN, chain);
        }
      } else if (equal_to(t->hi, seg[t->rseg].v1) &&
                 equal_to(t->lo, seg[t->lseg].v1)) {
        const size_t v0 = seg[t->rseg].next;
        const size_t v1 = seg[t->lseg].next;
 
        if (dir == TR_FROM_UP) {
          mnew = make_new_monotone_poly(vert, chain, mcur, v1, v0);
          traverse_polygon(vert, visited, decomp, seg, tr, mcur, t->u0, trnum,
                           flip, TR_FROM_DN, chain);
          traverse_polygon(vert, visited, decomp, seg, tr, mcur, t->u1, trnum,
                           flip, TR_FROM_DN, chain);
          traverse_polygon(vert, visited, decomp, seg, tr, mnew, t->d1, trnum,
                           flip, TR_FROM_UP, chain);
          traverse_polygon(vert, visited, decomp, seg, tr, mnew, t->d0, trnum,
                           flip, TR_FROM_UP, chain);
        } else {
          mnew = make_new_monotone_poly(vert, chain, mcur, v0, v1);
          traverse_polygon(vert, visited, decomp, seg, tr, mcur, t->d1, trnum,
                           flip, TR_FROM_UP, chain);
          traverse_polygon(vert, visited, decomp, seg, tr, mcur, t->d0, trnum,
                           flip, TR_FROM_UP, chain);
          traverse_polygon(vert, visited, decomp, seg, tr, mnew, t->u0, trnum,
                           flip, TR_FROM_DN, chain);
          traverse_polygon(vert, visited, decomp, seg, tr, mnew, t->u1, trnum,
                           flip, TR_FROM_DN, chain);
        }
      } else /* no split possible */
      {
        traverse_polygon(vert, visited, decomp, seg, tr, mcur, t->u0, trnum,
                         flip, TR_FROM_DN, chain);
        traverse_polygon(vert, visited, decomp, seg, tr, mcur, t->d0, trnum,
                         flip, TR_FROM_UP, chain);
        traverse_polygon(vert, visited, decomp, seg, tr, mcur, t->u1, trnum,
                         flip, TR_FROM_DN, chain);
        traverse_polygon(vert, visited, decomp, seg, tr, mcur, t->d1, trnum,
                         flip, TR_FROM_UP, chain);
      }
    }
  }
}
 
static void monotonate_trapezoids(size_t nsegs, segment_t *seg, traps_t *tr,
                                  int flip, boxes_t *decomp) {
  bitarray_t visited = bitarray_new(LIST_SIZE(tr));
 
  // Table to hold all the monotone polygons. Each monotone polygon is a
  // circularly linked list
  monchain_t *const mchain = gv_calloc(LIST_SIZE(tr), sizeof(monchain_t));
 
  // Chain initial information. This is used to decide which monotone polygon
  // to split if there are several other polygons touching the same vertex.
  vertexchain_t *const vert = gv_calloc(nsegs + 1, sizeof(vertexchain_t));
 
  mon = gv_calloc(nsegs, sizeof(size_t));
 
  /* First locate a trapezoid which lies inside the polygon */
  /* and which is triangular */
  size_t j;
  for (j = 0; j < LIST_SIZE(tr); j++)
    if (inside_polygon(LIST_AT(tr, j), seg))
      break;
  const size_t tr_start = j;
 
  /* Initialise the mon data-structure and start spanning all the */
  /* trapezoids within the polygon */
 
  for (size_t i = 1; i <= nsegs; i++) {
    mchain[i].prev = seg[i].prev;
    mchain[i].next = seg[i].next;
    mchain[i].vnum = i;
    vert[i].pt = seg[i].v0;
    LIST_APPEND(&vert[i].vertex, (vertex_t){
                                     .next = seg[i].next, // next vertex
                                     .pos = i // locn. of next vertex
                                 });
  }
 
  chain_idx = nsegs;
  mon_idx = 0;
  mon[0] = 1; /* position of any vertex in the first */
  /* chain  */
 
  /* traverse the polygon */
  if (is_valid_trap(LIST_GET(tr, tr_start).u0))
    traverse_polygon(vert, &visited, decomp, seg, tr, 0, tr_start,
                     LIST_GET(tr, tr_start).u0, flip, TR_FROM_UP, mchain);
  else if (is_valid_trap(LIST_GET(tr, tr_start).d0))
    traverse_polygon(vert, &visited, decomp, seg, tr, 0, tr_start,
                     LIST_GET(tr, tr_start).d0, flip, TR_FROM_DN, mchain);
 
  bitarray_reset(&visited);
  free(mchain);
  for (size_t i = 0; i < nsegs + 1; ++i) {
    LIST_FREE(&vert[i].vertex);
  }
  free(vert);
  free(mon);
}
 
static bool rectIntersect(boxf *d, const boxf r0, const boxf r1) {
  double t = fmax(r0.LL.x, r1.LL.x);
  d->UR.x = fmin(r0.UR.x, r1.UR.x);
  d->LL.x = t;
 
  t = fmax(r0.LL.y, r1.LL.y);
  d->UR.y = fmin(r0.UR.y, r1.UR.y);
  d->LL.y = t;
 
  return !(d->LL.x >= d->UR.x || d->LL.y >= d->UR.y);
}
 
boxes_t partition(cell *cells, size_t ncells, boxf bb) {
  const size_t nsegs = 4 * (ncells + 1);
  segment_t *segs = gv_calloc(nsegs + 1, sizeof(segment_t));
 
  if (DEBUG) {
    fprintf(stderr,
            "cells = %" PRISIZE_T " segs = %" PRISIZE_T " traps = dynamic\n",
            ncells, nsegs);
  }
  genSegments(cells, ncells, bb, segs, 0);
  if (DEBUG) {
    fprintf(stderr, "%" PRISIZE_T "\n\n", ncells + 1);
    for (size_t i = 1; i <= nsegs; i++) {
      if (i % 4 == 1)
        fprintf(stderr, "4\n");
      fprintf(stderr, "%f %f\n", segs[i].v0.x, segs[i].v0.y);
      if (i % 4 == 0)
        fprintf(stderr, "\n");
    }
  }
  size_t *permute = generateRandomOrdering(nsegs);
  traps_t hor_traps = construct_trapezoids(nsegs, segs, permute);
  free(permute);
  if (DEBUG) {
    fprintf(stderr, "hor traps = %" PRISIZE_T "\n", LIST_SIZE(&hor_traps));
  }
  if (LIST_SIZE(&hor_traps) == 0) {
    free(segs);
    fprintf(stderr, "horizontal trapezoid construction failed\n");
    return (boxes_t){0};
  }
  boxes_t hor_decomp = {0};
  monotonate_trapezoids(nsegs, segs, &hor_traps, 0, &hor_decomp);
  LIST_FREE(&hor_traps);
 
  genSegments(cells, ncells, bb, segs, 1);
  permute = generateRandomOrdering(nsegs);
  traps_t ver_traps = construct_trapezoids(nsegs, segs, permute);
  free(permute);
  if (DEBUG) {
    fprintf(stderr, "ver traps = %" PRISIZE_T "\n", LIST_SIZE(&ver_traps));
  }
  if (LIST_SIZE(&ver_traps) == 0) {
    LIST_FREE(&hor_decomp);
    free(segs);
    fprintf(stderr, "vertical trapezoid construction failed\n");
    return (boxes_t){0};
  }
  boxes_t vert_decomp = {0};
  monotonate_trapezoids(nsegs, segs, &ver_traps, 1, &vert_decomp);
  LIST_FREE(&ver_traps);
 
  boxes_t rs = {0};
  for (size_t i = 0; i < LIST_SIZE(&vert_decomp); ++i)
    for (size_t j = 0; j < LIST_SIZE(&hor_decomp); ++j) {
      boxf newbox = {0};
      if (rectIntersect(&newbox, LIST_GET(&vert_decomp, i),
                        LIST_GET(&hor_decomp, j)))
        LIST_APPEND(&rs, newbox);
    }
 
  free(segs);
  LIST_FREE(&hor_decomp);
  LIST_FREE(&vert_decomp);
  return rs;
}