diff --git a/Ball_merge_surface_reconstruction/attic/AutomaticBallMerge.cpp b/Ball_merge_surface_reconstruction/attic/AutomaticBallMerge.cpp new file mode 100644 index 000000000000..bb4232d5d997 --- /dev/null +++ b/Ball_merge_surface_reconstruction/attic/AutomaticBallMerge.cpp @@ -0,0 +1,225 @@ +#include "happly.h" +#include +#include +#include +#include +#include +#include +#include +#include +#include +#include +#include +#include +#include +using namespace CGAL; +using namespace std; +typedef CGAL::Exact_predicates_inexact_constructions_kernel K; +typedef CGAL::Triangulation_vertex_base_with_info_3 Vb; +typedef CGAL::Triangulation_cell_base_with_info_3 Cb; +#ifdef CGAL_LINKED_WITH_TBB +typedef CGAL::Triangulation_data_structure_3 Tds; +#else +typedef CGAL::Triangulation_data_structure_3 Tds; +#endif +typedef CGAL::Delaunay_triangulation_3 Delaunay; +typedef Delaunay::Point Point; +struct node { + float priority; + Delaunay::Cell_handle ch; + struct node *link; +} tempn; +bool Comparen(node n1, node n2) { return n1.priority < n2.priority; } +std::priority_queue, std::function> pq(Comparen); +std::ofstream myfile; +clock_t start, _end; +int group = 1, gcount = 0, option, maxg = 0, maingroup, max1 = 0, secondgroup; +float minx = 9999, miny = 9999, minz = 9999, maxx = -9999, maxy = -9999, maxz = -9999; +double par, bbdiaglen, tlen; +Delaunay T; +double sq(double a) { return a * a; } +Point cc(Delaunay::Cell_handle v) { + CGAL::Tetrahedron_3 t1(v->vertex(0)->point(), v->vertex(1)->point(),v->vertex(2)->point(), v->vertex(3)->point()); + return CGAL::circumcenter(t1); +} +double distance(Point p1, Point p2) { + return sqrt(sq(p2.x() - p1.x()) + sq(p2.y() - p1.y()) + sq(p2.z() - p1.z())); +} +double ir(Delaunay::Cell_handle fh1, Delaunay::Cell_handle fh2) { + Point p1 = cc(fh1); + Point p2 = cc(fh2); + double d = distance(p1, p2); + double a = distance(fh1->vertex(0)->point(), p1); + double b = distance(fh2->vertex(0)->point(), p2); + if (std::isnan(d)) + return 0; + if ((a + b - d) / a > (a + b - d) / b) + return (a + b - d) / a; + return (a + b - d) / b; +} + +int main(int argc, char **argv) { + const char *fname = argv[1]; + const char *inFilename = argv[1]; + std::ifstream inStream(inFilename); + std::vector> meshVertexColors; + Point p; + float minx = 9999, miny = 9999, minz = 9999, maxx = -9999, maxy = -9999, maxz = -9999; + double min[3], max[3]; + int vIndex = 0; + std::vector> points; + for (int i = 0; i < 3; i++) { + min[i] = std::numeric_limits::max(); + max[i] = -std::numeric_limits::max(); + } + std::cout<<"Reading the file\n"; + while (!inStream.eof()) { + string line; + getline(inStream, line); + istringstream str(line); + if (str >> p) { + int r, g, b; + if (str >> r) { + str >> g; + str >> b; + meshVertexColors.push_back({(unsigned char)r, (unsigned char)g, (unsigned char)b}); + } + if (p.x() < minx) + minx = p.x(); + if (p.x() < miny) + miny = p.y(); + if (p.x() < minz) + minz = p.z(); + if (p.x() > maxx) + maxx = p.x(); + if (p.x() > maxy) + maxy = p.y(); + if (p.x() > maxz) + maxz = p.z(); + points.push_back(std::make_pair(Point(p.x(), p.y(), p.z()), vIndex++)); + for (int i = 0; i < 3; i++) { + if (p[i] < min[i]) + min[i] = p[i]; + + if (p[i] > max[i]) + max[i] = p[i]; + } + } + } + std::cout<<"Computing Delaunay\n"; +#ifdef CGAL_LINKED_WITH_TBB + Bbox_3 bbox = Bbox_3(min[0], min[1], min[2], max[0], max[1], max[2]); + Delaunay::Lock_data_structure locking_ds(bbox, 50); + Delaunay T(points.begin(), points.end(), &locking_ds); +#else + Delaunay T(points.begin(), points.end()); +#endif +std::cout<<"Running BallMerge\n"; + double maxr = -9999; + Delaunay::Cell_handle maxcell; + Delaunay::Finite_cells_iterator vit; + for (vit = T.finite_cells_begin(); vit != T.finite_cells_end(); vit++) {//Loop to pick the seed tetrahedron to grow from + CGAL::Tetrahedron_3 t1(vit->vertex(0)->point(), vit->vertex(1)->point(),vit->vertex(2)->point(), vit->vertex(3)->point()); + Point tp1 = CGAL::circumcenter(t1); + double dis1 = distance(tp1, vit->vertex(0)->point()); + if (dis1 > maxr) { + Delaunay::Cell_handle ch = T.locate(tp1); + Delaunay::Cell_handle ch1 = vit; + if (!T.is_infinite(ch) && (T.is_infinite(vit->neighbor(0)) || T.is_infinite(vit->neighbor(2)) || T.is_infinite(vit->neighbor(1)) || T.is_infinite(vit->neighbor(3)))){//Tetrahedron with largest circumradius which is part of the convex hull + maxcell = ch; + maxr = dis1; + } + } + vit->info() = 0; + } + maxcell->info() = 1; + for (int i = 0; i < 4; i++) + if (!T.is_infinite(maxcell->neighbor(i))) { + tempn.priority = ir(maxcell, maxcell->neighbor(i)); + tempn.ch = maxcell->neighbor(i); + pq.push(tempn); + } + double maxir = ir(maxcell, maxcell->neighbor(0)); + for (int i = 0; i < 4; i++) + if (maxir < ir(maxcell, maxcell->neighbor(i))) + maxir = ir(maxcell, maxcell->neighbor(i)); + int loop = 0; + maxir = 2.0; + while (maxir > 1.0) {//Loop for 200 iterations maximum since maxir decreases in 0.01 steps + loop++; + while (!pq.empty()) {//Create groups from the seed tetrahedron for a particular parameter + node tch = pq.top(); + pq.pop(); + if (tch.ch->info() == 0) { + if (tch.priority > maxir) { + tch.ch->info() = loop; + for (int i = 0; i < 4; i++) + if ((!T.is_infinite(tch.ch->neighbor(i))) && (tch.ch->neighbor(i)->info() == 0)) { + tempn.priority = ir(tch.ch, tch.ch->neighbor(i)); + tempn.ch = tch.ch->neighbor(i); + pq.push(tempn); + } + } else { + pq.push(tch); + break; + } + } + } + int good = 0, masked = 0, separate = 0; + for (Delaunay::Finite_vertices_iterator vIter = T.finite_vertices_begin(); + vIter != T.finite_vertices_end(); vIter++) { + std::vector vif; + T.incident_cells(vIter, back_inserter(vif)); + int convfl = 0, naddedfl = 0, addedfl = 0; + for (int i = 0; i < vif.size(); i++) { + if (T.is_infinite(vif.at(i))) + convfl = 1; + else if (vif.at(i)->info() == 0) + naddedfl = 1; + else + addedfl = 1; + } + if (addedfl == 1 && (naddedfl == 1 || convfl == 1)) + good++; + if (addedfl == 1 && (naddedfl == 0 && convfl == 0)) + masked++; + if (addedfl == 0) + separate++; + } + if (good < points.size() * 0.9999999 && masked < points.size() * 0.05)//Check how many points got removed compared to the result generated using the previous parameter + maxir -= 0.01;//Step size for parameter change + else + break; + } + std::cout << "Automatic Parameter is " << maxir + 0.01 << "\n"; + std::cout<<"Writing to the file\n"; + std::vector> meshVertexPositions(points.size()); + std::vector> meshFaceIndices; + for (Delaunay::Finite_vertices_iterator vIter = T.finite_vertices_begin(); + vIter != T.finite_vertices_end(); vIter++) { + Point point = vIter->point(); + int vIndex = vIter->info(); + meshVertexPositions[vIndex][0] = point.x(); + meshVertexPositions[vIndex][1] = point.y(); + meshVertexPositions[vIndex][2] = point.z(); + } + std::string st = argv[1]; + st = st + "out.ply"; + char outname[100]; + strcpy(outname, st.c_str()); + bbdiaglen = distance(Point(minx, miny, minz), Point(maxx, maxy, maxz)); + for (vit = T.finite_cells_begin(); vit != T.finite_cells_end(); vit++) + for (int i = 0; i < 4; i++) + if ((vit->info() != 0) && (vit->info() != loop) && ((vit->neighbor(i)->info() == 0 || vit->neighbor(i)->info() == loop) || T.is_infinite(vit->neighbor(i)))) { + std::vector indices(3); + for (int j = 0; j < 3; j++) + indices[j] = vit->vertex((i + 1 + j) % 4)->info(); + meshFaceIndices.push_back(indices); + } + happly::PLYData plyOut; + plyOut.addVertexPositions(meshVertexPositions); + if (meshVertexColors.size() > 0) + plyOut.addVertexColors(meshVertexColors); + plyOut.addFaceIndices(meshFaceIndices); + plyOut.write(outname, happly::DataFormat::Binary); +} diff --git a/Ball_merge_surface_reconstruction/attic/BallMergeParallel.cpp b/Ball_merge_surface_reconstruction/attic/BallMergeParallel.cpp index a2105cebad50..a71b23420022 100644 --- a/Ball_merge_surface_reconstruction/attic/BallMergeParallel.cpp +++ b/Ball_merge_surface_reconstruction/attic/BallMergeParallel.cpp @@ -46,7 +46,7 @@ int _function(Delaunay::Cell_handle fh1, Delaunay::Cell_handle fh2){ //Function return 1; return 0; } -void Group(Delaunay::Cell_handle fh){ //Function to recursively group all the cells that is mergable from the cell fh with user specified parameter +void Group(Delaunay &T, Delaunay::Cell_handle fh){ //Function to recursively group all the cells that is mergable from the cell fh with user specified parameter std::queue q; //A queue to do the grouping q.push(fh); //The cell itself is mergable while (!q.empty()){ //A traversal @@ -55,7 +55,7 @@ void Group(Delaunay::Cell_handle fh){ //Function to recursively group all the ce fh->info() = group; //A label gcount++; for (int i = 0; i < 4; i++) //For each neighbors - if (fh->neighbor(i)->info() == 0 && _function(fh, fh->neighbor(i)) == 1){//If it is unlabeled and can be merged - using the function that computes IR + if ((!T.is_infinite(fh->neighbor(i)))&&(fh->neighbor(i)->info() == 0 && _function(fh, fh->neighbor(i)) == 1)){//If it is unlabeled and can be merged - using the function that computes IR fh->neighbor(i)->info() = group;//If mergable, then change the label q.push(fh->neighbor(i));//Push it into the traversal list } @@ -119,12 +119,12 @@ int main(int argc, char **argv){ } } #ifdef CGAL_LINKED_WITH_TBB //Parallel Delaunay computation - Bbox_3 bbox = Bbox_3(min[0], min[1], min[2], max[0], max[1], max[2]); - Delaunay::Lock_data_structure locking_ds(bbox, 50); + Bbox_3 bbox = Bbox_3(min[0], min[1], min[2], max[0], max[1], max[2]); + Delaunay::Lock_data_structure locking_ds(bbox, 50); Delaunay T(points.begin(), points.end(), &locking_ds); -#else + #else Delaunay T(points.begin(), points.end());//General Delaunay computation -#endif + #endif Delaunay::Finite_cells_iterator vit; for (vit = T.finite_cells_begin(); vit != T.finite_cells_end(); vit++)//Initialize the labels of all tetrahedrons vit->info() = 0; @@ -132,7 +132,7 @@ int main(int argc, char **argv){ for (vit = T.finite_cells_begin(); vit != T.finite_cells_end(); vit++){//For each cell if (vit->info() == 0){//If the cell label is not altered vit->info() = group;//Assign a label - Group(vit);//Group all the mergable cells + Group(T,vit);//Group all the mergable cells if (maxg < gcount){//Remember the largest group - based on the number of tetrahedrons in the group max1 = maxg; secondgroup = maingroup;//To remember the second largest group @@ -158,7 +158,7 @@ int main(int argc, char **argv){ meshVertexPositions[vIndex][2] = point.z(); } std::string st = argv[1]; - st = st + std::to_string(par) + "out.ply"; + st = st + std::to_string(par) + "out"+std::to_string(tlen)+".ply"; char outname[100]; strcpy(outname, st.c_str()); bbdiaglen = distance(Point(minx, miny, minz), Point(maxx, maxy, maxz)); @@ -174,7 +174,7 @@ int main(int argc, char **argv){ } else if (vit->neighbor(i)->info() != 9999)//If local if (bblen(vit->vertex((i + 1) % 4)->point(), vit->vertex((i + 2) % 4)->point(), vit->vertex((i + 3) % 4)->point()))//If the triangle crosses our bbdiagonal based criteria - if (!_function(vit, vit->neighbor(i)) == 1){//If the cells cannot be merged, then wirte the triangle between these two cells to the PLY file + if (!_function(vit, vit->neighbor(i)) == 1||T.is_infinite(vit->neighbor(i))){//If the cells cannot be merged, then wirte the triangle between these two cells to the PLY file vit->info() = 9999; std::vector indices(3); for (int j = 0; j < 3; j++) @@ -206,4 +206,4 @@ int main(int argc, char **argv){ plyOut2.addFaceIndices(meshFaceIndices); plyOut2.write(outname, happly::DataFormat::Binary); } -} \ No newline at end of file +} diff --git a/Ball_merge_surface_reconstruction/examples/Ball_merge_surface_reconstruction/ball_merge_reconstruction.cpp b/Ball_merge_surface_reconstruction/examples/Ball_merge_surface_reconstruction/ball_merge_reconstruction.cpp index e09e85ee8e51..8ef0463bc886 100644 --- a/Ball_merge_surface_reconstruction/examples/Ball_merge_surface_reconstruction/ball_merge_reconstruction.cpp +++ b/Ball_merge_surface_reconstruction/examples/Ball_merge_surface_reconstruction/ball_merge_reconstruction.cpp @@ -57,7 +57,7 @@ int main(int argc, char **argv) bmsr.result(meshVertexPositions, meshFaceIndices); std::string st = argv[1]; - st = st + std::to_string(par) + "out.ply"; + st = st + std::to_string(par) + "out"+std::to_string(tlen)+".ply"; happly::PLYData plyOut; plyOut.addVertexPositions(meshVertexPositions); diff --git a/Ball_merge_surface_reconstruction/include/CGAL/Ball_merge_surface_reconstruction.h b/Ball_merge_surface_reconstruction/include/CGAL/Ball_merge_surface_reconstruction.h index ee19ebaa88ee..c14c23394fa8 100644 --- a/Ball_merge_surface_reconstruction/include/CGAL/Ball_merge_surface_reconstruction.h +++ b/Ball_merge_surface_reconstruction/include/CGAL/Ball_merge_surface_reconstruction.h @@ -65,7 +65,7 @@ class Ball_merge_surface_reconstruction { } //Function to recursively group all the cells that are mergeable from the cell fh with user specified parameter - void regroup(Cell_handle fh) + void regroup(const Delaunay& dt3, Cell_handle fh) { std::queue q; //A queue to do the grouping q.push(fh); //The cell itself is mergeable @@ -75,8 +75,7 @@ class Ball_merge_surface_reconstruction { fh->info() = group; //A label ++gcount; for (int i = 0; i < 4; ++i){ //For each neighbor - if (fh->neighbor(i)->info() == 0 && _function(fh, fh->neighbor(i)) == 1){//If it is unlabeled and can be merged - using the function that computes IR - fh->neighbor(i)->info() = group;//If mergeable, then change the label + if ((!dt3.is_infinite(fh->neighbor(i)))&&(fh->neighbor(i)->info() == 0 && _function(fh, fh->neighbor(i)) == 1)){//If it is unlabeled and can be merged - using the function that computes IR fh->neighbor(i)->info() = group;//If mergeable, then change the label q.push(fh->neighbor(i));//Push it into the traversal list } } @@ -120,7 +119,7 @@ class Ball_merge_surface_reconstruction { for (vit = dt3.finite_cells_begin(); vit != dt3.finite_cells_end(); ++vit){//For each cell if (vit->info() == 0){//If the cell label is not altered vit->info() = group;//Assign a label - regroup(vit);//Group all the mergeable cells + regroup(dt3, vit);//Group all the mergeable cells if (maxg < gcount){ //Remember the largest group - based on the number of tetrahedra in the group max1 = maxg; @@ -168,7 +167,7 @@ class Ball_merge_surface_reconstruction { //If local AF: replace 9999. -1, numeric_limits, number_of_vertices ?? if (bblen(vit->vertex((i + 1) % 4)->point(), vit->vertex((i + 2) % 4)->point(), vit->vertex((i + 3) % 4)->point())) //If the triangle crosses our bbdiagonal based criteria - if (!_function(vit, vit->neighbor(i)) == 1){ + if (!_function(vit, vit->neighbor(i)) == 1||dt3.is_infinite(vit->neighbor(i))){ //If the cells cannot be merged, then write the triangle between these two cells to the PLY file vit->info() = 9999; std::vector indices(3);