diff options
| -rw-r--r-- | source/blender/blenlib/BLI_math_geom.h | 8 | ||||
| -rw-r--r-- | source/blender/blenlib/intern/math_geom.c | 75 | ||||
| -rw-r--r-- | source/blender/python/mathutils/mathutils_geometry.c | 93 |
3 files changed, 114 insertions, 62 deletions
diff --git a/source/blender/blenlib/BLI_math_geom.h b/source/blender/blenlib/BLI_math_geom.h index a9a55b10a9e..c0a9ea91e75 100644 --- a/source/blender/blenlib/BLI_math_geom.h +++ b/source/blender/blenlib/BLI_math_geom.h @@ -358,6 +358,14 @@ bool isect_plane_plane_v3(const float plane_a[4], float r_isect_co[3], float r_isect_no[3]) ATTR_WARN_UNUSED_RESULT; +bool isect_planes_v3_fn( + const float planes[][4], + const int planes_len, + const float eps_coplanar, + const float eps_isect, + void (*callback_fn)(const float co[3], int i, int j, int k, void *user_data), + void *user_data); + /* line/ray triangle */ bool isect_line_segment_tri_v3(const float p1[3], const float p2[3], diff --git a/source/blender/blenlib/intern/math_geom.c b/source/blender/blenlib/intern/math_geom.c index 1d2480f4d62..2b0018e7662 100644 --- a/source/blender/blenlib/intern/math_geom.c +++ b/source/blender/blenlib/intern/math_geom.c @@ -2298,6 +2298,81 @@ bool isect_plane_plane_v3(const float plane_a[4], } /** + * Intersect all planes, calling `callback_fn` for each point that intersects + * 3 of the planes that isn't outside any of the other planes. + * + * This can be thought of as calculating a convex-hull from an array of planes. + * + * \param eps_coplanar: Epsilon for testing if two planes are aligned (co-planar). + * \param eps_isect: Epsilon for testing of a point is behind any of the planes. + * + * \warning As complexity is a little under `O(N^3)`, this is only suitable for small arrays. + * + * \note This function could be optimized by some spatial structure. + */ +bool isect_planes_v3_fn( + const float planes[][4], + const int planes_len, + const float eps_coplanar, + const float eps_isect, + void (*callback_fn)(const float co[3], int i, int j, int k, void *user_data), + void *user_data) +{ + bool found = false; + + float n1n2[3], n2n3[3], n3n1[3]; + + for (int i = 0; i < planes_len; i++) { + const float *n1 = planes[i]; + for (int j = i + 1; j < planes_len; j++) { + const float *n2 = planes[j]; + cross_v3_v3v3(n1n2, n1, n2); + if (len_squared_v3(n1n2) <= eps_coplanar) { + continue; + } + for (int k = j + 1; k < planes_len; k++) { + const float *n3 = planes[k]; + cross_v3_v3v3(n2n3, n2, n3); + if (len_squared_v3(n2n3) <= eps_coplanar) { + continue; + } + + cross_v3_v3v3(n3n1, n3, n1); + if (len_squared_v3(n3n1) <= eps_coplanar) { + continue; + } + const float quotient = -dot_v3v3(n1, n2n3); + if (fabsf(quotient) < eps_coplanar) { + continue; + } + const float co_test[3] = { + ((n2n3[0] * n1[3]) + (n3n1[0] * n2[3]) + (n1n2[0] * n3[3])) / quotient, + ((n2n3[1] * n1[3]) + (n3n1[1] * n2[3]) + (n1n2[1] * n3[3])) / quotient, + ((n2n3[2] * n1[3]) + (n3n1[2] * n2[3]) + (n1n2[2] * n3[3])) / quotient, + }; + int i_test; + for (i_test = 0; i_test < planes_len; i_test++) { + const float *np_test = planes[i_test]; + if (((dot_v3v3(np_test, co_test) + np_test[3]) > eps_isect)) { + /* For low epsilon values the point could intersect it's own plane. */ + if (!ELEM(i_test, i, j, k)) { + break; + } + } + } + + if (i_test == planes_len) { /* ok */ + callback_fn(co_test, i, j, k, user_data); + found = true; + } + } + } + } + + return found; +} + +/** * Intersect two triangles. * * \param r_i1, r_i2: Retrieve the overlapping edge between the 2 triangles. diff --git a/source/blender/python/mathutils/mathutils_geometry.c b/source/blender/python/mathutils/mathutils_geometry.c index ad5f1a486b4..77ced169dab 100644 --- a/source/blender/python/mathutils/mathutils_geometry.c +++ b/source/blender/python/mathutils/mathutils_geometry.c @@ -1062,6 +1062,20 @@ static PyObject *M_Geometry_barycentric_transform(PyObject *UNUSED(self), PyObje return Vector_CreatePyObject(pt_dst, 3, NULL); } +struct PointsInPlanes_UserData { + PyObject *py_verts; + char *planes_used; +}; + +static void points_in_planes_fn(const float co[3], int i, int j, int k, void *user_data_p) +{ + struct PointsInPlanes_UserData *user_data = user_data_p; + PyList_APPEND(user_data->py_verts, Vector_CreatePyObject(co, 3, NULL)); + user_data->planes_used[i] = true; + user_data->planes_used[j] = true; + user_data->planes_used[k] = true; +} + PyDoc_STRVAR(M_Geometry_points_in_planes_doc, ".. function:: points_in_planes(planes)\n" "\n" @@ -1073,7 +1087,6 @@ PyDoc_STRVAR(M_Geometry_points_in_planes_doc, " :return: two lists, once containing the vertices inside the planes, another " "containing the plane indices used\n" " :rtype: pair of lists\n"); -/* note: this function could be optimized by some spatial structure */ static PyObject *M_Geometry_points_in_planes(PyObject *UNUSED(self), PyObject *args) { PyObject *py_planes; @@ -1090,81 +1103,37 @@ static PyObject *M_Geometry_points_in_planes(PyObject *UNUSED(self), PyObject *a } /* note, this could be refactored into plain C easy - py bits are noted */ - const float eps = 0.0001f; - const uint len = (uint)planes_len; - uint i, j, k, l; - float n1n2[3], n2n3[3], n3n1[3]; - float potentialVertex[3]; - char *planes_used = PyMem_Malloc(sizeof(char) * len); + struct PointsInPlanes_UserData user_data = { + .py_verts = PyList_New(0), + .planes_used = PyMem_Malloc(sizeof(char) * planes_len), + }; /* python */ - PyObject *py_verts = PyList_New(0); PyObject *py_plane_index = PyList_New(0); - memset(planes_used, 0, sizeof(char) * len); + memset(user_data.planes_used, 0, sizeof(char) * planes_len); - for (i = 0; i < len; i++) { - const float *N1 = planes[i]; - for (j = i + 1; j < len; j++) { - const float *N2 = planes[j]; - cross_v3_v3v3(n1n2, N1, N2); - if (len_squared_v3(n1n2) > eps) { - for (k = j + 1; k < len; k++) { - const float *N3 = planes[k]; - cross_v3_v3v3(n2n3, N2, N3); - if (len_squared_v3(n2n3) > eps) { - cross_v3_v3v3(n3n1, N3, N1); - if (len_squared_v3(n3n1) > eps) { - const float quotient = dot_v3v3(N1, n2n3); - if (fabsf(quotient) > eps) { - /** - * <pre> - * potentialVertex = ( - * (n2n3 * N1[3] + n3n1 * N2[3] + n1n2 * N3[3]) * - * (-1.0 / quotient)); - * </pre> - */ - const float quotient_ninv = -1.0f / quotient; - potentialVertex[0] = ((n2n3[0] * N1[3]) + (n3n1[0] * N2[3]) + (n1n2[0] * N3[3])) * - quotient_ninv; - potentialVertex[1] = ((n2n3[1] * N1[3]) + (n3n1[1] * N2[3]) + (n1n2[1] * N3[3])) * - quotient_ninv; - potentialVertex[2] = ((n2n3[2] * N1[3]) + (n3n1[2] * N2[3]) + (n1n2[2] * N3[3])) * - quotient_ninv; - for (l = 0; l < len; l++) { - const float *NP = planes[l]; - if ((dot_v3v3(NP, potentialVertex) + NP[3]) > 0.000001f) { - break; - } - } - - if (l == len) { /* ok */ - /* python */ - PyList_APPEND(py_verts, Vector_CreatePyObject(potentialVertex, 3, NULL)); - planes_used[i] = planes_used[j] = planes_used[k] = true; - } - } - } - } - } - } - } - } + const float eps_coplanar = 1e-4f; + const float eps_isect = 1e-6f; + const bool has_isect = isect_planes_v3_fn( + planes, planes_len, eps_coplanar, eps_isect, points_in_planes_fn, &user_data); PyMem_Free(planes); - /* now make a list of used planes */ - for (i = 0; i < len; i++) { - if (planes_used[i]) { - PyList_APPEND(py_plane_index, PyLong_FromLong(i)); + /* Now make user_data list of used planes. */ + if (has_isect) { + for (int i = 0; i < planes_len; i++) { + if (user_data.planes_used[i]) { + PyList_APPEND(py_plane_index, PyLong_FromLong(i)); + } } } - PyMem_Free(planes_used); + PyMem_Free(user_data.planes_used); { PyObject *ret = PyTuple_New(2); - PyTuple_SET_ITEMS(ret, py_verts, py_plane_index); + PyTuple_SET_ITEMS(ret, user_data.py_verts, py_plane_index); return ret; } } |