box.c

#include <math.h>
#include <stdlib.h>
#include <stdio.h>

#include "box.h"

// --------------------------------------------- Helper function declarations

int min(int a, int b);
int max(int a, int b);
double random_split_ratio(); 

void box_heap_insert(box_t *heap, int *heap_size, const box_t *box);

void box_heap_take_first(box_t *heap, int *heap_size, box_t *first);

// ==========================================================================

void box_set_volume(box_t *box)
{
   int *extents = box->extents;
   long volume = 1;
   int i;
   for(i = 0; i < N_DIMS; i++) {
      volume *= extents[i];
   }
   box->volume = volume;
}

void box_intersect(const box_t *a, const box_t *b, box_t *result)
{
   const int *ac = a->coords, *ae = a->extents;
   const int *bc = b->coords, *be = b->extents;

   int i, coordinate, extent;

   for(i = 0; i < N_DIMS; i++) {
      coordinate = max(ac[i], bc[i]);
      extent = min(ac[i] + ae[i], bc[i] + be[i]) - coordinate;

      result->coords[i] = coordinate;
      result->extents[i] = extent;
   }
   box_set_volume(result);
}

int box_is_empty(const box_t *b)
{
   const int *extents = b->extents;
   return extents[X] <= 0 || extents[Y] <= 0 || extents[Z] <= 0;
}

void box_split(const box_t *original, box_t *a, box_t *b, int randomize)
{
   const int *extents = original->extents;

   static int split_dim = 0;
   int i, j;
   int max_extent = -1;

   for(i = 0; i < N_DIMS; i++) {
      j = (split_dim + i) % N_DIMS;
      if(max_extent < extents[j]) {
         max_extent = extents[j];
         split_dim  = j;
      }
   }

   double split_ratio = randomize ? random_split_ratio() : 0.5;
   int ea = (int)round(max_extent * split_ratio);
   int eb = max_extent - ea;

   *a = *original;
   *b = *original;

   a->extents[split_dim] = ea;
   b->coords[split_dim] += ea;
   b->extents[split_dim] = eb;

   box_set_volume(a);
   box_set_volume(b);
}

void box_decompose(const box_t *original, int n_boxes, box_t *boxes)
{
   box_t largest_box, part1, part2;
   int n, randomize;

   boxes[0] = *original;
   n = 1;
   while(n <  n_boxes) {
      randomize = n < n_boxes/2;

      // take out box with maximum volume and split in two
      box_heap_take_first(boxes, &n, &largest_box);
      box_split(&largest_box, &part1, &part2, randomize);

      // insert boxes
      box_heap_insert(boxes, &n, &part1);
      box_heap_insert(boxes, &n, &part2);
   }
}

void box_grow(box_t *box, int extra_width)
{
   int *coords = box->coords;
   int *extents = box->extents;
   int i;
   for(i = 0; i < N_DIMS; i++) {
      coords[i] -= extra_width;
      extents[i] += 2 * extra_width;
   }
   box_set_volume(box);
}

// --------------------------------------------------------- Helper functions

int min(int a, int b)
{
   return b < a ? b : a;
}

int max(int a, int b)
{
   return b > a ? b : a;
}

double random_split_ratio()
{
   const double jitter = 0.085;

   double ratio = 0.5 - jitter + 2 * jitter * random() / RAND_MAX;
   return ratio;
}

void box_heap_insert(box_t *heap, int *heap_size, const box_t *box)
{
   int current = *heap_size;
   int parent = (current - 1) / 2;
   box_t tmp = *box;

   while(current > 0 && heap[parent].volume < tmp.volume) {
      heap[current] = heap[parent];
      current = parent;
      parent = (current - 1) / 2;
   }
   heap[current] = tmp;
   (*heap_size)++;
}

void box_heap_take_first(box_t *heap, int *heap_size, box_t *first)
{
   int current = 0;
   int child = current * 2 + 1;
   int n = *heap_size - 1;
   box_t *tail = &heap[n];

   *first = heap[current];
   *heap_size = n;
   while(child < n) {
      if(child+1 < n && heap[child+1].volume > heap[child].volume) child++;
      if(tail->volume >= heap[child].volume) break;

      heap[current] = heap[child];
      current = child;
      child = current * 2 + 1;
   }
   heap[current] = *tail;
}

char *box_to_string(const box_t *box, char *buf, int size)
{
   const int *c = box->coords;
   const int *e = box->extents;
   snprintf(buf, size, "x: [%5d,%5d); y: [%5d,%5d); z: [%5d,%5d); cells: %ld",
            c[X], c[X]+e[X], c[Y], c[Y]+e[Y], c[Z], c[Z]+e[Z], box->volume);
   return buf;
}