Imported bsphere.c, mostly ifdef'd out for now
[blender.git] / source / blender / blenkernel / intern / bsphere.c
1 /*
2  * $Id$
3  *
4  * ***** BEGIN GPL LICENSE BLOCK *****
5  *
6  * This program is free software; you can redistribute it and/or
7  * modify it under the terms of the GNU General Public License
8  * as published by the Free Software Foundation; either version 2
9  * of the License, or (at your option) any later version.
10  *
11  * This program is distributed in the hope that it will be useful,
12  * but WITHOUT ANY WARRANTY; without even the implied warranty of
13  * MERCHANTABILITY or FITNESS FOR A PARTICULAR PURPOSE.  See the
14  * GNU General Public License for more details.
15  *
16  * You should have received a copy of the GNU General Public License
17  * along with this program; if not, write to the Free Software  Foundation,
18  * Inc., 51 Franklin Street, Fifth Floor, Boston, MA 02110-1301, USA.
19  *
20  * The Original Code is Copyright (C) 2011 by Nicholas Bishop
21  * All rights reserved.
22  *
23  * The Original Code is: all of this file.
24  *
25  * Contributor(s): none yet.
26  *
27  * ***** END GPL LICENSE BLOCK *****
28  */
29
30 #include "MEM_guardedalloc.h"
31
32 #include "DNA_mesh_types.h"
33 #include "DNA_meshdata_types.h"
34 #include "DNA_modifier_types.h"
35
36 #include "BLI_utildefines.h"
37 #include "BLI_heap.h"
38 #include "BLI_ghash.h"
39 #include "BLI_listbase.h"
40 #include "BLI_math.h"
41
42 #include "BKE_cdderivedmesh.h"
43 #include "BKE_DerivedMesh.h"
44 #include "BKE_global.h"
45 #include "BKE_library.h"
46 #include "BKE_main.h"
47 #include "BKE_subsurf.h"
48
49 #include <assert.h>
50
51 typedef struct Tri {
52         struct Tri *next, *prev;
53         unsigned int v[3];
54         float no[3];
55         float area;
56         int marked;
57 } Tri;
58
59 typedef struct Edge {
60         struct Edge *next, *prev;
61         unsigned int v[2];
62         Tri *tri[2];
63 } Edge;
64
65 /* note: a `frame' is four vertices (contiguous within the MVert
66    array), stored simply as the index of the first vertex */
67
68 static void create_frame(float frame[4][3], const float co[3], float radius,
69                          const float mat[3][3], float x_offset)
70 {
71         float rx[3], ry[3], rz[3];
72         int i;
73
74         mul_v3_v3fl(ry, mat[1], radius);
75         mul_v3_v3fl(rz, mat[2], radius);
76         
77         add_v3_v3v3(frame[3], co, ry);
78         add_v3_v3v3(frame[3], frame[3], rz);
79
80         sub_v3_v3v3(frame[2], co, ry);
81         add_v3_v3v3(frame[2], frame[2], rz);
82
83         sub_v3_v3v3(frame[1], co, ry);
84         sub_v3_v3v3(frame[1], frame[1], rz);
85
86         add_v3_v3v3(frame[0], co, ry);
87         sub_v3_v3v3(frame[0], frame[0], rz);
88
89         mul_v3_v3fl(rx, mat[0], x_offset);
90         for(i = 0; i < 4; i++)
91                 add_v3_v3v3(frame[i], frame[i], rx);
92 }
93
94 static void create_mid_frame(float frame[4][3],
95                              const float co1[3], const float co2[3],
96                              const SkinNode *n1, const SkinNode *n2,
97                              const float mat[3][3])
98 {
99         create_frame(frame, co1, (n1->radius + n2->radius) / 2,
100                      mat, len_v3v3(co1, co2) / 2);
101 }
102
103 static void add_face(MFace *mface, int *totface,
104                      int v1, int v2, int v3, int v4)
105 {
106         MFace *f;
107
108         f = &mface[*totface];
109         f->v1 = v1;
110         f->v2 = v2;
111         f->v3 = v3;
112         f->v4 = v4 == -1 ? 0 : v4;
113         if(!v4) {
114                 f->v1 = v3;
115                 f->v2 = v4;
116                 f->v3 = v1;
117                 f->v4 = v2;
118         }
119         (*totface)++;
120 }
121
122 static void connect_frames(MFace *mface, int *totface, int frame1, int frame2)
123 {
124         int i;
125
126         for(i = 0; i < 4; i++) {
127                 add_face(mface, totface,
128                          frame2 + i,
129                          frame2 + (i+1) % 4,
130                          frame1 + (i+1) % 4,
131                          frame1 + i);
132         }
133 }
134
135 static Tri *add_triangle(ListBase *tris, MVert *mvert, int v1, int v2, int v3)
136 {
137         Tri *t;
138
139         t = MEM_callocN(sizeof(Tri), "add_triangle.tri");
140         t->v[0] = v1;
141         t->v[1] = v2;
142         t->v[2] = v3;
143         BLI_addtail(tris, t);
144
145         normal_tri_v3(t->no, mvert[t->v[0]].co, mvert[t->v[1]].co, mvert[t->v[2]].co);
146         t->area = area_tri_v3(mvert[t->v[0]].co, mvert[t->v[1]].co, mvert[t->v[2]].co);
147
148         return t;
149 }
150
151 static int point_tri_side(const Tri *t, MVert *mvert, int v)
152 {
153         float p[3], d;
154         sub_v3_v3v3(p, mvert[v].co, mvert[t->v[0]].co);
155         d = dot_v3v3(t->no, p);
156         if(d < 0) return -1;
157         else if(d > 0) return 1;
158         else return 0;
159 }
160
161 static int mark_v_outside_tris(ListBase *tris, MVert *mvert, int v)
162 {
163         Tri *t;
164         int outside = 0;
165
166         /* probably there's a much better way to do this */
167         for(t = tris->first; t; t = t->next) {
168                 if((t->marked = point_tri_side(t, mvert, v) > 0))
169                         outside = 1;
170         }
171         return outside;
172 }
173
174 static int edge_match(int e1_0, int e1_1, const unsigned int e2[2])
175 {
176         /* XXX: maybe isn't necesseary to check both directions? */
177         return (e1_0 == e2[0] && e1_1 == e2[1]) ||
178                (e1_0 == e2[1] && e1_1 == e2[0]);
179 }
180
181 /* returns true if the edge (e1, e2) is already in edges; that edge is
182    deleted here as well. if not found just returns 0 */
183 static int check_for_dup(ListBase *edges, unsigned int e1, unsigned int e2)
184 {
185         Edge *e, *next;
186
187         for(e = edges->first; e; e = next) {
188                 next = e->next;
189
190                 if(edge_match(e1, e2, e->v)) {
191                         /* remove the interior edge */
192                         BLI_freelinkN(edges, e);
193                         return 1;
194                 }
195         }
196
197         return 0;
198 }
199
200 static void expand_boundary_edges(ListBase *edges, Tri *t)
201 {
202         Edge *new;
203         int i;
204
205         /* insert each triangle edge into the boundary list; if any of
206            its edges are already in there, remove the edge entirely */
207         for(i = 0; i < 3; i++) {
208                 if(!check_for_dup(edges, t->v[i], t->v[(i+1)%3])) {
209                         new = MEM_callocN(sizeof(Edge), "Edge");
210                         new->v[0] = t->v[i];
211                         new->v[1] = t->v[(i+1)%3];
212                         BLI_addtail(edges, new);
213                 }
214         }
215 }
216
217 static int tri_matches_frame(const Tri *t, int frame)
218 {
219         int i, j, match = 0;
220         for(i = 0; i < 3; i++) {
221                 for(j = 0; j < 4; j++) {
222                         if(t->v[i] == frame + j) {
223                                 match++;
224                                 if(match >= 3)
225                                         return 1;
226                         }
227                 }
228         }
229         return 0;
230 }
231
232 static void quad_from_tris(const Edge *e, int ndx[4])
233 {
234         int i, j, opp = -1;
235
236         /* find what the second tri has that the first doesn't */
237         for(i = 0; i < 3; i++) {
238                 if(e->tri[1]->v[i] != e->tri[0]->v[0] &&
239                    e->tri[1]->v[i] != e->tri[0]->v[1] &&
240                    e->tri[1]->v[i] != e->tri[0]->v[2]) {
241                         opp = e->tri[1]->v[i];
242                         break;
243                 }
244         }
245         assert(opp != -1);
246
247         for(i = 0, j = 0; i < 3; i++, j++) {
248                 ndx[j] = e->tri[0]->v[i];
249                 /* when the triangle edge cuts across our quad-to-be,
250                    throw in the second triangle's vertex */
251                 if((e->tri[0]->v[i] == e->v[0] || e->tri[0]->v[i] == e->v[1]) &&
252                    (e->tri[0]->v[(i+1)%3] == e->v[0] || e->tri[0]->v[(i+1)%3] == e->v[1])) {
253                         j++;
254                         ndx[j] = opp;
255                 }
256         }
257 }
258
259 static void add_quad_from_tris(MFace *mface, int *totface, const Edge *e)
260 {
261         int ndx[4];
262
263         quad_from_tris(e, ndx);
264
265         add_face(mface, totface, ndx[0], ndx[1], ndx[2], ndx[3]);
266 }
267
268 static GHash *calc_edges(ListBase *tris)
269 {
270         GHash *adj;
271         ListBase *edges;
272         Edge *e;
273         Tri *t;
274         int i, e1, e2;
275
276         /* XXX: vertex range might be a little funky? so using a
277            hash here */
278         adj = BLI_ghash_new(BLI_ghashutil_inthash,
279                             BLI_ghashutil_intcmp,
280                             "calc_edges adj");
281
282         for(t = tris->first; t; t = t->next) {
283                 for(i = 0; i < 3; i++) {
284                         e1 = t->v[i];
285                         e2 = t->v[(i+1)%3];
286                         assert(e1 != e2);
287                         if(e1 > e2)
288                                 SWAP(int, e1, e2);
289
290                         edges = BLI_ghash_lookup(adj, SET_INT_IN_POINTER(e1));
291                         if(!edges) {
292                                 edges = MEM_callocN(sizeof(ListBase),
293                                                     "calc_edges ListBase");
294                                 BLI_ghash_insert(adj, SET_INT_IN_POINTER(e1), edges);
295                         }
296
297                         /* find the edge in the adjacency list */
298                         for(e = edges->first; e; e = e->next) {
299                                 assert(e->v[0] == e1);
300                                 if(e->v[1] == e2)
301                                         break;
302                         }
303
304                         /* if not found, create the edge */
305                         if(!e) {
306                                 e = MEM_callocN(sizeof(Edge), "calc_edges Edge");
307                                 e->v[0] = e1;
308                                 e->v[1] = e2;
309                                 BLI_addtail(edges, e);
310                         }
311
312                         /* should never be more than two faces
313                            attached to an edge here */
314                         assert(!e->tri[0] || !e->tri[1]);
315
316                         if(!e->tri[0])
317                                 e->tri[0] = t;
318                         else if(!e->tri[1])
319                                 e->tri[1] = t;
320                 }
321         }
322
323         return adj;
324 }
325
326 static int is_quad_symmetric(const MVert *mvert, const Edge *e)
327 {
328         int ndx[4];
329         float a[3];
330
331         quad_from_tris(e, ndx);
332
333         copy_v3_v3(a, mvert[ndx[0]].co);
334         a[0] = -a[0];
335
336         if(len_v3v3(a, mvert[ndx[1]].co) < FLT_EPSILON) {
337                 copy_v3_v3(a, mvert[ndx[2]].co);
338                 a[0] = -a[0];
339                 if(len_v3v3(a, mvert[ndx[3]].co) < FLT_EPSILON)
340                         return 1;
341         }
342         else if(len_v3v3(a, mvert[ndx[3]].co) < FLT_EPSILON) {
343                 copy_v3_v3(a, mvert[ndx[2]].co);
344                 a[0] = -a[0];
345                 if(len_v3v3(a, mvert[ndx[1]].co) < FLT_EPSILON)
346                         return 1;
347         }
348
349         return 0;
350 }
351
352 static void output_hull(const MVert *mvert, MFace *mface, int *totface, ListBase *tris)
353 {
354         Heap *heap;
355         GHash *adj;
356         GHashIterator *iter;
357         ListBase *edges;
358         Edge *e;
359         Tri *t;
360         float score;
361
362         heap = BLI_heap_new();
363         adj = calc_edges(tris);
364
365         /* unmark all triangles */
366         for(t = tris->first; t; t = t->next)
367                 t->marked = 0;
368
369         /* build heap */
370         iter = BLI_ghashIterator_new(adj);
371         while(!BLI_ghashIterator_isDone(iter)) {
372                 edges = BLI_ghashIterator_getValue(iter);
373                 for(e = edges->first; e; e = e->next) {
374                         /* only care if the edge is used by more than
375                            one triangle */
376                         if(e->tri[0] && e->tri[1]) {
377                                 score = (e->tri[0]->area + e->tri[1]->area) *
378                                         dot_v3v3(e->tri[0]->no, e->tri[1]->no);
379
380                                 /* increase score if the triangles
381                                    form a symmetric quad */
382                                 if(is_quad_symmetric(mvert, e))
383                                         score *= 10;
384
385                                 BLI_heap_insert(heap, -score, e);
386                         }
387                 }
388
389                 BLI_ghashIterator_step(iter);
390         }
391         BLI_ghashIterator_free(iter);
392
393         while(!BLI_heap_empty(heap)) {
394                 e = BLI_heap_popmin(heap);
395                 
396                 /* if both triangles still free, outupt as a quad */
397                 if(!e->tri[0]->marked && !e->tri[1]->marked) {
398                         add_quad_from_tris(mface, totface, e);
399                         e->tri[0]->marked = 1;
400                         e->tri[1]->marked = 1;
401                 }
402         }
403
404         /* free edge list */
405         iter = BLI_ghashIterator_new(adj);
406         while(!BLI_ghashIterator_isDone(iter)) {
407                 edges = BLI_ghashIterator_getValue(iter);
408                 BLI_freelistN(edges);
409                 MEM_freeN(edges);       
410                 BLI_ghashIterator_step(iter);
411         }
412         BLI_ghashIterator_free(iter);
413         BLI_ghash_free(adj, NULL, NULL);
414         BLI_heap_free(heap, NULL);
415
416         /* write out any remaining triangles */
417         for(t = tris->first; t; t = t->next) {
418                 if(!t->marked)
419                         add_face(mface, totface, t->v[0], t->v[1], t->v[2], -1);
420         }
421 }       
422
423 /* for vertex `v', find which triangles must be deleted to extend the
424    hull; find the boundary edges of that hole so that it can be filled
425    with connections to the new vertex, and update the `tri' list to
426    delete the marked triangles */
427 static void add_point(ListBase *tris, MVert *mvert, int v)
428 {
429         ListBase edges = {NULL, NULL};
430         Tri *t, *next;
431         Edge *e;
432
433         for(t = tris->first; t; t = next) {
434                 next = t->next;
435
436                 /* check if triangle is `visible' to v */
437                 if(t->marked) {
438                         expand_boundary_edges(&edges, t);
439                         /* remove the triangle */
440                         BLI_freelinkN(tris, t);
441                 }
442         }
443
444         /* fill hole boundary with triangles to new point */
445         for(e = edges.first; e; e = e->next)
446                 add_triangle(tris, mvert, e->v[0], e->v[1], v);
447         BLI_freelistN(&edges);
448 }
449
450 static void build_hull(MFace *mface, int *totface,
451                        MVert *mvert, int *frames, int totframe)
452 {
453         ListBase tris = {NULL, NULL};
454         Tri *t, *next;
455         int i, j;
456
457         /* use first frame to make initial degenerate pyramid */
458         add_triangle(&tris, mvert, frames[0], frames[0] + 1, frames[0] + 2);
459         add_triangle(&tris, mvert, frames[0] + 1, frames[0], frames[0] + 3);
460         add_triangle(&tris, mvert, frames[0] + 2, frames[0] + 1, frames[0] + 3);
461         add_triangle(&tris, mvert, frames[0], frames[0] + 2, frames[0] + 3);
462
463         for(i = 1; i < totframe; i++) {
464                 for(j = 0; j < 4; j++) {
465                         int v = frames[i] + j;
466
467                         /* ignore point already inside hull */
468                         if(mark_v_outside_tris(&tris, mvert, v)) {
469                                 /* expand hull and delete interior triangles */
470                                 add_point(&tris, mvert, v);
471                         }
472                 }
473         }
474
475         /* remove triangles that would fill the original frames */
476         for(t = tris.first; t; t = next) {
477                 next = t->next;
478
479                 for(i = 0; i < totframe; i++) {
480                         if(tri_matches_frame(t, frames[i])) {
481                                 BLI_freelinkN(&tris, t);
482                                 break;
483                         }
484                 }
485         }
486
487         output_hull(mvert, mface, totface, &tris);
488
489         BLI_freelistN(&tris);
490 }
491
492 /* test: build hull on origdm */
493 #if 0
494 static DerivedMesh *test_hull(DerivedMesh *origdm)
495 {
496         DerivedMesh *dm;
497         MVert *mvert, *origmvert;
498         MFace *mface;
499         int *frames;
500         int totvert = 0, totface = 0, totframe, i;
501
502         /* TODO */
503         totface = 2000;
504         totvert = origdm->getNumVerts(origdm);
505         dm = CDDM_new(totvert, 0, totface);
506
507         mvert = CDDM_get_verts(dm);
508         mface = CDDM_get_faces(dm);
509
510         origmvert = CDDM_get_verts(origdm);
511         for(i = 0; i < totvert; i++)
512                 mvert[i] = origmvert[i];
513
514         assert(totvert % 4 == 0);
515         totframe = totvert / 4;
516
517         frames = MEM_callocN(sizeof(int) * totframe, "frames");
518         for(i = 0; i < totframe; i++)
519                 frames[i] = i*4;
520
521         totface = 0;
522         build_hull(mface, &totface, mvert, frames, totframe);
523
524         CDDM_calc_edges(dm);
525         //CDDM_calc_normals(dm);
526
527         MEM_freeN(frames);
528
529         return dm;
530 }
531 #endif
532
533 /* TODO */
534 static void calc_edge_mat(float mat[3][3], const float a[3], const float b[3])
535 {
536         float Z[3] = {0, 0, 1};
537         float dot;
538
539         /* x = edge direction */
540         sub_v3_v3v3(mat[0], b, a);
541         normalize_v3(mat[0]);
542
543         dot = dot_v3v3(mat[0], Z);
544         if(dot > -1 + FLT_EPSILON && dot < 1 - FLT_EPSILON) {
545                 /* y = Z cross x */
546                 cross_v3_v3v3(mat[1], Z, mat[0]);
547                 normalize_v3(mat[1]);
548
549                 /* z = x cross y */
550                 cross_v3_v3v3(mat[2], mat[0], mat[1]);
551                 normalize_v3(mat[2]);
552         }
553         else {
554                 mat[1][0] = 1;
555                 mat[1][1] = 0;
556                 mat[1][2] = 0;
557                 mat[2][0] = 0;
558                 mat[2][1] = 1;
559                 mat[2][2] = 0;
560         }
561 }
562
563 /* BMesh paper */
564 #if 0
565 static float calc_R_sub_i_squared(const BSphereNode *n)
566 {
567         const float alpha = 1.5;
568         //const float alpha = 1;
569         const float R_sub_i = n->size[0] * alpha;
570         return R_sub_i * R_sub_i;
571 }
572
573 /* equation (2) */
574 static float calc_r_squared(const float v[3], const BSphereNode *n)
575 {
576         return len_squared_v3v3(v, n->co);
577 }
578
579 /* equation (1) */
580 static float calc_f_sub_i(const float v[3], const BSphereNode *n)
581 {
582         float r_squared;
583         float R_sub_i_squared;
584
585         r_squared = calc_r_squared(v, n);
586         R_sub_i_squared = calc_R_sub_i_squared(n);
587
588         if(r_squared <= R_sub_i_squared)
589                 return powf(1.0f - (r_squared / R_sub_i_squared), 2);
590         else
591                 return 0;
592 }
593
594 /* equation (3) */
595 static float calc_I(const float x[3], float T, const ListBase *base)
596 {
597         BSphereNode *n;
598         float a = -T;
599
600         for(n = base->first; n; n = n->next) {
601                 a += calc_I(x, 0, &n->children);
602                 a += calc_f_sub_i(x, n);
603         }
604
605         //printf("calc_I: %f, %f\n", -T, a);
606
607         return a;
608 }
609
610 /* kinda my own guess here */
611 static void calc_gradient(float g[3], const float x[3], const ListBase *base)
612 {
613         BSphereNode *n;
614
615         for(n = base->first; n; n = n->next) {
616                 float R_sub_i_squared = calc_R_sub_i_squared(n);
617                 float dist = len_v3v3(x, n->co);
618                 float f = calc_f_sub_i(x, n);
619                 float d = 1.0f / R_sub_i_squared;
620                 float tmp[3];
621
622                 /* XXX: other possibilities... */
623                 if(f > 0) {
624                         sub_v3_v3v3(tmp, x, n->co);
625                         mul_v3_fl(tmp, -4*d * (1 - d * dist));
626
627                         /* not sure if I'm doing this right; intent is
628                            to flip direction when x is `beneath' the
629                            level set */
630                         if(len_v3v3(x, n->co) > n->size[0])
631                                 negate_v3(tmp);
632
633                         add_v3_v3(g, tmp);
634                 }
635
636                 calc_gradient(g, x, &n->children);
637         }
638 }
639
640 /* equation (5) */
641 static float calc_F(const float x[3], const float K[2], float T, const ListBase *base)
642 {
643         float f_of_K;
644         float I_target;
645
646         /* TODO */
647         f_of_K = 1.0f / (1 + fabs(K[0]) + fabs(K[1]));
648         //f_of_K = 1;
649
650         /* TODO */
651         I_target = 0;
652
653         return (calc_I(x, T, base) - I_target) /*TODO: * f_of_K*/;
654 }
655
656 static float smallest_radius(const ListBase *base)
657 {
658         BSphereNode *n;
659         float s = FLT_MAX, t;
660
661         for(n = base->first; n; n = n->next) {
662                 if(n->size[0] < s)
663                         s = n->size[0];
664                 if((t = smallest_radius(&n->children)) < s)
665                         s = t;
666         }
667
668         return s;
669 }
670
671 /* equation (7) */
672 static float calc_delta_t(float F_max, float k, const ListBase *base, int print)
673 {
674         float min_r_sub_i;
675         float step;
676
677         min_r_sub_i = smallest_radius(base);
678         step = min_r_sub_i / powf(2, k);
679
680         if(print) {
681                 printf("min_r: %f, 2^k=%f, step=%f\n", min_r_sub_i, powf(2, k), step);
682         }
683
684         return step / F_max;
685 }
686
687 /* equation (6) */
688 static float evolve_x(float x[3], const float K[2], float T,
689                       int k, const ListBase *base, int i, float F_max)
690 {
691         float tmp[3], no[3];
692         float F, dt;
693
694         F = calc_F(x, K, T, base);
695
696         if(F < FLT_EPSILON) {
697                 //printf("stability reached for ndx=%d\n", i);
698                 return 0;
699         }
700
701         /* TODO ;-) */
702         if(F > F_max)
703                 F_max = F;
704
705         dt = calc_delta_t(F_max, k, base, 0);
706
707         dt = F / 2;
708
709         if(i == 0)
710                 ;//printf("F=%.3f, dt=%.3f\n", F, calc_delta_t(F_max, k, base, 1));
711         
712         zero_v3(no);
713         calc_gradient(no, x, base);
714         normalize_v3(no);
715         negate_v3(no);
716
717         mul_v3_v3fl(tmp, no, F * dt);
718         add_v3_v3(x, tmp);
719
720         return F_max;
721 }
722
723 static void bsphere_evolve(MVert *mvert, int totvert, const ListBase *base,
724                            float T, int steps, float subdiv_level) 
725 {
726         int i, s;
727
728         for(i = 0; i < totvert; i++) {
729                 float F_max = 0;
730                 /* TODO: for now just doing a constant number of steps */
731                 for(s = 0; s < steps; s++) {
732                         /* TODO */
733                         float K[2] = {0, 0};
734
735                         F_max = evolve_x(mvert[i].co, K, T, subdiv_level, base, i, F_max);
736                 }
737         }
738 }
739 #endif
740
741 typedef enum {
742         START_CAP = 1,
743         END_CAP = 2,
744         NODE_BISECT = 4,
745         CHILD_EDGE_BISECT = 8,
746         
747         IN_FRAME = 16,
748         OUT_FRAME = 32,
749 } FrameFlag;
750
751 typedef struct {
752         int totframe;
753         int outbase;
754         FrameFlag flag;
755         float co[0][4][3];
756 } Frames;
757
758 static Frames *frames_create(FrameFlag flag)
759 {
760         Frames *frames;
761         int totframe;
762         totframe = (!!(flag & START_CAP) + !!(flag & END_CAP) +
763                     !!(flag & NODE_BISECT) + !!(flag & CHILD_EDGE_BISECT));
764         assert(totframe);
765         frames = MEM_callocN(sizeof(Frames) + sizeof(float)*totframe*4*3,
766                              "frames_create.frames");
767         frames->totframe = totframe;
768         frames->flag = flag;
769         assert((flag & (END_CAP|CHILD_EDGE_BISECT)) != (END_CAP|CHILD_EDGE_BISECT));
770         return frames;
771 }
772
773 static int frame_offset(Frames *frames, FrameFlag flag)
774 {
775         int n = 0;
776         
777         assert(flag == IN_FRAME || flag == OUT_FRAME ||
778                !(flag & (IN_FRAME|OUT_FRAME)));
779
780         if(flag == IN_FRAME)
781                 return 0;
782         else if(flag == OUT_FRAME)
783                 return frames->totframe-1;
784         
785         assert((flag & frames->flag) == flag);
786
787         if(flag == START_CAP) return n;
788         if(frames->flag & START_CAP) n++;
789
790         if(flag == NODE_BISECT) return n;
791         if(frames->flag & NODE_BISECT) n++;
792         
793         if(flag == END_CAP) return n;
794         if(frames->flag & END_CAP) n++;
795
796         if(flag == CHILD_EDGE_BISECT) return n;
797         if(frames->flag & CHILD_EDGE_BISECT) n++;
798
799         assert(!"bad frame offset flag");
800 }
801
802 static int frame_base(Frames *frames, FrameFlag flag)
803 {
804         return frames->outbase + 4*frame_offset(frames, flag);
805 }
806
807 static void *frame_co(Frames *frames, FrameFlag flag)
808 {
809         return frames->co[frame_offset(frames, flag)];
810 }
811
812 static int frames_totvert(Frames *frames)
813 {
814         return frames->totframe*4;
815 }
816
817 #if 0
818 static int bsphere_build_frames(BSphereNode *n, GHash *ghash,
819                                 float (*parent_mat)[3], FrameFlag flag);
820
821 /* builds the `cap' of polygons for the end of a BSphere chain */
822 static int bsphere_end_node_frames(BSphereNode *n,
823                                    GHash *ghash,
824                                    float parent_mat[3][3],
825                                    FrameFlag flag)
826 {
827         Frames *frames;
828         float rad;
829
830         assert(flag == 0 || flag == START_CAP);
831         frames = frames_create(NODE_BISECT|END_CAP|flag, n, ghash);
832
833         /* TODO */
834         rad = n->size[0];
835
836         if(flag & START_CAP)
837                 create_frame(frame_co(frames, START_CAP), n, parent_mat, -rad, rad);
838
839         /* frame to bisect the node */
840         create_frame(frame_co(frames, NODE_BISECT), n, parent_mat, 0, rad);
841
842         /* frame to cap end node */
843         create_frame(frame_co(frames, END_CAP), n, parent_mat, rad, rad);
844
845         return frames_totvert(frames);
846 }
847
848 static int bsphere_connection_node_frames(BSphereNode *n,
849                                           GHash *ghash,
850                                           float parent_mat[3][3],
851                                           FrameFlag flag)
852 {
853         Frames *frames;
854         BSphereNode *child;
855         float childmat[3][3], nodemat[3][3], axis[3], angle, rad, child_rad;
856         int totvert;
857
858         assert(flag == 0 || flag == START_CAP);
859
860         child = n->children.first;
861
862         /* if child is not a branch node, bisect the child edge */
863         if(BLI_countlist(&child->children) <= 1)
864                 flag |= CHILD_EDGE_BISECT;
865
866         frames = frames_create(NODE_BISECT|flag, n, ghash);
867         totvert = frames_totvert(frames);
868
869         /* TODO */
870         rad = n->size[0];
871         child_rad = child->size[0];
872
873         /* build matrix to give to child */
874         sub_v3_v3v3(childmat[0], child->co, n->co);
875         normalize_v3(childmat[0]);
876         angle = angle_normalized_v3v3(parent_mat[0], childmat[0]);
877         cross_v3_v3v3(axis, parent_mat[0], childmat[0]);
878         normalize_v3(axis);
879         rotate_normalized_v3_v3v3fl(childmat[1], parent_mat[1], axis, angle);
880         rotate_normalized_v3_v3v3fl(childmat[2], parent_mat[2], axis, angle);
881
882         /* build child */
883         totvert += bsphere_build_frames(child, ghash, childmat, 0);
884
885         /* frame that bisects the node */
886         angle /= 2;
887         copy_v3_v3(nodemat[0], parent_mat[0]);
888         rotate_normalized_v3_v3v3fl(nodemat[1], parent_mat[1], axis, angle);
889         rotate_normalized_v3_v3v3fl(nodemat[2], parent_mat[2], axis, angle);
890         create_frame(frame_co(frames, NODE_BISECT), n, nodemat, 0, rad);
891
892         if(flag & START_CAP)
893                 create_frame(frame_co(frames, START_CAP), n, nodemat, -rad, rad);
894
895         if(flag & CHILD_EDGE_BISECT)
896                 create_mid_frame(frame_co(frames, CHILD_EDGE_BISECT), n, child, childmat);
897         
898         return totvert;
899 }
900
901 static int bsphere_branch_node_frames(BSphereNode *n, GHash *ghash)
902 {
903         BSphereNode *child;
904         int totvert = 0;
905         float rad;
906
907         rad = n->size[0];
908
909         /* build children */
910         for(child = n->children.first; child; child = child->next) {
911                 float childmat[3][3];
912
913                 calc_edge_mat(childmat, n->co, child->co);
914
915                 totvert += bsphere_build_frames(child, ghash, childmat, 0);
916         }
917
918         /* TODO: for now, builds no frames of its own */
919
920         return totvert;
921 }
922
923 static int bsphere_build_frames(BSphereNode *n, GHash *ghash,
924                                 float parent_mat[3][3], FrameFlag flag)
925 {
926         BSphereNode *child;
927
928         child = n->children.first;
929
930         if(!child)
931                 return bsphere_end_node_frames(n, ghash, parent_mat, flag);
932         else if(!child->next)
933                 return bsphere_connection_node_frames(n, ghash, parent_mat, flag);
934         else
935                 return bsphere_branch_node_frames(n, ghash);
936 }
937
938 static GHash *bsphere_frames(BSphere *bs, int *totvert)
939 {
940         GHash *ghash;
941         BSphereNode *n, *child;
942
943         ghash = BLI_ghash_new(BLI_ghashutil_ptrhash, BLI_ghashutil_ptrcmp,
944                               "bsphere_frames.ghash");
945         (*totvert) = 0;
946
947         for(n = bs->rootbase.first; n; n = n->next) {
948                 float mat[3][3];
949
950                 child = n->children.first;
951                 if(!child || (child && child->next)) {
952                         set_v3(mat[0], 0, 0, 1);
953                         set_v3(mat[1], 1, 0, 0);
954                         set_v3(mat[2], 0, 1, 0);
955                 }
956                 else if(!child->next) {
957                         calc_edge_mat(mat, n->co, child->co);
958                 }
959
960                 (*totvert) += bsphere_build_frames(n, ghash, mat, START_CAP);
961         }
962
963         return ghash;
964 }
965
966 static void output_frames(Frames *frames, MVert *mvert, MFace *mface,
967                           int *totvert, int *totface)
968 {
969         int b, i, j;
970
971         frames->outbase = (*totvert);
972         for(i = 0; i < frames->totframe; i++) {
973                 for(j = 0; j < 4; j++) {
974                         copy_v3_v3(mvert[(*totvert)].co, frames->co[i][j]);
975                         (*totvert)++;
976                 }
977
978                 if(i > 0) {
979                         connect_frames(mface, totface,
980                                        frames->outbase + 4*(i-1),
981                                        frames->outbase + 4*i);
982                 }
983         }
984
985         if(frames->flag & START_CAP) {
986                 b = frame_base(frames, START_CAP);
987                 add_face(mface, totface, b, b+1, b+2, b+3);
988         }
989         if(frames->flag & END_CAP) {
990                 b = frame_base(frames, END_CAP);
991                 add_face(mface, totface, b+3, b+2, b+1, b+0);
992         }
993 }
994
995 static void bsphere_build_skin(BSphereNode *parent, ListBase *base,
996                                GHash *frames, MVert *mvert,
997                                MFace *mface, int *totvert, int *totface)
998 {
999         BSphereNode *n, *child;
1000         Frames *node_frames, *child_frames, *parent_frames;
1001         int *hull_frames, totchild, i;
1002
1003         for(n = base->first; n; n = n->next) {
1004                 if((node_frames = BLI_ghash_lookup(frames, n)))
1005                         output_frames(node_frames, mvert, mface, totvert, totface);
1006
1007                 bsphere_build_skin(n, &n->children, frames, mvert, mface, totvert, totface);
1008
1009                 child = n->children.first;
1010                 if(child && !child->next) {
1011                         child_frames = BLI_ghash_lookup(frames, child);
1012                         if(child_frames) {
1013                                 connect_frames(mface, totface,
1014                                                frame_base(node_frames, OUT_FRAME),
1015                                                frame_base(child_frames, IN_FRAME));
1016                         }
1017                 }
1018                 else if(child && child->next) {
1019                         totchild = BLI_countlist(&n->children);
1020                         hull_frames = MEM_callocN(sizeof(int) * (1+totchild),
1021                                                   "bsphere_build_skin.hull_frames");
1022                         i = 0;
1023                         for(child = n->children.first; child; child = child->next) {
1024                                 child_frames = BLI_ghash_lookup(frames, child);
1025                                 assert(child_frames);
1026                                 hull_frames[i++] = frame_base(child_frames, IN_FRAME);
1027                         }
1028                         parent_frames = BLI_ghash_lookup(frames, parent);
1029                         if(parent_frames)
1030                                 hull_frames[i] = frame_base(parent_frames, OUT_FRAME);
1031                         else
1032                                 totchild--;
1033                         build_hull(mface, totface, mvert, hull_frames, totchild+1);
1034                         MEM_freeN(hull_frames);
1035                 }
1036         }
1037 }
1038
1039 static DerivedMesh *bsphere_base_skin(BSphere *bs)
1040 {
1041         DerivedMesh *dm;
1042         GHash *ghash;
1043         MVert *mvert;
1044         MFace *mface;
1045         int totvert, totface = 0;
1046
1047         /* no nodes, nothing to do */
1048         if(!bs->rootbase.first)
1049                 return NULL;
1050
1051         /* could probably do all this in one fell [recursive] swoop
1052            using some tricksies and clevers, but multiple passes is
1053            easier :) */
1054         ghash = bsphere_frames(bs, &totvert);
1055
1056         totface = totvert * 1.5;
1057         /* XXX: sizes OK? */
1058         dm = CDDM_new(totvert, 0, totface);
1059
1060         mvert = CDDM_get_verts(dm);
1061         mface = CDDM_get_faces(dm);
1062
1063         totvert = totface = 0;
1064         bsphere_build_skin(NULL, &bs->rootbase, ghash, mvert, mface, &totvert, &totface);
1065         BLI_ghash_free(ghash, NULL, (void*)MEM_freeN);
1066
1067         CDDM_lower_num_verts(dm, totvert);
1068         CDDM_lower_num_faces(dm, totface);
1069
1070         CDDM_calc_edges(dm);
1071         CDDM_calc_normals(dm);
1072
1073         return dm;
1074 }
1075
1076 /* generate random points (just for testing) */
1077 DerivedMesh *bsphere_random_points(BSphere *UNUSED(bs))
1078 {
1079         DerivedMesh *dm;
1080         MVert *mvert;
1081         const int totvert = 1;//50000;
1082         int i;
1083
1084         dm = CDDM_new(totvert, 0, 0);
1085
1086         mvert = CDDM_get_verts(dm);
1087         
1088         srand(1);
1089         for(i = 0; i < totvert; i++) {
1090                 mvert[i].co[0] = rand() / (float)RAND_MAX - 0.5;
1091                 mvert[i].co[1] = rand() / (float)RAND_MAX - 0.5;
1092                 mvert[i].co[2] = rand() / (float)RAND_MAX - 0.5;
1093         }
1094         
1095         return dm;
1096 }
1097
1098 DerivedMesh *bsphere_test_surface(BSphere *bs, DerivedMesh *input)
1099 {
1100         DerivedMesh *dm;
1101         MVert *mvert;
1102
1103         if(!input)
1104                 return NULL;
1105
1106         dm = CDDM_copy(input);
1107         mvert = CDDM_get_verts(dm);
1108         bsphere_evolve(mvert, input->getNumVerts(input), &bs->rootbase,
1109                        bs->skin_threshold, bs->evolve, 2);
1110
1111         return dm;
1112 }
1113
1114 DerivedMesh *bsphere_test_gradient(BSphere *bs, DerivedMesh *input)
1115 {
1116         DerivedMesh *dm;
1117         MVert *mvert;
1118         MEdge *medge;
1119         int totvert;
1120         int i;
1121
1122         if(!input)
1123                 return NULL;
1124
1125         totvert = input->getNumVerts(input);
1126         dm = CDDM_new(totvert*2, totvert, 0);
1127
1128         mvert = CDDM_get_verts(dm);
1129         medge = CDDM_get_edges(dm);
1130
1131         memcpy(mvert, CDDM_get_verts(input), sizeof(MVert) * totvert);
1132         
1133         for(i = 0; i < totvert; i++) {
1134                 const float *b = mvert[i].co;
1135                 float *g = mvert[totvert+i].co;
1136                 float K[2] = {0, 0};
1137
1138                 zero_v3(g);
1139                 calc_gradient(g, b, &bs->rootbase);
1140                 normalize_v3(g);
1141                 
1142                 mul_v3_fl(g, -calc_F(b, K, bs->skin_threshold, &bs->rootbase));
1143                 add_v3_v3(g, b);
1144
1145                 medge[i].v1 = i;
1146                 medge[i].v2 = totvert+i;
1147         }
1148         
1149         return dm;
1150 }
1151
1152 DerivedMesh *bsphere_get_skin_dm(BSphere *bs)
1153 {
1154         DerivedMesh *dm, *subdm;
1155
1156         dm = bsphere_base_skin(bs);
1157
1158         /* subdivide once */
1159         if(dm && bs->subdiv_level > 0) {
1160                 SubsurfModifierData smd;
1161
1162                 memset(&smd, 0, sizeof(smd));
1163                 smd.subdivType = ME_CC_SUBSURF;
1164                 smd.levels = bs->subdiv_level;
1165                 subdm = subsurf_make_derived_from_derived(dm, &smd, 0, 0, 0, 0);
1166                 dm->release(dm);
1167                 /* for convenience, convert back to cddm (could use ccgdm
1168                    later for better performance though? (TODO)) */
1169                 dm = CDDM_copy(subdm);
1170                 subdm->release(subdm);
1171
1172                 bsphere_evolve(CDDM_get_verts(dm), dm->getNumVerts(dm), &bs->rootbase,
1173                                bs->skin_threshold, bs->evolve, bs->subdiv_level);
1174         }
1175
1176         return dm;
1177 }
1178 #endif