blob_set.c

Go to the documentation of this file.

// mightymandel -- GPU-based Mandelbrot Set explorer
// Copyright (C) 2012,2013,2014,2015 Claude Heiland-Allen
// License GPL3+ http://www.gnu.org/licenses/gpl.html

#include <assert.h>
#include <math.h>
#include <stdio.h>
#include <stdint.h>
#include <stdlib.h>
#include <string.h>

#include "blob_set.h"
#include "utility.h"

struct blob_set_node {
  struct blob_set_node *next; 
  struct blob blob; 
};

struct blob_set {
  struct blob_set_node *set; 
};

struct blob_set *blob_set_new() {
  return (struct blob_set *) calloc(1, sizeof(struct blob_set));
}

void blob_set_delete(struct blob_set *s) {
  struct blob_set_node *node = s->set;
  s->set = 0;
  while (node) {
    struct blob_set_node *next = node->next;
    node->next = 0;
    free(node);
    node = next;
  }
  free(s);
}

void blob_set_insert(struct blob_set *s, const struct blob *blob) {
  struct blob_set_node *node = (struct blob_set_node *) calloc(1, sizeof(struct blob_set_node));
  node->next = s->set;
  node->blob.label = 0;
  node->blob.count = blob->count;
  node->blob.i = blob->i;
  node->blob.j = blob->j;
  s->set = node;
}

bool blob_set_contains(const struct blob_set *s, const struct blob *blob) {
  struct blob_set_node *node = s->set;
  while (node) {
    if (node->blob.count == blob->count && node->blob.i == blob->i && node->blob.j == blob->j) {
      return true;
    }
    node = node->next;
  }
  return false;
}

int cmp_blob_count_desc(const void *a, const void *b) {
  const struct blob *x = (const struct blob *) a;
  const struct blob *y = (const struct blob *) b;
  return y->count - x->count;
}

int cmp_blob_error_desc(const void *a, const void *b) {
  const struct blob *x = (const struct blob *) a;
  const struct blob *y = (const struct blob *) b;
  float d = y->error - x->error;
  if (d > 0) return  1;
  if (d < 0) return -1;
  return 0;
}

struct label {
  int label;
  int i;
  int j;
  float error;
};

int cmp_label(const void *a, const void *b) {
  const struct label *x = (const struct label *) a;
  const struct label *y = (const struct label *) b;
  return x->label - y->label;
}

struct uf {
  int parent;
  int rank;
};

struct ufs {
  int width;
  int height;
  struct uf *nodes;
  struct label *labels;
};

void uf_singleton(struct ufs *ufs, int x) {
  ufs->nodes[x].parent = x;
  ufs->nodes[x].rank = 0;
}

int uf_find(struct ufs *ufs, int x) {
  assert(0 <= x);
  if (ufs->nodes[x].parent != x) {
    ufs->nodes[x].parent = uf_find(ufs, ufs->nodes[x].parent);
  }
  return ufs->nodes[x].parent;
}

void uf_union(struct ufs *ufs, int x, int y) {
  assert(0 <= x);
  assert(0 <= y);
  int xroot = uf_find(ufs, x);
  int yroot = uf_find(ufs, y);
  assert(0 <= xroot);
  assert(0 <= yroot);
  if (xroot == yroot) { return; }
  if (ufs->nodes[xroot].rank < ufs->nodes[yroot].rank) {
    ufs->nodes[xroot].parent = yroot;
  } else if (ufs->nodes[xroot].rank < ufs->nodes[yroot].rank) {
    ufs->nodes[yroot].parent = xroot;
  } else {
    ufs->nodes[yroot].parent = xroot;
    ufs->nodes[xroot].rank++;
  }
}

