diff options
Diffstat (limited to 'source/blender/modifiers/intern/MOD_surfacedeform.c')
-rw-r--r-- | source/blender/modifiers/intern/MOD_surfacedeform.c | 1189 |
1 files changed, 1189 insertions, 0 deletions
diff --git a/source/blender/modifiers/intern/MOD_surfacedeform.c b/source/blender/modifiers/intern/MOD_surfacedeform.c new file mode 100644 index 00000000000..5e852e84516 --- /dev/null +++ b/source/blender/modifiers/intern/MOD_surfacedeform.c @@ -0,0 +1,1189 @@ +#include "DNA_object_types.h" +#include "DNA_scene_types.h" + +#include "BLI_alloca.h" +#include "BLI_math.h" +#include "BLI_math_geom.h" +#include "BLI_task.h" + +#include "BKE_cdderivedmesh.h" +#include "BKE_editmesh.h" +#include "BKE_library_query.h" +#include "BKE_modifier.h" + +#include "depsgraph_private.h" + +#include "MEM_guardedalloc.h" + +#include "MOD_util.h" + +typedef struct SDefAdjacency { + struct SDefAdjacency *next; + unsigned int index; +} SDefAdjacency; + +typedef struct SDefAdjacencyArray { + SDefAdjacency *first; + unsigned int num; /* Careful, this is twice the number of polygons (avoids an extra loop) */ +} SDefAdjacencyArray; + +typedef struct SDefEdgePolys { + unsigned int polys[2], num; +} SDefEdgePolys; + +typedef struct SDefBindCalcData { + BVHTreeFromMesh * const treeData; + const SDefAdjacencyArray * const vert_edges; + const SDefEdgePolys * const edge_polys; + SDefVert * const bind_verts; + const MLoopTri * const looptri; + const MPoly * const mpoly; + const MEdge * const medge; + const MLoop * const mloop; + const MVert * const mvert; + float (* const vertexCos)[3]; + const float falloff; + int success; +} SDefBindCalcData; + +typedef struct SDefBindPoly { + float (*coords)[3]; + float (*coords_v2)[2]; + float point_v2[2]; + float weight_angular; + float weight_dist_proj; + float weight_dist; + float weight; + float scales[2]; + float centroid[3]; + float centroid_v2[2]; + float normal[3]; + float cent_edgemid_vecs_v2[2][2]; + float edgemid_angle; + float point_edgemid_angles[2]; + float corner_edgemid_angles[2]; + float dominant_angle_weight; + unsigned int index; + unsigned int numverts; + unsigned int loopstart; + unsigned int edge_inds[2]; + unsigned int edge_vert_inds[2]; + unsigned int corner_ind; + unsigned int dominant_edge; + bool inside; +} SDefBindPoly; + +typedef struct SDefBindWeightData { + SDefBindPoly *bind_polys; + unsigned int numpoly; + unsigned int numbinds; +} SDefBindWeightData; + +typedef struct SDefDeformData { + const SDefVert * const bind_verts; + const MVert * const mvert; + float (* const vertexCos)[3]; +} SDefDeformData; + +/* Bind result values */ +enum { + MOD_SDEF_BIND_RESULT_SUCCESS = 1, + MOD_SDEF_BIND_RESULT_GENERIC_ERR = 0, + MOD_SDEF_BIND_RESULT_MEM_ERR = -1, + MOD_SDEF_BIND_RESULT_NONMANY_ERR = -2, + MOD_SDEF_BIND_RESULT_CONCAVE_ERR = -3, + MOD_SDEF_BIND_RESULT_OVERLAP_ERR = -4, +}; + +/* Infinite weight flags */ +enum { + MOD_SDEF_INFINITE_WEIGHT_ANGULAR = (1 << 0), + MOD_SDEF_INFINITE_WEIGHT_DIST_PROJ = (1 << 1), + MOD_SDEF_INFINITE_WEIGHT_DIST = (1 << 2), +}; + +static void initData(ModifierData *md) +{ + SurfaceDeformModifierData *smd = (SurfaceDeformModifierData *)md; + smd->target = NULL; + smd->verts = NULL; + smd->flags = 0; + smd->falloff = 4.0f; +} + +static void freeData(ModifierData *md) +{ + SurfaceDeformModifierData *smd = (SurfaceDeformModifierData *)md; + + if (smd->verts) { + for (int i = 0; i < smd->numverts; i++) { + if (smd->verts[i].binds) { + for (int j = 0; j < smd->verts[i].numbinds; j++) { + MEM_SAFE_FREE(smd->verts[i].binds[j].vert_inds); + MEM_SAFE_FREE(smd->verts[i].binds[j].vert_weights); + } + + MEM_freeN(smd->verts[i].binds); + } + } + + MEM_freeN(smd->verts); + smd->verts = NULL; + } +} + +static void copyData(ModifierData *md, ModifierData *target) +{ + SurfaceDeformModifierData *smd = (SurfaceDeformModifierData *)md; + SurfaceDeformModifierData *tsmd = (SurfaceDeformModifierData *)target; + + *tsmd = *smd; + + if (smd->verts) { + tsmd->verts = MEM_dupallocN(smd->verts); + + for (int i = 0; i < smd->numverts; i++) { + if (smd->verts[i].binds) { + tsmd->verts[i].binds = MEM_dupallocN(smd->verts[i].binds); + + for (int j = 0; j < smd->verts[i].numbinds; j++) { + if (smd->verts[i].binds[j].vert_inds) { + tsmd->verts[i].binds[j].vert_inds = MEM_dupallocN(smd->verts[i].binds[j].vert_inds); + } + + if (smd->verts[i].binds[j].vert_weights) { + tsmd->verts[i].binds[j].vert_weights = MEM_dupallocN(smd->verts[i].binds[j].vert_weights); + } + } + } + } + } +} + +static void foreachObjectLink(ModifierData *md, Object *ob, ObjectWalkFunc walk, void *userData) +{ + SurfaceDeformModifierData *smd = (SurfaceDeformModifierData *)md; + + walk(userData, ob, &smd->target, IDWALK_NOP); +} + +static void updateDepgraph(ModifierData *md, DagForest *forest, + struct Main *UNUSED(bmain), + struct Scene *UNUSED(scene), + Object *UNUSED(ob), + DagNode *obNode) +{ + SurfaceDeformModifierData *smd = (SurfaceDeformModifierData *)md; + + if (smd->target) { + DagNode *curNode = dag_get_node(forest, smd->target); + + dag_add_relation(forest, curNode, obNode, DAG_RL_DATA_DATA, "Surface Deform Modifier"); + } +} + +static void updateDepsgraph(ModifierData *md, + struct Main *UNUSED(bmain), + struct Scene *UNUSED(scene), + Object *UNUSED(ob), + struct DepsNodeHandle *node) +{ + SurfaceDeformModifierData *smd = (SurfaceDeformModifierData *)md; + if (smd->target != NULL) { + DEG_add_object_relation(node, smd->target, DEG_OB_COMP_GEOMETRY, "Surface Deform Modifier"); + } +} + +static void freeAdjacencyMap(SDefAdjacencyArray * const vert_edges, SDefAdjacency * const adj_ref, SDefEdgePolys * const edge_polys) +{ + MEM_freeN(edge_polys); + + MEM_freeN(adj_ref); + + MEM_freeN(vert_edges); +} + +static int buildAdjacencyMap(const MPoly *poly, const MEdge *edge, const MLoop * const mloop, const unsigned int numpoly, const unsigned int numedges, + SDefAdjacencyArray * const vert_edges, SDefAdjacency *adj, SDefEdgePolys * const edge_polys) +{ + const MLoop *loop; + + /* Fing polygons adjacent to edges */ + for (int i = 0; i < numpoly; i++, poly++) { + loop = &mloop[poly->loopstart]; + + for (int j = 0; j < poly->totloop; j++, loop++) { + if (edge_polys[loop->e].num == 0) { + edge_polys[loop->e].polys[0] = i; + edge_polys[loop->e].polys[1] = -1; + edge_polys[loop->e].num++; + } + else if (edge_polys[loop->e].num == 1) { + edge_polys[loop->e].polys[1] = i; + edge_polys[loop->e].num++; + } + else { + return MOD_SDEF_BIND_RESULT_NONMANY_ERR; + } + } + } + + /* Find edges adjacent to vertices */ + for (int i = 0; i < numedges; i++, edge++) { + adj->next = vert_edges[edge->v1].first; + adj->index = i; + vert_edges[edge->v1].first = adj; + vert_edges[edge->v1].num += edge_polys[i].num; + adj++; + + adj->next = vert_edges[edge->v2].first; + adj->index = i; + vert_edges[edge->v2].first = adj; + vert_edges[edge->v2].num += edge_polys[i].num; + adj++; + } + + return MOD_SDEF_BIND_RESULT_SUCCESS; +} + +BLI_INLINE void sortPolyVertsEdge(unsigned int *indices, const MLoop * const mloop, const unsigned int edge, const unsigned int num) +{ + bool found = false; + + for (int i = 0; i < num; i++) { + if (mloop[i].e == edge) { + found = true; + } + if (found) { + *indices = mloop[i].v; + indices++; + } + } + + /* Fill in remaining vertex indices that occur before the edge */ + for (int i = 0; mloop[i].e != edge; i++) { + *indices = mloop[i].v; + indices++; + } +} + +BLI_INLINE void sortPolyVertsTri(unsigned int *indices, const MLoop * const mloop, const unsigned int loopstart, const unsigned int num) +{ + for (int i = loopstart; i < num; i++) { + *indices = mloop[i].v; + indices++; + } + + for (int i = 0; i < loopstart; i++) { + *indices = mloop[i].v; + indices++; + } +} + +BLI_INLINE unsigned int nearestVert(SDefBindCalcData * const data, const float point_co[3]) +{ + const MVert * const mvert = data->mvert; + BVHTreeNearest nearest = {.dist_sq = FLT_MAX, .index = -1}; + const MPoly *poly; + const MEdge *edge; + const MLoop *loop; + float max_dist = FLT_MAX; + float dist; + unsigned int index = 0; + + BLI_bvhtree_find_nearest(data->treeData->tree, point_co, &nearest, data->treeData->nearest_callback, data->treeData); + + poly = &data->mpoly[data->looptri[nearest.index].poly]; + loop = &data->mloop[poly->loopstart]; + + for (int i = 0; i < poly->totloop; i++, loop++) { + edge = &data->medge[loop->e]; + dist = dist_squared_to_line_segment_v3(point_co, mvert[edge->v1].co, mvert[edge->v2].co); + + if (dist < max_dist) { + max_dist = dist; + index = loop->e; + } + } + + edge = &data->medge[index]; + if (len_squared_v3v3(point_co, mvert[edge->v1].co) < len_squared_v3v3(point_co, mvert[edge->v2].co)) { + return edge->v1; + } + else { + return edge->v2; + } +} + +BLI_INLINE int isPolyValid(const float coords[][2], const unsigned int nr) +{ + float prev_co[2]; + float curr_vec[2], prev_vec[2]; + + if (!is_poly_convex_v2(coords, nr)) { + return MOD_SDEF_BIND_RESULT_CONCAVE_ERR; + } + + copy_v2_v2(prev_co, coords[nr - 1]); + sub_v2_v2v2(prev_vec, prev_co, coords[nr - 2]); + + for (int i = 0; i < nr; i++) { + sub_v2_v2v2(curr_vec, coords[i], prev_co); + + if (len_squared_v2(curr_vec) < FLT_EPSILON) { + return MOD_SDEF_BIND_RESULT_OVERLAP_ERR; + } + + if (1.0f - dot_v2v2(prev_vec, curr_vec) < FLT_EPSILON) { + return MOD_SDEF_BIND_RESULT_CONCAVE_ERR; + } + + copy_v2_v2(prev_co, coords[i]); + copy_v2_v2(prev_vec, curr_vec); + } + + return MOD_SDEF_BIND_RESULT_SUCCESS; +} + +static void freeBindData(SDefBindWeightData * const bwdata) +{ + SDefBindPoly *bpoly = bwdata->bind_polys; + + if (bwdata->bind_polys) { + for (int i = 0; i < bwdata->numpoly; bpoly++, i++) { + MEM_SAFE_FREE(bpoly->coords); + MEM_SAFE_FREE(bpoly->coords_v2); + } + + MEM_freeN(bwdata->bind_polys); + } + + MEM_freeN(bwdata); +} + +BLI_INLINE float computeAngularWeight(const float point_angle, const float edgemid_angle) +{ + float weight; + + weight = point_angle; + weight /= edgemid_angle; + weight *= M_PI_2; + + return sinf(weight); +} + +BLI_INLINE SDefBindWeightData *computeBindWeights(SDefBindCalcData * const data, const float point_co[3]) +{ + const unsigned int nearest = nearestVert(data, point_co); + const SDefAdjacency * const vert_edges = data->vert_edges[nearest].first; + const SDefEdgePolys * const edge_polys = data->edge_polys; + + const SDefAdjacency *vedge; + const MPoly *poly; + const MLoop *loop; + + SDefBindWeightData *bwdata; + SDefBindPoly *bpoly; + + float world[3] = {0.0f, 0.0f, 1.0f}; + float avg_point_dist = 0.0f; + float tot_weight = 0.0f; + int inf_weight_flags = 0; + + bwdata = MEM_callocN(sizeof(*bwdata), "SDefBindWeightData"); + if (bwdata == NULL) { + data->success = MOD_SDEF_BIND_RESULT_MEM_ERR; + return NULL; + } + + bwdata->numpoly = data->vert_edges[nearest].num / 2; + + bpoly = MEM_callocN(sizeof(*bpoly) * bwdata->numpoly, "SDefBindPoly"); + if (bpoly == NULL) { + freeBindData(bwdata); + data->success = MOD_SDEF_BIND_RESULT_MEM_ERR; + return NULL; + } + + bwdata->bind_polys = bpoly; + + /* Loop over all adjacent edges, and build the SDefBindPoly data for each poly adjacent to those */ + for (vedge = vert_edges; vedge; vedge = vedge->next) { + unsigned int edge_ind = vedge->index; + + for (int i = 0; i < edge_polys[edge_ind].num; i++) { + { + bpoly = bwdata->bind_polys; + + for (int j = 0; j < bwdata->numpoly; bpoly++, j++) { + /* If coords isn't allocated, we have reached the first uninitialized bpoly */ + if ((bpoly->index == edge_polys[edge_ind].polys[i]) || (!bpoly->coords)) { + break; + } + } + } + + /* Check if poly was already created by another edge or still has to be initialized */ + if (!bpoly->coords) { + float angle; + float axis[3]; + float tmp_vec_v2[2]; + int is_poly_valid; + + bpoly->index = edge_polys[edge_ind].polys[i]; + bpoly->coords = NULL; + bpoly->coords_v2 = NULL; + + /* Copy poly data */ + poly = &data->mpoly[bpoly->index]; + loop = &data->mloop[poly->loopstart]; + + bpoly->numverts = poly->totloop; + bpoly->loopstart = poly->loopstart; + + bpoly->coords = MEM_mallocN(sizeof(*bpoly->coords) * poly->totloop, "SDefBindPolyCoords"); + if (bpoly->coords == NULL) { + freeBindData(bwdata); + data->success = MOD_SDEF_BIND_RESULT_MEM_ERR; + return NULL; + } + + bpoly->coords_v2 = MEM_mallocN(sizeof(*bpoly->coords_v2) * poly->totloop, "SDefBindPolyCoords_v2"); + if (bpoly->coords_v2 == NULL) { + freeBindData(bwdata); + data->success = MOD_SDEF_BIND_RESULT_MEM_ERR; + return NULL; + } + + for (int j = 0; j < poly->totloop; j++, loop++) { + copy_v3_v3(bpoly->coords[j], data->mvert[loop->v].co); + + /* Find corner and edge indices within poly loop array */ + if (loop->v == nearest) { + bpoly->corner_ind = j; + bpoly->edge_vert_inds[0] = (j == 0) ? (poly->totloop - 1) : (j - 1); + bpoly->edge_vert_inds[1] = (j == poly->totloop - 1) ? (0) : (j + 1); + + bpoly->edge_inds[0] = data->mloop[poly->loopstart + bpoly->edge_vert_inds[0]].e; + bpoly->edge_inds[1] = loop->e; + } + } + + /* Compute poly's parametric data */ + mid_v3_v3_array(bpoly->centroid, bpoly->coords, poly->totloop); + normal_poly_v3(bpoly->normal, bpoly->coords, poly->totloop); + + /* Compute poly skew angle and axis */ + angle = angle_normalized_v3v3(bpoly->normal, world); + + cross_v3_v3v3(axis, bpoly->normal, world); + normalize_v3(axis); + + /* Map coords onto 2d normal plane */ + map_to_plane_axis_angle_v2_v3v3fl(bpoly->point_v2, point_co, axis, angle); + + zero_v2(bpoly->centroid_v2); + for (int j = 0; j < poly->totloop; j++) { + map_to_plane_axis_angle_v2_v3v3fl(bpoly->coords_v2[j], bpoly->coords[j], axis, angle); + madd_v2_v2fl(bpoly->centroid_v2, bpoly->coords_v2[j], 1.0f / poly->totloop); + } + + is_poly_valid = isPolyValid(bpoly->coords_v2, poly->totloop); + + if (is_poly_valid != MOD_SDEF_BIND_RESULT_SUCCESS) { + freeBindData(bwdata); + data->success = is_poly_valid; + return NULL; + } + + bpoly->inside = isect_point_poly_v2(bpoly->point_v2, bpoly->coords_v2, poly->totloop, false); + + /* Initialize weight components */ + bpoly->weight_angular = 1.0f; + bpoly->weight_dist_proj = len_v2v2(bpoly->centroid_v2, bpoly->point_v2); + bpoly->weight_dist = len_v3v3(bpoly->centroid, point_co); + + avg_point_dist += bpoly->weight_dist; + + /* Compute centroid to mid-edge vectors */ + mid_v2_v2v2(bpoly->cent_edgemid_vecs_v2[0], + bpoly->coords_v2[bpoly->edge_vert_inds[0]], + bpoly->coords_v2[bpoly->corner_ind]); + + mid_v2_v2v2(bpoly->cent_edgemid_vecs_v2[1], + bpoly->coords_v2[bpoly->edge_vert_inds[1]], + bpoly->coords_v2[bpoly->corner_ind]); + + sub_v2_v2(bpoly->cent_edgemid_vecs_v2[0], bpoly->centroid_v2); + sub_v2_v2(bpoly->cent_edgemid_vecs_v2[1], bpoly->centroid_v2); + + /* Compute poly scales with respect to mid-edges, and normalize the vectors */ + bpoly->scales[0] = normalize_v2(bpoly->cent_edgemid_vecs_v2[0]); + bpoly->scales[1] = normalize_v2(bpoly->cent_edgemid_vecs_v2[1]); + + /* Compute the required polygon angles */ + bpoly->edgemid_angle = angle_normalized_v2v2(bpoly->cent_edgemid_vecs_v2[0], bpoly->cent_edgemid_vecs_v2[1]); + + sub_v2_v2v2(tmp_vec_v2, bpoly->coords_v2[bpoly->corner_ind], bpoly->centroid_v2); + normalize_v2(tmp_vec_v2); + + bpoly->corner_edgemid_angles[0] = angle_normalized_v2v2(tmp_vec_v2, bpoly->cent_edgemid_vecs_v2[0]); + bpoly->corner_edgemid_angles[1] = angle_normalized_v2v2(tmp_vec_v2, bpoly->cent_edgemid_vecs_v2[1]); + + /* Check for inifnite weights, and compute angular data otherwise */ + if (bpoly->weight_dist < FLT_EPSILON) { + inf_weight_flags |= MOD_SDEF_INFINITE_WEIGHT_DIST_PROJ; + inf_weight_flags |= MOD_SDEF_INFINITE_WEIGHT_DIST; + } + else if (bpoly->weight_dist_proj < FLT_EPSILON) { + inf_weight_flags |= MOD_SDEF_INFINITE_WEIGHT_DIST_PROJ; + } + else { + float cent_point_vec[2]; + + sub_v2_v2v2(cent_point_vec, bpoly->point_v2, bpoly->centroid_v2); + normalize_v2(cent_point_vec); + + bpoly->point_edgemid_angles[0] = angle_normalized_v2v2(cent_point_vec, bpoly->cent_edgemid_vecs_v2[0]); + bpoly->point_edgemid_angles[1] = angle_normalized_v2v2(cent_point_vec, bpoly->cent_edgemid_vecs_v2[1]); + } + } + } + } + + avg_point_dist /= bwdata->numpoly; + + /* If weights 1 and 2 are not infinite, loop over all adjacent edges again, + * and build adjacency dependent angle data (depends on all polygons having been computed) */ + if (!inf_weight_flags) { + for (vedge = vert_edges; vedge; vedge = vedge->next) { + SDefBindPoly *bpolys[2]; + const SDefEdgePolys *epolys; + float ang_weights[2]; + unsigned int edge_ind = vedge->index; + unsigned int edge_on_poly[2]; + + epolys = &edge_polys[edge_ind]; + + /* Find bind polys corresponding to the edge's adjacent polys */ + bpoly = bwdata->bind_polys; + + for (int i = 0, j = 0; (i < bwdata->numpoly) && (j < epolys->num); bpoly++, i++) { + if (ELEM(bpoly->index, epolys->polys[0], epolys->polys[1])) { + bpolys[j] = bpoly; + + if (bpoly->edge_inds[0] == edge_ind) { + edge_on_poly[j] = 0; + } + else { + edge_on_poly[j] = 1; + } + + j++; + } + } + + /* Compute angular weight component */ + if (epolys->num == 1) { + ang_weights[0] = computeAngularWeight(bpolys[0]->point_edgemid_angles[edge_on_poly[0]], bpolys[0]->edgemid_angle); + bpolys[0]->weight_angular *= ang_weights[0] * ang_weights[0]; + } + else if (epolys->num == 2) { + ang_weights[0] = computeAngularWeight(bpolys[0]->point_edgemid_angles[edge_on_poly[0]], bpolys[0]->edgemid_angle); + ang_weights[1] = computeAngularWeight(bpolys[1]->point_edgemid_angles[edge_on_poly[1]], bpolys[1]->edgemid_angle); + + bpolys[0]->weight_angular *= ang_weights[0] * ang_weights[1]; + bpolys[1]->weight_angular *= ang_weights[0] * ang_weights[1]; + } + } + } + + /* Compute scalings and falloff. + * Scale all weights if no infinite weight is found, + * scale only unprojected weight if projected weight is infinite, + * scale none if both are infinite. */ + if (!inf_weight_flags) { + bpoly = bwdata->bind_polys; + + for (int i = 0; i < bwdata->numpoly; bpoly++, i++) { + float corner_angle_weights[2]; + float scale_weight, sqr, inv_sqr; + + corner_angle_weights[0] = bpoly->point_edgemid_angles[0] / bpoly->corner_edgemid_angles[0]; + corner_angle_weights[1] = bpoly->point_edgemid_angles[1] / bpoly->corner_edgemid_angles[1]; + + if (isnan(corner_angle_weights[0]) || isnan(corner_angle_weights[1])) { + freeBindData(bwdata); + data->success = MOD_SDEF_BIND_RESULT_GENERIC_ERR; + return NULL; + } + + /* Find which edge the point is closer to */ + if (corner_angle_weights[0] < corner_angle_weights[1]) { + bpoly->dominant_edge = 0; + bpoly->dominant_angle_weight = corner_angle_weights[0]; + } + else { + bpoly->dominant_edge = 1; + bpoly->dominant_angle_weight = corner_angle_weights[1]; + } + + bpoly->dominant_angle_weight = sinf(bpoly->dominant_angle_weight * M_PI_2); + + /* Compute quadratic angular scale interpolation weight */ + scale_weight = bpoly->point_edgemid_angles[bpoly->dominant_edge] / bpoly->edgemid_angle; + scale_weight /= scale_weight + (bpoly->point_edgemid_angles[!bpoly->dominant_edge] / bpoly->edgemid_angle); + + sqr = scale_weight * scale_weight; + inv_sqr = 1.0f - scale_weight; + inv_sqr *= inv_sqr; + scale_weight = sqr / (sqr + inv_sqr); + + /* Compute interpolated scale (no longer need the individual scales, + * so simply storing the result over the scale in index zero) */ + bpoly->scales[0] = bpoly->scales[bpoly->dominant_edge] * (1.0f - scale_weight) + + bpoly->scales[!bpoly->dominant_edge] * scale_weight; + + /* Scale the point distance weights, and introduce falloff */ + bpoly->weight_dist_proj /= bpoly->scales[0]; + bpoly->weight_dist_proj = powf(bpoly->weight_dist_proj, data->falloff); + + bpoly->weight_dist /= avg_point_dist; + bpoly->weight_dist = powf(bpoly->weight_dist, data->falloff); + + /* Re-check for infinite weights, now that all scalings and interpolations are computed */ + if (bpoly->weight_dist < FLT_EPSILON) { + inf_weight_flags |= MOD_SDEF_INFINITE_WEIGHT_DIST_PROJ; + inf_weight_flags |= MOD_SDEF_INFINITE_WEIGHT_DIST; + } + else if (bpoly->weight_dist_proj < FLT_EPSILON) { + inf_weight_flags |= MOD_SDEF_INFINITE_WEIGHT_DIST_PROJ; + } + else if (bpoly->weight_angular < FLT_EPSILON) { + inf_weight_flags |= MOD_SDEF_INFINITE_WEIGHT_ANGULAR; + } + } + } + else if (!(inf_weight_flags & MOD_SDEF_INFINITE_WEIGHT_DIST)) { + bpoly = bwdata->bind_polys; + + for (int i = 0; i < bwdata->numpoly; bpoly++, i++) { + /* Scale the point distance weight by average point distance, and introduce falloff */ + bpoly->weight_dist /= avg_point_dist; + bpoly->weight_dist = powf(bpoly->weight_dist, data->falloff); + + /* Re-check for infinite weights, now that all scalings and interpolations are computed */ + if (bpoly->weight_dist < FLT_EPSILON) { + inf_weight_flags |= MOD_SDEF_INFINITE_WEIGHT_DIST; + } + } + } + + /* Final loop, to compute actual weights */ + bpoly = bwdata->bind_polys; + + for (int i = 0; i < bwdata->numpoly; bpoly++, i++) { + /* Weight computation from components */ + if (inf_weight_flags & MOD_SDEF_INFINITE_WEIGHT_DIST) { + bpoly->weight = bpoly->weight_dist < FLT_EPSILON ? 1.0f : 0.0f; + } + else if (inf_weight_flags & MOD_SDEF_INFINITE_WEIGHT_DIST_PROJ) { + bpoly->weight = bpoly->weight_dist_proj < FLT_EPSILON ? + 1.0f / bpoly->weight_dist : 0.0f; + } + else if (inf_weight_flags & MOD_SDEF_INFINITE_WEIGHT_ANGULAR) { + bpoly->weight = bpoly->weight_angular < FLT_EPSILON ? + 1.0f / bpoly->weight_dist_proj / bpoly->weight_dist : 0.0f; + } + else { + bpoly->weight = 1.0f / bpoly->weight_angular / + bpoly->weight_dist_proj / + bpoly->weight_dist; + } + + tot_weight += bpoly->weight; + } + + bpoly = bwdata->bind_polys; + + for (int i = 0; i < bwdata->numpoly; bpoly++, i++) { + bpoly->weight /= tot_weight; + + /* Evaluate if this poly is relevant to bind */ + /* Even though the weights should add up to 1.0, + * the losses of weights smaller than epsilon here + * should be negligible... */ + if (bpoly->weight >= FLT_EPSILON) { + if (bpoly->inside) { + bwdata->numbinds += 1; + } + else { + if (bpoly->dominant_angle_weight < FLT_EPSILON || 1.0f - bpoly->dominant_angle_weight < FLT_EPSILON) { + bwdata->numbinds += 1; + } + else { + bwdata->numbinds += 2; + } + } + } + } + + return bwdata; +} + +BLI_INLINE float computeNormalDisplacement(const float point_co[3], const float point_co_proj[3], const float normal[3]) +{ + float disp_vec[3]; + float normal_dist; + + sub_v3_v3v3(disp_vec, point_co, point_co_proj); + normal_dist = len_v3(disp_vec); + + if (dot_v3v3(disp_vec, normal) < 0) { + normal_dist *= -1; + } + + return normal_dist; +} + +static void bindVert(void *userdata, void *UNUSED(userdata_chunk), const int index, const int UNUSED(threadid)) +{ + SDefBindCalcData * const data = (SDefBindCalcData *)userdata; + float point_co[3]; + float point_co_proj[3]; + + SDefBindWeightData *bwdata; + SDefVert *sdvert = data->bind_verts + index; + SDefBindPoly *bpoly; + SDefBind *sdbind; + + if (data->success != MOD_SDEF_BIND_RESULT_SUCCESS) { + sdvert->binds = NULL; + sdvert->numbinds = 0; + return; + } + + copy_v3_v3(point_co, data->vertexCos[index]); + bwdata = computeBindWeights(data, point_co); + + if (bwdata == NULL) { + sdvert->binds = NULL; + sdvert->numbinds = 0; + return; + } + + sdvert->binds = MEM_callocN(sizeof(*sdvert->binds) * bwdata->numbinds, "SDefVertBindData"); + if (sdvert->binds == NULL) { + data->success = MOD_SDEF_BIND_RESULT_MEM_ERR; + sdvert->numbinds = 0; + return; + } + + sdvert->numbinds = bwdata->numbinds; + + sdbind = sdvert->binds; + + bpoly = bwdata->bind_polys; + + for (int i = 0; i < bwdata->numbinds; bpoly++) { + if (bpoly->weight >= FLT_EPSILON) { + if (bpoly->inside) { + const MLoop *loop = &data->mloop[bpoly->loopstart]; + + sdbind->influence = bpoly->weight; + sdbind->numverts = bpoly->numverts; + + sdbind->mode = MOD_SDEF_MODE_NGON; + sdbind->vert_weights = MEM_mallocN(sizeof(*sdbind->vert_weights) * bpoly->numverts, "SDefNgonVertWeights"); + if (sdbind->vert_weights == NULL) { + data->success = MOD_SDEF_BIND_RESULT_MEM_ERR; + return; + } + + sdbind->vert_inds = MEM_mallocN(sizeof(*sdbind->vert_inds) * bpoly->numverts, "SDefNgonVertInds"); + if (sdbind->vert_inds == NULL) { + data->success = MOD_SDEF_BIND_RESULT_MEM_ERR; + return; + } + + interp_weights_poly_v2(sdbind->vert_weights, bpoly->coords_v2, bpoly->numverts, bpoly->point_v2); + + /* Reproject vert based on weights and original poly verts, to reintroduce poly non-planarity */ + zero_v3(point_co_proj); + for (int j = 0; j < bpoly->numverts; j++, loop++) { + madd_v3_v3fl(point_co_proj, bpoly->coords[j], sdbind->vert_weights[j]); + sdbind->vert_inds[j] = loop->v; + } + + sdbind->normal_dist = computeNormalDisplacement(point_co, point_co_proj, bpoly->normal); + + sdbind++; + i++; + } + else { + float tmp_vec[3]; + float cent[3], norm[3]; + float v1[3], v2[3], v3[3]; + + if (1.0f - bpoly->dominant_angle_weight >= FLT_EPSILON) { + sdbind->influence = bpoly->weight * (1.0f - bpoly->dominant_angle_weight); + sdbind->numverts = bpoly->numverts; + + sdbind->mode = MOD_SDEF_MODE_CENTROID; + sdbind->vert_weights = MEM_mallocN(sizeof(*sdbind->vert_weights) * 3, "SDefCentVertWeights"); + if (sdbind->vert_weights == NULL) { + data->success = MOD_SDEF_BIND_RESULT_MEM_ERR; + return; + } + + sdbind->vert_inds = MEM_mallocN(sizeof(*sdbind->vert_inds) * bpoly->numverts, "SDefCentVertInds"); + if (sdbind->vert_inds == NULL) { + data->success = MOD_SDEF_BIND_RESULT_MEM_ERR; + return; + } + + sortPolyVertsEdge(sdbind->vert_inds, &data->mloop[bpoly->loopstart], + bpoly->edge_inds[bpoly->dominant_edge], bpoly->numverts); + + copy_v3_v3(v1, data->mvert[sdbind->vert_inds[0]].co); + copy_v3_v3(v2, data->mvert[sdbind->vert_inds[1]].co); + copy_v3_v3(v3, bpoly->centroid); + + mid_v3_v3v3v3(cent, v1, v2, v3); + normal_tri_v3(norm, v1, v2, v3); + + add_v3_v3v3(tmp_vec, point_co, bpoly->normal); + + /* We are sure the line is not parallel to the plane. + * Checking return value just to avoid warning... */ + if (!isect_line_plane_v3(point_co_proj, point_co, tmp_vec, cent, norm)) { + BLI_assert(false); + } + + interp_weights_tri_v3(sdbind->vert_weights, v1, v2, v3, point_co_proj); + + sdbind->normal_dist = computeNormalDisplacement(point_co, point_co_proj, bpoly->normal); + + sdbind++; + i++; + } + + if (bpoly->dominant_angle_weight >= FLT_EPSILON) { + sdbind->influence = bpoly->weight * bpoly->dominant_angle_weight; + sdbind->numverts = bpoly->numverts; + + sdbind->mode = MOD_SDEF_MODE_LOOPTRI; + sdbind->vert_weights = MEM_mallocN(sizeof(*sdbind->vert_weights) * 3, "SDefTriVertWeights"); + if (sdbind->vert_weights == NULL) { + data->success = MOD_SDEF_BIND_RESULT_MEM_ERR; + return; + } + + sdbind->vert_inds = MEM_mallocN(sizeof(*sdbind->vert_inds) * bpoly->numverts, "SDefTriVertInds"); + if (sdbind->vert_inds == NULL) { + data->success = MOD_SDEF_BIND_RESULT_MEM_ERR; + return; + } + + sortPolyVertsTri(sdbind->vert_inds, &data->mloop[bpoly->loopstart], bpoly->edge_vert_inds[0], bpoly->numverts); + + copy_v3_v3(v1, data->mvert[sdbind->vert_inds[0]].co); + copy_v3_v3(v2, data->mvert[sdbind->vert_inds[1]].co); + copy_v3_v3(v3, data->mvert[sdbind->vert_inds[2]].co); + + mid_v3_v3v3v3(cent, v1, v2, v3); + normal_tri_v3(norm, v1, v2, v3); + + add_v3_v3v3(tmp_vec, point_co, bpoly->normal); + + /* We are sure the line is not parallel to the plane. + * Checking return value just to avoid warning... */ + if (!isect_line_plane_v3(point_co_proj, point_co, tmp_vec, cent, norm)) { + BLI_assert(false); + } + + interp_weights_tri_v3(sdbind->vert_weights, v1, v2, v3, point_co_proj); + + sdbind->normal_dist = computeNormalDisplacement(point_co, point_co_proj, bpoly->normal); + + sdbind++; + i++; + } + } + } + } + + freeBindData(bwdata); +} + +static bool surfacedeformBind(SurfaceDeformModifierData *smd, float (*vertexCos)[3], + unsigned int numverts, unsigned int tnumpoly, DerivedMesh *tdm) +{ + BVHTreeFromMesh treeData = {NULL}; + const MPoly *mpoly = tdm->getPolyArray(tdm); + const MEdge *medge = tdm->getEdgeArray(tdm); + const MLoop *mloop = tdm->getLoopArray(tdm); + unsigned int tnumedges = tdm->getNumEdges(tdm); + unsigned int tnumverts = tdm->getNumVerts(tdm); + int adj_result; + SDefAdjacencyArray *vert_edges; + SDefAdjacency *adj_array; + SDefEdgePolys *edge_polys; + + vert_edges = MEM_callocN(sizeof(*vert_edges) * tnumverts, "SDefVertEdgeMap"); + if (vert_edges == NULL) { + modifier_setError((ModifierData *)smd, "Out of memory"); + return false; + } + + adj_array = MEM_mallocN(sizeof(*adj_array) * tnumedges * 2, "SDefVertEdge"); + if (adj_array == NULL) { + modifier_setError((ModifierData *)smd, "Out of memory"); + MEM_freeN(vert_edges); + return false; + } + + edge_polys = MEM_callocN(sizeof(*edge_polys) * tnumedges, "SDefEdgeFaceMap"); + if (edge_polys == NULL) { + modifier_setError((ModifierData *)smd, "Out of memory"); + MEM_freeN(vert_edges); + MEM_freeN(adj_array); + return false; + } + + smd->verts = MEM_mallocN(sizeof(*smd->verts) * numverts, "SDefBindVerts"); + if (smd->verts == NULL) { + modifier_setError((ModifierData *)smd, "Out of memory"); + freeAdjacencyMap(vert_edges, adj_array, edge_polys); + return false; + } + + bvhtree_from_mesh_looptri(&treeData, tdm, 0.0, 2, 6); + if (treeData.tree == NULL) { + modifier_setError((ModifierData *)smd, "Out of memory"); + freeAdjacencyMap(vert_edges, adj_array, edge_polys); + MEM_freeN(smd->verts); + smd->verts = NULL; + return false; + } + + adj_result = buildAdjacencyMap(mpoly, medge, mloop, tnumpoly, tnumedges, vert_edges, adj_array, edge_polys); + + if (adj_result == MOD_SDEF_BIND_RESULT_NONMANY_ERR) { + modifier_setError((ModifierData *)smd, "Target has edges with more than two polys"); + freeAdjacencyMap(vert_edges, adj_array, edge_polys); + free_bvhtree_from_mesh(&treeData); + MEM_freeN(smd->verts); + smd->verts = NULL; + return false; + } + + smd->numverts = numverts; + smd->numpoly = tnumpoly; + + SDefBindCalcData data = {.treeData = &treeData, + .vert_edges = vert_edges, + .edge_polys = edge_polys, + .mpoly = mpoly, + .medge = medge, + .mloop = mloop, + .looptri = tdm->getLoopTriArray(tdm), + .mvert = tdm->getVertArray(tdm), + .bind_verts = smd->verts, + .vertexCos = vertexCos, + .falloff = smd->falloff, + .success = MOD_SDEF_BIND_RESULT_SUCCESS}; + + BLI_task_parallel_range_ex(0, numverts, &data, NULL, 0, bindVert, + numverts > 10000, false); + + if (data.success == MOD_SDEF_BIND_RESULT_MEM_ERR) { + modifier_setError((ModifierData *)smd, "Out of memory"); + freeData((ModifierData *)smd); + } + else if (data.success == MOD_SDEF_BIND_RESULT_NONMANY_ERR) { + modifier_setError((ModifierData *)smd, "Target has edges with more than two polys"); + freeData((ModifierData *)smd); + } + else if (data.success == MOD_SDEF_BIND_RESULT_CONCAVE_ERR) { + modifier_setError((ModifierData *)smd, "Target contains concave polys"); + freeData((ModifierData *)smd); + } + else if (data.success == MOD_SDEF_BIND_RESULT_OVERLAP_ERR) { + modifier_setError((ModifierData *)smd, "Target contains overlapping verts"); + freeData((ModifierData *)smd); + } + else if (data.success == MOD_SDEF_BIND_RESULT_GENERIC_ERR) { + /* I know this message is vague, but I could not think of a way + * to explain this whith a reasonably sized message. + * Though it shouldn't really matter all that much, + * because this is very unlikely to occur */ + modifier_setError((ModifierData *)smd, "Target contains invalid polys"); + freeData((ModifierData *)smd); + } + + freeAdjacencyMap(vert_edges, adj_array, edge_polys); + free_bvhtree_from_mesh(&treeData); + + return data.success == 1; +} + +static void deformVert(void *userdata, void *UNUSED(userdata_chunk), const int index, const int UNUSED(threadid)) +{ + const SDefDeformData * const data = (SDefDeformData *)userdata; + const SDefBind *sdbind = data->bind_verts[index].binds; + const MVert * const mvert = data->mvert; + float * const vertexCos = data->vertexCos[index]; + float norm[3], temp[3]; + + zero_v3(vertexCos); + + for (int j = 0; j < data->bind_verts[index].numbinds; j++, sdbind++) { + /* Mode-generic operations (allocate poly coordinates) */ + float (*coords)[3] = MEM_mallocN(sizeof(*coords) * sdbind->numverts, "SDefDoPolyCoords"); + + for (int k = 0; k < sdbind->numverts; k++) { + copy_v3_v3(coords[k], mvert[sdbind->vert_inds[k]].co); + } + + normal_poly_v3(norm, coords, sdbind->numverts); + zero_v3(temp); + + /* ---------- looptri mode ---------- */ + if (sdbind->mode == MOD_SDEF_MODE_LOOPTRI) { + madd_v3_v3fl(temp, mvert[sdbind->vert_inds[0]].co, sdbind->vert_weights[0]); + madd_v3_v3fl(temp, mvert[sdbind->vert_inds[1]].co, sdbind->vert_weights[1]); + madd_v3_v3fl(temp, mvert[sdbind->vert_inds[2]].co, sdbind->vert_weights[2]); + } + else { + /* ---------- ngon mode ---------- */ + if (sdbind->mode == MOD_SDEF_MODE_NGON) { + for (int k = 0; k < sdbind->numverts; k++) { + madd_v3_v3fl(temp, coords[k], sdbind->vert_weights[k]); + } + } + + /* ---------- centroid mode ---------- */ + else if (sdbind->mode == MOD_SDEF_MODE_CENTROID) { + float cent[3]; + mid_v3_v3_array(cent, coords, sdbind->numverts); + + madd_v3_v3fl(temp, mvert[sdbind->vert_inds[0]].co, sdbind->vert_weights[0]); + madd_v3_v3fl(temp, mvert[sdbind->vert_inds[1]].co, sdbind->vert_weights[1]); + madd_v3_v3fl(temp, cent, sdbind->vert_weights[2]); + } + } + + MEM_freeN(coords); + + /* Apply normal offset (generic for all modes) */ + madd_v3_v3fl(temp, norm, sdbind->normal_dist); + + madd_v3_v3fl(vertexCos, temp, sdbind->influence); + } +} + +static void surfacedeformModifier_do(ModifierData *md, float (*vertexCos)[3], unsigned int numverts) +{ + SurfaceDeformModifierData *smd = (SurfaceDeformModifierData *)md; + DerivedMesh *tdm; + unsigned int tnumpoly; + + /* Exit function if bind flag is not set (free bind data if any) */ + if (!(smd->flags & MOD_SDEF_BIND)) { + freeData(md); + return; + } + + /* Handle target mesh both in and out of edit mode */ + if (smd->target == md->scene->obedit) { + BMEditMesh *em = BKE_editmesh_from_object(smd->target); + tdm = em->derivedFinal; + } + else { + tdm = smd->target->derivedFinal; + } + + tnumpoly = tdm->getNumPolys(tdm); + + /* If not bound, execute bind */ + if (!(smd->verts)) { + if (!surfacedeformBind(smd, vertexCos, numverts, tnumpoly, tdm)) { + smd->flags &= ~MOD_SDEF_BIND; + return; + } + } + + /* Poly count checks */ + if (smd->numverts != numverts) { + modifier_setError(md, "Verts changed from %u to %u", smd->numverts, numverts); + tdm->release(tdm); + return; + } + else if (smd->numpoly != tnumpoly) { + modifier_setError(md, "Target polygons changed from %u to %u", smd->numpoly, tnumpoly); + tdm->release(tdm); + return; + } + + /* Actual vertex location update starts here */ + SDefDeformData data = {.bind_verts = smd->verts, + .mvert = tdm->getVertArray(tdm), + .vertexCos = vertexCos}; + + BLI_task_parallel_range_ex(0, numverts, &data, NULL, 0, deformVert, + numverts > 10000, false); + + tdm->release(tdm); +} + +static void deformVerts(ModifierData *md, Object *UNUSED(ob), + DerivedMesh *UNUSED(derivedData), + float (*vertexCos)[3], int numVerts, + ModifierApplyFlag UNUSED(flag)) +{ + surfacedeformModifier_do(md, vertexCos, numVerts); +} + +static void deformVertsEM(ModifierData *md, Object *UNUSED(ob), + struct BMEditMesh *UNUSED(editData), + DerivedMesh *UNUSED(derivedData), + float (*vertexCos)[3], int numVerts) +{ + surfacedeformModifier_do(md, vertexCos, numVerts); +} + +static bool isDisabled(ModifierData *md, int UNUSED(useRenderParams)) +{ + SurfaceDeformModifierData *smd = (SurfaceDeformModifierData *)md; + + return !smd->target; +} + +ModifierTypeInfo modifierType_SurfaceDeform = { + /* name */ "Surface Deform", + /* structName */ "SurfaceDeformModifierData", + /* structSize */ sizeof(SurfaceDeformModifierData), + /* type */ eModifierTypeType_OnlyDeform, + /* flags */ eModifierTypeFlag_AcceptsMesh | + eModifierTypeFlag_SupportsEditmode, + + /* copyData */ copyData, + /* deformVerts */ deformVerts, + /* deformMatrices */ NULL, + /* deformVertsEM */ deformVertsEM, + /* deformMatricesEM */ NULL, + /* applyModifier */ NULL, + /* applyModifierEM */ NULL, + /* initData */ initData, + /* requiredDataMask */ NULL, + /* freeData */ freeData, + /* isDisabled */ isDisabled, + /* updateDepgraph */ updateDepgraph, + /* updateDepsgraph */ updateDepsgraph, + /* dependsOnTime */ NULL, + /* dependsOnNormals */ NULL, + /* foreachObjectLink */ foreachObjectLink, + /* foreachIDLink */ NULL, + /* foreachTexLink */ NULL, +}; |