Hello all,
I've been working on the Automatic polygonal mesh healing project this
summer.
As it was decided that the mesh healing module be made into a portable
module, I started my work with that.
Here is a list of my milestones this summer.
1. Created a portable module which will work as a storage as well as
computation structure.
2. Implemented and tested the zippering gaps algorithm to heal gaps in
meshes.
The involvement of a lot of geometry in my project introduced a lot of
special scenarios that forced me to change my implementations in a quite a
few places, which slowed down my work.
I'm attaching a patch to my code till now.
My development logs are at:
http://brlcad.org/wiki/User:Tandoorichick/GSoC2016/Logs#Development_Logs
Week wise progress up to week 7 at: tandoorichick.github.io
Right now, I'm adding support for skewed lines, after modifying the initial
method of checking if the 2-D projection of the lines intersect.
Regards,
Rakshika.
Index: include/analyze.h
===================================================================
--- include/analyze.h (revision 68040)
+++ include/analyze.h (working copy)
@@ -212,7 +212,10 @@
const char *pbrep, struct rt_gen_worker_vars *pbrep_rtvars,
const char *curr_comb, struct bu_ptbl *candidates, void *curr_union_data, int ncpus);
+ANALYZE_EXPORT void
+analyze_heal_bot(struct rt_bot_internal *bot);
+
__END_DECLS
#endif /* ANALYZE_H */
Index: include/ged/objects.h
===================================================================
--- include/ged/objects.h (revision 68040)
+++ include/ged/objects.h (working copy)
@@ -215,6 +215,11 @@
GED_EXPORT extern int ged_group(struct ged *gedp, int argc, const char *argv[]);
/**
+ * Heal command to heal the defects in bots
+ */
+GED_EXPORT extern int ged_heal(struct ged *gedp, int argc, const char *argv[]);
+
+/**
* Set the "hidden" flag for the specified objects so they do not
* appear in an "ls" command output
*/
Index: src/libanalyze/CMakeLists.txt
===================================================================
--- src/libanalyze/CMakeLists.txt (revision 68040)
+++ src/libanalyze/CMakeLists.txt (working copy)
@@ -20,6 +20,12 @@
util.cpp
volume.c
voxels.c
+ MeshHealing/MeshConversion_brlcad.cpp
+ MeshHealing/MeshConversion.cpp
+ MeshHealing/Zipper.cpp
+ MeshHealing/Geometry.cpp
+ MeshHealing/Stitch.cpp
+ heal_mesh.cpp
)
add_subdirectory(tests)
Index: src/libanalyze/MeshHealing/DCEL.h
===================================================================
--- src/libanalyze/MeshHealing/DCEL.h (revision 0)
+++ src/libanalyze/MeshHealing/DCEL.h (working copy)
@@ -0,0 +1,55 @@
+#ifndef SRC_LIBANALYZE_MESHHEALING_DCEL_H_
+#define SRC_LIBANALYZE_MESHHEALING_DCEL_H_
+
+#include <utility>
+
+struct DCEL_Edge;
+
+/* DCEL_Vertex record of a doubly-connected DCEL_Edge list includes:
+ 1. DCEL_Vertex ID
+ 2. Coordinates
+ 3. An arbitrary incident DCEL_Edge - one with the particular DCEL_Vertex as its origin
+*/
+struct DCEL_Vertex {
+ int vertex_id;
+ double coordinates[3];
+ DCEL_Edge* incident_edge;
+};
+
+/* DCEL_Face record of a doubly-connected DCEL_Edge list includes:
+ 1. DCEL_Face ID,
+ 2. ID of any DCEL_Edge that can be used to traverse the DCEL_Face in counter-clockwise order
+ */
+struct DCEL_Face {
+ int face_id;
+ DCEL_Edge* start_edge;
+};
+
+/* The half-DCEL_Edge record of a doubly-connected DCEL_Edge list includes:
+ 1. DCEL_Edge ID
+ 2. ID of the origin DCEL_Vertex
+ 3. ID of the twin DCEL_Edge
+ 4. ID of the incident DCEL_Face - DCEL_Face to its left (since we orient the faces in counter-clockwise order)
+ 5. ID of the previous DCEL_Edge in the incident DCEL_Face
+ 6. ID of the next DCEL_Edge in the incident DCEL_Face
+ */
+struct DCEL_Edge {
+ std::pair<int, int> edge_id;
+ DCEL_Vertex* origin;
+ DCEL_Edge* twin;
+ DCEL_Face* incident_face;
+ DCEL_Edge* next;
+ DCEL_Edge* previous;
+};
+
+#endif /* SRC_LIBANALYZE_MESHHEALING_DCEL_H_ */
+
+/*
+ * Local Variables:
+ * mode: C
+ * tab-width: 8
+ * indent-tabs-mode: t
+ * c-file-style: "stroustrup"
+ * End:
+ * ex: shiftwidth=4 tabstop=8
+ */
Index: src/libanalyze/MeshHealing/Geometry.cpp
===================================================================
--- src/libanalyze/MeshHealing/Geometry.cpp (revision 0)
+++ src/libanalyze/MeshHealing/Geometry.cpp (working copy)
@@ -0,0 +1,185 @@
+#include "Geometry.h"
+
+#include <cmath>
+
+/* Calculates the determinant of a 3x3 matrix */
+double
+calcDet(const double m[3][3])
+{
+
+ double d;
+
+ d = m[0][0] * (m[1][1]*m[2][2] - m[1][2] * m[2][1]);
+ d -= m[0][1] * (m[1][0]*m[2][2] - m[1][2] * m[2][0]);
+ d += m[0][2] * (m[1][0]*m[2][1] - m[1][1] * m[2][0]);
+
+ return d;
+}
+
+void
+findOrthogonalProjection(const double A[3], const double B[3], const double C[3], double D[3])
+{
+ /* Parameter in the line equation of AB */
+ double t;
+
+ double BA[3];
+ double CA[3];
+
+ for (int i = 0; i < 3; ++i) {
+
+ BA[i] = B[i] - A[i];
+ CA[i] = C[i] - A[i];
+ }
+
+ /* Finding parameter t for a point D on line AB such that AB and DC are perpendicular to each other using the formula:
+ * t = dot(B-A, C-A)/dot(B-A, B-A) */
+ double num = DOT_PROD(BA, CA);
+ double den = DOT_PROD(BA, BA);
+
+ t = num / den;
+
+ /* t=0 and t=1 correspond to the end points */
+ if (t < 0)
+ t = 0;
+ else if (t > 1)
+ t = 1;
+
+ for (int i = 0; i < 3; ++i) {
+ D[i] = A[i] + t * BA[i];
+ }
+}
+
+bool
+isOrthoProjPossible(const double A[3], const double B[3], const double C[3])
+{
+ /* Parameter in the line equation of AB */
+ double t;
+
+ double BA[3];
+ double CA[3];
+
+ for (int i = 0; i < 3; ++i) {
+
+ BA[i] = B[i] - A[i];
+ CA[i] = C[i] - A[i];
+ }
+
+ /* Finding parameter t for a point D on line AB such that AB and DC are perpendicular to each other using this formula:
+ * t = dot(B-A, C-A)/dot(B-A, B-A)
+ */
+ double num = DOT_PROD(BA, CA);
+ double den = DOT_PROD(BA, BA);
+
+ t = num / den;
+
+ if (t > 0 && t < 1)
+ return true;
+
+ return false;
+}
+
+double
+shortestDistBwPoints(const double A[3], const double B[3])
+{
+ double dist = 0;
+ for (int i = 0; i < 3; i++) {
+ dist += (A[i] - B[i]) * (A[i] - B[i]);
+ }
+ dist = std::sqrt(dist);
+
+ return dist;
+}
+
+double
+shortestDistToLine(const double line_point1[3], const double line_point2[3], const double point[3])
+{
+ double dist;
+ double ortho[3];
+ findOrthogonalProjection(line_point1, line_point2, point, ortho);
+
+ dist = shortestDistBwPoints(point, ortho);
+
+ return dist;
+}
+
+bool
+checkIfIntersectsInterior(const double line1_point1[3], const double line1_point2[3], const double line2_point1[3], const double line2_point2[3])
+{
+ int o1, o2, o3, o4;
+ o1 = orientation(line1_point1, line1_point2, line2_point1);
+ o2 = orientation(line1_point1, line1_point2, line2_point2);
+ o3 = orientation(line2_point1, line2_point2, line1_point1);
+ o4 = orientation(line2_point1, line2_point2, line1_point2);
+
+ if (o1 != o2 && o3 != o4 && o1 != COLL && o2 != COLL && o3 != COLL && o4 != COLL)
+ return true;
+
+ return false;
+}
+bool
+checkIfCollinear(const double line1_point1[3], const double line1_point2[3], const double line2_point1[3], const double line2_point2[3])
+{
+ int o1, o2, o3, o4;
+ o1 = orientation(line1_point1, line1_point2, line2_point1);
+ o2 = orientation(line1_point1, line1_point2, line2_point2);
+ o3 = orientation(line2_point1, line2_point2, line1_point1);
+ o4 = orientation(line2_point1, line2_point2, line1_point2);
+ if(o1 == COLL && o2 == COLL && o3 == COLL && o4 == COLL)
+ return true;
+ return false;
+}
+
+int
+orientation(const double P[3], const double Q[3], const double R[3])
+{
+ double m[3][3];
+ double det;
+ for (int i = 0; i < 3; i++) {
+ m[0][i] = P[i];
+ m[1][i] = Q[i];
+ m[2][i] = R[i];
+ }
+
+ if(NEAR_0(m[0][2]) && NEAR_0(m[1][2]) && NEAR_0(m[2][2])) {
+ m[0][2] = 1;
+ m[1][2] = 1;
+ m[2][2] = 1;
+ }
+ det = calcDet(m);
+
+ if(NEAR_0(det))
+ return COLL;
+ else if (det > 0)
+ return CCW;
+ else
+ return CW;
+}
+
+void
+twoDimCoords(const double coordinates[3], double twod_coords[3])
+{
+ twod_coords[0] = coordinates[0];
+ twod_coords[1] = coordinates[1];
+ twod_coords[2] = 0;
+}
+
+double
+perimeter(const double A[3], const double B[3], const double C[3])
+{
+ double dist;
+
+ dist = shortestDistBwPoints(A, B);
+ dist += shortestDistBwPoints(B, C);
+ dist += shortestDistBwPoints(C, A);
+
+ return dist;
+}
+/*
+ * Local Variables:
+ * mode: C
+ * tab-width: 8
+ * indent-tabs-mode: t
+ * c-file-style: "stroustrup"
+ * End:
+ * ex: shiftwidth=4 tabstop=8
+ */
Index: src/libanalyze/MeshHealing/Geometry.h
===================================================================
--- src/libanalyze/MeshHealing/Geometry.h (revision 0)
+++ src/libanalyze/MeshHealing/Geometry.h (working copy)
@@ -0,0 +1,60 @@
+#ifndef SRC_LIBANALYZE_MESHHEALING_GEOMETRY_H_
+#define SRC_LIBANALYZE_MESHHEALING_GEOMETRY_H_
+
+#include <cstdlib>
+
+#define COLL 0
+#define CW 1
+#define CCW 2
+#define DOT_PROD(A, B) \
+ A[0] * B[0] \
+ + A[1] * B[1] \
+ + A[2] * B[2]
+
+#define TOLERANCE (1.0e-20)
+#define NEAR_0(val) (((val) > -TOLERANCE) && ((val) < TOLERANCE))
+
+/* Finds the orthogonal projection of a point C on a line segment formed by points A and B */
+void findOrthogonalProjection(const double A[3], const double B[3], const double C[3], double D[3]);
+
+/* Returns true if orthogonal projection of C on the line segment AB is possible, false otherwise. */
+bool isOrthoProjPossible(const double A[3], const double B[3], const double C[3]);
+
+/* Returns the shortest distance between the two points */
+double shortestDistBwPoints(const double A[3], const double B[3]);
+
+/* Returns the shortest distance between a point and a line segment */
+double shortestDistToLine(const double line_point1[3], const double line_point2[3], const double point[3]);
+
+/* Checks if line segment AB and CD intersect */
+bool checkIfIntersectsInterior(const double line1_point1[3], const double line1_point2[3], const double line2_point1[3], const double line2_point2[3]);
+
+/* Checks if the two lines are made up of collinear points */
+bool checkIfCollinear(const double line1_point1[3], const double line1_point2[3], const double line2_point1[3], const double line2_point2[3]);
+
+/* To check the orientation of the triangle formed by the three points P, Q, and R
+ * 0 - Collinear
+ * 1 - Clockwise
+ * 2 - Counter-clockwise*/
+int orientation(const double P[3], const double Q[3], const double R[3]);
+
+/* Calculates the determinant of the 3x3 matrix m */
+double calcDet(const double m[3][3]);
+
+/* Returns the 2D coordinates of the 3D point */
+void twoDimCoords(const double coordinates[3], double twod_coords[3]);
+
+/* Returns the perimeter of the triangle formed by the vertices A, B, and C */
+double perimeter(const double A[3], const double B[3], const double C[3]);
+
+#endif /* SRC_LIBANALYZE_MESHHEALING_GEOMETRY_H_ */
+
+/*
+ * Local Variables:
+ * mode: C
+ * tab-width: 8
+ * indent-tabs-mode: t
+ * c-file-style: "stroustrup"
+ * End:
+ * ex: shiftwidth=4 tabstop=8
+ */
Index: src/libanalyze/MeshHealing/MeshConversion.cpp
===================================================================
--- src/libanalyze/MeshHealing/MeshConversion.cpp (revision 0)
+++ src/libanalyze/MeshHealing/MeshConversion.cpp (working copy)
@@ -0,0 +1,881 @@
+#include "MeshConversion.h"
+
+#include <cstddef>
+#include <iterator>
+#include <limits>
+
+#include "Geometry.h"
+
+/* Vertex record function */
+
+void
+PolygonalMesh::initIncidentEdge()
+{
+
+ unsigned int num_vertices = this->getNumVertices();
+
+ for (unsigned int i = 0; i < num_vertices; ++i) {
+ for (unsigned int j = 0; j < edgelist.size(); ++j) {
+
+ if (vertexlist[i].vertex_id == edgelist[j].origin->vertex_id) {
+ vertexlist[i].incident_edge = &(edgelist[j]);
+ break;
+ }
+ }
+ }
+}
+
+/* Face record function */
+
+void
+PolygonalMesh::initFaces()
+{
+ unsigned int num_faces = this->getNumFaces();
+ DCEL_Face face;
+ for (unsigned int i = 0; i <= num_faces; i++) {
+ face.face_id = i;
+ facelist.push_back(face);
+ }
+}
+
+/* Edge record functions */
+
+void
+PolygonalMesh::initTwinEdge()
+{
+
+ for (unsigned int i = 0; i < edgelist.size(); ++i) {
+
+ for (unsigned int j = 0; j < edgelist.size(); j++) {
+ if (edgelist[i].edge_id.first == edgelist[j].edge_id.second && \
+ edgelist[j].edge_id.first == edgelist[i].edge_id.second) {
+
+ edgelist[i].twin = &(edgelist[j]);
+ edgelist[j].twin = &(edgelist[i]);
+ break;
+ }
+ }
+ }
+}
+
+/* Getters and setters */
+
+DCEL_Vertex*
+PolygonalMesh::getVertex(int ID)
+{
+ if (ID < (signed int)vertexlist.size())
+ return &vertexlist[ID];
+ return NULL;
+}
+
+void
+PolygonalMesh::setVertexCoordsInRecord(int ID, double coordinates[3])
+{
+ if (ID < (signed int)vertexlist.size()) {
+ vertexlist[ID].coordinates[0] = coordinates[0];
+ vertexlist[ID].coordinates[1] = coordinates[1];
+ vertexlist[ID].coordinates[2] = coordinates[2];
+
+ setVertexCoords(ID);
+ }
+}
+
+void
+PolygonalMesh::setVertexRecord(DCEL_Vertex *vertex, DCEL_Vertex *vertex_to_be_replaced)
+{
+ for (unsigned int i = 0; i < edgelist.size(); ++i) {
+ /* Edges whose origin is the vertex to be replaced */
+ if (edgelist[i].origin->vertex_id == vertex_to_be_replaced->vertex_id) {
+ edgelist[i].origin = getVertex(vertex->vertex_id);
+ edgelist[i].edge_id.first = vertex->vertex_id;
+ }
+
+ /* Edges that end at the vertex to be replaced */
+ if (edgelist[i].edge_id.second == vertex_to_be_replaced->vertex_id) {
+ edgelist[i].edge_id.second = vertex->vertex_id;
+ }
+ }
+ setVertex(vertex->vertex_id, vertex_to_be_replaced->vertex_id);
+}
+
+void
+PolygonalMesh::deleteVertexRecord(int ID)
+{
+ if ((unsigned)ID >= vertexlist.size())
+ return;
+
+ /* Setting the references for the origin of edges whose origin comes after the vertex being deleted */
+ for (unsigned int i = 0; i < edgelist.size(); i++) {
+ if (edgelist[i].origin->vertex_id > ID) {
+ edgelist[i].origin = &vertexlist[edgelist[i].origin->vertex_id - 1];
+ }
+ }
+
+ /* Delete the vertex record with vertex id=ID. */
+
+ vertexlist.erase(vertexlist.begin() + ID);
+ setNumVertices();
+
+ for (unsigned int i = 0; i < edgelist.size(); i++) {
+ if (edgelist[i].edge_id.first > ID)
+ edgelist[i].edge_id.first--;
+
+ if (edgelist[i].edge_id.second > ID)
+ edgelist[i].edge_id.second--;
+ }
+
+ /* Reassigning vertex IDs to the vertices */
+ for (unsigned int i = ID; i < vertexlist.size(); i++) {
+ vertexlist[i].vertex_id = i;
+ }
+}
+
+void
+PolygonalMesh::insertVertexOnEdge(DCEL_Vertex *vertex, DCEL_Edge *edge)
+{
+ /* edge11 and edge12 make up the edge passed as argument. edge21 and edge 22 make up the twin_edge */
+ DCEL_Edge *edge11, *edge12, *edge21, *edge22;
+
+ DCEL_Edge *twin_edge;
+
+ edge11 = addEdge(std::make_pair(edge->edge_id.first, vertex->vertex_id), edge->origin, NULL, edge->incident_face, NULL, edge->previous);
+
+ edge12 = addEdge(std::make_pair(vertex->vertex_id, edge->edge_id.second), vertex, NULL, edge->incident_face, edge->next, edge11);
+
+ edge11->next = edge12;
+
+ vertex->incident_edge = edge12;
+
+ /* If the edge was the starting edge of its incident face, change it to either one of the two new edges */
+ if (edge->incident_face->face_id != UNBOUNDED_FACE) {
+ if (edge->incident_face->start_edge->edge_id.first == edge->edge_id.first)
+ edge->incident_face->start_edge = edge11;
+ }
+
+ twin_edge = edge->twin;
+
+ edge21 = addEdge(std::make_pair(twin_edge->edge_id.first, vertex->vertex_id), twin_edge->origin, NULL, twin_edge->incident_face, NULL, twin_edge->previous);
+
+ edge22 = addEdge(std::make_pair(vertex->vertex_id, twin_edge->edge_id.second), vertex, NULL, twin_edge->incident_face, twin_edge->next, edge21);
+
+ edge21->next = edge22;
+
+
+ /* If the edge was the starting edge of its incident face, change it to either one of the two new edges */
+ if (twin_edge->incident_face->face_id != UNBOUNDED_FACE) {
+ if (twin_edge->incident_face->start_edge->edge_id.first == twin_edge->edge_id.first)
+ twin_edge->incident_face->start_edge = edge21;
+ }
+
+ edge11->twin = edge22;
+ edge12->twin = edge21;
+ edge21->twin = edge12;
+ edge22->twin = edge11;
+
+ checkAndSetTwins(vertex);
+
+ /* Delete the argument edge and its twin */
+ int v1 = edge->edge_id.first;
+ int v2 = edge->edge_id.second;
+ int face = edge->incident_face->face_id;
+ int twin_face = edge->twin->incident_face->face_id;
+ deleteEdge(std::make_pair(v1, v2), face);
+ deleteEdge(std::make_pair(v2, v1), twin_face);
+}
+
+void
+PolygonalMesh::deleteEdge(std::pair<int, int> ID, int face_id)
+{
+ int index = getEdgeIndex(ID, face_id);
+
+ int e_id;
+
+ for (unsigned int i = 0; i < vertexlist.size(); i++) {
+ e_id = getEdgeIndex(vertexlist[i].incident_edge->edge_id, vertexlist[i].incident_edge->incident_face->face_id);
+ if (e_id > index) {
+ vertexlist[i].incident_edge = &edgelist[e_id - 1];
+ }
+
+ else if (e_id == index) {
+ for (unsigned int j = 0; j < edgelist.size(); j++) {
+ if (vertexlist[i].vertex_id == edgelist[j].edge_id.first && j != (unsigned)index) {
+ vertexlist[i]. incident_edge = &edgelist[j];
+ break;
+ }
+ }
+ }
+ }
+
+ for (unsigned int i = 1; i < facelist.size(); i++) {
+ e_id = getEdgeIndex(facelist[i].start_edge->edge_id, i);
+ if (e_id > index) {
+ facelist[i].start_edge = &edgelist[e_id - 1];
+ }
+ }
+
+ for (unsigned int i =0; i < edgelist.size(); i++) {
+ e_id = getEdgeIndex(edgelist[i].twin->edge_id, edgelist[i].twin->incident_face->face_id);
+ if (e_id > index) {
+ edgelist[i].twin = &edgelist[e_id - 1];
+ }
+
+ e_id = getEdgeIndex(edgelist[i].next->edge_id, edgelist[i].next->incident_face->face_id);
+ if(e_id > index) {
+ edgelist[i].next = &edgelist[e_id - 1];
+ }
+
+ e_id = getEdgeIndex(edgelist[i].previous->edge_id, edgelist[i].previous->incident_face->face_id);
+ if ( e_id > index) {
+ edgelist[i].previous = &edgelist[e_id - 1];
+ }
+ }
+
+
+ if(index != -1) {
+ edgelist.erase(edgelist.begin() + index);
+ is_edge_checked.erase(is_edge_checked.begin() + index);
+ }
+}
+
+void
+PolygonalMesh::splitFace(int face_id, DCEL_Vertex *vertex_on_edge, DCEL_Vertex *opp_vertex)
+{
+ DCEL_Edge *start_edge1 = NULL, *start_edge2 = NULL;
+ DCEL_Edge *new_edge1, *new_edge2;
+ /* Indices of the newly added faces */
+ int face1, face2;
+
+ /* Find the two starting edges of the resultant faces */
+ for (unsigned int i = 0; i < edgelist.size(); ++i) {
+
+ if (edgelist[i].incident_face->face_id == face_id) {
+
+ if (edgelist[i].edge_id.first == vertex_on_edge->vertex_id)
+ start_edge1 = &edgelist[i];
+ else if (edgelist[i].edge_id.second== vertex_on_edge->vertex_id)
+ start_edge2 = &edgelist[i];
+ }
+
+ if (start_edge1 != NULL && start_edge2 != NULL)
+ break;
+ }
+
+ start_edge1->next->previous = start_edge1;
+ start_edge2->previous->next = start_edge2;
+
+ /* Create two half edges from the vertex passed as argument and the vertex opposite to the edge it is present on */
+ new_edge1 = addEdge(std::make_pair(opp_vertex->vertex_id, vertex_on_edge->vertex_id), opp_vertex, NULL, NULL, start_edge1, start_edge1->next);
+ new_edge2 = addEdge(std::make_pair(vertex_on_edge->vertex_id, opp_vertex->vertex_id), vertex_on_edge, NULL, NULL, start_edge2->previous, start_edge2);
+
+ new_edge1->twin = new_edge2;
+ new_edge2->twin = new_edge1;
+
+ face1 = addFaceRecord(start_edge1);
+ face2 = addFaceRecord(start_edge2);
+
+ start_edge1->incident_face = &facelist[face1];
+ new_edge1->incident_face = &facelist[face1];
+ start_edge1->next->incident_face = &facelist[face1];
+
+ start_edge2->incident_face = &facelist[face2];
+ new_edge2->incident_face = &facelist[face2];
+ start_edge2->previous->incident_face = &facelist[face2];
+
+ deleteFaceRecord(face_id);
+}
+
+int
+PolygonalMesh::addFaceRecord(DCEL_Edge *start_edge)
+{
+ int last_face_id = facelist.back().face_id;
+ DCEL_Face new_face;
+ new_face.face_id = last_face_id + 1;
+ new_face.start_edge = start_edge;
+
+ facelist.push_back(new_face);
+
+ /*addFace();*/
+
+ return last_face_id + 1;
+}
+
+DCEL_Edge*
+PolygonalMesh::addEdge(std::pair<int, int> edge_id, DCEL_Vertex *origin, DCEL_Edge *twin, DCEL_Face *inc_face, DCEL_Edge *next, DCEL_Edge *prev)
+{
+ DCEL_Edge *new_edge = new DCEL_Edge;
+
+ new_edge->edge_id.first = edge_id.first;
+ new_edge->edge_id.second = edge_id.second;
+
+ new_edge->origin = origin;
+
+ new_edge->twin = twin;
+ if (twin != NULL)
+ twin->twin = new_edge;
+
+ new_edge->incident_face = inc_face;
+
+ new_edge->next = next;
+ if (next != NULL)
+ new_edge->next->previous = new_edge;
+
+ new_edge->previous = prev;
+ if (prev != NULL)
+ new_edge->previous->next = new_edge;
+
+ edgelist.push_back(*new_edge);
+ is_edge_checked.push_back(false);
+
+ delete new_edge;
+
+ if (twin != NULL)
+ twin->twin = &edgelist[edgelist.size() - 1];
+ if (next != NULL)
+ next->previous = &edgelist[edgelist.size() - 1];
+ if (prev != NULL)
+ prev->next = &edgelist[edgelist.size() - 1];
+
+ return &edgelist[edgelist.size() - 1];
+}
+
+void
+PolygonalMesh::deleteFaceRecord(int face_id)
+{
+ /*deleteFace(face_id - 1);*/
+
+ /* Setting incident face references */
+ for (unsigned int i = 0; i < edgelist.size(); i++) {
+ if (edgelist[i].incident_face->face_id > face_id) {
+ edgelist[i].incident_face = &facelist[edgelist[i].incident_face->face_id - 1];
+ }
+ }
+ facelist.erase(facelist.begin() + face_id);
+
+ /* Reassigning face IDs */
+ for (unsigned int i = face_id; i < facelist.size(); i++) {
+ facelist[i].face_id = i;
+ }
+}
+
+DCEL_Edge*
+PolygonalMesh::getEdge(std::pair<int, int> ID)
+{
+ for (unsigned int i = 0; i < edgelist.size(); i++) {
+ if (edgelist[i].edge_id.first == ID.first && edgelist[i].edge_id.second == ID.second)
+ return &edgelist[i];
+ }
+ return NULL;
+}
+
+int
+PolygonalMesh::addVertexRecord(double coordinates[3])
+{
+ int last_vertex_id = vertexlist[vertexlist.size() - 1].vertex_id;
+ DCEL_Vertex new_vertex;
+ new_vertex.vertex_id = last_vertex_id + 1;
+
+ new_vertex.coordinates[0] = coordinates[0];
+ new_vertex.coordinates[1] = coordinates[1];
+ new_vertex.coordinates[2] = coordinates[2];
+
+ vertexlist.push_back(new_vertex);
+
+ addVertex(new_vertex.vertex_id);
+
+ return last_vertex_id + 1;
+}
+
+DCEL_Edge*
+PolygonalMesh::findClosestEdge(DCEL_Vertex *vertex)
+{
+ double min_dist = std::numeric_limits<int>::max();
+ double dist;
+ double *ortho_proj = new double[3];
+ DCEL_Edge *edge = NULL;
+ DCEL_Vertex *closer_vertex;
+ bool is_eligible;
+ DCEL_Face *face;
+
+ for (unsigned int i = 0; i < edgelist.size(); ++i) {
+ is_eligible = true;
+ /* The edge should have the unbounded face incident on it */
+ if (edgelist[i].incident_face->face_id != UNBOUNDED_FACE)
+ continue;
+
+ /* The edge should not be incident on the vertex itself */
+ if (edgelist[i].edge_id.first == vertex->vertex_id || edgelist[i].edge_id.second == vertex->vertex_id)
+ continue;
+
+ /* Line joining the vertex to the edge should not cross any bounded face */
+
+ /* The feature edge is an edge if an orthogonal projection is possible */
+ if(isOrthoProjPossible(getVertex(edgelist[i].edge_id.first)->coordinates, \
+ getVertex(edgelist[i].edge_id.second)->coordinates, vertex->coordinates)){
+ /* If the face incident on this edge's twin contains the vertex, it will cross a bounded face */
+ face = edgelist[i].twin->incident_face;
+ if (isVertexInFace(vertex->vertex_id, face->face_id))
+ continue;
+
+ findOrthogonalProjection(getVertex(edgelist[i].edge_id.first)->coordinates, \
+ getVertex(edgelist[i].edge_id.second)->coordinates, vertex->coordinates, ortho_proj);
+
+ /* Check if the line from the vertex to the edge intersects the face incident on the edge's twin,
+ * or any of the faces around the vertex
+ */
+ if(!checkEligibleEdge(vertex, ortho_proj))
+ continue;
+ if (doesLineIntersectFace(edgelist[i].twin->incident_face, vertex, ortho_proj))
+ continue;
+ }
+
+ /* The feature is a vertex */
+ else {
+
+ closer_vertex = findCloserVertex(&edgelist[i], vertex);
+
+ /* If there exists an edge between the closer_vertex and vertex, make the current edge in-eligible */
+ for (unsigned int j = 0; j < edgelist.size(); j++) {
+ if (edgelist[j].edge_id.first == vertex->vertex_id && edgelist[j].edge_id.second == closer_vertex->vertex_id) {
+ is_eligible = false;
+ break;
+ }
+ }
+ if(!is_eligible)
+ continue;
+
+ /* Check if the line joining vertex to closer_vertex crosses any other bounded face. Precisely:
+ * 1. Any face around the vertex
+ * 2. Any face around the closer vertex
+ */
+ if (!checkEligibleEdge(vertex, closer_vertex->coordinates))
+ continue;
+ if(!checkEligibleEdge(closer_vertex, vertex->coordinates))
+ continue;
+
+ for (int j = 0; j < 3; j++) {
+ ortho_proj[j] = closer_vertex->coordinates[j];
+ }
+
+ }
+
+ /* The edge is now eligible to be the closest edge */
+ dist = shortestDistBwPoints(vertex->coordinates, ortho_proj);
+
+ if (dist < min_dist) {
+ min_dist = dist;
+ edge = &edgelist[i];
+ }
+ }
+ delete[] ortho_proj;
+ return edge;
+}
+
+DCEL_Vertex*
+PolygonalMesh::findCloserVertex(DCEL_Edge *edge, DCEL_Vertex *vertex)
+{
+ double dist1, dist2;
+ DCEL_Vertex *vertex1 = getVertex(edge->edge_id.first);
+ DCEL_Vertex *vertex2 = getVertex(edge->edge_id.second);
+
+ dist1 = shortestDistBwPoints(vertex->coordinates, vertex1->coordinates);
+ dist2 = shortestDistBwPoints(vertex->coordinates, vertex2->coordinates);
+
+ if(dist1 < dist2)
+ return vertex1;
+ return vertex2;
+}
+
+DCEL_Face*
+PolygonalMesh::getFace(int ID)
+{
+ return &facelist[ID];
+}
+
+bool
+PolygonalMesh::isVertexInFace(int vertex_id, int face_id)
+{
+ DCEL_Edge *edge = (getFace(face_id))->start_edge;
+ DCEL_Edge *trav_edge = edge;
+
+ if (face_id == UNBOUNDED_FACE)
+ return false;
+
+ while (true) {
+ if (trav_edge->edge_id.first == vertex_id)
+ return true;
+ trav_edge = trav_edge->next;
+ if (trav_edge == edge)
+ break;
+ }
+ return false;
+}
+
+bool
+PolygonalMesh::isFreeEdge(DCEL_Edge* edge)
+{
+ if (edge->incident_face->face_id == UNBOUNDED_FACE)
+ return true;
+
+ return false;
+}
+
+int
+PolygonalMesh::getEdgeIndex(std::pair<int, int> ID, int face_id)
+{
+ for (unsigned int i = 0; i < edgelist.size(); ++i) {
+ if (edgelist[i].edge_id.first == ID.first && edgelist[i].edge_id.second == ID.second && \
+ edgelist[i].incident_face->face_id == face_id) {
+ return i;
+ }
+ }
+ return -1;
+}
+
+int
+PolygonalMesh::getNumEdges()
+{
+ return edgelist.size();
+}
+
+DCEL_Edge*
+PolygonalMesh::findFreeEdgeChain()
+{
+ DCEL_Edge *edge, *trav_edge;
+
+ for (unsigned int i = 0; i < edgelist.size(); ++i) {
+
+ if (!isFreeEdge(&edgelist[i]) || is_edge_checked[i])
+ continue;
+
+ edge = &edgelist[i];
+ trav_edge = edge;
+
+ while (isFreeEdge(trav_edge)) {
+ is_edge_checked[getEdgeIndex(trav_edge->edge_id, trav_edge->incident_face->face_id)] = true;
+ trav_edge = trav_edge->next;
+ if (trav_edge == edge) {
+ return trav_edge;
+ }
+ }
+ }
+ return NULL;
+}
+
+bool
+PolygonalMesh::doesLineIntersectFaceWithVertex(DCEL_Face *face, DCEL_Vertex *vertex, double other_vertex[3])
+{
+ if(face->face_id == UNBOUNDED_FACE)
+ return false;
+
+ double *twod_vertex = new double[3];
+ double *twod_other_vertex = new double[3];
+ double *twod_v1 = new double[3];
+ double *twod_v2 = new double[3];
+
+ DCEL_Edge *edge_to_be_checked = face->start_edge;
+
+ /* Vertices of the end points of the opposite edge */
+ DCEL_Vertex *v1, *v2;
+
+ bool flag = false;
+
+ for (int i = 0; i < VERTICES_PER_FACE; i++) {
+ if (vertex->vertex_id != edge_to_be_checked->edge_id.first && vertex->vertex_id != edge_to_be_checked->edge_id.second) {
+ flag = true;
+ break;
+ }
+ edge_to_be_checked = edge_to_be_checked->next;
+ }
+
+ if(!flag) {
+ delete[] twod_vertex;
+ delete[] twod_other_vertex;
+ delete[] twod_v1;
+ delete[] twod_v2;
+
+ return false;
+ }
+
+ v1 = getVertex(edge_to_be_checked->edge_id.first);
+ v2 = getVertex(edge_to_be_checked->edge_id.second);
+
+ twoDimCoords(vertex->coordinates, twod_vertex);
+ twoDimCoords(other_vertex, twod_other_vertex);
+ twoDimCoords(v1->coordinates, twod_v1);
+ twoDimCoords(v2->coordinates, twod_v2);
+
+ if (checkIfIntersectsInterior(vertex->coordinates, other_vertex, v1->coordinates, v2->coordinates)) {
+ delete[] twod_vertex;
+ delete[] twod_other_vertex;
+ delete[] twod_v1;
+ delete[] twod_v2;
+
+ return true;
+ }
+
+ /* Check if it is collinear with the other two edges */
+ edge_to_be_checked = edge_to_be_checked->next;
+ v1 = getVertex(edge_to_be_checked->edge_id.first);
+ v2 = getVertex(edge_to_be_checked->edge_id.second);
+
+ twoDimCoords(v1->coordinates, twod_v1);
+ twoDimCoords(v2->coordinates, twod_v2);
+
+ if (checkIfCollinear(vertex->coordinates, other_vertex, v1->coordinates, v2->coordinates)) {
+ delete[] twod_vertex;
+ delete[] twod_other_vertex;
+ delete[] twod_v1;
+ delete[] twod_v2;
+
+ return true;
+ }
+
+ edge_to_be_checked = edge_to_be_checked->next;
+ v1 = getVertex(edge_to_be_checked->edge_id.first);
+ v2 = getVertex(edge_to_be_checked->edge_id.second);
+
+ twoDimCoords(v1->coordinates, twod_v1);
+ twoDimCoords(v2->coordinates, twod_v2);
+
+ bool result = checkIfCollinear(vertex->coordinates, other_vertex, v1->coordinates, v2->coordinates);
+
+ delete[] twod_vertex;
+ delete[] twod_other_vertex;
+ delete[] twod_v1;
+ delete[] twod_v2;
+
+ return result;
+}
+
+DCEL_Face*
+PolygonalMesh::getNextFace(DCEL_Face *face, DCEL_Vertex *vertex)
+{
+ if(face->face_id == UNBOUNDED_FACE) {
+ for (unsigned int i = 0; i < edgelist.size(); i++) {
+ if(edgelist[i].edge_id.first == vertex->vertex_id && edgelist[i].incident_face->face_id == UNBOUNDED_FACE)
+ return edgelist[i].twin->incident_face;
+ }
+ }
+ DCEL_Edge *edge = face->start_edge;
+ DCEL_Edge *trav_edge = edge;
+ do {
+ if (trav_edge->edge_id.first == vertex->vertex_id) {
+ return trav_edge->twin->incident_face;
+ }
+ trav_edge = trav_edge->next;
+ }
+ while(trav_edge != edge);
+ return NULL;
+}
+
+bool
+PolygonalMesh::checkEligibleEdge(DCEL_Vertex *vertex, double closer_vertex[3])
+{
+ DCEL_Face *face, *trav_face;
+ bool is_elgible = true;
+ face = vertex->incident_edge->incident_face;
+ trav_face = face;
+
+ DCEL_Edge *edge;
+ do {
+ if (doesLineIntersectFaceWithVertex(trav_face, vertex, closer_vertex)) {
+ is_elgible = false;
+ break;
+ }
+
+
+ /* Check if orientation of the face is changing - i.e is the face getting turned over */
+ if (trav_face->face_id != UNBOUNDED_FACE) {
+
+ edge = trav_face->start_edge;
+ while(edge->edge_id.first != vertex->vertex_id) {
+ edge = edge->next;
+ }
+
+ if (orientation(closer_vertex, edge->next->origin->coordinates, edge->next->next->origin->coordinates) == CW) {
+ is_elgible = false;
+ break;
+ }
+ }
+
+ trav_face = getNextFace(trav_face, vertex);
+ }
+ while(trav_face != face);
+
+ return is_elgible;
+}
+
+bool
+PolygonalMesh::isEdgeOnChain(DCEL_Edge *edge_to_be_checked, DCEL_Edge *chain_edge)
+{
+ if(edge_to_be_checked == NULL)
+ return false;
+ DCEL_Edge *trav_edge = chain_edge;
+ do {
+ if(edge_to_be_checked == trav_edge)
+ return true;
+ trav_edge = trav_edge->next;
+ }
+ while(trav_edge != chain_edge);
+ return false;
+}
+
+void
+PolygonalMesh::squeezeOutEdges(DCEL_Edge* edge, DCEL_Edge *next_edge)
+{
+ if (edge->incident_face->face_id != UNBOUNDED_FACE)
+ return;
+
+ edge->previous->next = next_edge->next;
+ next_edge->next->previous = edge->previous;
+
+ /* Delete the edge records */
+ int v1 = next_edge->edge_id.first, v2 = next_edge->edge_id.second, face = next_edge->incident_face->face_id;
+
+ deleteEdge(std::make_pair(edge->edge_id.first, edge->edge_id.second), edge->incident_face->face_id);
+ deleteEdge(std::make_pair(v1, v2), face);
+
+ /* If this edge is the incident edge of any vertex, change the reference to some other edge whose origin is this edge's origin */
+ for (unsigned int i = 0; i < vertexlist.size(); i++) {
+ if (vertexlist[i].vertex_id != vertexlist[i].incident_edge->edge_id.first) {
+ for (unsigned int j = 0; j < edgelist.size(); j++) {
+ if (edgelist[j].edge_id.first == vertexlist[i].vertex_id && edgelist[j].incident_face->face_id != UNBOUNDED_FACE) {
+ vertexlist[i].incident_edge = &edgelist[j];
+ break;
+ }
+ }
+ }
+ }
+}
+
+bool
+PolygonalMesh::doesLineIntersectFace(DCEL_Face* face, DCEL_Vertex* vertex, double ortho_proj[3])
+{
+ if(face->face_id == UNBOUNDED_FACE)
+ return false;
+
+ double *twod_vertex, *twod_ortho_proj, *twod_v1, *twod_v2;
+ twod_vertex = new double[3];
+ twod_ortho_proj = new double[3];
+ twod_v1 = new double[3];
+ twod_v2 = new double[3];
+
+ DCEL_Edge *edge_to_be_checked = face->start_edge;
+ DCEL_Edge *trav_edge = edge_to_be_checked;
+
+ /* Vertices of the end points of the opposite edge */
+ DCEL_Vertex *v1, *v2;
+
+
+ twoDimCoords(vertex->coordinates, twod_vertex);
+
+ twoDimCoords(ortho_proj, twod_ortho_proj);
+
+ do {
+ v1 = getVertex(trav_edge->edge_id.first);
+ v2 = getVertex(trav_edge->edge_id.second);
+
+ twoDimCoords(v1->coordinates, twod_v1);
+
+ twoDimCoords(v2->coordinates, twod_v2);
+
+ if (checkIfIntersectsInterior(vertex->coordinates, ortho_proj, v1->coordinates, v2->coordinates)) {
+ delete[] twod_vertex;
+ delete[] twod_ortho_proj;
+ delete[] twod_v1;
+ delete[] twod_v2;
+ return true;
+ }
+ trav_edge = trav_edge->next;
+ }
+ while(trav_edge != edge_to_be_checked);
+
+ delete[] twod_vertex;
+ delete[] twod_ortho_proj;
+ delete[] twod_v1;
+ delete[] twod_v2;
+ return false;
+}
+
+void
+PolygonalMesh::checkAndSetTwins(DCEL_Vertex* vertex)
+{
+ DCEL_Edge *edge = vertex->incident_edge;
+ DCEL_Edge *trav_edge = edge;
+ DCEL_Edge *other_edge;
+ bool flag1 = false, flag2 = false;
+ std::pair<int, int> edge1_id, edge2_id;
+
+ DCEL_Edge *edge1, *edge2;
+ do {
+ trav_edge = trav_edge->twin->next;
+ }
+ while (trav_edge->incident_face->face_id != UNBOUNDED_FACE && trav_edge != edge);
+
+ if(trav_edge == edge && trav_edge->incident_face->face_id != UNBOUNDED_FACE)
+ return;
+
+ edge1 = trav_edge;
+ edge1_id.first = edge1->edge_id.first;
+ edge1_id.second = edge1->edge_id.second;
+
+ edge2 = edge1->previous;
+ edge2_id.first = edge2->edge_id.first;
+ edge2_id.second = edge2->edge_id.second;
+
+
+ if (edge1->next->edge_id.second == edge1->edge_id.first)
+ flag1 = true;
+
+ if (edge2->previous->edge_id.first == edge2->edge_id.second)
+ flag2 = true;
+
+ /* Look for twins, for the twins of these edges */
+ for (unsigned int i = 0; i < edgelist.size(); i++) {
+ if (flag1) {
+ if (edgelist[i].edge_id.first == edge1->edge_id.first && edgelist[i].edge_id.second == edge1->edge_id.second &&\
+ edgelist[i].incident_face->face_id != UNBOUNDED_FACE) {
+
+ edge1->twin->twin = &edgelist[i];
+ edgelist[i].twin = edge1->twin;
+ other_edge = edge1->next;
+ squeezeOutEdges(edge1, other_edge);
+ edge2 = &edgelist[getEdgeIndex(edge2_id, UNBOUNDED_FACE)];
+ flag1 = false;
+ }
+ }
+
+ if (flag2) {
+ if (edgelist[i].edge_id.first == edge2->edge_id.first && edgelist[i].edge_id.second == edge2->edge_id.second &&\
+ edgelist[i].incident_face->face_id != UNBOUNDED_FACE) {
+
+ edge2->twin->twin = &edgelist[i];
+ edgelist[i].twin = edge2->twin;
+ other_edge = edge2->previous;
+ squeezeOutEdges(other_edge, edge2);
+ edge2 = &edgelist[getEdgeIndex(edge1_id, UNBOUNDED_FACE)];
+ flag2 = false;
+ }
+ }
+
+ if (!flag1 && !flag2)
+ break;
+ }
+}
+
+DCEL_Edge*
+PolygonalMesh::getEdgeWithIndex(int index)
+{
+ if ((unsigned)index < edgelist.size())
+ return &edgelist[index];
+ return NULL;
+}
+/*
+ * Local Variables:
+ * mode: C
+ * tab-width: 8
+ * indent-tabs-mode: t
+ * c-file-style: "stroustrup"
+ * End:
+ * ex: shiftwidth=4 tabstop=8
+ */
Index: src/libanalyze/MeshHealing/MeshConversion.h
===================================================================
--- src/libanalyze/MeshHealing/MeshConversion.h (revision 0)
+++ src/libanalyze/MeshHealing/MeshConversion.h (working copy)
@@ -0,0 +1,156 @@
+#ifndef SRC_LIBANALYZE_MESHHEALING_MESHCONVERSION_H_
+#define SRC_LIBANALYZE_MESHHEALING_MESHCONVERSION_H_
+
+#include <utility>
+#include <vector>
+
+#include "DCEL.h"
+
+#define VERTICES_PER_FACE 3
+#define UNBOUNDED_FACE 0
+
+
+class PolygonalMesh {
+
+protected:
+
+ /* DCEL records */
+ std::vector<DCEL_Vertex> vertexlist;
+ std::vector<DCEL_Face> facelist;
+ std::vector<DCEL_Edge> edgelist;
+
+ /* Virtual functions to update the native structures */
+
+ /* DCEL_Vertex record functions */
+
+ /* Returns number of vertices */
+ virtual int getNumVertices() = 0;
+ /* Sets the DCEL_Vertex IDs and coordinates for all the vertices */
+ virtual void initVertices() = 0;
+ /* For every DCEL_Vertex, sets an arbitrary DCEL_Vertex incident on it */
+ void initIncidentEdge();
+
+ /* DCEL_Face record functions
+ */
+
+ /* Returns number of faces */
+ virtual int getNumFaces() = 0;
+
+ /* Sets a DCEL_Face ID for each DCEL_Face */
+ void initFaces();
+
+ /* Sets an arbitrary DCEL_Edge that is a prt of the DCEL_Face, so that that it can be traversed in counter-clockwise order */
+ virtual void initStartEdge() = 0;
+
+ /* DCEL_Edge record functions */
+
+ /* Sets DCEL_Edge ID and the origin of the DCEL_Edge */
+ virtual void initEdges() = 0;
+
+ /* Sets the DCEL_Face to the half-DCEL_Edge's left */
+ virtual void initIncidentFace() = 0;
+
+ /* Sets the twin DCEL_Edge of the half-DCEL_Edge */
+ void initTwinEdge();
+
+ /* Sets the DCEL_Edge previous to this DCEL_Edge in the incident DCEL_Face */
+ virtual void initPrevEdge() = 0;
+
+ /* Sets the DCEL_Edge next to this DCEL_Edge in the incident DCEL_Face */
+ virtual void initNextEdge() = 0;
+
+ virtual void setNumVertices() = 0;
+ virtual void setNumFaces() = 0;
+ virtual void setVertexCoords(int ID) = 0;
+ virtual void setVertex(int v1, int v2) = 0;
+ virtual void deleteVertex(int ID) = 0;
+ virtual void deleteFace(int ID) = 0;
+ virtual void addFace() = 0;
+ virtual void addVertex(int ID) = 0;
+
+ /* Test function to create bot from scratch
+ virtual void createBot();*/
+
+
+public:
+ virtual void setVertices() = 0;
+ virtual void setFaces() = 0;
+
+ std::vector<bool> is_edge_checked;
+
+ /* Getters and setters */
+
+ /* Get vertex with vertex id=ID */
+ DCEL_Vertex* getVertex(int ID);
+ /* Get edge with edge id=ID */
+ DCEL_Edge* getEdge(std::pair<int, int> ID);
+ /* Returns index of the edge given it's ID pair */
+ int getEdgeIndex(std::pair<int, int> ID, int face_id);
+ /* Get face with face id=ID */
+ DCEL_Face* getFace(int ID);
+ /* Set the coordinates of the vertex with vertex id=ID */
+ void setVertexCoordsInRecord(int ID, double coordinates[3]);
+ /* Replace all references of a vertex to be replaced with the other vertex */
+ void setVertexRecord(DCEL_Vertex *vertex, DCEL_Vertex *vertex_to_be_replaced);
+ /* Deletes a vertex with vertex id=ID */
+ void deleteVertexRecord(int ID);
+ /* Insert a vertex on the half edge */
+ void insertVertexOnEdge(DCEL_Vertex *vertex, DCEL_Edge *edge);
+ /* Deletes edge with edge id=ID */
+ void deleteEdge(std::pair<int, int> ID, int face_id);
+ /* Split the face through the line joining the given vertex (on an edge) to the opposite vertex */
+ void splitFace(int face_id, DCEL_Vertex *vertex_on_edge, DCEL_Vertex *opp_vertex);
+ /* Adds new face and sets the starting edge and returns the ID */
+ int addFaceRecord(DCEL_Edge *start_edge);
+ /* Adds an edge and sets the DCEL edge record attributes */
+ DCEL_Edge* addEdge(std::pair<int, int> edge_id, DCEL_Vertex *origin, DCEL_Edge *twin, DCEL_Face *inc_face, DCEL_Edge *next, DCEL_Edge *prev);
+ /* Deletes face with given face_id */
+ void deleteFaceRecord(int face_id);
+ /* Distance to an edge is the perpendicular distance if an orthogonal projection is possible,
+ * else it is the distance to the closer end point.
+ * Set the feature edge or the vertex accordingly. Set the other one to NULL.
+ * The edge if set, it is the half edge with the unbounded face.
+ * Line joining the vertex to the edge should not intersect any existing edge/vertex */
+ DCEL_Edge* findClosestEdge(DCEL_Vertex *vertex);
+ /* Returns the end point of the edge closer to the vertex */
+ DCEL_Vertex* findCloserVertex(DCEL_Edge *edge, DCEL_Vertex *vertex);
+ /*Adds new vertex with the given coordinates and returns ID*/
+ int addVertexRecord(double coordinates[3]);
+ /* Checks if the vertex with the given ID is present in the face with the given ID) */
+ bool isVertexInFace(int vertex_id, int face_id);
+ /* Checks if an edge is a free edge */
+ bool isFreeEdge(DCEL_Edge *edge);
+ /* Returns the start and end edges of a free-edge chain */
+ DCEL_Edge* findFreeEdgeChain();
+ /* Returns number of half edges in the DCEL */
+ int getNumEdges();
+ /* Checks whether the edge opposite to vertex in face intersects the line from vertex to other_vertx */
+ bool doesLineIntersectFaceWithVertex(DCEL_Face *face, DCEL_Vertex *vertex, double other_vertex[3]);
+ /* Checks whether line from vertex to the ortho_proj intersects the face */
+ bool doesLineIntersectFace(DCEL_Face *face, DCEL_Vertex *vertex, double ortho_proj[3]);
+ /* Returns the face adjacent to 'face' around the vertex */
+ DCEL_Face* getNextFace(DCEL_Face *face, DCEL_Vertex *vertex);
+ /* Checks whether the line from vertex to closer_vertex crosses any faces around vertex */
+ bool checkEligibleEdge(DCEL_Vertex *vertex, double closer_vertex[3]);
+ /* Sets the boundary edge of the mesh, and hence the boundary of the mesh can be traversed through this one edge */
+ void squeezeOutEdges(DCEL_Edge *edge, DCEL_Edge *next_edge);
+ /* Sets twin reference after vertex/edge contraction */
+ void checkAndSetTwins(DCEL_Vertex *vertex);
+ /* Checking if one edge lies on the chain of another edge */
+ bool isEdgeOnChain(DCEL_Edge *edge_to_be_checked, DCEL_Edge *chain_edge);
+ /* Returns edge in the given index in the edge list */
+ DCEL_Edge* getEdgeWithIndex(int index);
+
+};
+
+#endif /* SRC_LIBANALYZE_MESHHEALING_MESHCONVERSION_H_ */
+
+/*
+ * Local Variables:
+ * mode: C
+ * tab-width: 8
+ * indent-tabs-mode: t
+ * c-file-style: "stroustrup"
+ * End:
+ * ex: shiftwidth=4 tabstop=8
+ */
Index: src/libanalyze/MeshHealing/MeshConversion_brlcad.cpp
===================================================================
--- src/libanalyze/MeshHealing/MeshConversion_brlcad.cpp (revision 0)
+++ src/libanalyze/MeshHealing/MeshConversion_brlcad.cpp (working copy)
@@ -0,0 +1,436 @@
+#include "MeshConversion_brlcad.h"
+
+#include <bu/defines.h>
+#include <bu/malloc.h>
+#include <malloc.h>
+#include <rt/geom.h>
+#include <rt/primitives/bot.h>
+#include <stddef.h>
+#include <utility>
+#include <vector>
+
+#include "DCEL.h"
+#include "Geometry.h"
+
+
+BrlcadMesh::BrlcadMesh(rt_bot_internal *bot_mesh)
+{
+ bot = bot_mesh;
+ initVertices();
+ initFaces();
+ initEdges();
+ initIncidentFace();
+ initIncidentEdge();
+ initStartEdge();
+ initTwinEdge();
+ initPrevEdge();
+ initNextEdge();
+}
+
+/* DCEL_Vertex record functions */
+
+int
+BrlcadMesh::getNumVertices()
+{
+ return bot->num_vertices;
+}
+
+void
+BrlcadMesh::initVertices()
+{
+
+ unsigned int num_vertices = this->getNumVertices();
+ DCEL_Vertex vertex;
+
+ for (unsigned int i = 0; i < num_vertices; i++) {
+
+ vertex.vertex_id = i;
+
+ vertex.coordinates[0] = bot->vertices[VERTICES_PER_FACE * i];
+ vertex.coordinates[1] = bot->vertices[VERTICES_PER_FACE * i + 1];
+ vertex.coordinates[2] = bot->vertices[VERTICES_PER_FACE * i + 2];
+ vertexlist.push_back(vertex);
+ }
+}
+
+
+/* DCEL_Face record functions
+ * Note: The unbounded DCEL_Face will have face_id = 0. Its start_edge will be set to NULL.
+ */
+
+int
+BrlcadMesh::getNumFaces()
+{
+ return bot->num_faces;
+}
+
+
+void
+BrlcadMesh::initStartEdge()
+{
+
+ int vertex1, vertex2, vertex3;
+ int v1, v2;
+ unsigned int num_faces = this->getNumFaces();
+ facelist[UNBOUNDED_FACE].start_edge = NULL;
+
+ for (unsigned int i = 0; i < num_faces; ++i) {
+ vertex1 = bot->faces[VERTICES_PER_FACE * i];
+ vertex2 = bot->faces[VERTICES_PER_FACE * i + 1];
+ vertex3 = bot->faces[VERTICES_PER_FACE * i + 2];
+
+ if (bot->orientation == RT_BOT_CCW)
+ vertex2 = bot->faces[VERTICES_PER_FACE * i + 1];
+ else if (bot->orientation == RT_BOT_CW)
+ vertex2 = bot->faces[VERTICES_PER_FACE * i + 2];
+ else {
+ /* Determine the CCW orientation and set the vertex2 accordingly */
+
+ /* If the orientation is clockwise, switch the second and third vertices */
+ if (orientation(getVertex(vertex1)->coordinates, getVertex(vertex2)->coordinates, \
+ getVertex(vertex3)->coordinates) != CCW)
+ vertex2 = vertex3;
+
+ }
+
+ for (unsigned int j = 0; j < edgelist.size(); j++) {
+ v1 = edgelist[j].edge_id.first;
+ v2 = edgelist[j].edge_id.second;
+
+ if (vertex1 == v1 && vertex2 == v2) {
+ /* i+1 due to the presence of the unbounded face record */
+ facelist[i + 1].start_edge = &(edgelist[j]);
+ break;
+ }
+ }
+ }
+}
+
+/* DCEL_Edge record functions */
+
+void
+BrlcadMesh::initEdges()
+{
+ unsigned int num_faces = this->getNumFaces();
+ int v1, v2;
+
+ /* To check if the half-edges are already present */
+ bool flag;
+
+ for (unsigned int i = 0; i < num_faces; ++i) {
+
+ DCEL_Edge temp;
+ std::pair<int, int> id;
+
+ for (unsigned int j = 0; j < VERTICES_PER_FACE; ++j){
+ v1 = bot->faces[VERTICES_PER_FACE * i + j];
+ v2 = bot->faces[VERTICES_PER_FACE * i + (j + 1) % VERTICES_PER_FACE];
+
+ flag = false;
+
+ /* Checking if the DCEL_Edge from v1 to v2 is already present */
+ for (unsigned int k = 0; k < edgelist.size(); k++) {
+ if (v1 == edgelist[k].edge_id.first && v2 == edgelist[k].edge_id.second)
+ flag = true;
+
+ if (flag)
+ break;
+ }
+
+ if (!flag) {
+ /* Setting the DCEL_Edge record for the half record from v1 to v2 */
+ temp.origin = &(vertexlist[v1]);
+ temp.edge_id = std::make_pair(v1, v2);
+ edgelist.push_back(temp);
+
+ /* Setting the DCEL_Edge record for the half record from v2 to v1 */
+ temp.origin = &(vertexlist[v2]);
+ temp.edge_id = std::make_pair(v2, v1);
+ edgelist.push_back(temp);
+ }
+ }
+ }
+}
+
+void
+BrlcadMesh::initIncidentFace()
+{
+
+ int vertex1, vertex2;
+ unsigned int num_faces = this->getNumFaces();
+ int v1, v2, v3, temp;
+ bool face_set ;
+
+ int *a = new int[2];
+ int *b = new int[3];
+
+ for (unsigned int i = 0; i < edgelist.size(); ++i) {
+ face_set = false;
+
+ vertex1 = edgelist[i].edge_id.first;
+ vertex2 = edgelist[i].edge_id.second;
+
+ a[0] = vertex1;
+ a[1] = vertex2;
+
+ for (unsigned int j = 0; j < num_faces; j++) {
+ v1 = bot->faces[VERTICES_PER_FACE * j];
+ v2 = bot->faces[VERTICES_PER_FACE * j + 1];
+ v3 = bot->faces[VERTICES_PER_FACE * j + 2];
+
+ if (bot->orientation == RT_BOT_CCW) {
+ v2 = bot->faces[VERTICES_PER_FACE * i + 1];
+ v3 = bot->faces[VERTICES_PER_FACE * i + 2];
+ }
+ else if (bot->orientation == RT_BOT_CW) {
+ v2 = bot->faces[VERTICES_PER_FACE * i + 2];
+ v3 = bot->faces[VERTICES_PER_FACE * i + 1];
+ }
+ else {
+
+ /* If the orientation is clockwise, switch the second and third vertices */
+ if(orientation(getVertex(v1)->coordinates, getVertex(v2)->coordinates, getVertex(v3)->coordinates) != CCW) {
+ temp = v2;
+ v2 = v3;
+ v3 = temp;
+ }
+ }
+
+ b[0] = v1;
+ b[1] = v2;
+ b[2] = v3;
+
+ /* Check if a (vertex1 and vertex2) is present in the same order as b (the triangle made of v1, v2, v3) */
+ if (rt_bot_same_orientation(a, b)) {
+ edgelist[i].incident_face = &(facelist[j + 1]);
+ face_set = true;
+ break;
+ }
+ }
+ /* If no DCEL_Face has been set till here, set the unbounded DCEL_Face to be the incident DCEL_Face */
+ if (!face_set) {
+ edgelist[i].incident_face = &(facelist[UNBOUNDED_FACE]);
+ }
+ }
+ delete[] a;
+ delete[] b;
+}
+
+void
+BrlcadMesh::initPrevEdge()
+{
+
+ int inc_face, third_vertex;
+ int vertex1, vertex2;
+
+ for (unsigned int i = 0; i < edgelist.size(); i++) {
+
+ /* inc_face has the face_id of the incident face for the current edge */
+ inc_face = edgelist[i].incident_face->face_id;
+
+ /* The vertex_id of the vertex list comprising the current half edge */
+ vertex1 = edgelist[i].edge_id.first;
+ vertex2 = edgelist[i].edge_id.second;
+
+ /* If the face is the unbounded face, check in the edge list for the vertex with vertex1 as its end vertex,
+ * and has the unbounded face incident on it. */
+ if(inc_face == UNBOUNDED_FACE) {
+ for(unsigned int j = 0; j < edgelist.size(); j++) {
+ if(edgelist[j].edge_id.second == vertex1 && edgelist[j].incident_face->face_id == UNBOUNDED_FACE) {
+ edgelist[i].previous = &edgelist[j];
+ break;
+ }
+ }
+ } else {
+ for (unsigned int j = 0; j < VERTICES_PER_FACE; j++) {
+
+ /*Finding the third vertex, that is one that is not vertex1 or vertex2 */
+ if (bot->faces[VERTICES_PER_FACE * (inc_face - 1) + j] != vertex1 && \
+ bot->faces[VERTICES_PER_FACE * (inc_face - 1) + j] != vertex2) {
+ third_vertex = bot->faces[VERTICES_PER_FACE * (inc_face - 1) + j];
+ break;
+ }
+ }
+
+ /* Finding the edge_id of the edge from third_vertex to vertex1 */
+ edgelist[i].previous = getEdge(std::make_pair(third_vertex, vertex1));
+ }
+ }
+}
+void
+BrlcadMesh::initNextEdge()
+{
+
+ int inc_face, third_vertex;
+ int vertex1, vertex2;
+
+ for (unsigned int i = 0; i < edgelist.size(); i++) {
+
+ /* inc_face has the face_id of the incident face for the current edge */
+ inc_face = edgelist[i].incident_face->face_id;
+
+ /* The vertex_id of the vertices comprising the current half edge */
+ vertex1 = edgelist[i].edge_id.first;
+ vertex2 = edgelist[i].edge_id.second;
+
+ /* If the face is the unbounded face, check in the edge list for the vertex with vertex 1 as its end vertex,
+ * and has the unbounded face incident on it. */
+ if(inc_face == UNBOUNDED_FACE) {
+ for(unsigned int j = 0; j < edgelist.size(); j++) {
+ if(edgelist[j].edge_id.first == vertex2 && edgelist[j].incident_face->face_id == UNBOUNDED_FACE) {
+ edgelist[i].next = &edgelist[j];
+ break;
+ }
+ }
+ } else {
+ for (unsigned int j = 0; j < VERTICES_PER_FACE; j++) {
+
+ /*Finding the third vertex, that is one that is not vertex1 or vertex2 */
+ if (bot->faces[VERTICES_PER_FACE * (inc_face - 1) + j] != vertex1 && bot->faces[VERTICES_PER_FACE * (inc_face - 1) + j] != vertex2) {
+ third_vertex = bot->faces[VERTICES_PER_FACE * (inc_face - 1) + j];
+ break;
+ }
+ }
+
+ /* Finding the edge_id of the edge from vertex2 to third_vertex */
+ edgelist[i].next = getEdge(std::make_pair(vertex2, third_vertex));
+ }
+ }
+}
+
+void
+BrlcadMesh::setNumVertices()
+{
+ bot->num_vertices = vertexlist.size();
+}
+
+void
+BrlcadMesh::setNumFaces()
+{
+ bot->num_faces = facelist.size();
+}
+
+void
+BrlcadMesh::setVertices()
+{
+ bu_free(bot->vertices, "vertices");
+ bot->vertices = (fastf_t*)bu_malloc(sizeof(fastf_t) * vertexlist.size() * 3, "vertices reallocation");
+ bot->num_vertices = vertexlist.size();
+
+ for (unsigned int i = 0; i < vertexlist.size(); i++) {
+ for (int k = 0; k < 3; k++)
+ bot->vertices[3 * i + k] = vertexlist[i].coordinates[k];
+ }
+}
+
+void
+BrlcadMesh::setFaces()
+{
+ DCEL_Edge *edge;
+
+ bot->faces = (int*)realloc(bot->faces, 3* (facelist.size() - 1) * sizeof(int));
+ bot->num_faces = facelist.size() - 1;
+
+ /* Unbounded face is the first record, so skip */
+ for (unsigned int i = 1; i < facelist.size(); ++i) {
+
+ edge = facelist[i].start_edge;
+ for (int j = 0; j < VERTICES_PER_FACE; ++j) {
+
+ bot->faces[VERTICES_PER_FACE * (i - 1) + j] = edge->edge_id.first;
+ edge = edge->next;
+ }
+ }
+}
+
+void
+BrlcadMesh::setVertexCoords(int ID)
+{
+ for (int i = 0; i < 3; i++) {
+ bot->vertices[3 * ID + i] = vertexlist[ID].coordinates[i];
+ }
+}
+
+void BrlcadMesh::setVertex(int vertex, int vertex_to_be_replaced)
+{
+ for (int i = 0; i < getNumFaces(); ++i) {
+ for (int j = 0; j < VERTICES_PER_FACE; ++j) {
+ if (bot->faces[VERTICES_PER_FACE * i + j] == vertex_to_be_replaced)
+ bot->faces[VERTICES_PER_FACE * i + j] = vertex;
+ }
+ }
+}
+
+void
+BrlcadMesh::deleteVertex(int ID)
+{
+ for (int i = ID ; i < getNumVertices() - 1; i++) {
+ for (int j = 0; j < 3; j++) {
+ bot->vertices[3 * i + j] = bot->vertices[3 * (i + 1) + j];
+ }
+ }
+
+ for (unsigned int i = 0; i < bot->num_faces; i++) {
+ for (int j = 0; j < VERTICES_PER_FACE; j++) {
+ if (bot->faces[3* i + j] > ID) {
+ bot->faces[3* i + j]--;
+ }
+ }
+ }
+
+ bot->num_vertices -= 1;
+ /*bot->vertices = (fastf_t*)realloc(bot->vertices, 3 * getNumVertices() * sizeof(fastf_t));*/
+}
+
+void
+BrlcadMesh::deleteFace(int ID)
+{
+ for (int i = ID ; i < (getNumFaces() - 1); i++) {
+ for (int j = 0; j < VERTICES_PER_FACE; j++) {
+ bot->faces[3 * i + j] = bot->faces[3 * (i + 1) + j];
+ }
+ }
+ bot->num_faces -= 1;
+
+ /*bot->faces = (int*)realloc(bot->faces, 3 * getNumFaces() * sizeof(int));*/
+}
+
+void
+BrlcadMesh::addFace()
+{
+ int *new_faces = (int*)realloc(bot->faces, 3 * (getNumFaces() + 1) * sizeof(int));
+ bot->faces = new_faces;
+ int ID = facelist.size() - 1;
+ DCEL_Edge *trav_edge = facelist[ID].start_edge;
+
+ for (int i = 0; i < VERTICES_PER_FACE; i++) {
+ bot->faces[3 * (ID - 1) + i] = trav_edge->edge_id.first;
+ trav_edge = trav_edge->next;
+ }
+
+ bot->num_faces += 1;
+}
+
+void
+BrlcadMesh::addVertex(int ID)
+{
+ fastf_t *new_vertices = (fastf_t*)realloc(bot->vertices, 3 * (getNumVertices() + 1) * sizeof(fastf_t));
+
+ bot->vertices = new_vertices;
+
+ for (int i = 0; i < 3; i++) {
+ bot->vertices[3 * ID + i] = vertexlist[ID].coordinates[i];
+ }
+ bot->num_vertices += 1;
+}
+
+
+/*
+ * Local Variables:
+ * mode: C
+ * tab-width: 8
+ * indent-tabs-mode: t
+ * c-file-style: "stroustrup"
+ * End:
+ * ex: shiftwidth=4 tabstop=8
+ */
Index: src/libanalyze/MeshHealing/MeshConversion_brlcad.h
===================================================================
--- src/libanalyze/MeshHealing/MeshConversion_brlcad.h (revision 0)
+++ src/libanalyze/MeshHealing/MeshConversion_brlcad.h (working copy)
@@ -0,0 +1,49 @@
+#ifndef SRC_LIBANALYZE_MESHHEALING_MESHCONVERSION_BRLCAD_H_
+#define SRC_LIBANALYZE_MESHHEALING_MESHCONVERSION_BRLCAD_H_
+
+#include "MeshConversion.h"
+
+struct rt_bot_internal;
+
+class BrlcadMesh: public PolygonalMesh {
+
+public:
+ BrlcadMesh(rt_bot_internal *bot_mesh);
+
+private:
+ rt_bot_internal *bot;
+
+ int getNumVertices();
+ void initVertices();
+ int getNumFaces();
+ void initStartEdge();
+ void initEdges();
+ void initIncidentFace();
+ void initPrevEdge();
+ void initNextEdge();
+
+ void setNumVertices();
+ void setNumFaces();
+ void setVertices();
+ void setFaces();
+ void setVertexCoords(int ID);
+ void setVertex(int v1, int v2);
+ void deleteVertex(int ID);
+ void deleteFace(int ID);
+ void addFace();
+ void addVertex(int ID);
+
+ /* Test function
+ void createBot();*/
+};
+
+#endif /* SRC_LIBANALYZE_MESHHEALING_MESHCONVERSION_BRLCAD_H_ */
+/*
+ * Local Variables:
+ * mode: C
+ * tab-width: 8
+ * indent-tabs-mode: t
+ * c-file-style: "stroustrup"
+ * End:
+ * ex: shiftwidth=4 tabstop=8
+ */
Index: src/libanalyze/MeshHealing/Stitch.cpp
===================================================================
--- src/libanalyze/MeshHealing/Stitch.cpp (revision 0)
+++ src/libanalyze/MeshHealing/Stitch.cpp (working copy)
@@ -0,0 +1,138 @@
+#include "Stitch.h"
+
+#include <stddef.h>
+#include <algorithm>
+#include <utility>
+
+#include "Geometry.h"
+
+void
+stitchGaps(PolygonalMesh *polymesh, DCEL_Edge* start_A, DCEL_Edge* start_B, double tolerance)
+{
+ /* The edges that will be used to traverse both meshes, according to the advancing rule */
+ DCEL_Edge *edge_A = start_A, *edge_B = start_B;
+
+ DCEL_Vertex *first_A, *first_B;
+ DCEL_Vertex *next_A, *next_B;
+ int flag_A = 0, flag_B = 0;
+
+ double perimeter_A, perimeter_B;
+
+ do {
+
+ first_A = polymesh->getVertex(edge_A->edge_id.first);
+ first_B = polymesh->getVertex(edge_B->edge_id.second);
+
+ next_A = polymesh->getVertex(edge_A->edge_id.second);
+ next_B = polymesh->getVertex(edge_B->edge_id.first);
+
+ perimeter_A = perimeter(first_A->coordinates, next_A->coordinates, first_B->coordinates);
+ perimeter_B = perimeter(first_B->coordinates, next_B->coordinates, first_A->coordinates);
+
+ if (perimeter_A < perimeter_B) {
+ /* Add triangle comprising the vertices first_A, next_A, and first_B, and
+ * Advance the edge pointer on chain A
+ */
+ addFace(polymesh, first_A, next_A, first_B, edge_A->previous, edge_A->next, edge_B->previous, edge_B, tolerance);
+ edge_A = edge_A->next;
+ if (edge_A == start_A)
+ flag_A = 1;
+
+ } else {
+ /* Add triangle comprising the vertices first_B, next_B, and first_A, and
+ * Advance the edge pointer on chain B
+ */
+ addFace(polymesh, next_B, first_B, first_A, edge_B->previous, edge_A->next, edge_A->previous, edge_A, tolerance);
+ edge_B = edge_B->previous;
+ if (edge_B == start_B)
+ flag_B = 1;
+ }
+ }
+ while (!(edge_A == start_A && flag_A) && !(edge_B == start_B && flag_B));
+
+ if (edge_A == start_A) {
+
+ while (edge_B != start_B) {
+ /* Make triangles with the the ending vertex in the A chain
+ * And advance the edge pointer on the B chain till it reaches the end of the chain
+ */
+ /*addFace(polymesh, next_B, first_B, edge_A->origin, edge_A->previous, );*/
+ edge_B = edge_B->next;
+ }
+ }
+
+ else if (edge_B == start_B) {
+ while (edge_A != start_A) {
+ /* Make triangles with the the ending vertex in the A chain
+ * And advance the edge pointer on the B chain till it reaches the end of the chain
+ */
+
+ edge_A = edge_A->previous;
+ }
+ }
+}
+
+void
+addFace(PolygonalMesh *polymesh, DCEL_Vertex* first_vertex_A, DCEL_Vertex* second_vertex_A, DCEL_Vertex* first_vertex_B, \
+ DCEL_Edge *prev_A, DCEL_Edge *next_A, DCEL_Edge *prev_B, DCEL_Edge *next_B, double tolerance)
+{
+ if (std::max(shortestDistBwPoints(first_vertex_A->coordinates, first_vertex_B->coordinates), \
+ shortestDistBwPoints(second_vertex_A->coordinates, first_vertex_B->coordinates)) > tolerance)
+ return;
+
+ if (!polymesh->checkEligibleEdge(first_vertex_A, first_vertex_B->coordinates))
+ return;
+ if (!polymesh->checkEligibleEdge(first_vertex_B, first_vertex_A->coordinates))
+ return;
+ if (!polymesh->checkEligibleEdge(second_vertex_A, first_vertex_B->coordinates))
+ return;
+ if (!polymesh->checkEligibleEdge(first_vertex_B, second_vertex_A->coordinates))
+ return;
+
+ DCEL_Edge *existing_edge = polymesh->getEdge(std::make_pair(first_vertex_A->vertex_id, second_vertex_A->vertex_id));
+
+ DCEL_Edge *b_edge_btoa = polymesh->getEdge(std::make_pair(first_vertex_B->vertex_id, first_vertex_A->vertex_id));
+ DCEL_Edge *ub_edge_atob = NULL;
+
+ /* Adding half-edges from the first vertex of chain A to the vertex on chain B, if they are not already present */
+ if (b_edge_btoa == NULL) {
+ ub_edge_atob = polymesh->addEdge(std::make_pair(first_vertex_A->vertex_id, first_vertex_B->vertex_id), first_vertex_A, \
+ NULL, polymesh->getFace(UNBOUNDED_FACE), prev_A, next_B);
+ b_edge_btoa = polymesh->addEdge(std::make_pair(first_vertex_B->vertex_id, first_vertex_A->vertex_id), first_vertex_B, \
+ ub_edge_atob, NULL, existing_edge, NULL);
+ }
+
+ DCEL_Edge *b_edge_atob = polymesh->getEdge(std::make_pair(second_vertex_A->vertex_id, first_vertex_B->vertex_id));
+ DCEL_Edge *ub_edge_btoa = NULL;
+
+ if (b_edge_atob == NULL) {
+ ub_edge_btoa = polymesh->addEdge(std::make_pair(first_vertex_B->vertex_id, second_vertex_A->vertex_id), first_vertex_B, \
+ NULL, polymesh->getFace(UNBOUNDED_FACE), next_A, prev_B);
+ b_edge_atob = polymesh->addEdge(std::make_pair(second_vertex_A->vertex_id, first_vertex_B->vertex_id), second_vertex_A, \
+ ub_edge_btoa, NULL, NULL, existing_edge);
+ }
+ if (ub_edge_atob != NULL)
+ ub_edge_atob->twin = b_edge_btoa;
+
+ if (ub_edge_btoa != NULL)
+ ub_edge_btoa->twin = b_edge_atob;
+
+ b_edge_atob->next = b_edge_btoa;
+ b_edge_btoa->previous = b_edge_atob;
+
+ /* Add face record for the new face */
+ int new_face = polymesh->addFaceRecord(b_edge_atob);
+
+ b_edge_atob->incident_face = polymesh->getFace(new_face);
+ b_edge_btoa->incident_face = polymesh->getFace(new_face);
+ existing_edge->incident_face = polymesh->getFace(new_face);
+}
+/*
+ * Local Variables:
+ * mode: C
+ * tab-width: 8
+ * indent-tabs-mode: t
+ * c-file-style: "stroustrup"
+ * End:
+ * ex: shiftwidth=4 tabstop=8
+ */
Index: src/libanalyze/MeshHealing/Zipper.cpp
===================================================================
--- src/libanalyze/MeshHealing/Zipper.cpp (revision 0)
+++ src/libanalyze/MeshHealing/Zipper.cpp (working copy)
@@ -0,0 +1,302 @@
+#include "Zipper.h"
+
+#include <stddef.h>
+#include <climits>
+#include <utility>
+#include <vector>
+
+#include "Geometry.h"
+
+void
+calcFeaturePair(PolygonalMesh *polymesh, queue_element *node)
+{
+ DCEL_Edge *closest_edge;
+ DCEL_Vertex *vertex = node->vertex;
+ /* Vertices of the closest edge */
+ DCEL_Vertex *v1, *v2;
+
+ /* Get the closest edge */
+ closest_edge = polymesh->findClosestEdge(vertex);
+ if(closest_edge == NULL) {
+ node->feature_edge = NULL;
+ node->feature_vertex = NULL;
+ node->dist = INT_MAX;
+ return;
+ }
+ v1 = (polymesh->getVertex(closest_edge->edge_id.first));
+ v2 = (polymesh->getVertex(closest_edge->edge_id.second));
+
+ /* Check if orthogonal projection is possible, if yes set the feature_edge, else set the feature_vertex */
+ if (isOrthoProjPossible(v1->coordinates, v2->coordinates, vertex->coordinates)) {
+ node->feature_edge = closest_edge;
+ node->dist = shortestDistToLine(v1->coordinates, v2->coordinates, vertex->coordinates);
+ node->feature_vertex = NULL;
+ }
+
+ else {
+ /* Returns the closest end point */
+ node->feature_vertex = polymesh->findCloserVertex(closest_edge, vertex);
+ node->dist = shortestDistBwPoints(vertex->coordinates, node->feature_vertex->coordinates);
+ node->feature_edge = NULL;
+ }
+}
+
+std::priority_queue<queue_element, std::vector<queue_element>, compareDist>
+initPriorityQueue(PolygonalMesh *polymesh, DCEL_Edge *start, double tolerance)
+{
+ queue_element node;
+ DCEL_Vertex *vertex;
+ DCEL_Edge *edge;
+ int v_id;
+
+ /* Priority queue ordered based on the error measure */
+ std::priority_queue<queue_element, std::vector<queue_element>, compareDist> PQ;
+
+ if(start == NULL)
+ return PQ;
+
+ /* Check if the gap is in the shape of a square.
+ * If yes pick a vertex and set it's feature vertex to the opposite vertex on the square.
+ */
+ int no_of_sides = 0;
+ DCEL_Edge *trav_edge = start;
+
+ edge = start;
+ do
+ {
+ /* Consider the origin v_id of each edge */
+ v_id = edge->origin->vertex_id;
+
+ vertex = polymesh->getVertex(v_id);
+ node.vertex = vertex;
+ calcFeaturePair(polymesh, &node);
+
+ /* Push the queue element to the priority queue */
+ if (node.dist < tolerance)
+ PQ.push(node);
+
+ edge = edge->next;
+ }
+ while (edge != start);
+
+ if (!PQ.empty()) {
+ do {
+ no_of_sides++;
+ trav_edge = trav_edge->next;
+ }
+ while (trav_edge != start);
+
+ double *vec1 = new double[3];
+ double *vec2 = new double[3];
+ double *vec3 = new double[3];
+ if (no_of_sides == 4) {
+ /* Check if there are two consecutive right angles */
+ for (int i = 0; i < 3; i++) {
+ vec1[i] = start->origin->coordinates[i] - start->next->origin->coordinates[i];
+ vec2[i] = start->next->origin->coordinates[i] - start->next->next->origin->coordinates[i];
+ vec3[i] = start->next->next->origin->coordinates[i] - start->next->next->next->origin->coordinates[i];
+ }
+ if (orientation(vec1, vec2, vec3) == CW) {
+ delete[] vec1;
+ delete[] vec2;
+ delete[] vec3;
+ return PQ;
+ }
+
+ if (NEAR_0(DOT_PROD(vec1, vec2)) && NEAR_0(DOT_PROD(vec2, vec3))) {
+ node.vertex = start->origin;
+ node.feature_vertex = start->next->next->origin;
+ node.dist = shortestDistBwPoints(node.vertex->coordinates, node.feature_vertex->coordinates);
+
+ while (!PQ.empty()) {
+ PQ.pop();
+ }
+
+ PQ.push(node);
+
+ delete[] vec1;
+ delete[] vec2;
+ delete[] vec3;
+
+ return PQ;
+ }
+ }
+ }
+
+ return PQ;
+}
+
+void
+zipperGaps(PolygonalMesh *polymesh, double tolerance, std::priority_queue<queue_element, std::vector<queue_element>, compareDist> PQ)
+{
+ double *A, *B, *C;
+ bool vertex_contract = false;
+ int deleted_vertex_id = INT_MAX;
+
+ /*queue_element *curr_element = new queue_element;*/
+ queue_element curr_element;
+
+ /* Temporary priority queue */
+ std::priority_queue<queue_element, std::vector<queue_element>, compareDist> temp_PQ;
+
+ if (PQ.empty()) {
+ return;
+ }
+
+ while (!PQ.empty()) {
+
+ curr_element = PQ.top();
+ PQ.pop();
+ /* Vertex has already been healed */
+ if(NEAR_0(curr_element.dist))
+ continue;
+
+ /* Check if the error measure is within the tolerance for zippering.
+ * If not break out of the loop, since all the error measures after this one are not going to be lesser than this one
+ */
+ if (curr_element.dist > tolerance)
+ break;
+
+ /* If the feature is a vertex, move both these vertices to the midpoint of the line connecting the two vertices */
+ if (curr_element.feature_vertex != NULL) {
+ vertex_contract = true;
+
+ double *coords = new double[3];
+
+ /* Find the mid-point */
+ coords[0] = (curr_element.vertex->coordinates[0] + curr_element.feature_vertex->coordinates[0]) / 2;
+ coords[1] = (curr_element.vertex->coordinates[1] + curr_element.feature_vertex->coordinates[1]) / 2;
+ coords[2] = (curr_element.vertex->coordinates[2] + curr_element.feature_vertex->coordinates[2]) / 2;
+
+ /* Changing the coordinates of the vertex to the midpoint of the vertex and its feature pair */
+ polymesh->setVertexCoordsInRecord(curr_element.vertex->vertex_id, coords);
+
+ /* Replace all references of the feature vertex with that of the vertex on the free edge chain */
+ polymesh->setVertexRecord(curr_element.vertex, curr_element.feature_vertex);
+
+ /* Delete the feature vertex's record */
+ deleted_vertex_id = curr_element.feature_vertex->vertex_id;
+
+ polymesh->deleteVertexRecord(curr_element.feature_vertex->vertex_id);
+
+ if (curr_element.vertex->vertex_id > curr_element.feature_vertex->vertex_id) {
+ polymesh->checkAndSetTwins(polymesh->getVertex(curr_element.vertex->vertex_id - 1));
+ } else {
+ polymesh->checkAndSetTwins(polymesh->getVertex(curr_element.vertex->vertex_id));
+ }
+
+ delete[] coords;
+ }
+
+ /* Else if the feature is an edge, project the vertex orthogonally onto the edge and split the edge and its incident face at the point */
+ else if (curr_element.feature_edge != NULL) {
+
+ double *D = new double[3];
+
+ vertex_contract = false;
+
+ A = polymesh->getVertex(curr_element.feature_edge->edge_id.first)->coordinates;
+ B = polymesh->getVertex(curr_element.feature_edge->edge_id.second)->coordinates;
+ C = polymesh->getVertex(curr_element.vertex->vertex_id)->coordinates;
+
+ findOrthogonalProjection(A, B, C, D);
+
+ polymesh->setVertexCoordsInRecord(curr_element.vertex->vertex_id, D);
+
+ std::pair<int, int> fe_id;
+ fe_id.first = curr_element.feature_edge->edge_id.first;
+ fe_id.second = curr_element.feature_edge->edge_id.second;
+
+ int face_id = curr_element.feature_edge->twin->incident_face->face_id;
+
+ /* Insert a vertex with coordinates D on the feature edge. */
+ polymesh->insertVertexOnEdge(polymesh->getVertex(curr_element.vertex->vertex_id), curr_element.feature_edge);
+
+ int index = polymesh->getEdgeIndex(std::make_pair(fe_id.second, curr_element.vertex->vertex_id), face_id);
+
+ curr_element.feature_edge = polymesh->getEdgeWithIndex(index);
+
+ /* Split the face incident on the half edge of the feature edge with a bounded face incident */
+ polymesh->splitFace(curr_element.feature_edge->incident_face->face_id, \
+ polymesh->getVertex(curr_element.vertex->vertex_id), curr_element.feature_edge->previous->origin);
+
+ delete[] D;
+ }
+
+
+ if (PQ.empty()) {
+ return;
+ }
+
+ /* To account for:
+ * 1. Those vertices with invalid elements
+ * 2. Those vertices whose coordinates have been changed, and hence the error measure - recalculate for all measures
+ */
+ queue_element *element = new queue_element;
+ do {
+ *element = PQ.top();
+ PQ.pop();
+
+ /* To account for the vertex deletion in the vertex contraction process */
+ if (vertex_contract && element->vertex->vertex_id == deleted_vertex_id) {
+ continue;
+ }
+ if(vertex_contract && element->vertex->vertex_id > deleted_vertex_id) {
+ element->vertex = polymesh->getVertex(element->vertex->vertex_id - 1);
+ }
+
+ calcFeaturePair(polymesh, element);
+
+ if (element->dist < tolerance)
+ temp_PQ.push(*element);
+ }
+ while (!PQ.empty());
+
+ delete element;
+
+ while (!temp_PQ.empty()) {
+ PQ.push(temp_PQ.top());
+ temp_PQ.pop();
+ }
+ }
+ /*delete curr_element;*/
+}
+
+int
+isValidFeature(PolygonalMesh *polymesh, queue_element *node)
+{
+ if (node->feature_vertex != NULL) {
+ if (polymesh->getVertex(node->feature_vertex->vertex_id) == NULL) {
+ return 0;
+ }
+ return 1;
+ }
+
+ else if(node->feature_edge != NULL){
+ if (polymesh->getEdge(std::make_pair(node->feature_edge->edge_id.first, \
+ node->feature_edge->edge_id.second)) == NULL)
+ return 0;
+ return 1;
+ }
+ return 1;
+}
+
+bool
+checkIfValidPQ(std::priority_queue<queue_element, std::vector<queue_element>, compareDist> PQ)
+{
+ while(!PQ.empty()) {
+ if (PQ.top().dist != INT_MAX)
+ return true;
+ PQ.pop();
+ }
+ return false;
+}
+/*
+ * Local Variables:
+ * mode: C
+ * tab-width: 8
+ * indent-tabs-mode: t
+ * c-file-style: "stroustrup"
+ * End:
+ * ex: shiftwidth=4 tabstop=8
+ */
Index: src/libanalyze/MeshHealing/Zipper.h
===================================================================
--- src/libanalyze/MeshHealing/Zipper.h (revision 0)
+++ src/libanalyze/MeshHealing/Zipper.h (working copy)
@@ -0,0 +1,58 @@
+#ifndef SRC_LIBANALYZE_MESHHEALING_ZIPPER_H_
+#define SRC_LIBANALYZE_MESHHEALING_ZIPPER_H_
+
+#include <queue>
+#include <vector>
+
+#include "MeshConversion.h"
+
+struct queue_element {
+ DCEL_Vertex *vertex;
+
+ /* Either vertex feature or the edge feature will be set.
+ * Feature is the closest edge.
+ * If an orthogonal projection of the vertex on the edge is possible, the edge feature will be set.
+ * Else, the vertex feature will be set to the closer end point of the edge.
+ */
+ DCEL_Vertex *feature_vertex;
+ DCEL_Edge *feature_edge;
+
+ /* Distance between the vertex and the feature - the error measure*/
+ double dist;
+};
+
+/* Operator to compare the error measures of two feature pairs */
+struct compareDist
+{
+ bool operator()(queue_element& lhs, queue_element& rhs)
+ {
+ return lhs.dist > rhs.dist;
+ }
+};
+
+
+/* Calculate error measures for the vertices on the free edge chain and push the queue_element variable to the priority queue */
+std::priority_queue<queue_element, std::vector<queue_element>, compareDist> initPriorityQueue(PolygonalMesh *polymesh, DCEL_Edge *start, double tol);
+
+/* Take the priority queue and zippers all those vertices whose error measures fall within the tolerance for zippering */
+void zipperGaps(PolygonalMesh *polymesh, double tolerance, std::priority_queue<queue_element, std::vector<queue_element>, compareDist> PQ);
+
+/* Check if the queue element has a valid feature */
+int isValidFeature(PolygonalMesh *polymesh, queue_element *node);
+
+/* Calculates the feature pair of a queue element */
+void calcFeaturePair(PolygonalMesh *polymesh, queue_element *node);
+
+/* Checks whether at least one node in the Priority queue has a valid feature pair */
+bool checkIfValidPQ(std::priority_queue<queue_element, std::vector<queue_element>, compareDist> PQ);
+
+#endif /* SRC_LIBANALYZE_MESHHEALING_ZIPPER_H_ */
+/*
+ * Local Variables:
+ * mode: C
+ * tab-width: 8
+ * indent-tabs-mode: t
+ * c-file-style: "stroustrup"
+ * End:
+ * ex: shiftwidth=4 tabstop=8
+ */
Index: src/libanalyze/heal_mesh.cpp
===================================================================
--- src/libanalyze/heal_mesh.cpp (revision 0)
+++ src/libanalyze/heal_mesh.cpp (working copy)
@@ -0,0 +1,65 @@
+#include "analyze.h"
+#include <rt/geom.h>
+
+#include <stddef.h>
+#include <queue>
+#include <vector>
+
+#include "MeshHealing/MeshConversion_brlcad.h"
+#include "MeshHealing/Stitch.h"
+#include "MeshHealing/Zipper.h"
+
+void
+analyze_heal_bot(struct rt_bot_internal *bot)
+{
+ PolygonalMesh *polymesh = new BrlcadMesh(bot);
+ double zipper_tol = 800;
+ /*double stitch_tol = 10;*/
+
+ /* Zippering the gaps */
+
+ polymesh->is_edge_checked.assign(polymesh->getNumEdges(), false);
+
+ DCEL_Edge *start;
+ do {
+ std::priority_queue<queue_element, std::vector<queue_element>, compareDist> PQ;
+ start = polymesh->findFreeEdgeChain();
+
+ PQ = initPriorityQueue(polymesh, start, zipper_tol);
+ zipperGaps(polymesh, zipper_tol, PQ);
+ }
+ while (start != NULL);
+
+ /* For testing purposes - creating the bot from scratch */
+ polymesh->setVertices();
+ polymesh->setFaces();
+
+
+ /*Stitching bigger defects
+
+ polymesh->is_edge_checked.assign(polymesh->getNumEdges(), false);
+
+ DCEL_Edge *closest_edge = NULL;
+
+ while (start != NULL) {
+ start = polymesh->findFreeEdgeChain();
+ if (start == NULL)
+ break;
+
+ closest_edge = polymesh->findClosestEdge(start->origin);
+
+ if (!polymesh->isEdgeOnChain(closest_edge, start)) {
+ stitchGaps(polymesh, start, closest_edge, stitch_tol);
+ }
+ }*/
+}
+
+/*
+ * Local Variables:
+ * tab-width: 8
+ * mode: C
+ * indent-tabs-mode: t
+ * c-file-style: "stroustrup"
+ * End:
+ * ex: shiftwidth=4 tabstop=8
+ */
Index: src/libged/CMakeLists.txt
===================================================================
--- src/libged/CMakeLists.txt (revision 68040)
+++ src/libged/CMakeLists.txt (working copy)
@@ -118,6 +118,7 @@
grid2model_lu.c
grid2view_lu.c
group.c
+ heal.c
hide.c
how.c
human.c
Index: src/libged/bot.c
===================================================================
--- src/libged/bot.c (revision 68040)
+++ src/libged/bot.c (working copy)
@@ -24,6 +24,7 @@
*/
#include "common.h"
+#include "analyze.h"
#include <stdlib.h>
#include <ctype.h>
Index: src/libtclcad/tclcad_obj.c
===================================================================
--- src/libtclcad/tclcad_obj.c (revision 68040)
+++ src/libtclcad/tclcad_obj.c (working copy)
@@ -1111,6 +1111,7 @@
{"grid2model_lu", "x y", 4, to_view_func_less, ged_grid2model_lu},
{"grid2view_lu", "x y", 4, to_view_func_less, ged_grid2view_lu},
{"handle_expose", "vname count", TO_UNLIMITED, to_handle_expose, GED_FUNC_PTR_NULL},
+ {"heal", (char *)0, TO_UNLIMITED, to_pass_through_func, ged_heal},
{"hide", (char *)0, TO_UNLIMITED, to_pass_through_func, ged_hide},
{"hide_view", "vname [0|1]", 3, to_hide_view, GED_FUNC_PTR_NULL},
{"how", (char *)0, TO_UNLIMITED, to_pass_through_func, ged_how},
Index: src/mged/setup.c
===================================================================
--- src/mged/setup.c (revision 68040)
+++ src/mged/setup.c (working copy)
@@ -170,6 +170,7 @@
{"grid2model_lu", cmd_ged_plain_wrapper, ged_grid2model_lu},
{"grid2view_lu", cmd_ged_plain_wrapper, ged_grid2view_lu},
{"has_embedded_fb", cmd_has_embedded_fb, GED_FUNC_PTR_NULL},
+ {"heal", cmd_ged_plain_wrapper, ged_heal},
{"hide", cmd_ged_plain_wrapper, ged_hide},
{"hist", cmd_hist, GED_FUNC_PTR_NULL},
{"history", f_history, GED_FUNC_PTR_NULL},
------------------------------------------------------------------------------
What NetFlow Analyzer can do for you? Monitors network bandwidth and traffic
patterns at an interface-level. Reveals which users, apps, and protocols are
consuming the most bandwidth. Provides multi-vendor support for NetFlow,
J-Flow, sFlow and other flows. Make informed decisions using capacity
planning reports. http://sdm.link/zohodev2dev
_______________________________________________
BRL-CAD Developer mailing list
brlcad-devel@lists.sourceforge.net
https://lists.sourceforge.net/lists/listinfo/brlcad-devel