struct blob *find_blobs0(int *blob_count, int width, int height, const float *glitched, enum blob_strategy strategy, int i0, int i1, int j0, int j1) {
  // union-find to label connected components
  struct ufs *ufs = (struct ufs *)calloc(1, sizeof(struct ufs));
  ufs->width = i1 - i0;
  ufs->height = j1 - j0;
  int pixels = ufs->width * ufs->height;
  ufs->nodes = (struct uf *) calloc(1, pixels * sizeof(struct uf));
  memset(ufs->nodes, -1, pixels * sizeof(struct uf));
  ufs->labels = (struct label *) calloc(1, pixels * sizeof(struct label));
  int maximum = pixels;
  int next = 1;
  for (int j = j0; j < j1; ++j) {
    for (int i = i0; i < i1; ++i) {
      int ki = j * width + i;
      int ko = (j - j0) * ufs->width + (i - i0);
      if (glitched[ki] != 0.0f) {
        int kileft = j * width + i - 1;
        int koleft = (j - j0) * ufs->width + (i - i0) - 1;
        int left = maximum;
        if (i > i0) {
          if (strategy == blob_boolean ? glitched[kileft] != 0.0f : glitched[kileft] == glitched[ki]) {
            left = ufs->labels[koleft].label;
          }
        }
        int kiabove = (j - 1) * width + i;
        int koabove = ((j - j0) - 1) * ufs->width + (i - i0);
        int above = maximum;
        if (j > j0) {
          if (strategy == blob_boolean ? glitched[kiabove] != 0.0f : glitched[kiabove] == glitched[ki]) {
            above = ufs->labels[koabove].label;
          }
        }
        ufs->labels[ko].i = i;
        ufs->labels[ko].j = j;
        ufs->labels[ko].error = glitched[ki];
        uf_singleton(ufs, ko);
        if (left == maximum && above == maximum) {
          ufs->labels[ko].label = next++;
        } else if (left == maximum) {
          ufs->labels[ko].label = above;
          uf_union(ufs, ko, koabove);
        } else if (above == maximum) {
          ufs->labels[ko].label = left;
          uf_union(ufs, ko, koleft);
        } else {
          ufs->labels[ko].label = left < above ? left : above;
          uf_union(ufs, ko, koleft);
          uf_union(ufs, ko, koabove);
        }
      }
    }
  }
  for (int j = j0; j < j1; ++j) {
    for (int i = i0; i < i1; ++i) {
      int ki = j * width + i;
      int ko = (j - j0) * ufs->width + (i - i0);
      if (glitched[ki] != 0.0f) {
        ufs->labels[ko].label = uf_find(ufs, ko);
      }
    }
  }
  // convert labels to blobs
  struct label *labels = (struct label *) malloc(pixels * sizeof(struct label));
  memcpy(labels, ufs->labels, pixels * sizeof(struct label));
  qsort(labels, pixels, sizeof(struct label), cmp_label);
  int last = 0;
  int count = 0;
  int64_t si = 0;
  int64_t sj = 0;
  int min_i = width;
  int max_i = -1;
  int min_j = height;
  int max_j = -1;
  float error = 0.0f;
  struct blob *blobs = (struct blob *) calloc(pixels, sizeof(struct blob));
  int index = 0;
  for (int k = 0; k < pixels; ++k) {
    if (labels[k].label == 0) {
      continue;
    }
    if (labels[k].label == last) {
      count++;
      si += labels[k].i;
      sj += labels[k].j;
      min_i = min(min_i, labels[k].i);
      max_i = max(max_i, labels[k].i);
      min_j = min(min_j, labels[k].j);
      max_j = max(max_j, labels[k].j);
    }
    if (labels[k].label != last || k == pixels - 1) {
      if (count > 0) {
        blobs[index].label = last;
        blobs[index].count = count;
        blobs[index].i = si;
        blobs[index].j = sj;
        blobs[index].min_i = min_i;
        blobs[index].max_i = max_i;
        blobs[index].min_j = min_j;
        blobs[index].max_j = max_j;
        blobs[index].error = error;
        index++;
      }
      last = labels[k].label;
      si = labels[k].i;
      sj = labels[k].j;
      min_i = labels[k].i;
      max_i = labels[k].i;
      min_j = labels[k].j;
      max_j = labels[k].j;
      error = labels[k].error;
      count = 1;
    }
  }
  free(labels); labels = 0;
  free(ufs->nodes); ufs->nodes = 0;
  free(ufs->labels); ufs->labels = 0;
  free(ufs); ufs = 0;
  // sort blobs, largest/erroneous first
  if (index) {
    qsort(blobs, index, sizeof(struct blob), strategy == blob_boolean ? cmp_blob_count_desc : cmp_blob_error_desc);
    *blob_count = index;
    return blobs;
  } else {
    free(blobs);
    *blob_count = 0;
    return 0;
  }
}

struct blob *find_blobs1(int *blob_count, int width, int height, const float *glitched) {
  return find_blobs0(blob_count, width, height, glitched, blob_boolean, 0, width, 0, height);
}

struct blob *find_blobs2(int *blob_count, int width, int height, const float *glitched, const struct blob *blob) {
  return find_blobs0(blob_count, width, height, glitched, blob_positive, blob->min_i, blob->max_i + 1, blob->min_j, blob->max_j + 1);
}