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