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