f33a1ef571ee958f526365e7e4738091e57b8ca3
[blender.git] / source / blender / bmesh / intern / bmesh_interp.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) 2007 Blender Foundation.
19  * All rights reserved.
20  *
21  * The Original Code is: all of this file.
22  *
23  * Contributor(s): Geoffrey Bantle.
24  *
25  * ***** END GPL LICENSE BLOCK *****
26  */
27
28 /** \file blender/bmesh/intern/bmesh_interp.c
29  *  \ingroup bmesh
30  *
31  * Functions for interpolating data across the surface of a mesh.
32  */
33
34 #include "MEM_guardedalloc.h"
35
36 #include "DNA_mesh_types.h"
37 #include "DNA_meshdata_types.h"
38
39 #include "BKE_customdata.h" 
40 #include "BKE_utildefines.h"
41 #include "BKE_multires.h"
42
43 #include "BLI_array.h"
44 #include "BLI_math.h"
45
46 #include "bmesh.h"
47 #include "bmesh_private.h"
48
49 /**
50  *                      bmesh_data_interp_from_verts
51  *
52  *  Interpolates per-vertex data from two sources to a target.
53  */
54 void BM_Data_Interp_From_Verts(BMesh *bm, BMVert *v1, BMVert *v2, BMVert *v, float fac)
55 {
56         if (v1->head.data && v2->head.data) {
57                 /* first see if we can avoid interpolation */
58                 if (fac <= 0.0f) {
59                         if (v1 == v) {
60                                 /* do nothing */
61                         }
62                         else {
63                                 CustomData_bmesh_free_block(&bm->vdata, &v->head.data);
64                                 CustomData_bmesh_copy_data(&bm->vdata, &bm->vdata, v1->head.data, &v->head.data);
65                         }
66                 }
67                 else if (fac >= 1.0f) {
68                         if (v2 == v) {
69                                 /* do nothing */
70                         }
71                         else {
72                                 CustomData_bmesh_free_block(&bm->vdata, &v->head.data);
73                                 CustomData_bmesh_copy_data(&bm->vdata, &bm->vdata, v2->head.data, &v->head.data);
74                         }
75                 }
76                 else {
77                         void *src[2];
78                         float w[2];
79
80                         src[0] = v1->head.data;
81                         src[1] = v2->head.data;
82                         w[0] = 1.0f-fac;
83                         w[1] = fac;
84                         CustomData_bmesh_interp(&bm->vdata, src, w, NULL, 2, v->head.data);
85                 }
86         }
87 }
88
89 /*
90  *  BM Data Vert Average
91  *
92  * Sets all the customdata (e.g. vert, loop) associated with a vert
93  * to the average of the face regions surrounding it.
94  */
95
96
97 static void UNUSED_FUNCTION(BM_Data_Vert_Average)(BMesh *UNUSED(bm), BMFace *UNUSED(f))
98 {
99         // BMIter iter;
100 }
101
102 /**
103  *                      bmesh_data_facevert_edgeinterp
104  *
105  *  Walks around the faces of an edge and interpolates the per-face-edge
106  *  data between two sources to a target.
107  * 
108  *  Returns -
109  *      Nothing
110  */
111  
112 void BM_Data_Facevert_Edgeinterp(BMesh *bm, BMVert *v1, BMVert *UNUSED(v2), BMVert *v, BMEdge *e1, float fac)
113 {
114         void *src[2];
115         float w[2];
116         BMLoop *l = NULL, *v1loop = NULL, *vloop = NULL, *v2loop = NULL;
117         
118         w[1] = 1.0f - fac;
119         w[0] = fac;
120
121         if (!e1->l) return;
122         l = e1->l;
123         do {
124                 if (l->v == v1) {
125                         v1loop = l;
126                         vloop = v1loop->next;
127                         v2loop = vloop->next;
128                 }
129                 else if (l->v == v) {
130                         v1loop = l->next;
131                         vloop = l;
132                         v2loop = l->prev;
133                 }
134                 
135                 if (!v1loop || !v2loop)
136                         return;
137                 
138                 src[0] = v1loop->head.data;
139                 src[1] = v2loop->head.data;                                     
140
141                 CustomData_bmesh_interp(&bm->ldata, src, w, NULL, 2, vloop->head.data);
142                 l = l->radial_next;
143         } while (l != e1->l);
144 }
145
146 void BM_loops_to_corners(BMesh *bm, Mesh *me, int findex,
147                          BMFace *f, int numTex, int numCol) 
148 {
149         BMLoop *l;
150         BMIter iter;
151         MTFace *texface;
152         MTexPoly *texpoly;
153         MCol *mcol;
154         MLoopCol *mloopcol;
155         MLoopUV *mloopuv;
156         int i, j;
157
158         for (i = 0; i < numTex; i++) {
159                 texface = CustomData_get_n(&me->fdata, CD_MTFACE, findex, i);
160                 texpoly = CustomData_bmesh_get_n(&bm->pdata, f->head.data, CD_MTEXPOLY, i);
161                 
162                 texface->tpage = texpoly->tpage;
163                 texface->flag = texpoly->flag;
164                 texface->transp = texpoly->transp;
165                 texface->mode = texpoly->mode;
166                 texface->tile = texpoly->tile;
167                 texface->unwrap = texpoly->unwrap;
168
169                 j = 0;
170                 BM_ITER(l, &iter, bm, BM_LOOPS_OF_FACE, f) {
171                         mloopuv = CustomData_bmesh_get_n(&bm->ldata, l->head.data, CD_MLOOPUV, i);
172                         copy_v2_v2(texface->uv[j], mloopuv->uv);
173
174                         j++;
175                 }
176
177         }
178
179         for (i = 0; i < numCol; i++) {
180                 mcol = CustomData_get_n(&me->fdata, CD_MCOL, findex, i);
181
182                 j = 0;
183                 BM_ITER(l, &iter, bm, BM_LOOPS_OF_FACE, f) {
184                         mloopcol = CustomData_bmesh_get_n(&bm->ldata, l->head.data, CD_MLOOPCOL, i);
185                         mcol[j].r = mloopcol->r;
186                         mcol[j].g = mloopcol->g;
187                         mcol[j].b = mloopcol->b;
188                         mcol[j].a = mloopcol->a;
189
190                         j++;
191                 }
192         }
193 }
194
195 /**
196  *                      BM_data_interp_from_face
197  *
198  *  projects target onto source, and pulls interpolated customdata from
199  *  source.
200  * 
201  *  Returns -
202  *      Nothing
203  */
204 void BM_face_interp_from_face(BMesh *bm, BMFace *target, BMFace *source)
205 {
206         BMLoop *l_iter;
207         BMLoop *l_first;
208
209         void **blocks = NULL;
210         float (*cos)[3] = NULL, *w = NULL;
211         BLI_array_fixedstack_declare(cos,     BM_NGON_STACK_SIZE, source->len, __func__);
212         BLI_array_fixedstack_declare(w,       BM_NGON_STACK_SIZE, source->len, __func__);
213         BLI_array_fixedstack_declare(blocks,  BM_NGON_STACK_SIZE, source->len, __func__);
214         int i;
215         
216         BM_Copy_Attributes(bm, bm, source, target);
217
218         i = 0;
219         l_iter = l_first = BM_FACE_FIRST_LOOP(source);
220         do {
221                 copy_v3_v3(cos[i], l_iter->v->co);
222                 blocks[i] = l_iter->head.data;
223                 i++;
224         } while ((l_iter = l_iter->next) != l_first);
225
226         i = 0;
227         l_iter = l_first = BM_FACE_FIRST_LOOP(target);
228         do {
229                 interp_weights_poly_v3(w, cos, source->len, l_iter->v->co);
230                 CustomData_bmesh_interp(&bm->ldata, blocks, w, NULL, source->len, l_iter->head.data);
231                 i++;
232         } while ((l_iter = l_iter->next) != l_first);
233
234         BLI_array_fixedstack_free(cos);
235         BLI_array_fixedstack_free(w);
236         BLI_array_fixedstack_free(blocks);
237 }
238
239 /* some math stuff for dealing with doubles, put here to
240  * avoid merge errors - joeedh */
241
242 #define VECMUL(a, b) (((a)[0] = (a)[0] * (b)), ((a)[1] = (a)[1] * (b)), ((a)[2] = (a)[2] * (b)))
243 #define VECADD2(a, b) (((a)[0] = (a)[0] + (b)[0]), ((a)[1] = (a)[1] + (b)[1]), ((a)[2] = (a)[2] + (b)[2]))
244 #define VECSUB2(a, b) (((a)[0] = (a)[0] - (b)[0]), ((a)[1] = (a)[1] - (b)[1]), ((a)[2] = (a)[2] - (b)[2]))
245
246 /* find closest point to p on line through l1, l2 and return lambda,
247  * where (0 <= lambda <= 1) when cp is in the line segement l1, l2
248  */
249 static double closest_to_line_v3_d(double cp[3], const double p[3], const double l1[3], const double l2[3])
250 {
251         double h[3], u[3], lambda;
252         VECSUB(u, l2, l1);
253         VECSUB(h, p, l1);
254         lambda = INPR(u, h) / INPR(u, u);
255         cp[0] = l1[0] + u[0] * lambda;
256         cp[1] = l1[1] + u[1] * lambda;
257         cp[2] = l1[2] + u[2] * lambda;
258         return lambda;
259 }
260
261 /* point closest to v1 on line v2-v3 in 3D */
262 static void UNUSED_FUNCTION(closest_to_line_segment_v3_d)(double *closest, double v1[3], double v2[3], double v3[3])
263 {
264         double lambda, cp[3];
265
266         lambda = closest_to_line_v3_d(cp, v1, v2, v3);
267
268         if (lambda <= 0.0) {
269                 VECCOPY(closest, v2);
270         }
271         else if (lambda >= 1.0) {
272                 VECCOPY(closest, v3);
273         }
274         else {
275                 VECCOPY(closest, cp);
276         }
277 }
278
279 static double UNUSED_FUNCTION(len_v3_d)(const double a[3])
280 {
281         return sqrt(INPR(a, a));
282 }
283
284 static double UNUSED_FUNCTION(len_v3v3_d)(const double a[3], const double b[3])
285 {
286         double d[3];
287
288         VECSUB(d, b, a);
289         return sqrt(INPR(d, d));
290 }
291
292 static void cent_quad_v3_d(double *cent, double *v1, double *v2, double *v3, double *v4)
293 {
294         cent[0] = 0.25 * (v1[0] + v2[0] + v3[0] + v4[0]);
295         cent[1] = 0.25 * (v1[1] + v2[1] + v3[1] + v4[1]);
296         cent[2] = 0.25 * (v1[2] + v2[2] + v3[2] + v4[2]);
297 }
298
299 static void UNUSED_FUNCTION(cent_tri_v3_d)(double *cent, double *v1, double *v2, double *v3)
300 {
301         cent[0] = 0.33333 * (v1[0] + v2[0] + v3[0]);
302         cent[1] = 0.33333 * (v1[1] + v2[1] + v3[1]);
303         cent[2] = 0.33333 * (v1[2] + v2[2] + v3[2]);
304 }
305
306 static void UNUSED_FUNCTION(cross_v3_v3v3_d)(double r[3], const double a[3], const double b[3])
307 {
308         r[0] = a[1] * b[2] - a[2] * b[1];
309         r[1] = a[2] * b[0] - a[0] * b[2];
310         r[2] = a[0] * b[1] - a[1] * b[0];
311 }
312
313 /* distance v1 to line-piece v2-v3 */
314 static double UNUSED_FUNCTION(dist_to_line_segment_v2_d)(double v1[3], double v2[3], double v3[3])
315 {
316         double labda, rc[2], pt[2], len;
317         
318         rc[0] = v3[0] - v2[0];
319         rc[1] = v3[1] - v2[1];
320         len = rc[0] * rc[0] + rc[1] * rc[1];
321         if (len == 0.0) {
322                 rc[0] = v1[0] - v2[0];
323                 rc[1] = v1[1] - v2[1];
324                 return sqrt(rc[0] * rc[0] + rc[1] * rc[1]);
325         }
326         
327         labda = (rc[0] * (v1[0] - v2[0]) + rc[1] * (v1[1] - v2[1])) / len;
328         if (labda <= 0.0) {
329                 pt[0] = v2[0];
330                 pt[1] = v2[1];
331         }
332         else if (labda >= 1.0) {
333                 pt[0] = v3[0];
334                 pt[1] = v3[1];
335         }
336         else {
337                 pt[0] = labda * rc[0] + v2[0];
338                 pt[1] = labda * rc[1] + v2[1];
339         }
340
341         rc[0] = pt[0] - v1[0];
342         rc[1] = pt[1] - v1[1];
343         return sqrt(rc[0] * rc[0] + rc[1] * rc[1]);
344 }
345
346
347 MINLINE double line_point_side_v2_d(const double *l1, const double *l2, const double *pt)
348 {
349         return  ((l1[0] - pt[0]) * (l2[1] - pt[1])) -
350                         ((l2[0] - pt[0]) * (l1[1] - pt[1]));
351 }
352
353 /* point in quad - only convex quads */
354 static int isect_point_quad_v2_d(double pt[2], double v1[2], double v2[2], double v3[2], double v4[2])
355 {
356         if (line_point_side_v2_d(v1, v2, pt) >= 0.0) {
357                 if (line_point_side_v2_d(v2, v3, pt) >= 0.0) {
358                         if (line_point_side_v2_d(v3, v4, pt) >= 0.0) {
359                                 if (line_point_side_v2_d(v4, v1, pt) >= 0.0) {
360                                         return 1;
361                                 }
362                         }
363                 }
364         }
365         else {
366                 if (! (line_point_side_v2_d(v2, v3, pt) >= 0.0)) {
367                         if (! (line_point_side_v2_d(v3, v4, pt) >= 0.0)) {
368                                 if (! (line_point_side_v2_d(v4, v1, pt) >= 0.0)) {
369                                         return -1;
370                                 }
371                         }
372                 }
373         }
374         
375         return 0;
376 }
377
378 /***** multires interpolation*****
379
380 mdisps is a grid of displacements, ordered thus:
381
382  v1/center----v4/next -> x
383      |           |
384      |           |
385   v2/prev------v3/cur
386      |
387      V
388      y
389  */
390
391 static int compute_mdisp_quad(BMLoop *l, double v1[3], double v2[3], double v3[3], double v4[3],
392                               double e1[3], double e2[3])
393 {
394         double cent[3] = {0.0, 0.0, 0.0}, n[3], p[3];
395         BMLoop *l_first;
396         BMLoop *l_iter;
397         
398         /* computer center */
399         l_iter = l_first = BM_FACE_FIRST_LOOP(l->f);
400         do {
401                 cent[0] += (double)l_iter->v->co[0];
402                 cent[1] += (double)l_iter->v->co[1];
403                 cent[2] += (double)l_iter->v->co[2];
404         } while ((l_iter = l_iter->next) != l_first);
405         
406         VECMUL(cent, (1.0 / (double)l->f->len));
407         
408         VECADD(p, l->prev->v->co, l->v->co);
409         VECMUL(p, 0.5);
410         VECADD(n, l->next->v->co, l->v->co);
411         VECMUL(n, 0.5);
412         
413         VECCOPY(v1, cent);
414         VECCOPY(v2, p);
415         VECCOPY(v3, l->v->co);
416         VECCOPY(v4, n);
417         
418         VECSUB(e1, v2, v1);
419         VECSUB(e2, v3, v4);
420         
421         return 1;
422 }
423
424 /* funnily enough, I think this is identical to face_to_crn_interp, heh */
425 static double quad_coord(double aa[3], double bb[3], double cc[3], double dd[3], int a1, int a2)
426 {
427         double x, y, z, f1;
428         
429         x = aa[a1] * cc[a2] - cc[a1] * aa[a2];
430         y = aa[a1] * dd[a2] + bb[a1] * cc[a2] - cc[a1] * bb[a2] - dd[a1] * aa[a2];
431         z = bb[a1] * dd[a2] - dd[a1] * bb[a2];
432         
433         if (fabs(2 * (x - y + z)) > DBL_EPSILON * 10.0) {
434                 double f2;
435
436                 f1 = (sqrt(y * y - 4.0 * x * z) - y + 2.0 * z) / (2.0 * (x - y + z));
437                 f2 = (-sqrt(y * y - 4.0 * x * z) - y + 2.0 * z) / (2.0 * (x - y + z));
438
439                 f1 = fabs(f1);
440                 f2 = fabs(f2);
441                 f1 = MIN2(f1, f2);
442                 CLAMP(f1, 0.0, 1.0 + DBL_EPSILON);
443         }
444         else {
445                 f1 = -z / (y - 2 * z);
446                 CLAMP(f1, 0.0, 1.0 + DBL_EPSILON);
447                 
448                 if (isnan(f1) || f1 > 1.0 || f1 < 0.0) {
449                         int i;
450                         
451                         for (i = 0; i < 2; i++) {
452                                 if (fabsf(aa[i]) < FLT_EPSILON * 100)
453                                         return aa[(i + 1) % 2] / fabs(bb[(i + 1) % 2] - aa[(i + 1) % 2]);
454                                 if (fabsf(cc[i]) < FLT_EPSILON * 100)
455                                         return cc[(i + 1) % 2] / fabs(dd[(i + 1) % 2] - cc[(i + 1) % 2]);
456                         }
457                 }
458         }
459
460         return f1;
461 }
462
463 static int quad_co(double *x, double *y, double v1[3], double v2[3], double v3[3], double v4[3],
464                    double p[3], float n[3])
465 {
466         float projverts[5][3], n2[3];
467         double dprojverts[4][3], origin[3] = {0.0f, 0.0f, 0.0f};
468         int i;
469
470         /* project points into 2d along normal */
471         VECCOPY(projverts[0], v1);
472         VECCOPY(projverts[1], v2);
473         VECCOPY(projverts[2], v3);
474         VECCOPY(projverts[3], v4);
475         VECCOPY(projverts[4], p);
476
477         normal_quad_v3(n2, projverts[0], projverts[1], projverts[2], projverts[3]);
478
479         if (INPR(n, n2) < -FLT_EPSILON) {
480                 return 0;
481         }
482
483         /* rotate */
484         poly_rotate_plane(n, projverts, 5);
485         
486         /* flatten */
487         for (i = 0; i < 5; i++) {
488                 projverts[i][2] = 0.0f;
489         }
490         
491         /* subtract origin */
492         for (i = 0; i < 4; i++) {
493                 VECSUB2(projverts[i], projverts[4]);
494         }
495         
496         VECCOPY(dprojverts[0], projverts[0]);
497         VECCOPY(dprojverts[1], projverts[1]);
498         VECCOPY(dprojverts[2], projverts[2]);
499         VECCOPY(dprojverts[3], projverts[3]);
500
501         if (!isect_point_quad_v2_d(origin, dprojverts[0], dprojverts[1], dprojverts[2], dprojverts[3])) {
502                 return 0;
503         }
504         
505         *y = quad_coord(dprojverts[1], dprojverts[0], dprojverts[2], dprojverts[3], 0, 1);
506         *x = quad_coord(dprojverts[2], dprojverts[1], dprojverts[3], dprojverts[0], 0, 1);
507
508         return 1;
509 }
510
511
512 /* tl is loop to project onto, l is loop whose internal displacement, co, is being
513  * projected.  x and y are location in loop's mdisps grid of point co. */
514 static int mdisp_in_mdispquad(BMesh *bm, BMLoop *l, BMLoop *tl, double p[3], double *x, double *y, int res)
515 {
516         double v1[3], v2[3], c[3], v3[3], v4[3], e1[3], e2[3];
517         double eps = FLT_EPSILON * 4000;
518         
519         if (len_v3(l->v->no) == 0.0f)
520                 BM_Vert_UpdateAllNormals(bm, l->v);
521         if (len_v3(tl->v->no) == 0.0f)
522                 BM_Vert_UpdateAllNormals(bm, tl->v);
523                 
524         compute_mdisp_quad(tl, v1, v2, v3, v4, e1, e2);
525
526         /* expand quad a bit */
527         cent_quad_v3_d(c, v1, v2, v3, v4);
528         
529         VECSUB2(v1, c); VECSUB2(v2, c);
530         VECSUB2(v3, c); VECSUB2(v4, c);
531         VECMUL(v1, 1.0 + eps); VECMUL(v2, 1.0 + eps);
532         VECMUL(v3, 1.0 + eps); VECMUL(v4, 1.0 + eps);
533         VECADD2(v1, c); VECADD2(v2, c);
534         VECADD2(v3, c); VECADD2(v4, c);
535         
536         if (!quad_co(x, y, v1, v2, v3, v4, p, l->v->no))
537                 return 0;
538         
539         *x *= res - 1;
540         *y *= res - 1;
541                   
542         return 1;
543 }
544
545 static void bmesh_loop_interp_mdisps(BMesh *bm, BMLoop *target, BMFace *source)
546 {
547         MDisps *mdisps;
548         BMLoop *l_iter;
549         BMLoop *l_first;
550         double x, y, d, v1[3], v2[3], v3[3], v4[3] = {0.0f, 0.0f, 0.0f}, e1[3], e2[3];
551         int ix, iy, res;
552         
553         /* ignore 2-edged faces */
554         if (target->f->len < 3)
555                 return;
556         
557         if (!CustomData_has_layer(&bm->ldata, CD_MDISPS))
558                 return;
559         
560         mdisps = CustomData_bmesh_get(&bm->ldata, target->head.data, CD_MDISPS);
561         compute_mdisp_quad(target, v1, v2, v3, v4, e1, e2);
562         
563         /* if no disps data allocate a new grid, the size of the first grid in source. */
564         if (!mdisps->totdisp) {
565                 MDisps *md2 = CustomData_bmesh_get(&bm->ldata, BM_FACE_FIRST_LOOP(source)->head.data, CD_MDISPS);
566                 
567                 mdisps->totdisp = md2->totdisp;
568                 if (mdisps->totdisp) {
569                         mdisps->disps = MEM_callocN(sizeof(float) * 3 * mdisps->totdisp, "mdisp->disps in bmesh_loop_intern_mdisps");
570                 }
571                 else {
572                         return;
573                 }
574         }
575         
576         res = (int)sqrt(mdisps->totdisp);
577         d = 1.0 / (double)(res - 1);
578         for (x = 0.0f, ix = 0; ix < res; x += d, ix++) {
579                 for (y = 0.0f, iy = 0; iy < res; y += d, iy++) {
580                         double co1[3], co2[3], co[3];
581                         /* double xx, yy; */ /* UNUSED */
582                         
583                         VECCOPY(co1, e1);
584                         
585                         /* if (!iy) yy = y + (double)FLT_EPSILON * 20; */
586                         /* else yy = y - (double)FLT_EPSILON * 20; */
587                         
588                         VECMUL(co1, y);
589                         VECADD2(co1, v1);
590                         
591                         VECCOPY(co2, e2);
592                         VECMUL(co2, y);
593                         VECADD2(co2, v4);
594                         
595                         /* if (!ix) xx = x + (double)FLT_EPSILON * 20; */ /* UNUSED */
596                         /* else xx = x - (double)FLT_EPSILON * 20; */ /* UNUSED */
597                         
598                         VECSUB(co, co2, co1);
599                         VECMUL(co, x);
600                         VECADD2(co, co1);
601                         
602                         l_iter = l_first = BM_FACE_FIRST_LOOP(source);
603                         do {
604                                 double x2, y2;
605                                 MDisps *md1, *md2;
606
607                                 md1 = CustomData_bmesh_get(&bm->ldata, target->head.data, CD_MDISPS);
608                                 md2 = CustomData_bmesh_get(&bm->ldata, l_iter->head.data, CD_MDISPS);
609                                 
610                                 if (mdisp_in_mdispquad(bm, target, l_iter, co, &x2, &y2, res)) {
611                                         /* int ix2 = (int)x2; */ /* UNUSED */
612                                         /* int iy2 = (int)y2; */ /* UNUSED */
613                                         
614                                         old_mdisps_bilinear(md1->disps[iy * res + ix], md2->disps, res, (float)x2, (float)y2);
615                                 }
616                         } while ((l_iter = l_iter->next) != l_first);
617                 }
618         }
619 }
620
621 void BM_multires_smooth_bounds(BMesh *bm, BMFace *f)
622 {
623         BMLoop *l;
624         BMIter liter;
625         
626         if (!CustomData_has_layer(&bm->ldata, CD_MDISPS))
627                 return;
628         
629         BM_ITER(l, &liter, bm, BM_LOOPS_OF_FACE, f) {
630                 MDisps *mdp = CustomData_bmesh_get(&bm->ldata, l->prev->head.data, CD_MDISPS);
631                 MDisps *mdl = CustomData_bmesh_get(&bm->ldata, l->head.data, CD_MDISPS);
632                 MDisps *mdn = CustomData_bmesh_get(&bm->ldata, l->next->head.data, CD_MDISPS);
633                 float co1[3];
634                 int sides;
635                 int y;
636                 
637                 /*
638                  *  mdisps is a grid of displacements, ordered thus:
639                  *  
640                  *                     v4/next
641                  *                       |                              
642                  *   |      v1/cent-----mid2 ---> x
643                  *   |         |         | 
644                  *   |         |         |
645                  *  v2/prev---mid1-----v3/cur
646                  *             |
647                  *             V
648                  *             y
649                  */
650                   
651                 sides = (int)sqrt(mdp->totdisp);
652                 for (y = 0; y < sides; y++) {
653                         add_v3_v3v3(co1, mdn->disps[y * sides], mdl->disps[y]);
654                         mul_v3_fl(co1, 0.5);
655
656                         copy_v3_v3(mdn->disps[y * sides], co1);
657                         copy_v3_v3(mdl->disps[y], co1);
658                 }
659         }
660         
661         BM_ITER(l, &liter, bm, BM_LOOPS_OF_FACE, f) {
662                 MDisps *mdl1 = CustomData_bmesh_get(&bm->ldata, l->head.data, CD_MDISPS);
663                 MDisps *mdl2;
664                 float co1[3], co2[3], co[3];
665                 int sides;
666                 int y;
667                 
668                 /*
669                  *  mdisps is a grid of displacements, ordered thus:
670                  *  
671                  *                     v4/next
672                  *                       |                              
673                  *   |      v1/cent-----mid2 ---> x
674                  *   |         |         | 
675                  *   |         |         |
676                  *  v2/prev---mid1-----v3/cur
677                  *             |
678                  *             V
679                  *             y
680                  */
681                  
682                 if (l->radial_next == l)
683                         continue;
684
685                 if (l->radial_next->v == l->v)
686                         mdl2 = CustomData_bmesh_get(&bm->ldata, l->radial_next->head.data, CD_MDISPS);
687                 else
688                         mdl2 = CustomData_bmesh_get(&bm->ldata, l->radial_next->next->head.data, CD_MDISPS);
689                         
690                 sides = (int)sqrt(mdl1->totdisp);
691                 for (y = 0; y < sides; y++) {
692                         int a1, a2, o1, o2;
693                         
694                         if (l->v != l->radial_next->v) {
695                                 a1 = sides * y + sides - 2;
696                                 a2 = (sides - 2) * sides + y;
697                                 
698                                 o1 = sides * y + sides - 1;
699                                 o2 = (sides - 1) * sides + y;
700                         }
701                         else {
702                                 a1 = sides * y + sides - 2;
703                                 a2 = sides * y + sides - 2;
704                                 o1 = sides * y + sides - 1;
705                                 o2 = sides * y + sides - 1;
706                         }
707                         
708                         /* magic blending numbers, hardcoded! */
709                         add_v3_v3v3(co1, mdl1->disps[a1], mdl2->disps[a2]);
710                         mul_v3_fl(co1, 0.18);
711                         
712                         add_v3_v3v3(co2, mdl1->disps[o1], mdl2->disps[o2]);
713                         mul_v3_fl(co2, 0.32);
714                         
715                         add_v3_v3v3(co, co1, co2);
716                         
717                         copy_v3_v3(mdl1->disps[o1], co);
718                         copy_v3_v3(mdl2->disps[o2], co);
719                 }
720         }
721 }
722
723 void BM_loop_interp_multires(BMesh *bm, BMLoop *target, BMFace *source)
724 {
725         bmesh_loop_interp_mdisps(bm, target, source);
726 }
727
728 void BM_loop_interp_from_face(BMesh *bm, BMLoop *target, BMFace *source, 
729                               int do_vertex, int do_multires)
730 {
731         BMLoop *l_iter;
732         BMLoop *l_first;
733         void **blocks = NULL;
734         void **vblocks = NULL;
735         float (*cos)[3] = NULL, co[3], *w = NULL, cent[3] = {0.0f, 0.0f, 0.0f};
736         BLI_array_fixedstack_declare(cos,      BM_NGON_STACK_SIZE, source->len, __func__);
737         BLI_array_fixedstack_declare(w,        BM_NGON_STACK_SIZE, source->len, __func__);
738         BLI_array_fixedstack_declare(blocks,   BM_NGON_STACK_SIZE, source->len, __func__);
739         BLI_array_fixedstack_declare(vblocks,  BM_NGON_STACK_SIZE, do_vertex ? source->len : 0, __func__);
740         int i, ax, ay;
741         
742         BM_Copy_Attributes(bm, bm, source, target->f);
743
744         i = 0;
745         l_iter = l_first = BM_FACE_FIRST_LOOP(source);
746         do {
747                 copy_v3_v3(cos[i], l_iter->v->co);
748                 add_v3_v3(cent, cos[i]);
749                 
750                 w[i] = 0.0f;
751                 blocks[i] = l_iter->head.data;
752         
753                 if (do_vertex) {
754                         vblocks[i] = l_iter->v->head.data;
755                 }
756                 i++;
757         
758         } while ((l_iter = l_iter->next) != l_first);
759
760         /* find best projection of face XY, XZ or YZ: barycentric weights of
761          * the 2d projected coords are the same and faster to compute */
762
763         axis_dominant_v3(&ax, &ay, source->no);
764
765         /* scale source face coordinates a bit, so points sitting directonly on an
766          * edge will work. */
767         mul_v3_fl(cent, 1.0f / (float)source->len);
768         for (i = 0; i < source->len; i++) {
769                 float vec[3], tmp[3];
770                 sub_v3_v3v3(vec, cent, cos[i]);
771                 mul_v3_fl(vec, 0.001f);
772                 add_v3_v3(cos[i], vec);
773                 
774                 copy_v3_v3(tmp, cos[i]);
775                 cos[i][0] = tmp[ax];
776                 cos[i][1] = tmp[ay];
777                 cos[i][2] = 0.0;
778         }
779         
780         
781         /* interpolate */
782         co[0] = target->v->co[ax];
783         co[1] = target->v->co[ay];
784         co[2] = 0.0f;
785         
786         interp_weights_poly_v3(w, cos, source->len, co);
787         CustomData_bmesh_interp(&bm->ldata, blocks, w, NULL, source->len, target->head.data);
788         if (do_vertex) {
789                 CustomData_bmesh_interp(&bm->vdata, vblocks, w, NULL, source->len, target->v->head.data);
790                 BLI_array_fixedstack_free(vblocks);
791         }
792
793         BLI_array_fixedstack_free(cos);
794         BLI_array_fixedstack_free(w);
795         BLI_array_fixedstack_free(blocks);
796         
797         if (do_multires) {
798                 if (CustomData_has_layer(&bm->ldata, CD_MDISPS)) {
799                         bmesh_loop_interp_mdisps(bm, target, source);
800                 }
801         }
802 }
803
804
805 void BM_vert_interp_from_face(BMesh *bm, BMVert *v, BMFace *source)
806 {
807         BMLoop *l_iter;
808         BMLoop *l_first;
809         void **blocks = NULL;
810         float (*cos)[3] = NULL, *w = NULL, cent[3] = {0.0f, 0.0f, 0.0f};
811         BLI_array_fixedstack_declare(cos,      BM_NGON_STACK_SIZE, source->len, __func__);
812         BLI_array_fixedstack_declare(w,        BM_NGON_STACK_SIZE, source->len, __func__);
813         BLI_array_fixedstack_declare(blocks,   BM_NGON_STACK_SIZE, source->len, __func__);
814         int i;
815
816         i = 0;
817         l_iter = l_first = BM_FACE_FIRST_LOOP(source);
818         do {
819                 copy_v3_v3(cos[i], l_iter->v->co);
820                 add_v3_v3(cent, cos[i]);
821
822                 w[i] = 0.0f;
823                 blocks[i] = l_iter->v->head.data;
824                 i++;
825         } while ((l_iter = l_iter->next) != l_first);
826
827         /* scale source face coordinates a bit, so points sitting directonly on an
828      * edge will work. */
829         mul_v3_fl(cent, 1.0f / (float)source->len);
830         for (i = 0; i < source->len; i++) {
831                 float vec[3];
832                 sub_v3_v3v3(vec, cent, cos[i]);
833                 mul_v3_fl(vec, 0.01);
834                 add_v3_v3(cos[i], vec);
835         }
836         
837         /* interpolate */
838         interp_weights_poly_v3(w, cos, source->len, v->co);
839         CustomData_bmesh_interp(&bm->vdata, blocks, w, NULL, source->len, v->head.data);
840         
841         BLI_array_fixedstack_free(cos);
842         BLI_array_fixedstack_free(w);
843         BLI_array_fixedstack_free(blocks);
844 }
845
846 static void update_data_blocks(BMesh *bm, CustomData *olddata, CustomData *data)
847 {
848         BMIter iter;
849         BLI_mempool *oldpool = olddata->pool;
850         void *block;
851
852         CustomData_bmesh_init_pool(data, data == &bm->ldata ? 2048 : 512);
853
854         if (data == &bm->vdata) {
855                 BMVert *eve;
856                 
857                 BM_ITER(eve, &iter, bm, BM_VERTS_OF_MESH, NULL) {
858                         block = NULL;
859                         CustomData_bmesh_set_default(data, &block);
860                         CustomData_bmesh_copy_data(olddata, data, eve->head.data, &block);
861                         CustomData_bmesh_free_block(olddata, &eve->head.data);
862                         eve->head.data = block;
863                 }
864         }
865         else if (data == &bm->edata) {
866                 BMEdge *eed;
867
868                 BM_ITER(eed, &iter, bm, BM_EDGES_OF_MESH, NULL) {
869                         block = NULL;
870                         CustomData_bmesh_set_default(data, &block);
871                         CustomData_bmesh_copy_data(olddata, data, eed->head.data, &block);
872                         CustomData_bmesh_free_block(olddata, &eed->head.data);
873                         eed->head.data = block;
874                 }
875         }
876         else if (data == &bm->pdata || data == &bm->ldata) {
877                 BMIter liter;
878                 BMFace *efa;
879                 BMLoop *l;
880
881                 BM_ITER(efa, &iter, bm, BM_FACES_OF_MESH, NULL) {
882                         if (data == &bm->pdata) {
883                                 block = NULL;
884                                 CustomData_bmesh_set_default(data, &block);
885                                 CustomData_bmesh_copy_data(olddata, data, efa->head.data, &block);
886                                 CustomData_bmesh_free_block(olddata, &efa->head.data);
887                                 efa->head.data = block;
888                         }
889
890                         if (data == &bm->ldata) {
891                                 BM_ITER(l, &liter, bm, BM_LOOPS_OF_FACE, efa) {
892                                         block = NULL;
893                                         CustomData_bmesh_set_default(data, &block);
894                                         CustomData_bmesh_copy_data(olddata, data, l->head.data, &block);
895                                         CustomData_bmesh_free_block(olddata, &l->head.data);
896                                         l->head.data = block;
897                                 }
898                         }
899                 }
900         }
901
902         if (oldpool) {
903                 /* this should never happen but can when dissolve fails - [#28960] */
904                 BLI_assert(data->pool != oldpool);
905
906                 BLI_mempool_destroy(oldpool);
907         }
908 }
909
910 void BM_add_data_layer(BMesh *bm, CustomData *data, int type)
911 {
912         CustomData olddata;
913
914         olddata = *data;
915         olddata.layers = (olddata.layers) ? MEM_dupallocN(olddata.layers): NULL;
916
917         /* the pool is now owned by olddata and must not be shared */
918         data->pool = NULL;
919
920         CustomData_add_layer(data, type, CD_DEFAULT, NULL, 0);
921
922         update_data_blocks(bm, &olddata, data);
923         if (olddata.layers) MEM_freeN(olddata.layers);
924 }
925
926 void BM_add_data_layer_named(BMesh *bm, CustomData *data, int type, const char *name)
927 {
928         CustomData olddata;
929
930         olddata = *data;
931         olddata.layers = (olddata.layers) ? MEM_dupallocN(olddata.layers): NULL;
932
933         /* the pool is now owned by olddata and must not be shared */
934         data->pool = NULL;
935
936         CustomData_add_layer_named(data, type, CD_DEFAULT, NULL, 0, name);
937
938         update_data_blocks(bm, &olddata, data);
939         if (olddata.layers) MEM_freeN(olddata.layers);
940 }
941
942 void BM_free_data_layer(BMesh *bm, CustomData *data, int type)
943 {
944         CustomData olddata;
945
946         olddata = *data;
947         olddata.layers = (olddata.layers) ? MEM_dupallocN(olddata.layers): NULL;
948
949         /* the pool is now owned by olddata and must not be shared */
950         data->pool = NULL;
951
952         CustomData_free_layer_active(data, type, 0);
953
954         update_data_blocks(bm, &olddata, data);
955         if (olddata.layers) MEM_freeN(olddata.layers);
956 }
957
958 void BM_free_data_layer_n(BMesh *bm, CustomData *data, int type, int n)
959 {
960         CustomData olddata;
961
962         olddata = *data;
963         olddata.layers = (olddata.layers) ? MEM_dupallocN(olddata.layers): NULL;
964
965         /* the pool is now owned by olddata and must not be shared */
966         data->pool = NULL;
967
968         CustomData_free_layer(data, type, 0, CustomData_get_layer_index_n(data, type, n));
969         
970         update_data_blocks(bm, &olddata, data);
971         if (olddata.layers) MEM_freeN(olddata.layers);
972 }
973
974 float BM_GetCDf(CustomData *cd, void *element, int type)
975 {
976         float *f = CustomData_bmesh_get(cd, ((BMHeader *)element)->data, type);
977         return f ? *f : 0.0f;
978 }
979
980 void BM_SetCDf(CustomData *cd, void *element, int type, float val)
981 {
982         float *f = CustomData_bmesh_get(cd, ((BMHeader *)element)->data, type);
983         if (f) *f = val;
984 }