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