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