mightymandel v16

GPU-based Mandelbrot set explorer

blob_set.c
Go to the documentation of this file.
1 // mightymandel -- GPU-based Mandelbrot Set explorer
2 // Copyright (C) 2012,2013,2014,2015 Claude Heiland-Allen
3 // License GPL3+ http://www.gnu.org/licenses/gpl.html
4 
10 #include <assert.h>
11 #include <math.h>
12 #include <stdio.h>
13 #include <stdint.h>
14 #include <stdlib.h>
15 #include <string.h>
16 
17 #include "blob_set.h"
18 #include "utility.h"
19 
23 struct blob_set_node {
24  struct blob_set_node *next;
25  struct blob blob;
26 };
27 
31 struct blob_set {
32  struct blob_set_node *set;
33 };
34 
36  return (struct blob_set *) calloc(1, sizeof(struct blob_set));
37 }
38 
39 void blob_set_delete(struct blob_set *s) {
40  struct blob_set_node *node = s->set;
41  s->set = 0;
42  while (node) {
43  struct blob_set_node *next = node->next;
44  node->next = 0;
45  free(node);
46  node = next;
47  }
48  free(s);
49 }
50 
51 void blob_set_insert(struct blob_set *s, const struct blob *blob) {
52  struct blob_set_node *node = (struct blob_set_node *) calloc(1, sizeof(struct blob_set_node));
53  node->next = s->set;
54  node->blob.label = 0;
55  node->blob.count = blob->count;
56  node->blob.i = blob->i;
57  node->blob.j = blob->j;
58  s->set = node;
59 }
60 
61 bool blob_set_contains(const struct blob_set *s, const struct blob *blob) {
62  struct blob_set_node *node = s->set;
63  while (node) {
64  if (node->blob.count == blob->count && node->blob.i == blob->i && node->blob.j == blob->j) {
65  return true;
66  }
67  node = node->next;
68  }
69  return false;
70 }
71 
75 int cmp_blob_count_desc(const void *a, const void *b) {
76  const struct blob *x = (const struct blob *) a;
77  const struct blob *y = (const struct blob *) b;
78  return y->count - x->count;
79 }
80 
84 int cmp_blob_error_desc(const void *a, const void *b) {
85  const struct blob *x = (const struct blob *) a;
86  const struct blob *y = (const struct blob *) b;
87  float d = y->error - x->error;
88  if (d > 0) return 1;
89  if (d < 0) return -1;
90  return 0;
91 }
92 
96 struct label {
97  int label;
98  int i;
99  int j;
100  float error;
101 };
102 
106 int cmp_label(const void *a, const void *b) {
107  const struct label *x = (const struct label *) a;
108  const struct label *y = (const struct label *) b;
109  return x->label - y->label;
110 }
111 
115 struct uf {
116  int parent;
117  int rank;
118 };
119 
123 struct ufs {
124  int width;
125  int height;
126  struct uf *nodes;
127  struct label *labels;
128 };
129 
133 void uf_singleton(struct ufs *ufs, int x) {
134  ufs->nodes[x].parent = x;
135  ufs->nodes[x].rank = 0;
136 }
137 
141 int uf_find(struct ufs *ufs, int x) {
142  assert(0 <= x);
143  if (ufs->nodes[x].parent != x) {
144  ufs->nodes[x].parent = uf_find(ufs, ufs->nodes[x].parent);
145  }
146  return ufs->nodes[x].parent;
147 }
148 
152 void uf_union(struct ufs *ufs, int x, int y) {
153  assert(0 <= x);
154  assert(0 <= y);
155  int xroot = uf_find(ufs, x);
156  int yroot = uf_find(ufs, y);
157  assert(0 <= xroot);
158  assert(0 <= yroot);
159  if (xroot == yroot) { return; }
160  if (ufs->nodes[xroot].rank < ufs->nodes[yroot].rank) {
161  ufs->nodes[xroot].parent = yroot;
162  } else if (ufs->nodes[xroot].rank < ufs->nodes[yroot].rank) {
163  ufs->nodes[yroot].parent = xroot;
164  } else {
165  ufs->nodes[yroot].parent = xroot;
166  ufs->nodes[xroot].rank++;
167  }
168 }
169 
170 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) {
171  // union-find to label connected components
172  struct ufs *ufs = (struct ufs *)calloc(1, sizeof(struct ufs));
173  ufs->width = i1 - i0;
174  ufs->height = j1 - j0;
175  int pixels = ufs->width * ufs->height;
176  ufs->nodes = (struct uf *) calloc(1, pixels * sizeof(struct uf));
177  memset(ufs->nodes, -1, pixels * sizeof(struct uf));
178  ufs->labels = (struct label *) calloc(1, pixels * sizeof(struct label));
179  int maximum = pixels;
180  int next = 1;
181  for (int j = j0; j < j1; ++j) {
182  for (int i = i0; i < i1; ++i) {
183  int ki = j * width + i;
184  int ko = (j - j0) * ufs->width + (i - i0);
185  if (glitched[ki] != 0.0f) {
186  int kileft = j * width + i - 1;
187  int koleft = (j - j0) * ufs->width + (i - i0) - 1;
188  int left = maximum;
189  if (i > i0) {
190  if (strategy == blob_boolean ? glitched[kileft] != 0.0f : glitched[kileft] == glitched[ki]) {
191  left = ufs->labels[koleft].label;
192  }
193  }
194  int kiabove = (j - 1) * width + i;
195  int koabove = ((j - j0) - 1) * ufs->width + (i - i0);
196  int above = maximum;
197  if (j > j0) {
198  if (strategy == blob_boolean ? glitched[kiabove] != 0.0f : glitched[kiabove] == glitched[ki]) {
199  above = ufs->labels[koabove].label;
200  }
201  }
202  ufs->labels[ko].i = i;
203  ufs->labels[ko].j = j;
204  ufs->labels[ko].error = glitched[ki];
205  uf_singleton(ufs, ko);
206  if (left == maximum && above == maximum) {
207  ufs->labels[ko].label = next++;
208  } else if (left == maximum) {
209  ufs->labels[ko].label = above;
210  uf_union(ufs, ko, koabove);
211  } else if (above == maximum) {
212  ufs->labels[ko].label = left;
213  uf_union(ufs, ko, koleft);
214  } else {
215  ufs->labels[ko].label = left < above ? left : above;
216  uf_union(ufs, ko, koleft);
217  uf_union(ufs, ko, koabove);
218  }
219  }
220  }
221  }
222  for (int j = j0; j < j1; ++j) {
223  for (int i = i0; i < i1; ++i) {
224  int ki = j * width + i;
225  int ko = (j - j0) * ufs->width + (i - i0);
226  if (glitched[ki] != 0.0f) {
227  ufs->labels[ko].label = uf_find(ufs, ko);
228  }
229  }
230  }
231  // convert labels to blobs
232  struct label *labels = (struct label *) malloc(pixels * sizeof(struct label));
233  memcpy(labels, ufs->labels, pixels * sizeof(struct label));
234  qsort(labels, pixels, sizeof(struct label), cmp_label);
235  int last = 0;
236  int count = 0;
237  int64_t si = 0;
238  int64_t sj = 0;
239  int min_i = width;
240  int max_i = -1;
241  int min_j = height;
242  int max_j = -1;
243  float error = 0.0f;
244  struct blob *blobs = (struct blob *) calloc(pixels, sizeof(struct blob));
245  int index = 0;
246  for (int k = 0; k < pixels; ++k) {
247  if (labels[k].label == 0) {
248  continue;
249  }
250  if (labels[k].label == last) {
251  count++;
252  si += labels[k].i;
253  sj += labels[k].j;
254  min_i = min(min_i, labels[k].i);
255  max_i = max(max_i, labels[k].i);
256  min_j = min(min_j, labels[k].j);
257  max_j = max(max_j, labels[k].j);
258  }
259  if (labels[k].label != last || k == pixels - 1) {
260  if (count > 0) {
261  blobs[index].label = last;
262  blobs[index].count = count;
263  blobs[index].i = si;
264  blobs[index].j = sj;
265  blobs[index].min_i = min_i;
266  blobs[index].max_i = max_i;
267  blobs[index].min_j = min_j;
268  blobs[index].max_j = max_j;
269  blobs[index].error = error;
270  index++;
271  }
272  last = labels[k].label;
273  si = labels[k].i;
274  sj = labels[k].j;
275  min_i = labels[k].i;
276  max_i = labels[k].i;
277  min_j = labels[k].j;
278  max_j = labels[k].j;
279  error = labels[k].error;
280  count = 1;
281  }
282  }
283  free(labels); labels = 0;
284  free(ufs->nodes); ufs->nodes = 0;
285  free(ufs->labels); ufs->labels = 0;
286  free(ufs); ufs = 0;
287  // sort blobs, largest/erroneous first
288  if (index) {
289  qsort(blobs, index, sizeof(struct blob), strategy == blob_boolean ? cmp_blob_count_desc : cmp_blob_error_desc);
290  *blob_count = index;
291  return blobs;
292  } else {
293  free(blobs);
294  *blob_count = 0;
295  return 0;
296  }
297 }
298 
299 struct blob *find_blobs1(int *blob_count, int width, int height, const float *glitched) {
300  return find_blobs0(blob_count, width, height, glitched, blob_boolean, 0, width, 0, height);
301 }
302 
303 struct blob *find_blobs2(int *blob_count, int width, int height, const float *glitched, const struct blob *blob) {
304  return find_blobs0(blob_count, width, height, glitched, blob_positive, blob->min_i, blob->max_i + 1, blob->min_j, blob->max_j + 1);
305 }