Cleanup: remove redundant doxygen \file argument
[blender.git] / source / blender / blenlib / intern / voronoi_2d.c
1 /*
2  * This program is free software; you can redistribute it and/or
3  * modify it under the terms of the GNU General Public License
4  * as published by the Free Software Foundation; either version 2
5  * of the License, or (at your option) any later version.
6  *
7  * This program is distributed in the hope that it will be useful,
8  * but WITHOUT ANY WARRANTY; without even the implied warranty of
9  * MERCHANTABILITY or FITNESS FOR A PARTICULAR PURPOSE.  See the
10  * GNU General Public License for more details.
11  *
12  * You should have received a copy of the GNU General Public License
13  * along with this program; if not, write to the Free Software Foundation,
14  * Inc., 51 Franklin Street, Fifth Floor, Boston, MA 02110-1301, USA.
15  *
16  * The Original Code is Copyright (C) 2012 Blender Foundation.
17  * All rights reserved.
18  */
19
20 /** \file \ingroup bli
21  *
22  * Fortune's algorithm implemented using explanation and some code snippets from
23  * http://blog.ivank.net/fortunes-algorithm-and-implementation.html
24  */
25
26 #include "MEM_guardedalloc.h"
27
28 #include "BLI_listbase.h"
29 #include "BLI_math.h"
30 #include "BLI_voronoi_2d.h"
31 #include "BLI_utildefines.h"
32
33 #define VORONOI_EPS 1e-2f
34
35 enum {
36         voronoiEventType_Site = 0,
37         voronoiEventType_Circle = 1,
38 };
39
40 typedef struct VoronoiEvent {
41         struct VoronoiEvent *next, *prev;
42
43         int type;       /* type of event (site or circle) */
44         float site[2];  /* site for which event was generated */
45
46         struct VoronoiParabola *parabola;   /* parabola for which event was generated */
47 } VoronoiEvent;
48
49 typedef struct VoronoiParabola {
50         struct VoronoiParabola *left, *right, *parent;
51         VoronoiEvent *event;
52         VoronoiEdge *edge;
53         float site[2];
54         bool is_leaf;
55 } VoronoiParabola;
56
57 typedef struct VoronoiProcess {
58         ListBase queue, edges;
59         VoronoiParabola *root;
60         int width, height;
61         float current_y;
62 } VoronoiProcess;
63
64 /* event */
65
66 static void voronoi_insertEvent(VoronoiProcess *process, VoronoiEvent *event)
67 {
68         VoronoiEvent *current_event = process->queue.first;
69
70         while (current_event) {
71                 if (current_event->site[1] < event->site[1]) {
72                         break;
73                 }
74                 if (current_event->site[1] == event->site[1]) {
75                         event->site[1] -= VORONOI_EPS;
76                 }
77
78                 current_event = current_event->next;
79         }
80
81         BLI_insertlinkbefore(&process->queue, current_event, event);
82 }
83
84 /* edge */
85 static VoronoiEdge *voronoiEdge_new(float start[2], float left[2], float right[2])
86 {
87         VoronoiEdge *edge = MEM_callocN(sizeof(VoronoiEdge), "voronoi edge");
88
89         copy_v2_v2(edge->start, start);
90         copy_v2_v2(edge->left, left);
91         copy_v2_v2(edge->right, right);
92
93         edge->neighbor = NULL;
94         edge->end[0] = 0;
95         edge->end[1] = 0;
96
97         edge->f = (right[0] - left[0]) / (left[1] - right[1]);
98         edge->g = start[1] - edge->f * start[0];
99
100         edge->direction[0] = right[1] - left[1];
101         edge->direction[1] = -(right[0] - left[0]);
102
103         return edge;
104 }
105
106 /* parabola */
107
108 static VoronoiParabola *voronoiParabola_new(void)
109 {
110         VoronoiParabola *parabola = MEM_callocN(sizeof(VoronoiParabola), "voronoi parabola");
111
112         parabola->is_leaf = false;
113         parabola->event = NULL;
114         parabola->edge = NULL;
115         parabola->parent = NULL;
116
117         return parabola;
118 }
119
120 static VoronoiParabola *voronoiParabola_newSite(float site[2])
121 {
122         VoronoiParabola *parabola = MEM_callocN(sizeof(VoronoiParabola), "voronoi parabola site");
123
124         copy_v2_v2(parabola->site, site);
125         parabola->is_leaf = true;
126         parabola->event = NULL;
127         parabola->edge = NULL;
128         parabola->parent = NULL;
129
130         return parabola;
131 }
132
133 /* returns the closest leave which is on the left of current node */
134 static VoronoiParabola *voronoiParabola_getLeftChild(VoronoiParabola *parabola)
135 {
136         VoronoiParabola *current_parabola;
137
138         if (!parabola)
139                 return NULL;
140
141         current_parabola = parabola->left;
142         while (!current_parabola->is_leaf) {
143                 current_parabola = current_parabola->right;
144         }
145
146         return current_parabola;
147 }
148
149 /* returns the closest leave which is on the right of current node */
150 static VoronoiParabola *voronoiParabola_getRightChild(VoronoiParabola *parabola)
151 {
152         VoronoiParabola *current_parabola;
153
154         if (!parabola)
155                 return NULL;
156
157         current_parabola = parabola->right;
158         while (!current_parabola->is_leaf) {
159                 current_parabola = current_parabola->left;
160         }
161
162         return current_parabola;
163 }
164
165 /* returns the closest parent which is on the left */
166 static VoronoiParabola *voronoiParabola_getLeftParent(VoronoiParabola *parabola)
167 {
168         VoronoiParabola *current_par = parabola->parent;
169         VoronoiParabola *last_parabola = parabola;
170
171         while (current_par->left == last_parabola) {
172                 if (!current_par->parent)
173                         return NULL;
174
175                 last_parabola = current_par;
176                 current_par = current_par->parent;
177         }
178
179         return current_par;
180 }
181
182 /* returns the closest parent which is on the right */
183 static VoronoiParabola *voronoiParabola_getRightParent(VoronoiParabola *parabola)
184 {
185         VoronoiParabola *current_parabola = parabola->parent;
186         VoronoiParabola *last_parabola = parabola;
187
188         while (current_parabola->right == last_parabola) {
189                 if (!current_parabola->parent)
190                         return NULL;
191
192                 last_parabola = current_parabola;
193                 current_parabola = current_parabola->parent;
194         }
195
196         return current_parabola;
197 }
198
199 static void voronoiParabola_setLeft(VoronoiParabola *parabola, VoronoiParabola *left)
200 {
201         parabola->left = left;
202         left->parent = parabola;
203 }
204
205 static void voronoiParabola_setRight(VoronoiParabola *parabola, VoronoiParabola *right)
206 {
207         parabola->right = right;
208         right->parent = parabola;
209 }
210
211 static float voronoi_getY(VoronoiProcess *process, float p[2], float x)
212 {
213         float ly = process->current_y;
214
215         float dp = 2 * (p[1] - ly);
216         float a1 = 1 / dp;
217         float b1 = -2 * p[0] / dp;
218         float c1 = ly + dp / 4 + p[0] * p[0] / dp;
219
220         return a1 * x * x + b1 * x + c1;
221 }
222
223 static float voronoi_getXOfEdge(VoronoiProcess *process, VoronoiParabola *par, float y)
224 {
225         VoronoiParabola *left = voronoiParabola_getLeftChild(par);
226         VoronoiParabola *right = voronoiParabola_getRightChild(par);
227         float p[2], r[2];
228         float dp, a1, b1, c1, a2, b2, c2, a, b, c, disc, ry, x1, x2;
229         float ly = process->current_y;
230
231         copy_v2_v2(p, left->site);
232         copy_v2_v2(r, right->site);
233
234         dp = 2.0f * (p[1] - y);
235         a1 = 1.0f / dp;
236         b1 = -2.0f * p[0] / dp;
237         c1 = y + dp / 4 + p[0] * p[0] / dp;
238
239         dp = 2.0f * (r[1] - y);
240         a2 = 1.0f / dp;
241         b2 = -2.0f * r[0] / dp;
242         c2 = ly + dp / 4 + r[0] * r[0] / dp;
243
244         a = a1 - a2;
245         b = b1 - b2;
246         c = c1 - c2;
247
248         disc = b * b - 4 * a * c;
249         x1 = (-b + sqrtf(disc)) / (2 * a);
250         x2 = (-b - sqrtf(disc)) / (2 * a);
251
252         if (p[1] < r[1])
253                 ry = max_ff(x1, x2);
254         else
255                 ry = min_ff(x1, x2);
256
257         return ry;
258 }
259
260 static VoronoiParabola *voronoi_getParabolaByX(VoronoiProcess *process, float xx)
261 {
262         VoronoiParabola *par = process->root;
263         float x = 0.0f;
264         float ly = process->current_y;
265
266         while (!par->is_leaf) {
267                 x = voronoi_getXOfEdge(process, par, ly);
268
269                 if (x > xx)
270                         par = par->left;
271                 else
272                         par = par->right;
273         }
274
275         return par;
276 }
277
278 static int voronoi_getEdgeIntersection(VoronoiEdge *a, VoronoiEdge *b, float p[2])
279 {
280         float x = (b->g - a->g) / (a->f - b->f);
281         float y = a->f * x + a->g;
282
283         if ((x - a->start[0]) / a->direction[0] < 0)
284                 return 0;
285
286         if ((y - a->start[1]) / a->direction[1] < 0)
287                 return 0;
288
289         if ((x - b->start[0]) / b->direction[0] < 0)
290                 return 0;
291
292         if ((y - b->start[1]) / b->direction[1] < 0)
293                 return 0;
294
295         p[0] = x;
296         p[1] = y;
297
298         return 1;
299 }
300
301 static void voronoi_checkCircle(VoronoiProcess *process, VoronoiParabola *b)
302 {
303         VoronoiParabola *lp = voronoiParabola_getLeftParent(b);
304         VoronoiParabola *rp = voronoiParabola_getRightParent(b);
305
306         VoronoiParabola *a  = voronoiParabola_getLeftChild(lp);
307         VoronoiParabola *c  = voronoiParabola_getRightChild(rp);
308
309         VoronoiEvent *event;
310
311         float ly = process->current_y;
312         float s[2], dx, dy, d;
313
314         if (!a || !c || len_squared_v2v2(a->site, c->site) < VORONOI_EPS)
315                 return;
316
317         if (!voronoi_getEdgeIntersection(lp->edge, rp->edge, s))
318                 return;
319
320         dx = a->site[0] - s[0];
321         dy = a->site[1] - s[1];
322
323         d = sqrtf((dx * dx) + (dy * dy));
324
325         if (s[1] - d >= ly)
326                 return;
327
328         event = MEM_callocN(sizeof(VoronoiEvent), "voronoi circle event");
329
330         event->type = voronoiEventType_Circle;
331
332         event->site[0] = s[0];
333         event->site[1] = s[1] - d;
334
335         b->event = event;
336         event->parabola = b;
337
338         voronoi_insertEvent(process, event);
339 }
340
341 static void voronoi_addParabola(VoronoiProcess *process, float site[2])
342 {
343         VoronoiParabola *root = process->root;
344         VoronoiParabola *par, *p0, *p1, *p2;
345         VoronoiEdge *el, *er;
346         float start[2];
347
348         if (!process->root) {
349                 process->root = voronoiParabola_newSite(site);
350
351                 return;
352         }
353
354         if (root->is_leaf && root->site[1] - site[1] < 0) {
355                 float *fp = root->site;
356                 float s[2];
357
358                 root->is_leaf = false;
359                 voronoiParabola_setLeft(root, voronoiParabola_newSite(fp));
360                 voronoiParabola_setRight(root, voronoiParabola_newSite(site));
361
362                 s[0] = (site[0] + fp[0]) / 2.0f;
363                 s[1] = process->height;
364
365                 if (site[0] > fp[0])
366                         root->edge = voronoiEdge_new(s, fp, site);
367                 else
368                         root->edge = voronoiEdge_new(s, site, fp);
369
370                 BLI_addtail(&process->edges, root->edge);
371
372                 return;
373         }
374
375         par = voronoi_getParabolaByX(process, site[0]);
376
377         if (par->event) {
378                 BLI_freelinkN(&process->queue, par->event);
379
380                 par->event = NULL;
381         }
382
383         start[0] = site[0];
384         start[1] = voronoi_getY(process, par->site, site[0]);
385
386         el = voronoiEdge_new(start, par->site, site);
387         er = voronoiEdge_new(start, site, par->site);
388
389         el->neighbor = er;
390         BLI_addtail(&process->edges, el);
391
392         par->edge = er;
393         par->is_leaf = false;
394
395         p0 = voronoiParabola_newSite(par->site);
396         p1 = voronoiParabola_newSite(site);
397         p2 = voronoiParabola_newSite(par->site);
398
399         voronoiParabola_setRight(par, p2);
400         voronoiParabola_setLeft(par, voronoiParabola_new());
401         par->left->edge = el;
402
403         voronoiParabola_setLeft(par->left, p0);
404         voronoiParabola_setRight(par->left, p1);
405
406         voronoi_checkCircle(process, p0);
407         voronoi_checkCircle(process, p2);
408 }
409
410 static void voronoi_removeParabola(VoronoiProcess *process, VoronoiEvent *event)
411 {
412         VoronoiParabola *p1 = event->parabola;
413
414         VoronoiParabola *xl = voronoiParabola_getLeftParent(p1);
415         VoronoiParabola *xr = voronoiParabola_getRightParent(p1);
416
417         VoronoiParabola *p0 = voronoiParabola_getLeftChild(xl);
418         VoronoiParabola *p2 = voronoiParabola_getRightChild(xr);
419
420         VoronoiParabola *higher = NULL, *par, *gparent;
421
422         float p[2];
423
424         if (p0->event) {
425                 BLI_freelinkN(&process->queue, p0->event);
426                 p0->event = NULL;
427         }
428
429         if (p2->event) {
430                 BLI_freelinkN(&process->queue, p2->event);
431                 p2->event = NULL;
432         }
433
434         p[0] = event->site[0];
435         p[1] = voronoi_getY(process, p1->site, event->site[0]);
436
437         copy_v2_v2(xl->edge->end, p);
438         copy_v2_v2(xr->edge->end, p);
439
440         par = p1;
441         while (par != process->root) {
442                 par = par->parent;
443
444                 if (par == xl)
445                         higher = xl;
446                 if (par == xr)
447                         higher = xr;
448         }
449
450         higher->edge = voronoiEdge_new(p, p0->site, p2->site);
451         BLI_addtail(&process->edges, higher->edge);
452
453         gparent = p1->parent->parent;
454         if (p1->parent->left == p1) {
455                 if (gparent->left == p1->parent)
456                         voronoiParabola_setLeft(gparent, p1->parent->right);
457                 if (gparent->right == p1->parent)
458                         voronoiParabola_setRight(gparent, p1->parent->right);
459         }
460         else {
461                 if (gparent->left == p1->parent)
462                         voronoiParabola_setLeft(gparent, p1->parent->left);
463                 if (gparent->right == p1->parent)
464                         voronoiParabola_setRight(gparent, p1->parent->left);
465         }
466
467         MEM_freeN(p1->parent);
468         MEM_freeN(p1);
469
470         voronoi_checkCircle(process, p0);
471         voronoi_checkCircle(process, p2);
472 }
473
474 static void voronoi_finishEdge(VoronoiProcess *process, VoronoiParabola *parabola)
475 {
476         float mx;
477
478         if (parabola->is_leaf) {
479                 MEM_freeN(parabola);
480                 return;
481         }
482
483         if (parabola->edge->direction[0] > 0.0f)
484                 mx = max_ff(process->width, parabola->edge->start[0] + 10);
485         else
486                 mx = min_ff(0.0f, parabola->edge->start[0] - 10.0f);
487
488         parabola->edge->end[0] = mx;
489         parabola->edge->end[1] = mx * parabola->edge->f + parabola->edge->g;
490
491         voronoi_finishEdge(process, parabola->left);
492         voronoi_finishEdge(process, parabola->right);
493
494         MEM_freeN(parabola);
495 }
496
497 static void voronoi_clampEdgeVertex(int width, int height, float *coord, float *other_coord)
498 {
499         const float corners[4][2] = {{0.0f, 0.0f},
500                                      {width - 1, 0.0f},
501                                      {width - 1, height - 1},
502                                      {0.0f, height - 1}};
503         int i;
504
505         if (IN_RANGE_INCL(coord[0], 0, width - 1) && IN_RANGE_INCL(coord[1], 0, height - 1)) {
506                 return;
507         }
508
509         for (i = 0; i < 4; i++) {
510                 float v1[2], v2[2];
511                 float p[2];
512
513                 copy_v2_v2(v1, corners[i]);
514
515                 if (i == 3)
516                         copy_v2_v2(v2, corners[0]);
517                 else
518                         copy_v2_v2(v2, corners[i + 1]);
519
520                 if (isect_seg_seg_v2_point(v1, v2, coord, other_coord, p) == 1) {
521                         if (i == 0 && coord[1] > p[1])
522                                 continue;
523                         if (i == 1 && coord[0] < p[0])
524                                 continue;
525                         if (i == 2 && coord[1] < p[1])
526                                 continue;
527                         if (i == 3 && coord[0] > p[0])
528                                 continue;
529
530                         copy_v2_v2(coord, p);
531                 }
532         }
533 }
534
535 static void voronoi_clampEdges(ListBase *edges, int width, int height, ListBase *clamped_edges)
536 {
537         VoronoiEdge *edge;
538
539         edge = edges->first;
540         while (edge) {
541                 VoronoiEdge *new_edge = MEM_callocN(sizeof(VoronoiEdge), "clamped edge");
542
543                 *new_edge = *edge;
544                 BLI_addtail(clamped_edges, new_edge);
545
546                 voronoi_clampEdgeVertex(width, height, new_edge->start, new_edge->end);
547                 voronoi_clampEdgeVertex(width, height, new_edge->end, new_edge->start);
548
549                 edge = edge->next;
550         }
551 }
552
553 static int voronoi_getNextSideCoord(ListBase *edges, float coord[2], int dim, int dir, float next_coord[2])
554 {
555         VoronoiEdge *edge = edges->first;
556         float distance = FLT_MAX;
557         int other_dim = dim ? 0 : 1;
558
559         while (edge) {
560                 bool ok = false;
561                 float co[2], cur_distance;
562
563                 if (fabsf(edge->start[other_dim] - coord[other_dim]) < VORONOI_EPS &&
564                     len_squared_v2v2(coord, edge->start) > VORONOI_EPS)
565                 {
566                         copy_v2_v2(co, edge->start);
567                         ok = true;
568                 }
569
570                 if (fabsf(edge->end[other_dim] - coord[other_dim]) < VORONOI_EPS &&
571                     len_squared_v2v2(coord, edge->end) > VORONOI_EPS)
572                 {
573                         copy_v2_v2(co, edge->end);
574                         ok = true;
575                 }
576
577                 if (ok) {
578                         if (dir > 0 && coord[dim] > co[dim]) {
579                                 ok = false;
580                         }
581                         else if (dir < 0 && coord[dim] < co[dim]) {
582                                 ok = false;
583                         }
584                 }
585
586                 if (ok) {
587                         cur_distance = len_squared_v2v2(coord, co);
588                         if (cur_distance < distance) {
589                                 copy_v2_v2(next_coord, co);
590                                 distance = cur_distance;
591                         }
592                 }
593
594                 edge = edge->next;
595         }
596
597         return distance < FLT_MAX;
598 }
599
600 static void voronoi_createBoundaryEdges(ListBase *edges, int width, int height)
601 {
602         const float corners[4][2] = {{width - 1, 0.0f},
603                                      {width - 1, height - 1},
604                                      {0.0f, height - 1},
605                                      {0.0f, 0.0f}};
606         int i, dim = 0, dir = 1;
607
608         float coord[2] = {0.0f, 0.0f};
609         float next_coord[2] = {0.0f, 0.0f};
610
611         for (i = 0; i < 4; i++) {
612                 while (voronoi_getNextSideCoord(edges, coord, dim, dir, next_coord)) {
613                         VoronoiEdge *edge = MEM_callocN(sizeof(VoronoiEdge), "boundary edge");
614
615                         copy_v2_v2(edge->start, coord);
616                         copy_v2_v2(edge->end, next_coord);
617                         BLI_addtail(edges, edge);
618
619                         copy_v2_v2(coord, next_coord);
620                 }
621
622                 if (len_squared_v2v2(coord, corners[i]) > VORONOI_EPS) {
623                         VoronoiEdge *edge = MEM_callocN(sizeof(VoronoiEdge), "boundary edge");
624
625                         copy_v2_v2(edge->start, coord);
626                         copy_v2_v2(edge->end, corners[i]);
627                         BLI_addtail(edges, edge);
628                         copy_v2_v2(coord, corners[i]);
629                 }
630
631                 dim = dim ? 0 : 1;
632                 if (i == 1)
633                         dir = -1;
634         }
635 }
636
637 void BLI_voronoi_compute(const VoronoiSite *sites, int sites_total, int width, int height, ListBase *edges)
638 {
639         VoronoiProcess process;
640         VoronoiEdge *edge;
641         int i;
642
643         memset(&process, 0, sizeof(VoronoiProcess));
644
645         process.width = width;
646         process.height = height;
647
648         for (i = 0; i < sites_total; i++) {
649                 VoronoiEvent *event = MEM_callocN(sizeof(VoronoiEvent), "voronoi site event");
650
651                 event->type = voronoiEventType_Site;
652                 copy_v2_v2(event->site, sites[i].co);
653
654                 voronoi_insertEvent(&process, event);
655         }
656
657         while (process.queue.first) {
658                 VoronoiEvent *event = process.queue.first;
659
660                 process.current_y = event->site[1];
661
662                 if (event->type == voronoiEventType_Site) {
663                         voronoi_addParabola(&process, event->site);
664                 }
665                 else {
666                         voronoi_removeParabola(&process, event);
667                 }
668
669                 BLI_freelinkN(&process.queue, event);
670         }
671
672         voronoi_finishEdge(&process, process.root);
673
674         edge = process.edges.first;
675         while (edge) {
676                 if (edge->neighbor) {
677                         copy_v2_v2(edge->start, edge->neighbor->end);
678                         MEM_freeN(edge->neighbor);
679                 }
680
681                 edge = edge->next;
682         }
683
684         BLI_movelisttolist(edges, &process.edges);
685 }
686
687 static bool testVoronoiEdge(const float site[2], const float point[2], const VoronoiEdge *edge)
688 {
689         float p[2];
690
691         if (isect_seg_seg_v2_point(site, point, edge->start, edge->end, p) == 1) {
692                 if (len_squared_v2v2(p, edge->start) > VORONOI_EPS &&
693                     len_squared_v2v2(p, edge->end) > VORONOI_EPS)
694                 {
695                         return false;
696                 }
697         }
698
699         return true;
700 }
701
702 static int voronoi_addTriangulationPoint(const float coord[2], const float color[3],
703                                          VoronoiTriangulationPoint **triangulated_points,
704                                          int *triangulated_points_total)
705 {
706         VoronoiTriangulationPoint *triangulation_point;
707         int i;
708
709         for (i = 0; i < *triangulated_points_total; i++) {
710                 if (equals_v2v2(coord, (*triangulated_points)[i].co)) {
711                         triangulation_point = &(*triangulated_points)[i];
712
713                         add_v3_v3(triangulation_point->color, color);
714                         triangulation_point->power++;
715
716                         return i;
717                 }
718         }
719
720         if (*triangulated_points) {
721                 *triangulated_points = MEM_reallocN(*triangulated_points,
722                                                     sizeof(VoronoiTriangulationPoint) * (*triangulated_points_total + 1));
723         }
724         else {
725                 *triangulated_points = MEM_callocN(sizeof(VoronoiTriangulationPoint), "triangulation points");
726         }
727
728         triangulation_point = &(*triangulated_points)[(*triangulated_points_total)];
729         copy_v2_v2(triangulation_point->co, coord);
730         copy_v3_v3(triangulation_point->color, color);
731
732         triangulation_point->power = 1;
733
734         (*triangulated_points_total)++;
735
736         return (*triangulated_points_total) - 1;
737 }
738
739 static void voronoi_addTriangle(int v1, int v2, int v3, int (**triangles)[3], int *triangles_total)
740 {
741         int *triangle;
742
743         if (*triangles) {
744                 *triangles = MEM_reallocN(*triangles, sizeof(int[3]) * (*triangles_total + 1));
745         }
746         else {
747                 *triangles = MEM_callocN(sizeof(int[3]), "trianglulation triangles");
748         }
749
750         triangle = (int *)&(*triangles)[(*triangles_total)];
751
752         triangle[0] = v1;
753         triangle[1] = v2;
754         triangle[2] = v3;
755
756         (*triangles_total)++;
757 }
758
759 void BLI_voronoi_triangulate(const VoronoiSite *sites, int sites_total, ListBase *edges, int width, int height,
760                              VoronoiTriangulationPoint **triangulated_points_r, int *triangulated_points_total_r,
761                              int (**triangles_r)[3], int *triangles_total_r)
762 {
763         VoronoiTriangulationPoint *triangulated_points = NULL;
764         int (*triangles)[3] = NULL;
765         int triangulated_points_total = 0, triangles_total = 0;
766         int i;
767         ListBase boundary_edges = {NULL, NULL};
768
769         voronoi_clampEdges(edges, width, height, &boundary_edges);
770         voronoi_createBoundaryEdges(&boundary_edges, width, height);
771
772         for (i = 0; i < sites_total; i++) {
773                 VoronoiEdge *edge;
774                 int v1;
775
776                 v1 = voronoi_addTriangulationPoint(sites[i].co, sites[i].color, &triangulated_points, &triangulated_points_total);
777
778                 edge = boundary_edges.first;
779                 while (edge) {
780                         VoronoiEdge *test_edge = boundary_edges.first;
781                         bool ok_start = true, ok_end = true;
782
783                         while (test_edge) {
784                                 if (ok_start && !testVoronoiEdge(sites[i].co, edge->start, test_edge)) {
785                                         ok_start = false;
786                                         break;
787                                 }
788
789                                 if (ok_end && !testVoronoiEdge(sites[i].co, edge->end, test_edge)) {
790                                         ok_end = false;
791                                         break;
792                                 }
793
794                                 test_edge = test_edge->next;
795                         }
796
797                         if (ok_start && ok_end) {
798                                 int v2, v3;
799
800                                 v2 = voronoi_addTriangulationPoint(edge->start, sites[i].color, &triangulated_points, &triangulated_points_total);
801                                 v3 = voronoi_addTriangulationPoint(edge->end, sites[i].color, &triangulated_points, &triangulated_points_total);
802
803                                 voronoi_addTriangle(v1, v2, v3, &triangles, &triangles_total);
804                         }
805
806                         edge = edge->next;
807                 }
808         }
809
810         for (i = 0; i < triangulated_points_total; i++) {
811                 VoronoiTriangulationPoint *triangulation_point = &triangulated_points[i];
812
813                 mul_v3_fl(triangulation_point->color, 1.0f / triangulation_point->power);
814         }
815
816         *triangulated_points_r = triangulated_points;
817         *triangulated_points_total_r = triangulated_points_total;
818
819         *triangles_r = triangles;
820         *triangles_total_r = triangles_total;
821
822         BLI_freelistN(&boundary_edges);
823 }