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