Merge branch 'master' into blender2.8
[blender.git] / source / blender / modifiers / intern / MOD_surfacedeform.c
1 #include "DNA_mesh_types.h"
2 #include "DNA_meshdata_types.h"
3 #include "DNA_object_types.h"
4 #include "DNA_scene_types.h"
5
6 #include "BLI_alloca.h"
7 #include "BLI_math.h"
8 #include "BLI_math_geom.h"
9 #include "BLI_task.h"
10
11 #include "BKE_bvhutils.h"
12 #include "BKE_mesh.h"
13 #include "BKE_editmesh.h"
14 #include "BKE_library.h"
15 #include "BKE_library_query.h"
16 #include "BKE_modifier.h"
17
18 #include "DEG_depsgraph.h"
19
20 #include "MEM_guardedalloc.h"
21
22 #include "MOD_util.h"
23
24 typedef struct SDefAdjacency {
25         struct SDefAdjacency *next;
26         unsigned int index;
27 } SDefAdjacency;
28
29 typedef struct SDefAdjacencyArray {
30         SDefAdjacency *first;
31         unsigned int num; /* Careful, this is twice the number of polygons (avoids an extra loop) */
32 } SDefAdjacencyArray;
33
34 typedef struct SDefEdgePolys {
35         unsigned int polys[2], num;
36 } SDefEdgePolys;
37
38 typedef struct SDefBindCalcData {
39         BVHTreeFromMesh * const treeData;
40         const SDefAdjacencyArray * const vert_edges;
41         const SDefEdgePolys * const edge_polys;
42         SDefVert * const bind_verts;
43         const MLoopTri * const looptri;
44         const MPoly * const mpoly;
45         const MEdge * const medge;
46         const MLoop * const mloop;
47         float (* const targetCos)[3];
48         float (* const vertexCos)[3];
49         float imat[4][4];
50         const float falloff;
51         int success;
52 } SDefBindCalcData;
53
54 typedef struct SDefBindPoly {
55         float (*coords)[3];
56         float (*coords_v2)[2];
57         float point_v2[2];
58         float weight_angular;
59         float weight_dist_proj;
60         float weight_dist;
61         float weight;
62         float scales[2];
63         float centroid[3];
64         float centroid_v2[2];
65         float normal[3];
66         float cent_edgemid_vecs_v2[2][2];
67         float edgemid_angle;
68         float point_edgemid_angles[2];
69         float corner_edgemid_angles[2];
70         float dominant_angle_weight;
71         unsigned int index;
72         unsigned int numverts;
73         unsigned int loopstart;
74         unsigned int edge_inds[2];
75         unsigned int edge_vert_inds[2];
76         unsigned int corner_ind;
77         unsigned int dominant_edge;
78         bool inside;
79 } SDefBindPoly;
80
81 typedef struct SDefBindWeightData {
82         SDefBindPoly *bind_polys;
83         unsigned int numpoly;
84         unsigned int numbinds;
85 } SDefBindWeightData;
86
87 typedef struct SDefDeformData {
88         const SDefVert * const bind_verts;
89         float (* const targetCos)[3];
90         float (* const vertexCos)[3];
91 } SDefDeformData;
92
93 /* Bind result values */
94 enum {
95         MOD_SDEF_BIND_RESULT_SUCCESS = 1,
96         MOD_SDEF_BIND_RESULT_GENERIC_ERR = 0,
97         MOD_SDEF_BIND_RESULT_MEM_ERR = -1,
98         MOD_SDEF_BIND_RESULT_NONMANY_ERR = -2,
99         MOD_SDEF_BIND_RESULT_CONCAVE_ERR = -3,
100         MOD_SDEF_BIND_RESULT_OVERLAP_ERR = -4,
101 };
102
103 /* Infinite weight flags */
104 enum {
105         MOD_SDEF_INFINITE_WEIGHT_ANGULAR = (1 << 0),
106         MOD_SDEF_INFINITE_WEIGHT_DIST_PROJ = (1 << 1),
107         MOD_SDEF_INFINITE_WEIGHT_DIST = (1 << 2),
108 };
109
110 static void initData(ModifierData *md)
111 {
112         SurfaceDeformModifierData *smd = (SurfaceDeformModifierData *)md;
113         smd->target = NULL;
114         smd->verts = NULL;
115         smd->flags = 0;
116         smd->falloff = 4.0f;
117 }
118
119 static void freeData(ModifierData *md)
120 {
121         SurfaceDeformModifierData *smd = (SurfaceDeformModifierData *)md;
122
123         if (smd->verts) {
124                 for (int i = 0; i < smd->numverts; i++) {
125                         if (smd->verts[i].binds) {
126                                 for (int j = 0; j < smd->verts[i].numbinds; j++) {
127                                         MEM_SAFE_FREE(smd->verts[i].binds[j].vert_inds);
128                                         MEM_SAFE_FREE(smd->verts[i].binds[j].vert_weights);
129                                 }
130
131                                 MEM_SAFE_FREE(smd->verts[i].binds);
132                         }
133                 }
134
135                 MEM_SAFE_FREE(smd->verts);
136         }
137 }
138
139 static void copyData(const ModifierData *md, ModifierData *target)
140 {
141         const SurfaceDeformModifierData *smd = (const SurfaceDeformModifierData *)md;
142         SurfaceDeformModifierData *tsmd = (SurfaceDeformModifierData *)target;
143
144         modifier_copyData_generic(md, target);
145
146         if (smd->verts) {
147                 tsmd->verts = MEM_dupallocN(smd->verts);
148
149                 for (int i = 0; i < smd->numverts; i++) {
150                         if (smd->verts[i].binds) {
151                                 tsmd->verts[i].binds = MEM_dupallocN(smd->verts[i].binds);
152
153                                 for (int j = 0; j < smd->verts[i].numbinds; j++) {
154                                         if (smd->verts[i].binds[j].vert_inds) {
155                                                 tsmd->verts[i].binds[j].vert_inds = MEM_dupallocN(smd->verts[i].binds[j].vert_inds);
156                                         }
157
158                                         if (smd->verts[i].binds[j].vert_weights) {
159                                                 tsmd->verts[i].binds[j].vert_weights = MEM_dupallocN(smd->verts[i].binds[j].vert_weights);
160                                         }
161                                 }
162                         }
163                 }
164         }
165 }
166
167 static void foreachObjectLink(ModifierData *md, Object *ob, ObjectWalkFunc walk, void *userData)
168 {
169         SurfaceDeformModifierData *smd = (SurfaceDeformModifierData *)md;
170
171         walk(userData, ob, &smd->target, IDWALK_NOP);
172 }
173
174 static void updateDepsgraph(ModifierData *md, const ModifierUpdateDepsgraphContext *ctx)
175 {
176         SurfaceDeformModifierData *smd = (SurfaceDeformModifierData *)md;
177         if (smd->target != NULL) {
178                 DEG_add_object_relation(ctx->node, smd->target, DEG_OB_COMP_GEOMETRY, "Surface Deform Modifier");
179         }
180 }
181
182 static void freeAdjacencyMap(SDefAdjacencyArray * const vert_edges, SDefAdjacency * const adj_ref, SDefEdgePolys * const edge_polys)
183 {
184         MEM_freeN(edge_polys);
185
186         MEM_freeN(adj_ref);
187
188         MEM_freeN(vert_edges);
189 }
190
191 static int buildAdjacencyMap(
192         const MPoly *poly, const MEdge *edge, const MLoop * const mloop, const unsigned int numpoly, const unsigned int numedges,
193         SDefAdjacencyArray * const vert_edges, SDefAdjacency *adj, SDefEdgePolys * const edge_polys)
194 {
195         const MLoop *loop;
196
197         /* Fing polygons adjacent to edges */
198         for (int i = 0; i < numpoly; i++, poly++) {
199                 loop = &mloop[poly->loopstart];
200
201                 for (int j = 0; j < poly->totloop; j++, loop++) {
202                         if (edge_polys[loop->e].num == 0) {
203                                 edge_polys[loop->e].polys[0] = i;
204                                 edge_polys[loop->e].polys[1] = -1;
205                                 edge_polys[loop->e].num++;
206                         }
207                         else if (edge_polys[loop->e].num == 1) {
208                                 edge_polys[loop->e].polys[1] = i;
209                                 edge_polys[loop->e].num++;
210                         }
211                         else {
212                                 return MOD_SDEF_BIND_RESULT_NONMANY_ERR;
213                         }
214                 }
215         }
216
217         /* Find edges adjacent to vertices */
218         for (int i = 0; i < numedges; i++, edge++) {
219                 adj->next = vert_edges[edge->v1].first;
220                 adj->index = i;
221                 vert_edges[edge->v1].first = adj;
222                 vert_edges[edge->v1].num += edge_polys[i].num;
223                 adj++;
224
225                 adj->next = vert_edges[edge->v2].first;
226                 adj->index = i;
227                 vert_edges[edge->v2].first = adj;
228                 vert_edges[edge->v2].num += edge_polys[i].num;
229                 adj++;
230         }
231
232         return MOD_SDEF_BIND_RESULT_SUCCESS;
233 }
234
235 BLI_INLINE void sortPolyVertsEdge(unsigned int *indices, const MLoop * const mloop, const unsigned int edge, const unsigned int num)
236 {
237         bool found = false;
238
239         for (int i = 0; i < num; i++) {
240                 if (mloop[i].e == edge) {
241                         found = true;
242                 }
243                 if (found) {
244                         *indices = mloop[i].v;
245                         indices++;
246                 }
247         }
248
249         /* Fill in remaining vertex indices that occur before the edge */
250         for (int i = 0; mloop[i].e != edge; i++) {
251                 *indices = mloop[i].v;
252                 indices++;
253         }
254 }
255
256 BLI_INLINE void sortPolyVertsTri(unsigned int *indices, const MLoop * const mloop, const unsigned int loopstart, const unsigned int num)
257 {
258         for (int i = loopstart; i < num; i++) {
259                 *indices = mloop[i].v;
260                 indices++;
261         }
262
263         for (int i = 0; i < loopstart; i++) {
264                 *indices = mloop[i].v;
265                 indices++;
266         }
267 }
268
269 BLI_INLINE unsigned int nearestVert(SDefBindCalcData * const data, const float point_co[3])
270 {
271         BVHTreeNearest nearest = {.dist_sq = FLT_MAX, .index = -1};
272         const MPoly *poly;
273         const MEdge *edge;
274         const MLoop *loop;
275         float t_point[3];
276         float max_dist = FLT_MAX;
277         float dist;
278         unsigned int index = 0;
279
280         mul_v3_m4v3(t_point, data->imat, point_co);
281
282         BLI_bvhtree_find_nearest(data->treeData->tree, t_point, &nearest, data->treeData->nearest_callback, data->treeData);
283
284         poly = &data->mpoly[data->looptri[nearest.index].poly];
285         loop = &data->mloop[poly->loopstart];
286
287         for (int i = 0; i < poly->totloop; i++, loop++) {
288                 edge = &data->medge[loop->e];
289                 dist = dist_squared_to_line_segment_v3(point_co, data->targetCos[edge->v1], data->targetCos[edge->v2]);
290
291                 if (dist < max_dist) {
292                         max_dist = dist;
293                         index = loop->e;
294                 }
295         }
296
297         edge = &data->medge[index];
298         if (len_squared_v3v3(point_co, data->targetCos[edge->v1]) < len_squared_v3v3(point_co, data->targetCos[edge->v2])) {
299                 return edge->v1;
300         }
301         else {
302                 return edge->v2;
303         }
304 }
305
306 BLI_INLINE int isPolyValid(const float coords[][2], const unsigned int nr)
307 {
308         float prev_co[2];
309         float curr_vec[2], prev_vec[2];
310
311         if (!is_poly_convex_v2(coords, nr)) {
312                 return MOD_SDEF_BIND_RESULT_CONCAVE_ERR;
313         }
314
315         copy_v2_v2(prev_co, coords[nr - 1]);
316         sub_v2_v2v2(prev_vec, prev_co, coords[nr - 2]);
317         normalize_v2(prev_vec);
318
319         for (int i = 0; i < nr; i++) {
320                 sub_v2_v2v2(curr_vec, coords[i], prev_co);
321
322                 const float curr_len = normalize_v2(curr_vec);
323                 if (curr_len < FLT_EPSILON) {
324                         return MOD_SDEF_BIND_RESULT_OVERLAP_ERR;
325                 }
326
327                 if (1.0f - dot_v2v2(prev_vec, curr_vec) < FLT_EPSILON) {
328                         return MOD_SDEF_BIND_RESULT_CONCAVE_ERR;
329                 }
330
331                 copy_v2_v2(prev_co, coords[i]);
332                 copy_v2_v2(prev_vec, curr_vec);
333         }
334
335         return MOD_SDEF_BIND_RESULT_SUCCESS;
336 }
337
338 static void freeBindData(SDefBindWeightData * const bwdata)
339 {
340         SDefBindPoly *bpoly = bwdata->bind_polys;
341
342         if (bwdata->bind_polys) {
343                 for (int i = 0; i < bwdata->numpoly; bpoly++, i++) {
344                         MEM_SAFE_FREE(bpoly->coords);
345                         MEM_SAFE_FREE(bpoly->coords_v2);
346                 }
347
348                 MEM_freeN(bwdata->bind_polys);
349         }
350
351         MEM_freeN(bwdata);
352 }
353
354 BLI_INLINE float computeAngularWeight(const float point_angle, const float edgemid_angle)
355 {
356         float weight;
357
358         weight = point_angle;
359         weight /= edgemid_angle;
360         weight *= M_PI_2;
361
362         return sinf(weight);
363 }
364
365 BLI_INLINE SDefBindWeightData *computeBindWeights(SDefBindCalcData * const data, const float point_co[3])
366 {
367         const unsigned int nearest = nearestVert(data, point_co);
368         const SDefAdjacency * const vert_edges = data->vert_edges[nearest].first;
369         const SDefEdgePolys * const edge_polys = data->edge_polys;
370
371         const SDefAdjacency *vedge;
372         const MPoly *poly;
373         const MLoop *loop;
374
375         SDefBindWeightData *bwdata;
376         SDefBindPoly *bpoly;
377
378         float world[3] = {0.0f, 0.0f, 1.0f};
379         float avg_point_dist = 0.0f;
380         float tot_weight = 0.0f;
381         int inf_weight_flags = 0;
382
383         bwdata = MEM_callocN(sizeof(*bwdata), "SDefBindWeightData");
384         if (bwdata == NULL) {
385                 data->success = MOD_SDEF_BIND_RESULT_MEM_ERR;
386                 return NULL;
387         }
388
389         bwdata->numpoly = data->vert_edges[nearest].num / 2;
390
391         bpoly = MEM_calloc_arrayN(bwdata->numpoly, sizeof(*bpoly), "SDefBindPoly");
392         if (bpoly == NULL) {
393                 freeBindData(bwdata);
394                 data->success = MOD_SDEF_BIND_RESULT_MEM_ERR;
395                 return NULL;
396         }
397
398         bwdata->bind_polys = bpoly;
399
400         /* Loop over all adjacent edges, and build the SDefBindPoly data for each poly adjacent to those */
401         for (vedge = vert_edges; vedge; vedge = vedge->next) {
402                 unsigned int edge_ind = vedge->index;
403
404                 for (int i = 0; i < edge_polys[edge_ind].num; i++) {
405                         {
406                                 bpoly = bwdata->bind_polys;
407
408                                 for (int j = 0; j < bwdata->numpoly; bpoly++, j++) {
409                                         /* If coords isn't allocated, we have reached the first uninitialized bpoly */
410                                         if ((bpoly->index == edge_polys[edge_ind].polys[i]) || (!bpoly->coords)) {
411                                                 break;
412                                         }
413                                 }
414                         }
415
416                         /* Check if poly was already created by another edge or still has to be initialized */
417                         if (!bpoly->coords) {
418                                 float angle;
419                                 float axis[3];
420                                 float tmp_vec_v2[2];
421                                 int is_poly_valid;
422
423                                 bpoly->index = edge_polys[edge_ind].polys[i];
424                                 bpoly->coords = NULL;
425                                 bpoly->coords_v2 = NULL;
426
427                                 /* Copy poly data */
428                                 poly = &data->mpoly[bpoly->index];
429                                 loop = &data->mloop[poly->loopstart];
430
431                                 bpoly->numverts = poly->totloop;
432                                 bpoly->loopstart = poly->loopstart;
433
434                                 bpoly->coords = MEM_malloc_arrayN(poly->totloop, sizeof(*bpoly->coords), "SDefBindPolyCoords");
435                                 if (bpoly->coords == NULL) {
436                                         freeBindData(bwdata);
437                                         data->success = MOD_SDEF_BIND_RESULT_MEM_ERR;
438                                         return NULL;
439                                 }
440
441                                 bpoly->coords_v2 = MEM_malloc_arrayN(poly->totloop, sizeof(*bpoly->coords_v2), "SDefBindPolyCoords_v2");
442                                 if (bpoly->coords_v2 == NULL) {
443                                         freeBindData(bwdata);
444                                         data->success = MOD_SDEF_BIND_RESULT_MEM_ERR;
445                                         return NULL;
446                                 }
447
448                                 for (int j = 0; j < poly->totloop; j++, loop++) {
449                                         copy_v3_v3(bpoly->coords[j], data->targetCos[loop->v]);
450
451                                         /* Find corner and edge indices within poly loop array */
452                                         if (loop->v == nearest) {
453                                                 bpoly->corner_ind = j;
454                                                 bpoly->edge_vert_inds[0] = (j == 0) ? (poly->totloop - 1) : (j - 1);
455                                                 bpoly->edge_vert_inds[1] = (j == poly->totloop - 1) ? (0) : (j + 1);
456
457                                                 bpoly->edge_inds[0] = data->mloop[poly->loopstart + bpoly->edge_vert_inds[0]].e;
458                                                 bpoly->edge_inds[1] = loop->e;
459                                         }
460                                 }
461
462                                 /* Compute poly's parametric data */
463                                 mid_v3_v3_array(bpoly->centroid, bpoly->coords, poly->totloop);
464                                 normal_poly_v3(bpoly->normal, bpoly->coords, poly->totloop);
465
466                                 /* Compute poly skew angle and axis */
467                                 angle = angle_normalized_v3v3(bpoly->normal, world);
468
469                                 cross_v3_v3v3(axis, bpoly->normal, world);
470                                 normalize_v3(axis);
471
472                                 /* Map coords onto 2d normal plane */
473                                 map_to_plane_axis_angle_v2_v3v3fl(bpoly->point_v2, point_co, axis, angle);
474
475                                 zero_v2(bpoly->centroid_v2);
476                                 for (int j = 0; j < poly->totloop; j++) {
477                                         map_to_plane_axis_angle_v2_v3v3fl(bpoly->coords_v2[j], bpoly->coords[j], axis, angle);
478                                         madd_v2_v2fl(bpoly->centroid_v2, bpoly->coords_v2[j], 1.0f / poly->totloop);
479                                 }
480
481                                 is_poly_valid = isPolyValid(bpoly->coords_v2, poly->totloop);
482
483                                 if (is_poly_valid != MOD_SDEF_BIND_RESULT_SUCCESS) {
484                                         freeBindData(bwdata);
485                                         data->success = is_poly_valid;
486                                         return NULL;
487                                 }
488
489                                 bpoly->inside = isect_point_poly_v2(bpoly->point_v2, bpoly->coords_v2, poly->totloop, false);
490
491                                 /* Initialize weight components */
492                                 bpoly->weight_angular = 1.0f;
493                                 bpoly->weight_dist_proj = len_v2v2(bpoly->centroid_v2, bpoly->point_v2);
494                                 bpoly->weight_dist = len_v3v3(bpoly->centroid, point_co);
495
496                                 avg_point_dist += bpoly->weight_dist;
497
498                                 /* Compute centroid to mid-edge vectors */
499                                 mid_v2_v2v2(bpoly->cent_edgemid_vecs_v2[0],
500                                             bpoly->coords_v2[bpoly->edge_vert_inds[0]],
501                                             bpoly->coords_v2[bpoly->corner_ind]);
502
503                                 mid_v2_v2v2(bpoly->cent_edgemid_vecs_v2[1],
504                                             bpoly->coords_v2[bpoly->edge_vert_inds[1]],
505                                             bpoly->coords_v2[bpoly->corner_ind]);
506
507                                 sub_v2_v2(bpoly->cent_edgemid_vecs_v2[0], bpoly->centroid_v2);
508                                 sub_v2_v2(bpoly->cent_edgemid_vecs_v2[1], bpoly->centroid_v2);
509
510                                 /* Compute poly scales with respect to mid-edges, and normalize the vectors */
511                                 bpoly->scales[0] = normalize_v2(bpoly->cent_edgemid_vecs_v2[0]);
512                                 bpoly->scales[1] = normalize_v2(bpoly->cent_edgemid_vecs_v2[1]);
513
514                                 /* Compute the required polygon angles */
515                                 bpoly->edgemid_angle = angle_normalized_v2v2(bpoly->cent_edgemid_vecs_v2[0], bpoly->cent_edgemid_vecs_v2[1]);
516
517                                 sub_v2_v2v2(tmp_vec_v2, bpoly->coords_v2[bpoly->corner_ind], bpoly->centroid_v2);
518                                 normalize_v2(tmp_vec_v2);
519
520                                 bpoly->corner_edgemid_angles[0] = angle_normalized_v2v2(tmp_vec_v2, bpoly->cent_edgemid_vecs_v2[0]);
521                                 bpoly->corner_edgemid_angles[1] = angle_normalized_v2v2(tmp_vec_v2, bpoly->cent_edgemid_vecs_v2[1]);
522
523                                 /* Check for inifnite weights, and compute angular data otherwise */
524                                 if (bpoly->weight_dist < FLT_EPSILON) {
525                                         inf_weight_flags |= MOD_SDEF_INFINITE_WEIGHT_DIST_PROJ;
526                                         inf_weight_flags |= MOD_SDEF_INFINITE_WEIGHT_DIST;
527                                 }
528                                 else if (bpoly->weight_dist_proj < FLT_EPSILON) {
529                                         inf_weight_flags |= MOD_SDEF_INFINITE_WEIGHT_DIST_PROJ;
530                                 }
531                                 else {
532                                         float cent_point_vec[2];
533
534                                         sub_v2_v2v2(cent_point_vec, bpoly->point_v2, bpoly->centroid_v2);
535                                         normalize_v2(cent_point_vec);
536
537                                         bpoly->point_edgemid_angles[0] = angle_normalized_v2v2(cent_point_vec, bpoly->cent_edgemid_vecs_v2[0]);
538                                         bpoly->point_edgemid_angles[1] = angle_normalized_v2v2(cent_point_vec, bpoly->cent_edgemid_vecs_v2[1]);
539                                 }
540                         }
541                 }
542         }
543
544         avg_point_dist /= bwdata->numpoly;
545
546         /* If weights 1 and 2 are not infinite, loop over all adjacent edges again,
547          * and build adjacency dependent angle data (depends on all polygons having been computed) */
548         if (!inf_weight_flags) {
549                 for (vedge = vert_edges; vedge; vedge = vedge->next) {
550                         SDefBindPoly *bpolys[2];
551                         const SDefEdgePolys *epolys;
552                         float ang_weights[2];
553                         unsigned int edge_ind = vedge->index;
554                         unsigned int edge_on_poly[2];
555
556                         epolys = &edge_polys[edge_ind];
557
558                         /* Find bind polys corresponding to the edge's adjacent polys */
559                         bpoly = bwdata->bind_polys;
560
561                         for (int i = 0, j = 0; (i < bwdata->numpoly) && (j < epolys->num); bpoly++, i++) {
562                                 if (ELEM(bpoly->index, epolys->polys[0], epolys->polys[1])) {
563                                         bpolys[j] = bpoly;
564
565                                         if (bpoly->edge_inds[0] == edge_ind) {
566                                                 edge_on_poly[j] = 0;
567                                         }
568                                         else {
569                                                 edge_on_poly[j] = 1;
570                                         }
571
572                                         j++;
573                                 }
574                         }
575
576                         /* Compute angular weight component */
577                         if (epolys->num == 1) {
578                                 ang_weights[0] = computeAngularWeight(bpolys[0]->point_edgemid_angles[edge_on_poly[0]], bpolys[0]->edgemid_angle);
579                                 bpolys[0]->weight_angular *= ang_weights[0] * ang_weights[0];
580                         }
581                         else if (epolys->num == 2) {
582                                 ang_weights[0] = computeAngularWeight(bpolys[0]->point_edgemid_angles[edge_on_poly[0]], bpolys[0]->edgemid_angle);
583                                 ang_weights[1] = computeAngularWeight(bpolys[1]->point_edgemid_angles[edge_on_poly[1]], bpolys[1]->edgemid_angle);
584
585                                 bpolys[0]->weight_angular *= ang_weights[0] * ang_weights[1];
586                                 bpolys[1]->weight_angular *= ang_weights[0] * ang_weights[1];
587                         }
588                 }
589         }
590
591         /* Compute scalings and falloff.
592          * Scale all weights if no infinite weight is found,
593          * scale only unprojected weight if projected weight is infinite,
594          * scale none if both are infinite. */
595         if (!inf_weight_flags) {
596                 bpoly = bwdata->bind_polys;
597
598                 for (int i = 0; i < bwdata->numpoly; bpoly++, i++) {
599                         float corner_angle_weights[2];
600                         float scale_weight, sqr, inv_sqr;
601
602                         corner_angle_weights[0] = bpoly->point_edgemid_angles[0] / bpoly->corner_edgemid_angles[0];
603                         corner_angle_weights[1] = bpoly->point_edgemid_angles[1] / bpoly->corner_edgemid_angles[1];
604
605                         if (isnan(corner_angle_weights[0]) || isnan(corner_angle_weights[1])) {
606                                 freeBindData(bwdata);
607                                 data->success = MOD_SDEF_BIND_RESULT_GENERIC_ERR;
608                                 return NULL;
609                         }
610
611                         /* Find which edge the point is closer to */
612                         if (corner_angle_weights[0] < corner_angle_weights[1]) {
613                                 bpoly->dominant_edge = 0;
614                                 bpoly->dominant_angle_weight = corner_angle_weights[0];
615                         }
616                         else {
617                                 bpoly->dominant_edge = 1;
618                                 bpoly->dominant_angle_weight = corner_angle_weights[1];
619                         }
620
621                         bpoly->dominant_angle_weight = sinf(bpoly->dominant_angle_weight * M_PI_2);
622
623                         /* Compute quadratic angular scale interpolation weight */
624                         scale_weight = bpoly->point_edgemid_angles[bpoly->dominant_edge] / bpoly->edgemid_angle;
625                         scale_weight /= scale_weight + (bpoly->point_edgemid_angles[!bpoly->dominant_edge] / bpoly->edgemid_angle);
626
627                         sqr = scale_weight * scale_weight;
628                         inv_sqr = 1.0f - scale_weight;
629                         inv_sqr *= inv_sqr;
630                         scale_weight = sqr / (sqr + inv_sqr);
631
632                         /* Compute interpolated scale (no longer need the individual scales,
633                          * so simply storing the result over the scale in index zero) */
634                         bpoly->scales[0] = bpoly->scales[bpoly->dominant_edge] * (1.0f - scale_weight) +
635                                            bpoly->scales[!bpoly->dominant_edge] * scale_weight;
636
637                         /* Scale the point distance weights, and introduce falloff */
638                         bpoly->weight_dist_proj /= bpoly->scales[0];
639                         bpoly->weight_dist_proj = powf(bpoly->weight_dist_proj, data->falloff);
640
641                         bpoly->weight_dist /= avg_point_dist;
642                         bpoly->weight_dist = powf(bpoly->weight_dist, data->falloff);
643
644                         /* Re-check for infinite weights, now that all scalings and interpolations are computed */
645                         if (bpoly->weight_dist < FLT_EPSILON) {
646                                 inf_weight_flags |= MOD_SDEF_INFINITE_WEIGHT_DIST_PROJ;
647                                 inf_weight_flags |= MOD_SDEF_INFINITE_WEIGHT_DIST;
648                         }
649                         else if (bpoly->weight_dist_proj < FLT_EPSILON) {
650                                 inf_weight_flags |= MOD_SDEF_INFINITE_WEIGHT_DIST_PROJ;
651                         }
652                         else if (bpoly->weight_angular < FLT_EPSILON) {
653                                 inf_weight_flags |= MOD_SDEF_INFINITE_WEIGHT_ANGULAR;
654                         }
655                 }
656         }
657         else if (!(inf_weight_flags & MOD_SDEF_INFINITE_WEIGHT_DIST)) {
658                 bpoly = bwdata->bind_polys;
659
660                 for (int i = 0; i < bwdata->numpoly; bpoly++, i++) {
661                         /* Scale the point distance weight by average point distance, and introduce falloff */
662                         bpoly->weight_dist /= avg_point_dist;
663                         bpoly->weight_dist = powf(bpoly->weight_dist, data->falloff);
664
665                         /* Re-check for infinite weights, now that all scalings and interpolations are computed */
666                         if (bpoly->weight_dist < FLT_EPSILON) {
667                                 inf_weight_flags |= MOD_SDEF_INFINITE_WEIGHT_DIST;
668                         }
669                 }
670         }
671
672         /* Final loop, to compute actual weights */
673         bpoly = bwdata->bind_polys;
674
675         for (int i = 0; i < bwdata->numpoly; bpoly++, i++) {
676                 /* Weight computation from components */
677                 if (inf_weight_flags & MOD_SDEF_INFINITE_WEIGHT_DIST) {
678                         bpoly->weight = bpoly->weight_dist < FLT_EPSILON ? 1.0f : 0.0f;
679                 }
680                 else if (inf_weight_flags & MOD_SDEF_INFINITE_WEIGHT_DIST_PROJ) {
681                         bpoly->weight = bpoly->weight_dist_proj < FLT_EPSILON ?
682                                         1.0f / bpoly->weight_dist : 0.0f;
683                 }
684                 else if (inf_weight_flags & MOD_SDEF_INFINITE_WEIGHT_ANGULAR) {
685                         bpoly->weight = bpoly->weight_angular < FLT_EPSILON ?
686                                         1.0f / bpoly->weight_dist_proj / bpoly->weight_dist : 0.0f;
687                 }
688                 else {
689                         bpoly->weight = 1.0f / bpoly->weight_angular /
690                                                bpoly->weight_dist_proj /
691                                                bpoly->weight_dist;
692                 }
693
694                 tot_weight += bpoly->weight;
695         }
696
697         bpoly = bwdata->bind_polys;
698
699         for (int i = 0; i < bwdata->numpoly; bpoly++, i++) {
700                 bpoly->weight /= tot_weight;
701
702                 /* Evaluate if this poly is relevant to bind */
703                 /* Even though the weights should add up to 1.0,
704                  * the losses of weights smaller than epsilon here
705                  * should be negligible... */
706                 if (bpoly->weight >= FLT_EPSILON) {
707                         if (bpoly->inside) {
708                                 bwdata->numbinds += 1;
709                         }
710                         else {
711                                 if (bpoly->dominant_angle_weight < FLT_EPSILON || 1.0f - bpoly->dominant_angle_weight < FLT_EPSILON) {
712                                         bwdata->numbinds += 1;
713                                 }
714                                 else {
715                                         bwdata->numbinds += 2;
716                                 }
717                         }
718                 }
719         }
720
721         return bwdata;
722 }
723
724 BLI_INLINE float computeNormalDisplacement(const float point_co[3], const float point_co_proj[3], const float normal[3])
725 {
726         float disp_vec[3];
727         float normal_dist;
728
729         sub_v3_v3v3(disp_vec, point_co, point_co_proj);
730         normal_dist = len_v3(disp_vec);
731
732         if (dot_v3v3(disp_vec, normal) < 0) {
733                 normal_dist *= -1;
734         }
735
736         return normal_dist;
737 }
738
739 static void bindVert(
740         void *__restrict userdata,
741         const int index,
742         const ParallelRangeTLS *__restrict UNUSED(tls))
743 {
744         SDefBindCalcData * const data = (SDefBindCalcData *)userdata;
745         float point_co[3];
746         float point_co_proj[3];
747
748         SDefBindWeightData *bwdata;
749         SDefVert *sdvert = data->bind_verts + index;
750         SDefBindPoly *bpoly;
751         SDefBind *sdbind;
752
753         if (data->success != MOD_SDEF_BIND_RESULT_SUCCESS) {
754                 sdvert->binds = NULL;
755                 sdvert->numbinds = 0;
756                 return;
757         }
758
759         copy_v3_v3(point_co, data->vertexCos[index]);
760         bwdata = computeBindWeights(data, point_co);
761
762         if (bwdata == NULL) {
763                 sdvert->binds = NULL;
764                 sdvert->numbinds = 0;
765                 return;
766         }
767
768         sdvert->binds = MEM_calloc_arrayN(bwdata->numbinds, sizeof(*sdvert->binds), "SDefVertBindData");
769         if (sdvert->binds == NULL) {
770                 data->success = MOD_SDEF_BIND_RESULT_MEM_ERR;
771                 sdvert->numbinds = 0;
772                 return;
773         }
774
775         sdvert->numbinds = bwdata->numbinds;
776
777         sdbind = sdvert->binds;
778
779         bpoly = bwdata->bind_polys;
780
781         for (int i = 0; i < bwdata->numbinds; bpoly++) {
782                 if (bpoly->weight >= FLT_EPSILON) {
783                         if (bpoly->inside) {
784                                 const MLoop *loop = &data->mloop[bpoly->loopstart];
785
786                                 sdbind->influence = bpoly->weight;
787                                 sdbind->numverts = bpoly->numverts;
788
789                                 sdbind->mode = MOD_SDEF_MODE_NGON;
790                                 sdbind->vert_weights = MEM_malloc_arrayN(bpoly->numverts, sizeof(*sdbind->vert_weights), "SDefNgonVertWeights");
791                                 if (sdbind->vert_weights == NULL) {
792                                         data->success = MOD_SDEF_BIND_RESULT_MEM_ERR;
793                                         return;
794                                 }
795
796                                 sdbind->vert_inds = MEM_malloc_arrayN(bpoly->numverts, sizeof(*sdbind->vert_inds), "SDefNgonVertInds");
797                                 if (sdbind->vert_inds == NULL) {
798                                         data->success = MOD_SDEF_BIND_RESULT_MEM_ERR;
799                                         return;
800                                 }
801
802                                 interp_weights_poly_v2(sdbind->vert_weights, bpoly->coords_v2, bpoly->numverts, bpoly->point_v2);
803
804                                 /* Reproject vert based on weights and original poly verts, to reintroduce poly non-planarity */
805                                 zero_v3(point_co_proj);
806                                 for (int j = 0; j < bpoly->numverts; j++, loop++) {
807                                         madd_v3_v3fl(point_co_proj, bpoly->coords[j], sdbind->vert_weights[j]);
808                                         sdbind->vert_inds[j] = loop->v;
809                                 }
810
811                                 sdbind->normal_dist = computeNormalDisplacement(point_co, point_co_proj, bpoly->normal);
812
813                                 sdbind++;
814                                 i++;
815                         }
816                         else {
817                                 float tmp_vec[3];
818                                 float cent[3], norm[3];
819                                 float v1[3], v2[3], v3[3];
820
821                                 if (1.0f - bpoly->dominant_angle_weight >= FLT_EPSILON) {
822                                         sdbind->influence = bpoly->weight * (1.0f - bpoly->dominant_angle_weight);
823                                         sdbind->numverts = bpoly->numverts;
824
825                                         sdbind->mode = MOD_SDEF_MODE_CENTROID;
826                                         sdbind->vert_weights = MEM_malloc_arrayN(3, sizeof(*sdbind->vert_weights), "SDefCentVertWeights");
827                                         if (sdbind->vert_weights == NULL) {
828                                                 data->success = MOD_SDEF_BIND_RESULT_MEM_ERR;
829                                                 return;
830                                         }
831
832                                         sdbind->vert_inds = MEM_malloc_arrayN(bpoly->numverts, sizeof(*sdbind->vert_inds), "SDefCentVertInds");
833                                         if (sdbind->vert_inds == NULL) {
834                                                 data->success = MOD_SDEF_BIND_RESULT_MEM_ERR;
835                                                 return;
836                                         }
837
838                                         sortPolyVertsEdge(sdbind->vert_inds, &data->mloop[bpoly->loopstart],
839                                                           bpoly->edge_inds[bpoly->dominant_edge], bpoly->numverts);
840
841                                         copy_v3_v3(v1, data->targetCos[sdbind->vert_inds[0]]);
842                                         copy_v3_v3(v2, data->targetCos[sdbind->vert_inds[1]]);
843                                         copy_v3_v3(v3, bpoly->centroid);
844
845                                         mid_v3_v3v3v3(cent, v1, v2, v3);
846                                         normal_tri_v3(norm, v1, v2, v3);
847
848                                         add_v3_v3v3(tmp_vec, point_co, bpoly->normal);
849
850                                         /* We are sure the line is not parallel to the plane.
851                                          * Checking return value just to avoid warning... */
852                                         if (!isect_line_plane_v3(point_co_proj, point_co, tmp_vec, cent, norm)) {
853                                                 BLI_assert(false);
854                                         }
855
856                                         interp_weights_tri_v3(sdbind->vert_weights, v1, v2, v3, point_co_proj);
857
858                                         sdbind->normal_dist = computeNormalDisplacement(point_co, point_co_proj, bpoly->normal);
859
860                                         sdbind++;
861                                         i++;
862                                 }
863
864                                 if (bpoly->dominant_angle_weight >= FLT_EPSILON) {
865                                         sdbind->influence = bpoly->weight * bpoly->dominant_angle_weight;
866                                         sdbind->numverts = bpoly->numverts;
867
868                                         sdbind->mode = MOD_SDEF_MODE_LOOPTRI;
869                                         sdbind->vert_weights = MEM_malloc_arrayN(3, sizeof(*sdbind->vert_weights), "SDefTriVertWeights");
870                                         if (sdbind->vert_weights == NULL) {
871                                                 data->success = MOD_SDEF_BIND_RESULT_MEM_ERR;
872                                                 return;
873                                         }
874
875                                         sdbind->vert_inds = MEM_malloc_arrayN(bpoly->numverts, sizeof(*sdbind->vert_inds), "SDefTriVertInds");
876                                         if (sdbind->vert_inds == NULL) {
877                                                 data->success = MOD_SDEF_BIND_RESULT_MEM_ERR;
878                                                 return;
879                                         }
880
881                                         sortPolyVertsTri(sdbind->vert_inds, &data->mloop[bpoly->loopstart], bpoly->edge_vert_inds[0], bpoly->numverts);
882
883                                         copy_v3_v3(v1, data->targetCos[sdbind->vert_inds[0]]);
884                                         copy_v3_v3(v2, data->targetCos[sdbind->vert_inds[1]]);
885                                         copy_v3_v3(v3, data->targetCos[sdbind->vert_inds[2]]);
886
887                                         mid_v3_v3v3v3(cent, v1, v2, v3);
888                                         normal_tri_v3(norm, v1, v2, v3);
889
890                                         add_v3_v3v3(tmp_vec, point_co, bpoly->normal);
891
892                                         /* We are sure the line is not parallel to the plane.
893                                          * Checking return value just to avoid warning... */
894                                         if (!isect_line_plane_v3(point_co_proj, point_co, tmp_vec, cent, norm)) {
895                                                 BLI_assert(false);
896                                         }
897
898                                         interp_weights_tri_v3(sdbind->vert_weights, v1, v2, v3, point_co_proj);
899
900                                         sdbind->normal_dist = computeNormalDisplacement(point_co, point_co_proj, bpoly->normal);
901
902                                         sdbind++;
903                                         i++;
904                                 }
905                         }
906                 }
907         }
908
909         freeBindData(bwdata);
910 }
911
912 static bool surfacedeformBind(
913         SurfaceDeformModifierData *smd, float (*vertexCos)[3],
914         unsigned int numverts, unsigned int tnumpoly, unsigned int tnumverts, Mesh *target)
915 {
916         BVHTreeFromMesh treeData = {NULL};
917         const MVert *mvert = target->mvert;
918         const MPoly *mpoly = target->mpoly;
919         const MEdge *medge = target->medge;
920         const MLoop *mloop = target->mloop;
921         unsigned int tnumedges = target->totedge;
922         int adj_result;
923         SDefAdjacencyArray *vert_edges;
924         SDefAdjacency *adj_array;
925         SDefEdgePolys *edge_polys;
926
927         vert_edges = MEM_calloc_arrayN(tnumverts, sizeof(*vert_edges), "SDefVertEdgeMap");
928         if (vert_edges == NULL) {
929                 modifier_setError((ModifierData *)smd, "Out of memory");
930                 return false;
931         }
932
933         adj_array = MEM_malloc_arrayN(tnumedges, 2 * sizeof(*adj_array), "SDefVertEdge");
934         if (adj_array == NULL) {
935                 modifier_setError((ModifierData *)smd, "Out of memory");
936                 MEM_freeN(vert_edges);
937                 return false;
938         }
939
940         edge_polys = MEM_calloc_arrayN(tnumedges, sizeof(*edge_polys), "SDefEdgeFaceMap");
941         if (edge_polys == NULL) {
942                 modifier_setError((ModifierData *)smd, "Out of memory");
943                 MEM_freeN(vert_edges);
944                 MEM_freeN(adj_array);
945                 return false;
946         }
947
948         smd->verts = MEM_malloc_arrayN(numverts, sizeof(*smd->verts), "SDefBindVerts");
949         if (smd->verts == NULL) {
950                 modifier_setError((ModifierData *)smd, "Out of memory");
951                 freeAdjacencyMap(vert_edges, adj_array, edge_polys);
952                 return false;
953         }
954
955         BKE_bvhtree_from_mesh_get(&treeData, target, BVHTREE_FROM_LOOPTRI, 2);
956         if (treeData.tree == NULL) {
957                 modifier_setError((ModifierData *)smd, "Out of memory");
958                 freeAdjacencyMap(vert_edges, adj_array, edge_polys);
959                 MEM_freeN(smd->verts);
960                 smd->verts = NULL;
961                 return false;
962         }
963
964         adj_result = buildAdjacencyMap(mpoly, medge, mloop, tnumpoly, tnumedges, vert_edges, adj_array, edge_polys);
965
966         if (adj_result == MOD_SDEF_BIND_RESULT_NONMANY_ERR) {
967                 modifier_setError((ModifierData *)smd, "Target has edges with more than two polygons");
968                 freeAdjacencyMap(vert_edges, adj_array, edge_polys);
969                 free_bvhtree_from_mesh(&treeData);
970                 MEM_freeN(smd->verts);
971                 smd->verts = NULL;
972                 return false;
973         }
974
975         smd->numverts = numverts;
976         smd->numpoly = tnumpoly;
977
978         SDefBindCalcData data = {.treeData = &treeData,
979                                      .vert_edges = vert_edges,
980                                      .edge_polys = edge_polys,
981                                      .mpoly = mpoly,
982                                      .medge = medge,
983                                      .mloop = mloop,
984                                      .looptri = BKE_mesh_runtime_looptri_ensure(target),
985                                      .targetCos = MEM_malloc_arrayN(tnumverts, sizeof(float[3]), "SDefTargetBindVertArray"),
986                                      .bind_verts = smd->verts,
987                                      .vertexCos = vertexCos,
988                                      .falloff = smd->falloff,
989                                      .success = MOD_SDEF_BIND_RESULT_SUCCESS};
990
991         if (data.targetCos == NULL) {
992                 modifier_setError((ModifierData *)smd, "Out of memory");
993                 freeData((ModifierData *)smd);
994                 return false;
995         }
996
997         invert_m4_m4(data.imat, smd->mat);
998
999         for (int i = 0; i < tnumverts; i++) {
1000                 mul_v3_m4v3(data.targetCos[i], smd->mat, mvert[i].co);
1001         }
1002
1003         ParallelRangeSettings settings;
1004         BLI_parallel_range_settings_defaults(&settings);
1005         settings.use_threading = (numverts > 10000);
1006         BLI_task_parallel_range(0, numverts,
1007                                 &data,
1008                                 bindVert,
1009                                 &settings);
1010
1011         MEM_freeN(data.targetCos);
1012
1013         if (data.success == MOD_SDEF_BIND_RESULT_MEM_ERR) {
1014                 modifier_setError((ModifierData *)smd, "Out of memory");
1015                 freeData((ModifierData *)smd);
1016         }
1017         else if (data.success == MOD_SDEF_BIND_RESULT_NONMANY_ERR) {
1018                 modifier_setError((ModifierData *)smd, "Target has edges with more than two polygons");
1019                 freeData((ModifierData *)smd);
1020         }
1021         else if (data.success == MOD_SDEF_BIND_RESULT_CONCAVE_ERR) {
1022                 modifier_setError((ModifierData *)smd, "Target contains concave polygons");
1023                 freeData((ModifierData *)smd);
1024         }
1025         else if (data.success == MOD_SDEF_BIND_RESULT_OVERLAP_ERR) {
1026                 modifier_setError((ModifierData *)smd, "Target contains overlapping verts");
1027                 freeData((ModifierData *)smd);
1028         }
1029         else if (data.success == MOD_SDEF_BIND_RESULT_GENERIC_ERR) {
1030                 /* I know this message is vague, but I could not think of a way
1031                  * to explain this whith a reasonably sized message.
1032                  * Though it shouldn't really matter all that much,
1033                  * because this is very unlikely to occur */
1034                 modifier_setError((ModifierData *)smd, "Target contains invalid polygons");
1035                 freeData((ModifierData *)smd);
1036         }
1037
1038         freeAdjacencyMap(vert_edges, adj_array, edge_polys);
1039         free_bvhtree_from_mesh(&treeData);
1040
1041         return data.success == 1;
1042 }
1043
1044 static Mesh *surfacedeform_get_mesh(SurfaceDeformModifierData *smd, bool *r_needsfree)
1045 {
1046         Mesh *mesh;
1047
1048         /* Handle target mesh both in and out of edit mode */
1049         if (smd->target->mode & OB_MODE_EDIT) {
1050                 BMEditMesh *em = BKE_editmesh_from_object(smd->target);
1051                 mesh = BKE_bmesh_to_mesh_nomain(em->bm, &(struct BMeshToMeshParams){0});
1052                 *r_needsfree = true;
1053         }
1054         else {
1055                 mesh = BKE_modifier_get_evaluated_mesh_from_object(
1056                            smd->target, smd->modifier.mode & eModifierMode_Render ? MOD_APPLY_RENDER : 0);
1057                 *r_needsfree = false;
1058         }
1059
1060         if (!mesh) {
1061                 mesh = get_mesh(smd->target, NULL, NULL, NULL, false, false);
1062                 *r_needsfree = true;
1063         }
1064
1065         return mesh;
1066 }
1067
1068 static void deformVert(
1069         void *__restrict userdata,
1070         const int index,
1071         const ParallelRangeTLS *__restrict UNUSED(tls))
1072 {
1073         const SDefDeformData * const data = (SDefDeformData *)userdata;
1074         const SDefBind *sdbind = data->bind_verts[index].binds;
1075         float * const vertexCos = data->vertexCos[index];
1076         float norm[3], temp[3];
1077
1078         zero_v3(vertexCos);
1079
1080         for (int j = 0; j < data->bind_verts[index].numbinds; j++, sdbind++) {
1081                 /* Mode-generic operations (allocate poly coordinates) */
1082                 float (*coords)[3] = MEM_malloc_arrayN(sdbind->numverts, sizeof(*coords), "SDefDoPolyCoords");
1083
1084                 for (int k = 0; k < sdbind->numverts; k++) {
1085                         copy_v3_v3(coords[k], data->targetCos[sdbind->vert_inds[k]]);
1086                 }
1087
1088                 normal_poly_v3(norm, coords, sdbind->numverts);
1089                 zero_v3(temp);
1090
1091                 /* ---------- looptri mode ---------- */
1092                 if (sdbind->mode == MOD_SDEF_MODE_LOOPTRI) {
1093                         madd_v3_v3fl(temp, data->targetCos[sdbind->vert_inds[0]], sdbind->vert_weights[0]);
1094                         madd_v3_v3fl(temp, data->targetCos[sdbind->vert_inds[1]], sdbind->vert_weights[1]);
1095                         madd_v3_v3fl(temp, data->targetCos[sdbind->vert_inds[2]], sdbind->vert_weights[2]);
1096                 }
1097                 else {
1098                         /* ---------- ngon mode ---------- */
1099                         if (sdbind->mode == MOD_SDEF_MODE_NGON) {
1100                                 for (int k = 0; k < sdbind->numverts; k++) {
1101                                         madd_v3_v3fl(temp, coords[k], sdbind->vert_weights[k]);
1102                                 }
1103                         }
1104
1105                         /* ---------- centroid mode ---------- */
1106                         else if (sdbind->mode == MOD_SDEF_MODE_CENTROID) {
1107                                 float cent[3];
1108                                 mid_v3_v3_array(cent, coords, sdbind->numverts);
1109
1110                                 madd_v3_v3fl(temp, data->targetCos[sdbind->vert_inds[0]], sdbind->vert_weights[0]);
1111                                 madd_v3_v3fl(temp, data->targetCos[sdbind->vert_inds[1]], sdbind->vert_weights[1]);
1112                                 madd_v3_v3fl(temp, cent, sdbind->vert_weights[2]);
1113                         }
1114                 }
1115
1116                 MEM_freeN(coords);
1117
1118                 /* Apply normal offset (generic for all modes) */
1119                 madd_v3_v3fl(temp, norm, sdbind->normal_dist);
1120
1121                 madd_v3_v3fl(vertexCos, temp, sdbind->influence);
1122         }
1123 }
1124
1125 static void surfacedeformModifier_do(
1126         ModifierData *md,
1127         float (*vertexCos)[3], unsigned int numverts, Object *ob)
1128 {
1129         SurfaceDeformModifierData *smd = (SurfaceDeformModifierData *)md;
1130         bool free_target;
1131         Mesh *target;
1132         unsigned int tnumverts, tnumpoly;
1133
1134         /* Exit function if bind flag is not set (free bind data if any) */
1135         if (!(smd->flags & MOD_SDEF_BIND)) {
1136                 freeData(md);
1137                 return;
1138         }
1139
1140         target = surfacedeform_get_mesh(smd, &free_target);
1141         if (!target) {
1142                 modifier_setError(md, "No valid target mesh");
1143                 return;
1144         }
1145
1146         tnumverts = target->totvert;
1147         tnumpoly = target->totpoly;
1148
1149         /* If not bound, execute bind */
1150         if (!(smd->verts)) {
1151                 float tmp_mat[4][4];
1152
1153                 invert_m4_m4(tmp_mat, ob->obmat);
1154                 mul_m4_m4m4(smd->mat, tmp_mat, smd->target->obmat);
1155
1156                 if (!surfacedeformBind(smd, vertexCos, numverts, tnumpoly, tnumverts, target)) {
1157                         smd->flags &= ~MOD_SDEF_BIND;
1158                         return;
1159                 }
1160         }
1161
1162         /* Poly count checks */
1163         if (smd->numverts != numverts) {
1164                 modifier_setError(md, "Verts changed from %u to %u", smd->numverts, numverts);
1165                 if (free_target) BKE_id_free(NULL, target);
1166                 return;
1167         }
1168         else if (smd->numpoly != tnumpoly) {
1169                 modifier_setError(md, "Target polygons changed from %u to %u", smd->numpoly, tnumpoly);
1170                 if (free_target) BKE_id_free(NULL, target);
1171                 return;
1172         }
1173
1174         /* Actual vertex location update starts here */
1175         SDefDeformData data = {
1176                 .bind_verts = smd->verts,
1177                 .targetCos = MEM_malloc_arrayN(tnumverts, sizeof(float[3]), "SDefTargetVertArray"),
1178                 .vertexCos = vertexCos,
1179         };
1180
1181         if (data.targetCos != NULL) {
1182                 const MVert * const mvert = target->mvert;
1183
1184                 for (int i = 0; i < tnumverts; i++) {
1185                         mul_v3_m4v3(data.targetCos[i], smd->mat, mvert[i].co);
1186                 }
1187
1188                 ParallelRangeSettings settings;
1189                 BLI_parallel_range_settings_defaults(&settings);
1190                 settings.use_threading = (numverts > 10000);
1191                 BLI_task_parallel_range(0, numverts,
1192                                         &data,
1193                                         deformVert,
1194                                         &settings);
1195
1196                 MEM_freeN(data.targetCos);
1197         }
1198
1199         if (free_target) BKE_id_free(NULL, target);
1200 }
1201
1202 static void deformVerts(
1203         ModifierData *md, const ModifierEvalContext *ctx,
1204         Mesh *UNUSED(mesh),
1205         float (*vertexCos)[3], int numVerts)
1206 {
1207         surfacedeformModifier_do(md, vertexCos, numVerts, ctx->object);
1208 }
1209
1210 static void deformVertsEM(
1211         ModifierData *md, const ModifierEvalContext *ctx,
1212         struct BMEditMesh *UNUSED(editData),
1213         Mesh *UNUSED(mesh),
1214         float (*vertexCos)[3], int numVerts)
1215 {
1216         surfacedeformModifier_do(md, vertexCos, numVerts, ctx->object);
1217 }
1218
1219 static bool isDisabled(ModifierData *md, int UNUSED(useRenderParams))
1220 {
1221         SurfaceDeformModifierData *smd = (SurfaceDeformModifierData *)md;
1222
1223         return !smd->target && !(smd->verts && !(smd->flags & MOD_SDEF_BIND));
1224 }
1225
1226 ModifierTypeInfo modifierType_SurfaceDeform = {
1227         /* name */              "Surface Deform",
1228         /* structName */        "SurfaceDeformModifierData",
1229         /* structSize */        sizeof(SurfaceDeformModifierData),
1230         /* type */              eModifierTypeType_OnlyDeform,
1231         /* flags */             eModifierTypeFlag_AcceptsMesh |
1232                                 eModifierTypeFlag_SupportsEditmode,
1233
1234         /* copyData */          copyData,
1235
1236         /* deformVerts_DM */    NULL,
1237         /* deformMatrices_DM */ NULL,
1238         /* deformVertsEM_DM */  NULL,
1239         /* deformMatricesEM_DM*/NULL,
1240         /* applyModifier_DM */  NULL,
1241         /* applyModifierEM_DM */NULL,
1242
1243         /* deformVerts */       deformVerts,
1244         /* deformMatrices */    NULL,
1245         /* deformVertsEM */     deformVertsEM,
1246         /* deformMatricesEM */  NULL,
1247         /* applyModifier */     NULL,
1248         /* applyModifierEM */   NULL,
1249
1250         /* initData */          initData,
1251         /* requiredDataMask */  NULL,
1252         /* freeData */          freeData,
1253         /* isDisabled */        isDisabled,
1254         /* updateDepsgraph */   updateDepsgraph,
1255         /* dependsOnTime */     NULL,
1256         /* dependsOnNormals */  NULL,
1257         /* foreachObjectLink */ foreachObjectLink,
1258         /* foreachIDLink */     NULL,
1259         /* foreachTexLink */    NULL,
1260 };