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