Daniel,
PFA the patch with all the changes in code.
Regrds,
Rakshika.
Index: libanalyze/CMakeLists.txt
===================================================================
--- libanalyze/CMakeLists.txt (revision 68040)
+++ libanalyze/CMakeLists.txt (working copy)
@@ -20,6 +20,14 @@
util.cpp
volume.c
voxels.c
+ MeshHealing/MeshConversion_brlcad.cpp
+ MeshHealing/MeshConversion_brlcad.h
+ MeshHealing/MeshConversion.cpp
+ MeshHealing/MeshConversion.h
+ MeshHealing/Zipper.cpp
+ MeshHealing/Zipper.h
+ MeshHealing/Geometry.cpp
+ MeshHealing/Geometry.h
)
add_subdirectory(tests)
Index: libanalyze/MeshHealing/DCEL.h
===================================================================
--- libanalyze/MeshHealing/DCEL.h (revision 0)
+++ libanalyze/MeshHealing/DCEL.h (working copy)
@@ -0,0 +1,50 @@
+#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;
+};
+
+/*
+ * Local Variables:
+ * mode: C
+ * tab-width: 8
+ * indent-tabs-mode: t
+ * c-file-style: "stroustrup"
+ * End:
+ * ex: shiftwidth=4 tabstop=8
+ */
Index: libanalyze/MeshHealing/Geometry.cpp
===================================================================
--- libanalyze/MeshHealing/Geometry.cpp (revision 0)
+++ libanalyze/MeshHealing/Geometry.cpp (working copy)
@@ -0,0 +1,154 @@
+#include "Geometry.h"
+
+#include <cmath>
+
+/* Calculates the determinant of a 3x3 matrix */
+double
+calcDet(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;
+}
+
+double*
+findOrthogonalProjection(double *A, double *B, double *C)
+{
+ /* Parameter in the line equation of AB */
+ double t;
+
+ double BA[3], CA[3];
+ double *D = new double[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) */
+
+ t = DOT_PROD(BA, CA) / DOT_PROD(BA, BA);
+
+ /* 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];
+ }
+
+ return D;
+}
+
+bool
+isOrthoProjPossible(double *A, double *B, double *C)
+{
+ /* Parameter in the line equation of AB */
+ double t;
+
+ double BA[3], 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) */
+
+ t = DOT_PROD(BA, CA) / DOT_PROD(BA, BA);
+
+ if (t > 0 && t < 1)
+ return true;
+
+ return false;
+}
+
+double
+shortestDistBwPoints(double *A, double *B)
+{
+ double dist = 0;
+ for (int i = 0; i < 3; i++) {
+ dist += (A[i] - B[i]) * (A[i] - B[i]);
+ }
+ dist = sqrt(dist);
+
+ return dist;
+}
+
+double
+shortestDistToLine(double *A, double *B, double *C)
+{
+ double dist;
+
+ dist = shortestDistBwPoints(C, findOrthogonalProjection(A, B, C));
+
+ return dist;
+}
+
+bool
+checkIfIntersects(double* A, double* B, double* C, double* D)
+{
+ int o1, o2, o3, o4;
+ o1 = orientation(A, B, C);
+ o2 = orientation(A, B, D);
+ o3 = orientation(C, D, A);
+ o4 = orientation(C, D, B);
+
+ if (o1 != o2 && o3 != o4)
+ return true;
+
+ return false;
+}
+
+int
+orientation(double* P, double* Q, double* R)
+{
+ 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];
+ }
+
+ det = calcDet(m);
+
+ if (det > 0)
+ return CW;
+ else if (det < 0)
+ return CCW;
+ else
+ return COLL;
+}
+
+double*
+twoDimCoords(double *coordinates)
+{
+ double *coords = new double[3];
+
+ coords[0] = coordinates[0];
+ coords[1] = coordinates[1];
+ coords[2] = 0;
+
+ return coords;
+}
+
+/*
+ * Local Variables:
+ * mode: C
+ * tab-width: 8
+ * indent-tabs-mode: t
+ * c-file-style: "stroustrup"
+ * End:
+ * ex: shiftwidth=4 tabstop=8
+ */
Index: libanalyze/MeshHealing/Geometry.h
===================================================================
--- libanalyze/MeshHealing/Geometry.h (revision 0)
+++ libanalyze/MeshHealing/Geometry.h (working copy)
@@ -0,0 +1,50 @@
+#ifndef SRC_LIBANALYZE_MESHHEALING_GEOMETRY_H_
+#define SRC_LIBANALYZE_MESHHEALING_GEOMETRY_H_
+
+
+#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]
+
+/* Finds the orthogonal projection of a point C on a line segment formed by points A and B */
+double* findOrthogonalProjection(double *A, double *B, double *C);
+
+/* Returns true if orthogonal projection of C on the line segment AB is possible, false otherwise. */
+bool isOrthoProjPossible(double *A, double *B, double *C);
+
+/* Returns the shortest distance between the two points */
+double shortestDistBwPoints(double *A, double *B);
+
+/* Returns the shortest distance between a point A and a line segment BC */
+double shortestDistToLine(double *A, double *B, double *C);
+
+/* Checks if line segment AB and CD intersect */
+bool checkIfIntersects(double *A, double *B, double *C, double *D);
+
+/* 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(double *P, double *Q, double *R);
+
+/* Calculates the determinant of the 3x3 matrix m */
+double calcDet(double m[3][3]);
+
+/* Returns the 2D coordinates of the 3D point */
+double* twoDimCoords(double *coordinates);
+
+#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: libanalyze/MeshHealing/MeshConversion.cpp
===================================================================
--- libanalyze/MeshHealing/MeshConversion.cpp (revision 0)
+++ libanalyze/MeshHealing/MeshConversion.cpp (working copy)
@@ -0,0 +1,490 @@
+#include "MeshConversion.h"
+#include "Geometry.h"
+
+#include <cstddef>
+#include <iterator>
+#include <limits>
+
+/* 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();
+ for (unsigned int i = 0; i <= num_faces; i++) {
+ facelist[i].face_id = i;
+ }
+}
+
+/* Edge record functions */
+
+void
+PolygonalMesh::initTwinEdge()
+{
+
+ for (unsigned int i = 0; i < edgelist.size(); ++i) {
+
+ /* If the edge's twin has already been set, move on to the next edge */
+ if (edgelist[i].twin != NULL)
+ continue;
+
+ 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)
+{
+ 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 *v1, DCEL_Vertex *v2)
+{
+ for (unsigned int i = 0; i < edgelist.size(); ++i) {
+ /* Edges whose origin is v2 */
+ if (edgelist[i].origin->vertex_id == v2->vertex_id) {
+ edgelist[i].origin = v1;
+ edgelist[i].edge_id.first = v1->vertex_id;
+ }
+
+ /* Edges that end at v2 */
+ if (edgelist[i].edge_id.second == v2->vertex_id) {
+ edgelist[i].edge_id.second = v1->vertex_id;
+ }
+ }
+ setVertex(v1->vertex_id, v2->vertex_id);
+}
+
+void
+PolygonalMesh::deleteVertexRecord(int ID)
+{
+ int inc_face_id;
+
+ deleteVertex(ID);
+ /* Delete the vertex record with vertex id=ID. */
+ if (ID < (signed int)vertexlist.size()) {
+ vertexlist.erase(vertexlist.begin() + ID);
+ setNumVertices();
+ }
+
+ /* Deleting the edges that contain the vertex */
+ for (unsigned int i = 0; i < edgelist.size(); ++i) {
+
+ if (edgelist[i].edge_id.first == ID || edgelist[i].edge_id.second == ID) {
+ edgelist.erase(edgelist.begin() + i);
+
+ /* Delete the face that this edge is incident on as well */
+ inc_face_id = edgelist[i].incident_face->face_id;
+ deleteFace(inc_face_id);
+ facelist.erase(facelist.begin() + inc_face_id);
+
+ setNumFaces();
+
+ /* Set the incident face of the edge to the unbounded face */
+ edgelist[i].incident_face = &facelist[getNumFaces()];
+ }
+ }
+}
+
+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, *twin_edge;
+
+ vertex->incident_edge = &edge12;
+
+ edge11.origin = edge->origin;
+ edge11.next = &edge12;
+ edge11.previous = edge->previous;
+ edge11.incident_face = edge->incident_face;
+
+ edge12.origin = vertex;
+ edge12.next = edge->next;
+ edge12.previous = &edge11;
+ edge12.incident_face = edge->incident_face;
+
+ edge->previous->next = &edge11;
+ edge->next->previous = &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->start_edge->edge_id.first == edge->edge_id.first)
+ edge->incident_face->start_edge = &edge11;
+
+
+ twin_edge = edge->twin;
+
+ edge21.origin = twin_edge->origin;
+ edge21.next = &edge22;
+ edge21.previous = twin_edge->previous;
+ edge21.incident_face = twin_edge->incident_face;
+
+ edge22.origin = vertex;
+ edge22.next = twin_edge->next;
+ edge22.previous = &edge21;
+ edge22.incident_face = twin_edge->incident_face;
+
+ twin_edge->previous->next = &edge21;
+ twin_edge->next->previous = &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->start_edge->edge_id.first == twin_edge->edge_id.first)
+ twin_edge->incident_face->start_edge = &edge11;
+
+ edge11.twin = &edge22;
+ edge21.twin = &edge12;
+ edge12.twin = &edge21;
+ edge22.twin = &edge11;
+
+
+ /* Delete the argument edge and its twin */
+ deleteEdge(std::make_pair(edge->edge_id.first, edge->edge_id.second));
+ deleteEdge(std::make_pair(edge->edge_id.second, edge->edge_id.first));
+}
+
+void
+PolygonalMesh::deleteEdge(std::pair<int, int> ID)
+{
+ for (unsigned int i = 0; i < edgelist.size(); ++i) {
+
+ if (ID.first == edgelist[i].edge_id.first && ID.second == edgelist[i].edge_id.second) {
+ edgelist.erase(edgelist.begin() + i);
+ }
+ }
+}
+
+void
+PolygonalMesh::splitFace(int face_id, DCEL_Vertex *vertex, DCEL_Vertex *opp_vertex)
+{
+ DCEL_Edge *start_edge1, *start_edge2;
+ 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].edge_id.first == vertex->vertex_id) {
+
+ if (edgelist[i].incident_face->face_id == face_id)
+ start_edge1 = &edgelist[i];
+ else
+ start_edge2 = &edgelist[i];
+ }
+
+ if (start_edge1 != NULL && start_edge2 != NULL)
+ break;
+ }
+
+ face1 = addFaceRecord(start_edge1);
+ face2 = addFaceRecord(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(vertex->vertex_id, opp_vertex->vertex_id), vertex, NULL, &facelist[face1], start_edge1->previous, start_edge1);
+ new_edge2 = addEdge(std::make_pair(opp_vertex->vertex_id, vertex->vertex_id), opp_vertex, NULL, &facelist[face2], start_edge2, start_edge2->next);
+
+ new_edge1->twin = new_edge2;
+ new_edge2->twin = new_edge1;
+
+ 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(new_face.face_id);
+
+ return last_face_id + 1;
+}
+
+void
+PolygonalMesh::deleteDuplicateEdges()
+{
+ for (unsigned int i = 0; i < edgelist.size(); ++i) {
+
+ if (edgelist[i].incident_face->face_id == UNBOUNDED_FACE)
+ continue;
+
+ for (unsigned int j = 0; j < edgelist.size(); ++j) {
+
+ /* Looking for the duplicate half edge whose incident face will be the unbounded face */
+ if (edgelist[i].edge_id.first == edgelist[j].edge_id.first && edgelist[i].edge_id.second == edgelist[j].edge_id.second) {
+ edgelist.erase(edgelist.begin() + j);
+ }
+ }
+ }
+}
+
+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 = edge_id;
+ new_edge->origin = origin;
+ new_edge->twin = twin;
+ new_edge->incident_face = inc_face;
+ new_edge->next = next;
+ new_edge->previous = prev;
+ return new_edge;
+}
+
+void
+PolygonalMesh::deleteFaceRecord(int face_id)
+{
+ facelist.erase(facelist.begin() + face_id);
+ deleteFace(face_id);
+}
+
+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)
+{
+ 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;
+ DCEL_Edge *edge = NULL;
+ DCEL_Face *face;
+ DCEL_Edge *prev, *next;
+ int v1, v2;
+
+ for (unsigned int i = 0; i < edgelist.size(); ++i) {
+ /* 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 */
+
+ /* If the face incident on this edge or its twin contains the vertex, it will cross a bounded face */
+ face = edgelist[i].incident_face;
+ if (isVertexInFace(vertex->vertex_id, face->face_id))
+ continue;
+
+ face = edgelist[i].twin->incident_face;
+ if (isVertexInFace(vertex->vertex_id, face->face_id))
+ continue;
+
+ ortho_proj = findOrthogonalProjection(getVertex(edgelist[i].edge_id.first)->coordinates, getVertex(edgelist[i].edge_id.second)->coordinates, vertex->coordinates);
+
+ /* Get the other two edges in the triangle.
+ * Check if the two dimensional projections of the line from the vertex to the edge intersects either of them
+ * If it does, get the edge incident on the twin and check for the above condition
+ */
+ prev = edgelist[i].previous;
+ next = edgelist[i].next;
+
+ v1 = prev->edge_id.first;
+ v2 = prev->edge_id.second;
+
+ if (checkIfIntersects(twoDimCoords(vertex->coordinates), twoDimCoords(ortho_proj), \
+ twoDimCoords(getVertex(v1)->coordinates), twoDimCoords(getVertex(v2)->coordinates))) {
+ continue;
+ }
+
+ v1 = next->edge_id.first;
+ v2 = next->edge_id.second;
+
+ if (checkIfIntersects(twoDimCoords(vertex->coordinates), twoDimCoords(ortho_proj), \
+ twoDimCoords(getVertex(v1)->coordinates), twoDimCoords(getVertex(v2)->coordinates))) {
+ continue;
+ }
+
+ /* 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];
+ }
+ }
+ return edge;
+}
+
+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;
+ if (face_id == UNBOUNDED_FACE)
+ return false;
+
+ for (int i = 0; i < 3; i++) {
+ if (edge->edge_id.first == vertex_id)
+ return true;
+ edge = edge->next;
+ }
+ 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)
+{
+ 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 i;
+ }
+ }
+ return -1;
+}
+
+DCEL_Edge**
+PolygonalMesh::findFreeEdgeChain()
+{
+ DCEL_Edge *edge, *trav_edge, *start_edge, *end_edge;
+ int num_edges = edgelist.size();
+ bool *is_cyclic = new bool[num_edges];
+
+ DCEL_Edge* *start_and_end = new DCEL_Edge*[2];
+
+ for (unsigned int i = 0; i < edgelist.size(); i++) {
+ is_cyclic[i] = true;
+ }
+
+ for (unsigned int i = 0; i < edgelist.size(); ++i) {
+
+ /* Edge is part of a loop */
+ if (is_cyclic[i])
+ continue;
+
+ edge = &edgelist[i];
+ trav_edge = edge;
+
+ while (isFreeEdge(trav_edge)) {
+ if (trav_edge == edge) {
+ end_edge = trav_edge->previous;
+ trav_edge = trav_edge->previous;
+ }
+ trav_edge = trav_edge->previous;
+ }
+ start_edge = trav_edge->next;
+
+ if (end_edge->next != start_edge) {
+ trav_edge = edge;
+ while (isFreeEdge(trav_edge)) {
+ trav_edge = trav_edge->next;
+ }
+ end_edge = trav_edge->previous;
+ }
+
+ /* Checking if there's a loop */
+ if (start_edge == end_edge->next) {
+ trav_edge = start_edge;
+ do {
+ is_cyclic[getEdgeIndex(std::make_pair(trav_edge->edge_id.first, trav_edge->edge_id.second))] = true;
+ trav_edge = trav_edge->next;
+ }
+ while (trav_edge != start_edge);
+ continue;
+ }
+
+ start_and_end[0] = start_edge;
+ start_and_end[1] = end_edge;
+ return start_and_end;
+ }
+ return NULL;
+}
+
+/*
+ * Local Variables:
+ * mode: C
+ * tab-width: 8
+ * indent-tabs-mode: t
+ * c-file-style: "stroustrup"
+ * End:
+ * ex: shiftwidth=4 tabstop=8
+ */
Index: libanalyze/MeshHealing/MeshConversion.h
===================================================================
--- libanalyze/MeshHealing/MeshConversion.h (revision 0)
+++ libanalyze/MeshHealing/MeshConversion.h (working copy)
@@ -0,0 +1,129 @@
+#include <utility>
+#include <vector>
+
+#include "DCEL.h"
+#include "Geometry.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 setVertices() = 0;
+ virtual void setVertexCoords(int ID) = 0;
+ virtual void setVertex(int v1, int v2) = 0;
+ virtual void deleteVertex(int ID) = 0;
+ virtual void setFaces() = 0;
+ virtual void deleteFace(int ID) = 0;
+ virtual void addFace(int ID) = 0;
+ virtual void addVertex(int ID) = 0;
+
+public:
+
+ /* 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);
+ /* 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);
+ /* Replace all references of vertex v2 with v1 */
+ void setVertexRecord(DCEL_Vertex *v1, DCEL_Vertex *v2);
+ /* 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);
+ /* 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, DCEL_Vertex *opp_vertex);
+ /* Adds new face and sets the starting edge and returns the ID */
+ int addFaceRecord(DCEL_Edge *start_edge);
+ /* Deletes duplicate unused edges-if there are edges with the same start point and end point, delete the one with the unbounded incident face */
+ void deleteDuplicateEdges();
+ /* 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);
+ /*Adds new vertex with the given coordinates and returns ID*/
+ int addVertexRecord(double *coordinates);
+ /* 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();
+};
+
+/*
+ * Local Variables:
+ * mode: C
+ * tab-width: 8
+ * indent-tabs-mode: t
+ * c-file-style: "stroustrup"
+ * End:
+ * ex: shiftwidth=4 tabstop=8
+ */
Index: libanalyze/MeshHealing/MeshConversion_brlcad.cpp
===================================================================
--- libanalyze/MeshHealing/MeshConversion_brlcad.cpp (revision 0)
+++ libanalyze/MeshHealing/MeshConversion_brlcad.cpp (working copy)
@@ -0,0 +1,390 @@
+#include "MeshConversion_brlcad.h"
+
+#include <malloc.h>
+#include <stddef.h>
+#include <utility>
+#include <vector>
+
+
+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();
+
+ for (unsigned int i = 0; i < num_vertices; i++) {
+
+ vertexlist[i].vertex_id = i;
+
+ vertexlist[i].coordinates[0] = bot->vertices[VERTICES_PER_FACE * i];
+ vertexlist[i].coordinates[1] = bot->vertices[VERTICES_PER_FACE * i + 1];
+ vertexlist[i].coordinates[2] = bot->vertices[VERTICES_PER_FACE * i + 2];
+ }
+}
+
+
+/* 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();
+
+ /* Will contain the coordinates of the three vertex list */
+ double coords[3][3];
+ double determinant;
+
+ 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 */
+
+ /* Finding the determinant formed by the coordinates of the three vertices to get the orientation */
+ for (unsigned int j = 0; j < 3; j++) {
+ coords[0][j] = bot->vertices[vertex1 * 3 + j];
+ coords[1][j] = bot->vertices[vertex2 * 3 + j];
+ coords[2][j] = bot->vertices[vertex3 * 3 + j];
+
+ }
+
+ determinant = calcDet(coords);
+
+ /* If the orientation is clockwise, switch the second and third vertices */
+ if (determinant > 0)
+ 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 && \
+ v2 == edgelist[k].edge_id.first && v1 == 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]);
+ id = std::make_pair(v1, v2);
+ temp.edge_id = id;
+ edgelist.push_back(temp);
+
+ /* Setting the DCEL_Edge record for the half record from v1 to v2 */
+ temp.origin = &(vertexlist[v2]);
+ id = std::make_pair(v2, v1);
+ temp.edge_id = id;
+ edgelist.push_back(temp);
+ }
+ }
+ }
+}
+
+void
+BrlcadMesh::initIncidentFace()
+{
+
+ int vertex1, vertex2;
+ unsigned int num_faces = this->getNumFaces();
+ int v1, v2, v3;
+ int *a, *b;
+
+ a = new int[2];
+ b = new int[3];
+
+ for (unsigned int i = 0; i < edgelist.size(); ++i) {
+
+ 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 + 1; j++) {
+ v1 = bot->faces[VERTICES_PER_FACE * j];
+ v2 = bot->faces[VERTICES_PER_FACE * j + 1];
+ v3 = bot->faces[VERTICES_PER_FACE * j + 2];
+
+ 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]);
+ break;
+ }
+ }
+ /* If no DCEL_Face has been set till here, set the unbounded DCEL_Face to be the incident DCEL_Face */
+ if (edgelist[i].incident_face == NULL) {
+ edgelist[i].incident_face = &(facelist[UNBOUNDED_FACE]);
+ }
+ }
+
+}
+
+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;
+
+ 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].next = 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;
+
+ 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].previous = 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()
+{
+ free(bot->vertices);
+ bot->vertices = (fastf_t*)malloc(sizeof(fastf_t) * vertexlist.size() * 3);
+
+ for (unsigned int i = 0; i < vertexlist.size(); i++) {
+
+ /* Check for edges incident on this vertex and update the ID to the current value of i */
+ for (unsigned int j = 0; j < edgelist.size(); ++j) {
+ if (edgelist[j].edge_id.first == vertexlist[i].vertex_id)
+ edgelist[j].edge_id.first = i;
+
+ else if (edgelist[i].edge_id.second == vertexlist[i].vertex_id)
+ edgelist[j].edge_id.second = i;
+ }
+
+ vertexlist[i].vertex_id = i;
+
+ for (int k = 0; k < 3; k++)
+ bot->vertices[3 * i + k] = vertexlist[i].coordinates[k];
+ }
+}
+
+void
+BrlcadMesh::setFaces()
+{
+ DCEL_Edge *edge;
+
+ /* Unbounded face is the first record, so skip */
+ for (unsigned int i = 1; i < facelist.size(); ++i) {
+ facelist[i].face_id = 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 v1, int v2)
+{
+ for (int i = 0; i < getNumFaces(); ++i) {
+ for (int j = 0; j < VERTICES_PER_FACE; ++j) {
+ if (bot->faces[VERTICES_PER_FACE * i + j] == v1)
+ bot->faces[VERTICES_PER_FACE * i + j] = v2;
+ }
+ }
+}
+
+void
+BrlcadMesh::deleteVertex(int ID)
+{
+ for (int i = 3 * ID ; i < 3 * getNumVertices(); i++) {
+ bot->vertices[i] = bot->vertices[i + 1];
+ }
+
+ bot->num_vertices -= 1;
+}
+
+void
+BrlcadMesh::deleteFace(int ID)
+{
+ for (int i = 3 * (ID - 1) ; i < 3 * getNumFaces(); i++) {
+ bot->faces[i] = bot->faces[i + 1];
+ }
+ bot->num_faces -= 1;
+}
+
+void
+BrlcadMesh::addFace(int ID)
+{
+ int *new_faces = (int*)realloc(bot->faces, getNumFaces() + 1);
+ bot->faces = new_faces;
+
+ DCEL_Edge *edge = facelist[ID].start_edge;
+
+ for (int i = 0; i < 3; i++) {
+ bot->faces[3 * (ID - 1) + i] = edge->edge_id.first;
+ edge = edge->next;
+ }
+ bot->num_faces += 1;
+}
+
+void
+BrlcadMesh::addVertex(int ID)
+{
+ fastf_t *new_vertices = (fastf_t*)realloc(bot->vertices, getNumVertices() + 1);
+
+ 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: libanalyze/MeshHealing/MeshConversion_brlcad.h
===================================================================
--- libanalyze/MeshHealing/MeshConversion_brlcad.h (revision 0)
+++ libanalyze/MeshHealing/MeshConversion_brlcad.h (working copy)
@@ -0,0 +1,43 @@
+#include "MeshConversion.h"
+
+#include "rt/geom.h"
+#include "rt/primitives/bot.h"
+#include "vmath.h"
+
+class BrlcadMesh: public PolygonalMesh {
+
+private:
+ rt_bot_internal *bot;
+
+ BrlcadMesh(rt_bot_internal *bot_mesh);
+
+ 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(int ID);
+ void addVertex(int ID);
+};
+
+/*
+ * Local Variables:
+ * mode: C
+ * tab-width: 8
+ * indent-tabs-mode: t
+ * c-file-style: "stroustrup"
+ * End:
+ * ex: shiftwidth=4 tabstop=8
+ */
Index: libanalyze/MeshHealing/Zipper.cpp
===================================================================
--- libanalyze/MeshHealing/Zipper.cpp (revision 0)
+++ libanalyze/MeshHealing/Zipper.cpp (working copy)
@@ -0,0 +1,175 @@
+#include "Zipper.h"
+#include "Geometry.h"
+#include <cstddef>
+
+void
+calcFeaturePair(PolygonalMesh *polymesh, queue_element *node)
+{
+ DCEL_Edge *closest_edge;
+ int new_vertex_id;
+ DCEL_Vertex *vertex = node->vertex;
+ /* IDs of the vertices of the closest edge */
+ DCEL_Vertex *v1, *v2;
+
+ /* Get the closest edge */
+ closest_edge = polymesh->findClosestEdge(vertex);
+ 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);
+ }
+
+ else {
+ /* Returns the closest end point */
+ new_vertex_id = polymesh->addVertexRecord(findOrthogonalProjection((polymesh->getVertex(closest_edge->edge_id.first))->coordinates, \
+ (polymesh->getVertex(closest_edge->edge_id.second))->coordinates, vertex->coordinates));
+
+ node->feature_vertex = polymesh->getVertex(new_vertex_id);
+ node->dist = shortestDistBwPoints(vertex->coordinates, node->feature_vertex->coordinates);
+ }
+}
+
+void
+initPriorityQueue(PolygonalMesh *polymesh, DCEL_Edge *start, DCEL_Edge *end)
+{
+ queue_element node;
+ DCEL_Vertex *vertex;
+ DCEL_Edge *edge;
+ int v_id;
+
+ edge = start;
+ do
+ {
+ /* Consider the starting v_id of each edge, the ending v_id of the end edge will be handled separately */
+ 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 */
+ PQ.push(node);
+
+ edge = edge->next;
+ }
+
+ while (end->next->next != edge);
+}
+
+void
+zipperGaps(PolygonalMesh *polymesh, double tolerance)
+{
+ double coords[3];
+
+ double *A, *B, *C, *D;
+
+ /* Temporary priority queue */
+ std::priority_queue<queue_element, std::vector<queue_element>, compareDist> temp_PQ;
+
+ /* IDs of the vertices of the closest edge */
+ DCEL_Vertex *v1, *v2;
+
+ while (!PQ.empty()) {
+
+ /* 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 curr_element;
+ do {
+ curr_element = PQ.top();
+ PQ.pop();
+
+ if (!isValidFeature(polymesh, &curr_element)) {
+ calcFeaturePair(polymesh, &curr_element);
+ }
+ else {
+ if (curr_element.feature_vertex != NULL)
+ curr_element.dist = shortestDistBwPoints(curr_element.vertex->coordinates, curr_element.feature_vertex->coordinates);
+ else {
+ v1 = (polymesh->getVertex(curr_element.feature_edge->edge_id.first));
+ v2 = (polymesh->getVertex(curr_element.feature_edge->edge_id.second));
+ curr_element.dist = shortestDistToLine(v1->coordinates, v2->coordinates, curr_element.vertex->coordinates);
+ }
+ }
+
+ temp_PQ.push(curr_element);
+ }
+ while (!PQ.empty());
+
+ PQ = temp_PQ;
+
+ /* 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) {
+
+ /* 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 */
+ polymesh->deleteVertexRecord(curr_element.feature_vertex->vertex_id);
+ }
+
+ /* 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) {
+
+ 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;
+
+ D = findOrthogonalProjection(A, B, C);
+
+ polymesh->setVertexCoordsInRecord(curr_element.vertex->vertex_id, D);
+
+ /* Insert a vertex with coordinates D on the feature edge. */
+ polymesh->insertVertexOnEdge(curr_element.vertex, curr_element.feature_edge);
+
+ /* Split the face incident on the half edge of the feature edge with a bounded face incident */
+ polymesh->splitFace(curr_element.feature_edge->twin->incident_face->face_id, curr_element.vertex, curr_element.feature_edge->next->origin);
+
+ }
+ }
+ polymesh->deleteDuplicateEdges();
+}
+
+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 (polymesh->getEdge(std::make_pair(node->feature_edge->edge_id.first, node->feature_edge->edge_id.second)) == NULL)
+ return 0;
+ return 1;
+ }
+}
+
+/*
+ * Local Variables:
+ * mode: C
+ * tab-width: 8
+ * indent-tabs-mode: t
+ * c-file-style: "stroustrup"
+ * End:
+ * ex: shiftwidth=4 tabstop=8
+ */
Index: libanalyze/MeshHealing/Zipper.h
===================================================================
--- libanalyze/MeshHealing/Zipper.h (revision 0)
+++ libanalyze/MeshHealing/Zipper.h (working copy)
@@ -0,0 +1,53 @@
+#include "MeshConversion.h"
+
+#include <queue>
+
+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;
+ }
+};
+
+
+/* Priority queue ordered based on the error measure */
+std::priority_queue<queue_element, std::vector<queue_element>, compareDist> PQ;
+
+/* Calculate error measures for the vertices on the free edge chain and push the queue_element variable to the priority queue */
+void initPriorityQueue(PolygonalMesh *polymesh, DCEL_Edge *start, DCEL_Edge *end);
+
+/* Take the priority queue and zippers all those vertices whose error measures fall within the tolerance for zippering */
+void zipperGaps(PolygonalMesh *polymesh, double tolerance);
+
+/* 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);
+
+/*
+ * Local Variables:
+ * mode: C
+ * tab-width: 8
+ * indent-tabs-mode: t
+ * c-file-style: "stroustrup"
+ * End:
+ * ex: shiftwidth=4 tabstop=8
+ */
Index: libged/bot.c
===================================================================
--- libged/bot.c (revision 68040)
+++ libged/bot.c (working copy)
@@ -34,7 +34,43 @@
#include "wdb.h"
#include "./ged_private.h"
+int
+ged_heal(struct ged *gedp, int argc, const char *argv[])
+{
+ struct directory *bot_dp;
+ struct rt_db_internal intern;
+ struct rt_bot_internal *bot;
+ const char *cmd = argv[0];
+ const char *primitive = NULL;
+ GED_CHECK_DATABASE_OPEN(gedp, GED_ERROR);
+ GED_CHECK_READ_ONLY(gedp, GED_ERROR);
+ GED_CHECK_ARGC_GT_0(gedp, argc, GED_ERROR);
+
+ /* initialize result */
+ bu_vls_trunc(gedp->ged_result_str, 0);
+
+ primitive = argv[argc - 1];
+
+ /* get bot */
+ GED_DB_LOOKUP(gedp, bot_dp, primitive, LOOKUP_NOISY, GED_ERROR & GED_QUIET);
+ GED_DB_GET_INTERNAL(gedp, &intern, bot_dp, bn_mat_identity, &rt_uniresource, GED_ERROR);
+
+ if (intern.idb_major_type != DB5_MAJORTYPE_BRLCAD || intern.idb_minor_type != DB5_MINORTYPE_BRLCAD_BOT) {
+ bu_vls_printf(gedp->ged_result_str, "%s: %s is not a BOT solid!", cmd, primitive);
+ return GED_ERROR;
+ }
+
+ bot = (struct rt_bot_internal *)intern.idb_ptr;
+ RT_BOT_CK_MAGIC(bot);
+
+ analyse_heal_bot(bot);
+
+ bu_vls_printf(gedp->ged_result_str, "Heal command");
+
+ return GED_OK;
+}
+
int
ged_bot(struct ged *gedp, int argc, const char *argv[])
{
Index: libtclcad/tclcad_obj.c
===================================================================
--- libtclcad/tclcad_obj.c (revision 68040)
+++ 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: libwdb/bot.c
===================================================================
--- libwdb/bot.c (revision 68040)
+++ libwdb/bot.c (working copy)
@@ -38,6 +38,8 @@
#include "raytrace.h"
#include "wdb.h"
+#include "../src/libanalyze/MeshHealing/MeshConversion.h"
+#include "../src/libanalyze/MeshHealing/Zipper.h"
int
mk_bot_w_normals(
@@ -142,7 +144,18 @@
faces, thickness, face_mode, 0, NULL, NULL));
}
+int
+analyse_heal_bot(struct rt_bot_internal *bot)
+{
+ PolygonalMesh polymesh = new BrlcadMesh(bot);
+ double tolerance = 5;
+ DCEL_Edge **start_and_end = polymesh.findFreeEdgeChain();
+ initPriorityQueue(polymesh, start_and_end[0], start_and_end[1]);
+ zipperGaps(polymesh, tolerance);
+}
+
+
/*
* Local Variables:
* mode: C
Index: mged/setup.c
===================================================================
--- mged/setup.c (revision 68040)
+++ 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},
------------------------------------------------------------------------------
Attend Shape: An AT&T Tech Expo July 15-16. Meet us at AT&T Park in San
Francisco, CA to explore cutting-edge tech and listen to tech luminaries
present their vision of the future. This family event has something for
everyone, including kids. Get more information and register today.
http://sdm.link/attshape
_______________________________________________
BRL-CAD Developer mailing list
brlcad-devel@lists.sourceforge.net
https://lists.sourceforge.net/lists/listinfo/brlcad-devel
