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