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