diff options
Diffstat (limited to 'source/blender/bmesh/operators/bmo_smooth_laplacian.c')
-rw-r--r-- | source/blender/bmesh/operators/bmo_smooth_laplacian.c | 910 |
1 files changed, 461 insertions, 449 deletions
diff --git a/source/blender/bmesh/operators/bmo_smooth_laplacian.c b/source/blender/bmesh/operators/bmo_smooth_laplacian.c index bf376a757ef..85f62241782 100644 --- a/source/blender/bmesh/operators/bmo_smooth_laplacian.c +++ b/source/blender/bmesh/operators/bmo_smooth_laplacian.c @@ -26,7 +26,6 @@ #include "eigen_capi.h" - #include "bmesh.h" #include "intern/bmesh_operators_private.h" /* own include */ @@ -37,23 +36,23 @@ #define SMOOTH_LAPLACIAN_MIN_EDGE_PERCENTAGE 0.15f struct BLaplacianSystem { - float *eweights; /* Length weights per Edge */ - float (*fweights)[3]; /* Cotangent weights per face */ - float *ring_areas; /* Total area per ring*/ - float *vlengths; /* Total sum of lengths(edges) per vertice*/ - float *vweights; /* Total sum of weights per vertice*/ - int numEdges; /* Number of edges*/ - int numFaces; /* Number of faces*/ - int numVerts; /* Number of verts*/ - short *zerola; /* Is zero area or length*/ - - /* Pointers to data*/ - BMesh *bm; - BMOperator *op; - LinearSolver *context; - - /*Data*/ - float min_area; + float *eweights; /* Length weights per Edge */ + float (*fweights)[3]; /* Cotangent weights per face */ + float *ring_areas; /* Total area per ring*/ + float *vlengths; /* Total sum of lengths(edges) per vertice*/ + float *vweights; /* Total sum of weights per vertice*/ + int numEdges; /* Number of edges*/ + int numFaces; /* Number of faces*/ + int numVerts; /* Number of verts*/ + short *zerola; /* Is zero area or length*/ + + /* Pointers to data*/ + BMesh *bm; + BMOperator *op; + LinearSolver *context; + + /*Data*/ + float min_area; }; typedef struct BLaplacianSystem LaplacianSystem; @@ -64,87 +63,89 @@ static void delete_laplacian_system(LaplacianSystem *sys); static void delete_void_pointer(void *data); static void fill_laplacian_matrix(LaplacianSystem *sys); static void memset_laplacian_system(LaplacianSystem *sys, int val); -static void validate_solution(LaplacianSystem *sys, int usex, int usey, int usez, int preserve_volume); -static void volume_preservation(BMOperator *op, float vini, float vend, int usex, int usey, int usez); +static void validate_solution( + LaplacianSystem *sys, int usex, int usey, int usez, int preserve_volume); +static void volume_preservation( + BMOperator *op, float vini, float vend, int usex, int usey, int usez); static void delete_void_pointer(void *data) { - if (data) { - MEM_freeN(data); - } + if (data) { + MEM_freeN(data); + } } static void delete_laplacian_system(LaplacianSystem *sys) { - delete_void_pointer(sys->eweights); - delete_void_pointer(sys->fweights); - delete_void_pointer(sys->ring_areas); - delete_void_pointer(sys->vlengths); - delete_void_pointer(sys->vweights); - delete_void_pointer(sys->zerola); - if (sys->context) { - EIG_linear_solver_delete(sys->context); - } - sys->bm = NULL; - sys->op = NULL; - MEM_freeN(sys); + delete_void_pointer(sys->eweights); + delete_void_pointer(sys->fweights); + delete_void_pointer(sys->ring_areas); + delete_void_pointer(sys->vlengths); + delete_void_pointer(sys->vweights); + delete_void_pointer(sys->zerola); + if (sys->context) { + EIG_linear_solver_delete(sys->context); + } + sys->bm = NULL; + sys->op = NULL; + MEM_freeN(sys); } static void memset_laplacian_system(LaplacianSystem *sys, int val) { - memset(sys->eweights, val, sizeof(float) * sys->numEdges); - memset(sys->fweights, val, sizeof(float) * sys->numFaces * 3); - memset(sys->ring_areas, val, sizeof(float) * sys->numVerts); - memset(sys->vlengths, val, sizeof(float) * sys->numVerts); - memset(sys->vweights, val, sizeof(float) * sys->numVerts); - memset(sys->zerola, val, sizeof(short) * sys->numVerts); + memset(sys->eweights, val, sizeof(float) * sys->numEdges); + memset(sys->fweights, val, sizeof(float) * sys->numFaces * 3); + memset(sys->ring_areas, val, sizeof(float) * sys->numVerts); + memset(sys->vlengths, val, sizeof(float) * sys->numVerts); + memset(sys->vweights, val, sizeof(float) * sys->numVerts); + memset(sys->zerola, val, sizeof(short) * sys->numVerts); } static LaplacianSystem *init_laplacian_system(int a_numEdges, int a_numFaces, int a_numVerts) { - LaplacianSystem *sys; - sys = MEM_callocN(sizeof(LaplacianSystem), "ModLaplSmoothSystem"); - sys->numEdges = a_numEdges; - sys->numFaces = a_numFaces; - sys->numVerts = a_numVerts; - - sys->eweights = MEM_callocN(sizeof(float) * sys->numEdges, "ModLaplSmoothEWeight"); - if (!sys->eweights) { - delete_laplacian_system(sys); - return NULL; - } - - sys->fweights = MEM_callocN(sizeof(float) * 3 * sys->numFaces, "ModLaplSmoothFWeight"); - if (!sys->fweights) { - delete_laplacian_system(sys); - return NULL; - } - - sys->ring_areas = MEM_callocN(sizeof(float) * sys->numVerts, "ModLaplSmoothRingAreas"); - if (!sys->ring_areas) { - delete_laplacian_system(sys); - return NULL; - } - - sys->vlengths = MEM_callocN(sizeof(float) * sys->numVerts, "ModLaplSmoothVlengths"); - if (!sys->vlengths) { - delete_laplacian_system(sys); - return NULL; - } - - sys->vweights = MEM_callocN(sizeof(float) * sys->numVerts, "ModLaplSmoothVweights"); - if (!sys->vweights) { - delete_laplacian_system(sys); - return NULL; - } - - sys->zerola = MEM_callocN(sizeof(short) * sys->numVerts, "ModLaplSmoothZeloa"); - if (!sys->zerola) { - delete_laplacian_system(sys); - return NULL; - } - - return sys; + LaplacianSystem *sys; + sys = MEM_callocN(sizeof(LaplacianSystem), "ModLaplSmoothSystem"); + sys->numEdges = a_numEdges; + sys->numFaces = a_numFaces; + sys->numVerts = a_numVerts; + + sys->eweights = MEM_callocN(sizeof(float) * sys->numEdges, "ModLaplSmoothEWeight"); + if (!sys->eweights) { + delete_laplacian_system(sys); + return NULL; + } + + sys->fweights = MEM_callocN(sizeof(float) * 3 * sys->numFaces, "ModLaplSmoothFWeight"); + if (!sys->fweights) { + delete_laplacian_system(sys); + return NULL; + } + + sys->ring_areas = MEM_callocN(sizeof(float) * sys->numVerts, "ModLaplSmoothRingAreas"); + if (!sys->ring_areas) { + delete_laplacian_system(sys); + return NULL; + } + + sys->vlengths = MEM_callocN(sizeof(float) * sys->numVerts, "ModLaplSmoothVlengths"); + if (!sys->vlengths) { + delete_laplacian_system(sys); + return NULL; + } + + sys->vweights = MEM_callocN(sizeof(float) * sys->numVerts, "ModLaplSmoothVweights"); + if (!sys->vweights) { + delete_laplacian_system(sys); + return NULL; + } + + sys->zerola = MEM_callocN(sizeof(short) * sys->numVerts, "ModLaplSmoothZeloa"); + if (!sys->zerola) { + delete_laplacian_system(sys); + return NULL; + } + + return sys; } /** @@ -165,387 +166,398 @@ static LaplacianSystem *init_laplacian_system(int a_numEdges, int a_numFaces, in static void init_laplacian_matrix(LaplacianSystem *sys) { - float areaf; - float *v1, *v2, *v3, *v4; - float w1, w2, w3, w4; - int i, j; - bool has_4_vert; - uint idv1, idv2, idv3, idv4, idv[4]; - BMEdge *e; - BMFace *f; - BMIter eiter; - BMIter fiter; - BMIter vi; - BMVert *vn; - BMVert *vf[4]; - - BM_ITER_MESH_INDEX (e, &eiter, sys->bm, BM_EDGES_OF_MESH, j) { - if (!BM_elem_flag_test(e, BM_ELEM_SELECT) && BM_edge_is_boundary(e)) { - v1 = e->v1->co; - v2 = e->v2->co; - idv1 = BM_elem_index_get(e->v1); - idv2 = BM_elem_index_get(e->v2); - - w1 = len_v3v3(v1, v2); - if (w1 > sys->min_area) { - w1 = 1.0f / w1; - i = BM_elem_index_get(e); - sys->eweights[i] = w1; - sys->vlengths[idv1] += w1; - sys->vlengths[idv2] += w1; - } - else { - sys->zerola[idv1] = 1; - sys->zerola[idv2] = 1; - } - } - } - - BM_ITER_MESH (f, &fiter, sys->bm, BM_FACES_OF_MESH) { - if (BM_elem_flag_test(f, BM_ELEM_SELECT)) { - - BM_ITER_ELEM_INDEX (vn, &vi, f, BM_VERTS_OF_FACE, i) { - vf[i] = vn; - } - has_4_vert = (i == 4) ? 1 : 0; - idv1 = BM_elem_index_get(vf[0]); - idv2 = BM_elem_index_get(vf[1]); - idv3 = BM_elem_index_get(vf[2]); - idv4 = has_4_vert ? BM_elem_index_get(vf[3]) : 0; - - v1 = vf[0]->co; - v2 = vf[1]->co; - v3 = vf[2]->co; - v4 = has_4_vert ? vf[3]->co : NULL; - - if (has_4_vert) { - areaf = area_quad_v3(v1, v2, v3, v4); - } - else { - areaf = area_tri_v3(v1, v2, v3); - } - - if (fabsf(areaf) < sys->min_area) { - sys->zerola[idv1] = 1; - sys->zerola[idv2] = 1; - sys->zerola[idv3] = 1; - if (has_4_vert) { - sys->zerola[idv4] = 1; - } - } - - sys->ring_areas[idv1] += areaf; - sys->ring_areas[idv2] += areaf; - sys->ring_areas[idv3] += areaf; - if (has_4_vert) { - sys->ring_areas[idv4] += areaf; - } - - if (has_4_vert) { - - idv[0] = idv1; - idv[1] = idv2; - idv[2] = idv3; - idv[3] = idv4; - - for (j = 0; j < 4; j++) { - idv1 = idv[j]; - idv2 = idv[(j + 1) % 4]; - idv3 = idv[(j + 2) % 4]; - idv4 = idv[(j + 3) % 4]; - - v1 = vf[j]->co; - v2 = vf[(j + 1) % 4]->co; - v3 = vf[(j + 2) % 4]->co; - v4 = vf[(j + 3) % 4]->co; - - w2 = cotangent_tri_weight_v3(v4, v1, v2) + cotangent_tri_weight_v3(v3, v1, v2); - w3 = cotangent_tri_weight_v3(v2, v3, v1) + cotangent_tri_weight_v3(v4, v1, v3); - w4 = cotangent_tri_weight_v3(v2, v4, v1) + cotangent_tri_weight_v3(v3, v4, v1); - - sys->vweights[idv1] += (w2 + w3 + w4) / 4.0f; - } - } - else { - i = BM_elem_index_get(f); - - w1 = cotangent_tri_weight_v3(v1, v2, v3); - w2 = cotangent_tri_weight_v3(v2, v3, v1); - w3 = cotangent_tri_weight_v3(v3, v1, v2); - - sys->fweights[i][0] += w1; - sys->fweights[i][1] += w2; - sys->fweights[i][2] += w3; - - sys->vweights[idv1] += w2 + w3; - sys->vweights[idv2] += w1 + w3; - sys->vweights[idv3] += w1 + w2; - } - } - } + float areaf; + float *v1, *v2, *v3, *v4; + float w1, w2, w3, w4; + int i, j; + bool has_4_vert; + uint idv1, idv2, idv3, idv4, idv[4]; + BMEdge *e; + BMFace *f; + BMIter eiter; + BMIter fiter; + BMIter vi; + BMVert *vn; + BMVert *vf[4]; + + BM_ITER_MESH_INDEX (e, &eiter, sys->bm, BM_EDGES_OF_MESH, j) { + if (!BM_elem_flag_test(e, BM_ELEM_SELECT) && BM_edge_is_boundary(e)) { + v1 = e->v1->co; + v2 = e->v2->co; + idv1 = BM_elem_index_get(e->v1); + idv2 = BM_elem_index_get(e->v2); + + w1 = len_v3v3(v1, v2); + if (w1 > sys->min_area) { + w1 = 1.0f / w1; + i = BM_elem_index_get(e); + sys->eweights[i] = w1; + sys->vlengths[idv1] += w1; + sys->vlengths[idv2] += w1; + } + else { + sys->zerola[idv1] = 1; + sys->zerola[idv2] = 1; + } + } + } + + BM_ITER_MESH (f, &fiter, sys->bm, BM_FACES_OF_MESH) { + if (BM_elem_flag_test(f, BM_ELEM_SELECT)) { + + BM_ITER_ELEM_INDEX(vn, &vi, f, BM_VERTS_OF_FACE, i) + { + vf[i] = vn; + } + has_4_vert = (i == 4) ? 1 : 0; + idv1 = BM_elem_index_get(vf[0]); + idv2 = BM_elem_index_get(vf[1]); + idv3 = BM_elem_index_get(vf[2]); + idv4 = has_4_vert ? BM_elem_index_get(vf[3]) : 0; + + v1 = vf[0]->co; + v2 = vf[1]->co; + v3 = vf[2]->co; + v4 = has_4_vert ? vf[3]->co : NULL; + + if (has_4_vert) { + areaf = area_quad_v3(v1, v2, v3, v4); + } + else { + areaf = area_tri_v3(v1, v2, v3); + } + + if (fabsf(areaf) < sys->min_area) { + sys->zerola[idv1] = 1; + sys->zerola[idv2] = 1; + sys->zerola[idv3] = 1; + if (has_4_vert) { + sys->zerola[idv4] = 1; + } + } + + sys->ring_areas[idv1] += areaf; + sys->ring_areas[idv2] += areaf; + sys->ring_areas[idv3] += areaf; + if (has_4_vert) { + sys->ring_areas[idv4] += areaf; + } + + if (has_4_vert) { + + idv[0] = idv1; + idv[1] = idv2; + idv[2] = idv3; + idv[3] = idv4; + + for (j = 0; j < 4; j++) { + idv1 = idv[j]; + idv2 = idv[(j + 1) % 4]; + idv3 = idv[(j + 2) % 4]; + idv4 = idv[(j + 3) % 4]; + + v1 = vf[j]->co; + v2 = vf[(j + 1) % 4]->co; + v3 = vf[(j + 2) % 4]->co; + v4 = vf[(j + 3) % 4]->co; + + w2 = cotangent_tri_weight_v3(v4, v1, v2) + cotangent_tri_weight_v3(v3, v1, v2); + w3 = cotangent_tri_weight_v3(v2, v3, v1) + cotangent_tri_weight_v3(v4, v1, v3); + w4 = cotangent_tri_weight_v3(v2, v4, v1) + cotangent_tri_weight_v3(v3, v4, v1); + + sys->vweights[idv1] += (w2 + w3 + w4) / 4.0f; + } + } + else { + i = BM_elem_index_get(f); + + w1 = cotangent_tri_weight_v3(v1, v2, v3); + w2 = cotangent_tri_weight_v3(v2, v3, v1); + w3 = cotangent_tri_weight_v3(v3, v1, v2); + + sys->fweights[i][0] += w1; + sys->fweights[i][1] += w2; + sys->fweights[i][2] += w3; + + sys->vweights[idv1] += w2 + w3; + sys->vweights[idv2] += w1 + w3; + sys->vweights[idv3] += w1 + w2; + } + } + } } static void fill_laplacian_matrix(LaplacianSystem *sys) { - float *v1, *v2, *v3, *v4; - float w2, w3, w4; - int i, j; - bool has_4_vert; - uint idv1, idv2, idv3, idv4, idv[4]; - - BMEdge *e; - BMFace *f; - BMIter eiter; - BMIter fiter; - BMIter vi; - BMVert *vn; - BMVert *vf[4]; - - BM_ITER_MESH (f, &fiter, sys->bm, BM_FACES_OF_MESH) { - if (BM_elem_flag_test(f, BM_ELEM_SELECT)) { - BM_ITER_ELEM_INDEX (vn, &vi, f, BM_VERTS_OF_FACE, i) { - vf[i] = vn; - } - has_4_vert = (i == 4) ? 1 : 0; - if (has_4_vert) { - idv[0] = BM_elem_index_get(vf[0]); - idv[1] = BM_elem_index_get(vf[1]); - idv[2] = BM_elem_index_get(vf[2]); - idv[3] = BM_elem_index_get(vf[3]); - for (j = 0; j < 4; j++) { - idv1 = idv[j]; - idv2 = idv[(j + 1) % 4]; - idv3 = idv[(j + 2) % 4]; - idv4 = idv[(j + 3) % 4]; - - v1 = vf[j]->co; - v2 = vf[(j + 1) % 4]->co; - v3 = vf[(j + 2) % 4]->co; - v4 = vf[(j + 3) % 4]->co; - - w2 = cotangent_tri_weight_v3(v4, v1, v2) + cotangent_tri_weight_v3(v3, v1, v2); - w3 = cotangent_tri_weight_v3(v2, v3, v1) + cotangent_tri_weight_v3(v4, v1, v3); - w4 = cotangent_tri_weight_v3(v2, v4, v1) + cotangent_tri_weight_v3(v3, v4, v1); - - w2 = w2 / 4.0f; - w3 = w3 / 4.0f; - w4 = w4 / 4.0f; - - if (!vert_is_boundary(vf[j]) && sys->zerola[idv1] == 0) { - EIG_linear_solver_matrix_add(sys->context, idv1, idv2, w2 * sys->vweights[idv1]); - EIG_linear_solver_matrix_add(sys->context, idv1, idv3, w3 * sys->vweights[idv1]); - EIG_linear_solver_matrix_add(sys->context, idv1, idv4, w4 * sys->vweights[idv1]); - } - } - } - else { - idv1 = BM_elem_index_get(vf[0]); - idv2 = BM_elem_index_get(vf[1]); - idv3 = BM_elem_index_get(vf[2]); - /* Is ring if number of faces == number of edges around vertice*/ - i = BM_elem_index_get(f); - if (!vert_is_boundary(vf[0]) && sys->zerola[idv1] == 0) { - EIG_linear_solver_matrix_add(sys->context, idv1, idv2, sys->fweights[i][2] * sys->vweights[idv1]); - EIG_linear_solver_matrix_add(sys->context, idv1, idv3, sys->fweights[i][1] * sys->vweights[idv1]); - } - if (!vert_is_boundary(vf[1]) && sys->zerola[idv2] == 0) { - EIG_linear_solver_matrix_add(sys->context, idv2, idv1, sys->fweights[i][2] * sys->vweights[idv2]); - EIG_linear_solver_matrix_add(sys->context, idv2, idv3, sys->fweights[i][0] * sys->vweights[idv2]); - } - if (!vert_is_boundary(vf[2]) && sys->zerola[idv3] == 0) { - EIG_linear_solver_matrix_add(sys->context, idv3, idv1, sys->fweights[i][1] * sys->vweights[idv3]); - EIG_linear_solver_matrix_add(sys->context, idv3, idv2, sys->fweights[i][0] * sys->vweights[idv3]); - } - } - } - } - BM_ITER_MESH (e, &eiter, sys->bm, BM_EDGES_OF_MESH) { - if (!BM_elem_flag_test(e, BM_ELEM_SELECT) && BM_edge_is_boundary(e)) { - v1 = e->v1->co; - v2 = e->v2->co; - idv1 = BM_elem_index_get(e->v1); - idv2 = BM_elem_index_get(e->v2); - if (sys->zerola[idv1] == 0 && sys->zerola[idv2] == 0) { - i = BM_elem_index_get(e); - EIG_linear_solver_matrix_add(sys->context, idv1, idv2, sys->eweights[i] * sys->vlengths[idv1]); - EIG_linear_solver_matrix_add(sys->context, idv2, idv1, sys->eweights[i] * sys->vlengths[idv2]); - } - } - } + float *v1, *v2, *v3, *v4; + float w2, w3, w4; + int i, j; + bool has_4_vert; + uint idv1, idv2, idv3, idv4, idv[4]; + + BMEdge *e; + BMFace *f; + BMIter eiter; + BMIter fiter; + BMIter vi; + BMVert *vn; + BMVert *vf[4]; + + BM_ITER_MESH (f, &fiter, sys->bm, BM_FACES_OF_MESH) { + if (BM_elem_flag_test(f, BM_ELEM_SELECT)) { + BM_ITER_ELEM_INDEX(vn, &vi, f, BM_VERTS_OF_FACE, i) + { + vf[i] = vn; + } + has_4_vert = (i == 4) ? 1 : 0; + if (has_4_vert) { + idv[0] = BM_elem_index_get(vf[0]); + idv[1] = BM_elem_index_get(vf[1]); + idv[2] = BM_elem_index_get(vf[2]); + idv[3] = BM_elem_index_get(vf[3]); + for (j = 0; j < 4; j++) { + idv1 = idv[j]; + idv2 = idv[(j + 1) % 4]; + idv3 = idv[(j + 2) % 4]; + idv4 = idv[(j + 3) % 4]; + + v1 = vf[j]->co; + v2 = vf[(j + 1) % 4]->co; + v3 = vf[(j + 2) % 4]->co; + v4 = vf[(j + 3) % 4]->co; + + w2 = cotangent_tri_weight_v3(v4, v1, v2) + cotangent_tri_weight_v3(v3, v1, v2); + w3 = cotangent_tri_weight_v3(v2, v3, v1) + cotangent_tri_weight_v3(v4, v1, v3); + w4 = cotangent_tri_weight_v3(v2, v4, v1) + cotangent_tri_weight_v3(v3, v4, v1); + + w2 = w2 / 4.0f; + w3 = w3 / 4.0f; + w4 = w4 / 4.0f; + + if (!vert_is_boundary(vf[j]) && sys->zerola[idv1] == 0) { + EIG_linear_solver_matrix_add(sys->context, idv1, idv2, w2 * sys->vweights[idv1]); + EIG_linear_solver_matrix_add(sys->context, idv1, idv3, w3 * sys->vweights[idv1]); + EIG_linear_solver_matrix_add(sys->context, idv1, idv4, w4 * sys->vweights[idv1]); + } + } + } + else { + idv1 = BM_elem_index_get(vf[0]); + idv2 = BM_elem_index_get(vf[1]); + idv3 = BM_elem_index_get(vf[2]); + /* Is ring if number of faces == number of edges around vertice*/ + i = BM_elem_index_get(f); + if (!vert_is_boundary(vf[0]) && sys->zerola[idv1] == 0) { + EIG_linear_solver_matrix_add( + sys->context, idv1, idv2, sys->fweights[i][2] * sys->vweights[idv1]); + EIG_linear_solver_matrix_add( + sys->context, idv1, idv3, sys->fweights[i][1] * sys->vweights[idv1]); + } + if (!vert_is_boundary(vf[1]) && sys->zerola[idv2] == 0) { + EIG_linear_solver_matrix_add( + sys->context, idv2, idv1, sys->fweights[i][2] * sys->vweights[idv2]); + EIG_linear_solver_matrix_add( + sys->context, idv2, idv3, sys->fweights[i][0] * sys->vweights[idv2]); + } + if (!vert_is_boundary(vf[2]) && sys->zerola[idv3] == 0) { + EIG_linear_solver_matrix_add( + sys->context, idv3, idv1, sys->fweights[i][1] * sys->vweights[idv3]); + EIG_linear_solver_matrix_add( + sys->context, idv3, idv2, sys->fweights[i][0] * sys->vweights[idv3]); + } + } + } + } + BM_ITER_MESH (e, &eiter, sys->bm, BM_EDGES_OF_MESH) { + if (!BM_elem_flag_test(e, BM_ELEM_SELECT) && BM_edge_is_boundary(e)) { + v1 = e->v1->co; + v2 = e->v2->co; + idv1 = BM_elem_index_get(e->v1); + idv2 = BM_elem_index_get(e->v2); + if (sys->zerola[idv1] == 0 && sys->zerola[idv2] == 0) { + i = BM_elem_index_get(e); + EIG_linear_solver_matrix_add( + sys->context, idv1, idv2, sys->eweights[i] * sys->vlengths[idv1]); + EIG_linear_solver_matrix_add( + sys->context, idv2, idv1, sys->eweights[i] * sys->vlengths[idv2]); + } + } + } } static bool vert_is_boundary(BMVert *v) { - BMEdge *ed; - BMFace *f; - BMIter ei; - BMIter fi; - BM_ITER_ELEM (ed, &ei, v, BM_EDGES_OF_VERT) { - if (BM_edge_is_boundary(ed)) { - return 1; - } - } - BM_ITER_ELEM (f, &fi, v, BM_FACES_OF_VERT) { - if (!BM_elem_flag_test(f, BM_ELEM_SELECT)) { - return 1; - } - } - return 0; + BMEdge *ed; + BMFace *f; + BMIter ei; + BMIter fi; + BM_ITER_ELEM (ed, &ei, v, BM_EDGES_OF_VERT) { + if (BM_edge_is_boundary(ed)) { + return 1; + } + } + BM_ITER_ELEM (f, &fi, v, BM_FACES_OF_VERT) { + if (!BM_elem_flag_test(f, BM_ELEM_SELECT)) { + return 1; + } + } + return 0; } -static void volume_preservation(BMOperator *op, float vini, float vend, int usex, int usey, int usez) +static void volume_preservation( + BMOperator *op, float vini, float vend, int usex, int usey, int usez) { - float beta; - BMOIter siter; - BMVert *v; - - if (vend != 0.0f) { - beta = pow(vini / vend, 1.0f / 3.0f); - BMO_ITER (v, &siter, op->slots_in, "verts", BM_VERT) { - if (usex) { - v->co[0] *= beta; - } - if (usey) { - v->co[1] *= beta; - } - if (usez) { - v->co[2] *= beta; - } - - } - } + float beta; + BMOIter siter; + BMVert *v; + + if (vend != 0.0f) { + beta = pow(vini / vend, 1.0f / 3.0f); + BMO_ITER (v, &siter, op->slots_in, "verts", BM_VERT) { + if (usex) { + v->co[0] *= beta; + } + if (usey) { + v->co[1] *= beta; + } + if (usez) { + v->co[2] *= beta; + } + } + } } -static void validate_solution(LaplacianSystem *sys, int usex, int usey, int usez, int preserve_volume) +static void validate_solution( + LaplacianSystem *sys, int usex, int usey, int usez, int preserve_volume) { - int m_vertex_id; - float leni, lene; - float vini, vend; - float *vi1, *vi2, ve1[3], ve2[3]; - uint idv1, idv2; - BMOIter siter; - BMVert *v; - BMEdge *e; - BMIter eiter; - - BM_ITER_MESH (e, &eiter, sys->bm, BM_EDGES_OF_MESH) { - idv1 = BM_elem_index_get(e->v1); - idv2 = BM_elem_index_get(e->v2); - vi1 = e->v1->co; - vi2 = e->v2->co; - ve1[0] = EIG_linear_solver_variable_get(sys->context, 0, idv1); - ve1[1] = EIG_linear_solver_variable_get(sys->context, 1, idv1); - ve1[2] = EIG_linear_solver_variable_get(sys->context, 2, idv1); - ve2[0] = EIG_linear_solver_variable_get(sys->context, 0, idv2); - ve2[1] = EIG_linear_solver_variable_get(sys->context, 1, idv2); - ve2[2] = EIG_linear_solver_variable_get(sys->context, 2, idv2); - leni = len_v3v3(vi1, vi2); - lene = len_v3v3(ve1, ve2); - if (lene > leni * SMOOTH_LAPLACIAN_MAX_EDGE_PERCENTAGE || lene < leni * SMOOTH_LAPLACIAN_MIN_EDGE_PERCENTAGE) { - sys->zerola[idv1] = 1; - sys->zerola[idv2] = 1; - } - } - - if (preserve_volume) { - vini = BM_mesh_calc_volume(sys->bm, false); - } - BMO_ITER (v, &siter, sys->op->slots_in, "verts", BM_VERT) { - m_vertex_id = BM_elem_index_get(v); - if (sys->zerola[m_vertex_id] == 0) { - if (usex) { - v->co[0] = EIG_linear_solver_variable_get(sys->context, 0, m_vertex_id); - } - if (usey) { - v->co[1] = EIG_linear_solver_variable_get(sys->context, 1, m_vertex_id); - } - if (usez) { - v->co[2] = EIG_linear_solver_variable_get(sys->context, 2, m_vertex_id); - } - } - } - if (preserve_volume) { - vend = BM_mesh_calc_volume(sys->bm, false); - volume_preservation(sys->op, vini, vend, usex, usey, usez); - } - + int m_vertex_id; + float leni, lene; + float vini, vend; + float *vi1, *vi2, ve1[3], ve2[3]; + uint idv1, idv2; + BMOIter siter; + BMVert *v; + BMEdge *e; + BMIter eiter; + + BM_ITER_MESH (e, &eiter, sys->bm, BM_EDGES_OF_MESH) { + idv1 = BM_elem_index_get(e->v1); + idv2 = BM_elem_index_get(e->v2); + vi1 = e->v1->co; + vi2 = e->v2->co; + ve1[0] = EIG_linear_solver_variable_get(sys->context, 0, idv1); + ve1[1] = EIG_linear_solver_variable_get(sys->context, 1, idv1); + ve1[2] = EIG_linear_solver_variable_get(sys->context, 2, idv1); + ve2[0] = EIG_linear_solver_variable_get(sys->context, 0, idv2); + ve2[1] = EIG_linear_solver_variable_get(sys->context, 1, idv2); + ve2[2] = EIG_linear_solver_variable_get(sys->context, 2, idv2); + leni = len_v3v3(vi1, vi2); + lene = len_v3v3(ve1, ve2); + if (lene > leni * SMOOTH_LAPLACIAN_MAX_EDGE_PERCENTAGE || + lene < leni * SMOOTH_LAPLACIAN_MIN_EDGE_PERCENTAGE) { + sys->zerola[idv1] = 1; + sys->zerola[idv2] = 1; + } + } + + if (preserve_volume) { + vini = BM_mesh_calc_volume(sys->bm, false); + } + BMO_ITER (v, &siter, sys->op->slots_in, "verts", BM_VERT) { + m_vertex_id = BM_elem_index_get(v); + if (sys->zerola[m_vertex_id] == 0) { + if (usex) { + v->co[0] = EIG_linear_solver_variable_get(sys->context, 0, m_vertex_id); + } + if (usey) { + v->co[1] = EIG_linear_solver_variable_get(sys->context, 1, m_vertex_id); + } + if (usez) { + v->co[2] = EIG_linear_solver_variable_get(sys->context, 2, m_vertex_id); + } + } + } + if (preserve_volume) { + vend = BM_mesh_calc_volume(sys->bm, false); + volume_preservation(sys->op, vini, vend, usex, usey, usez); + } } void bmo_smooth_laplacian_vert_exec(BMesh *bm, BMOperator *op) { - int i; - int m_vertex_id; - bool usex, usey, usez, preserve_volume; - float lambda_factor, lambda_border; - float w; - BMOIter siter; - BMVert *v; - LaplacianSystem *sys; - - if (bm->totface == 0) { - return; - } - sys = init_laplacian_system(bm->totedge, bm->totface, bm->totvert); - if (!sys) { - return; - } - sys->bm = bm; - sys->op = op; - - memset_laplacian_system(sys, 0); - - BM_mesh_elem_index_ensure(bm, BM_VERT); - lambda_factor = BMO_slot_float_get(op->slots_in, "lambda_factor"); - lambda_border = BMO_slot_float_get(op->slots_in, "lambda_border"); - sys->min_area = 0.00001f; - usex = BMO_slot_bool_get(op->slots_in, "use_x"); - usey = BMO_slot_bool_get(op->slots_in, "use_y"); - usez = BMO_slot_bool_get(op->slots_in, "use_z"); - preserve_volume = BMO_slot_bool_get(op->slots_in, "preserve_volume"); - - - sys->context = EIG_linear_least_squares_solver_new(bm->totvert, bm->totvert, 3); - - for (i = 0; i < bm->totvert; i++) { - EIG_linear_solver_variable_lock(sys->context, i); - } - BMO_ITER (v, &siter, op->slots_in, "verts", BM_VERT) { - m_vertex_id = BM_elem_index_get(v); - EIG_linear_solver_variable_unlock(sys->context, m_vertex_id); - EIG_linear_solver_variable_set(sys->context, 0, m_vertex_id, v->co[0]); - EIG_linear_solver_variable_set(sys->context, 1, m_vertex_id, v->co[1]); - EIG_linear_solver_variable_set(sys->context, 2, m_vertex_id, v->co[2]); - } - - init_laplacian_matrix(sys); - BMO_ITER (v, &siter, op->slots_in, "verts", BM_VERT) { - m_vertex_id = BM_elem_index_get(v); - EIG_linear_solver_right_hand_side_add(sys->context, 0, m_vertex_id, v->co[0]); - EIG_linear_solver_right_hand_side_add(sys->context, 1, m_vertex_id, v->co[1]); - EIG_linear_solver_right_hand_side_add(sys->context, 2, m_vertex_id, v->co[2]); - i = m_vertex_id; - if (sys->zerola[i] == 0) { - w = sys->vweights[i] * sys->ring_areas[i]; - sys->vweights[i] = (w == 0.0f) ? 0.0f : -lambda_factor / (4.0f * w); - w = sys->vlengths[i]; - sys->vlengths[i] = (w == 0.0f) ? 0.0f : -lambda_border * 2.0f / w; - - if (!vert_is_boundary(v)) { - EIG_linear_solver_matrix_add(sys->context, i, i, 1.0f + lambda_factor / (4.0f * sys->ring_areas[i])); - } - else { - EIG_linear_solver_matrix_add(sys->context, i, i, 1.0f + lambda_border * 2.0f); - } - } - else { - EIG_linear_solver_matrix_add(sys->context, i, i, 1.0f); - } - } - fill_laplacian_matrix(sys); - - if (EIG_linear_solver_solve(sys->context) ) { - validate_solution(sys, usex, usey, usez, preserve_volume); - } - - delete_laplacian_system(sys); + int i; + int m_vertex_id; + bool usex, usey, usez, preserve_volume; + float lambda_factor, lambda_border; + float w; + BMOIter siter; + BMVert *v; + LaplacianSystem *sys; + + if (bm->totface == 0) { + return; + } + sys = init_laplacian_system(bm->totedge, bm->totface, bm->totvert); + if (!sys) { + return; + } + sys->bm = bm; + sys->op = op; + + memset_laplacian_system(sys, 0); + + BM_mesh_elem_index_ensure(bm, BM_VERT); + lambda_factor = BMO_slot_float_get(op->slots_in, "lambda_factor"); + lambda_border = BMO_slot_float_get(op->slots_in, "lambda_border"); + sys->min_area = 0.00001f; + usex = BMO_slot_bool_get(op->slots_in, "use_x"); + usey = BMO_slot_bool_get(op->slots_in, "use_y"); + usez = BMO_slot_bool_get(op->slots_in, "use_z"); + preserve_volume = BMO_slot_bool_get(op->slots_in, "preserve_volume"); + + sys->context = EIG_linear_least_squares_solver_new(bm->totvert, bm->totvert, 3); + + for (i = 0; i < bm->totvert; i++) { + EIG_linear_solver_variable_lock(sys->context, i); + } + BMO_ITER (v, &siter, op->slots_in, "verts", BM_VERT) { + m_vertex_id = BM_elem_index_get(v); + EIG_linear_solver_variable_unlock(sys->context, m_vertex_id); + EIG_linear_solver_variable_set(sys->context, 0, m_vertex_id, v->co[0]); + EIG_linear_solver_variable_set(sys->context, 1, m_vertex_id, v->co[1]); + EIG_linear_solver_variable_set(sys->context, 2, m_vertex_id, v->co[2]); + } + + init_laplacian_matrix(sys); + BMO_ITER (v, &siter, op->slots_in, "verts", BM_VERT) { + m_vertex_id = BM_elem_index_get(v); + EIG_linear_solver_right_hand_side_add(sys->context, 0, m_vertex_id, v->co[0]); + EIG_linear_solver_right_hand_side_add(sys->context, 1, m_vertex_id, v->co[1]); + EIG_linear_solver_right_hand_side_add(sys->context, 2, m_vertex_id, v->co[2]); + i = m_vertex_id; + if (sys->zerola[i] == 0) { + w = sys->vweights[i] * sys->ring_areas[i]; + sys->vweights[i] = (w == 0.0f) ? 0.0f : -lambda_factor / (4.0f * w); + w = sys->vlengths[i]; + sys->vlengths[i] = (w == 0.0f) ? 0.0f : -lambda_border * 2.0f / w; + + if (!vert_is_boundary(v)) { + EIG_linear_solver_matrix_add( + sys->context, i, i, 1.0f + lambda_factor / (4.0f * sys->ring_areas[i])); + } + else { + EIG_linear_solver_matrix_add(sys->context, i, i, 1.0f + lambda_border * 2.0f); + } + } + else { + EIG_linear_solver_matrix_add(sys->context, i, i, 1.0f); + } + } + fill_laplacian_matrix(sys); + + if (EIG_linear_solver_solve(sys->context)) { + validate_solution(sys, usex, usey, usez, preserve_volume); + } + + delete_laplacian_system(sys); } |