Skip to content

Commit

Permalink
update from Amal's 24/01/31 email
Browse files Browse the repository at this point in the history
  • Loading branch information
sloriot committed Feb 22, 2024
1 parent abcfd74 commit 3cea793
Show file tree
Hide file tree
Showing 4 changed files with 240 additions and 16 deletions.
225 changes: 225 additions & 0 deletions Ball_merge_surface_reconstruction/attic/AutomaticBallMerge.cpp
Original file line number Diff line number Diff line change
@@ -0,0 +1,225 @@
#include "happly.h"
#include <CGAL/Delaunay_triangulation_3.h>
#include <CGAL/Delaunay_triangulation_cell_base_3.h>
#include <CGAL/Exact_predicates_inexact_constructions_kernel.h>
#include <CGAL/IO/Color.h>
#include <CGAL/Triangulation_cell_base_with_info_3.h>
#include <CGAL/Triangulation_vertex_base_with_info_3.h>
#include <bits/stdc++.h>
#include <fstream>
#include <iostream>
#include <queue>
#include <sstream>
#include <string>
#include <sys/time.h>
using namespace CGAL;
using namespace std;
typedef CGAL::Exact_predicates_inexact_constructions_kernel K;
typedef CGAL::Triangulation_vertex_base_with_info_3<int, K> Vb;
typedef CGAL::Triangulation_cell_base_with_info_3<int, K> Cb;
#ifdef CGAL_LINKED_WITH_TBB
typedef CGAL::Triangulation_data_structure_3<Vb, Cb, Parallel_tag> Tds;
#else
typedef CGAL::Triangulation_data_structure_3<Vb, Cb> Tds;
#endif
typedef CGAL::Delaunay_triangulation_3<K, Tds> 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<node, std::vector<node>, std::function<bool(node, node)>> 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<K> 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<std::array<unsigned char, 3>> 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<std::pair<Point, unsigned>> points;
for (int i = 0; i < 3; i++) {
min[i] = std::numeric_limits<double>::max();
max[i] = -std::numeric_limits<double>::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<K> 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<Delaunay::Cell_handle> 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<std::array<double, 3>> meshVertexPositions(points.size());
std::vector<std::vector<int>> 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<int> 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);
}
20 changes: 10 additions & 10 deletions Ball_merge_surface_reconstruction/attic/BallMergeParallel.cpp
Original file line number Diff line number Diff line change
Expand Up @@ -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<Delaunay::Cell_handle> q; //A queue to do the grouping
q.push(fh); //The cell itself is mergable
while (!q.empty()){ //A traversal
Expand All @@ -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
}
Expand Down Expand Up @@ -119,20 +119,20 @@ 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;
if (option == 1){//If the user opted for global algorithm
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
Expand All @@ -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));
Expand All @@ -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<int> indices(3);
for (int j = 0; j < 3; j++)
Expand Down Expand Up @@ -206,4 +206,4 @@ int main(int argc, char **argv){
plyOut2.addFaceIndices(meshFaceIndices);
plyOut2.write(outname, happly::DataFormat::Binary);
}
}
}
Original file line number Diff line number Diff line change
Expand Up @@ -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);
Expand Down
Original file line number Diff line number Diff line change
Expand Up @@ -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<Cell_handle> q; //A queue to do the grouping
q.push(fh); //The cell itself is mergeable
Expand All @@ -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
}
}
Expand Down Expand Up @@ -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;
Expand Down Expand Up @@ -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<int> indices(3);
Expand Down

0 comments on commit 3cea793

Please sign in to comment.