Subject: CGAL users discussion list
List archive
- From: Kwok Jasper <>
- To: <>
- Subject: RE: [cgal-discuss] advice about CGAL Delaunay Triangulation
- Date: Mon, 2 Feb 2009 23:05:08 +0800
- Importance: Normal
#include <string.h>
/********************************************************SamplePoint.h******************************************************/
#include <fstream>
#include <list>
#include <CGAL/Delaunay_triangulation_3.h>
#include <CGAL/Exact_predicates_inexact_constructions_kernel.h>
#include <CGAL/Fuzzy_sphere.h>
#include <CGAL/Cartesian.h>
#include <CGAL/Triangulation_hierarchy_3.h>
typedef CGAL::Exact_predicates_inexact_constructions_kernel K1;
//typedef CGAL::Cartesian<double> K1;
/**
#include <CGAL/Simple_cartesian.h>
#include <CGAL/Filtered_kernel.h>
typedef CGAL::Simple_cartesian<double> CK;
typedef CGAL::Filtered_kernel<CK> K1;
**/
/**
typedef CGAL::Exact_predicates_inexact_constructions_kernel epick;
typedef CGAL::Lazy_exact_nt<CGAL::Cartesian> K1;
**/
/**
#include <CGAL/MP_Float.h>
#include <CGAL/Lazy_exact_nt.h>
#include <CGAL/Quotient.h>
#include <CGAL/Exact_predicates_inexact_constructions_kernel.h>
typedef CGAL::Lazy_exact_nt<CGAL::Exact_predicates_inexact_constructions_kernel<double> > NT;
typedef CGAL::Cartesian<NT> K1;
**/
/**
typedef CGAL::Triangulation_vertex_base_3<K1> Vb;
typedef CGAL::Triangulation_hierarchy_vertex_base_3<Vb> Vbh;
typedef CGAL::Triangulation_data_structure_3<Vbh> Tds;
typedef CGAL::Delaunay_triangulation_3<K1,Tds> Triangulation;
typedef CGAL::Triangulation_hierarchy_3<Triangulation> Triangulations;
**/
/**
typedef CGAL::Delaunay_triangulation_3<K, Tds> Triangulation;
typedef CGAL::Triangulation_hierarchy_3<Triangulation> Triangulations;
**/
typedef CGAL::Delaunay_triangulation_3<K1> Triangulations;
typedef Triangulations::Point Point_tr;
typedef Triangulations::Vertex Vertex;
typedef Triangulations::Vertex_handle Vertex_handle;
typedef Triangulations::Finite_vertices_iterator finite_vertices_itr;
typedef Triangulations::Edge Edge;
typedef Triangulations::Facet Facet;
typedef Triangulations::Finite_facets_iterator finite_facets_itr;
typedef Triangulations::Triangle Triangle;
typedef Triangulations::Segment Segment;
typedef CGAL::Cartesian<double> Rep_class;
typedef Rep_class::Point_3 Point;
typedef CGAL::Vector_3<Rep_class > Vector;
//k-neighbour
#include <CGAL/Orthogonal_k_neighbor_search.h>
#include <CGAL/Search_traits_3.h>
#include <math.h>
#include <iostream>
#include "mex.h"
#include "find_surface_normal.h"
using namespace std;
class PriorityQueueForNeighborSearchNode
{
public:
double key_dist;
Point point;
PriorityQueueForNeighborSearchNode(){}
PriorityQueueForNeighborSearchNode(const double dist, const double x_i, const double y_i, const double z_i);
PriorityQueueForNeighborSearchNode(const PriorityQueueForNeighborSearchNode& node_i);
friend class PriorityQueueForNeighborSearchTree;
};
PriorityQueueForNeighborSearchNode::PriorityQueueForNeighborSearchNode(const PriorityQueueForNeighborSearchNode& node_i)
{
point = Point(node_i.point.x(), node_i.point.y(), node_i.point.z());
key_dist = node_i.key_dist;
}
PriorityQueueForNeighborSearchNode::PriorityQueueForNeighborSearchNode(const double dist, const double x_i, const double y_i, const double z_i)
{
point = Point(x_i, y_i, z_i);
key_dist = dist;
}
class PriorityQueueForNeighborSearchTree
{
private:
vector<PriorityQueueForNeighborSearchNode> tree_content;
PriorityQueueForNeighborSearchTree();
void insert(const double dist, const double x_i, const double y_i, const double z_i);
Point delete_min();
};
PriorityQueueForNeighborSearchTree::PriorityQueueForNeighborSearchTree()
{
PriorityQueueForNeighborSearchNode node(-1, -1, -1, -1);
tree_content.push_back(node);
}
void PriorityQueueForNeighborSearchTree::insert(const double dist, const double x_i, const double y_i, const double z_i)
{
//tree_content.resize(tree_content.size());
PriorityQueueForNeighborSearchNode node(dist, x_i, y_i, z_i);
tree_content.push_back(node);
PriorityQueueForNeighborSearchNode temp;
for(int hole = tree_content.size()-1; ((hole > 1) && (tree_content.at(hole/2).key_dist > dist)); hole /= 2)
{
temp = tree_content.at(hole);
tree_content.at(hole) = tree_content.at(hole/2);
tree_content.at(hole/2) = temp;
}
}
Point PriorityQueueForNeighborSearchTree::delete_min()
{
double tempX = tree_content.at(1).point.x();
double tempY = tree_content.at(1).point.y();
double tempZ = tree_content.at(1).point.z();
Point result(tempX, tempY, tempZ);
tree_content.at(1) = tree_content.at(tree_content.size() - 1);
tree_content.pop_back();
int child;
PriorityQueueForNeighborSearchNode temp = tree_content.at(1);
int hole = 1;
for(hole = 1; ((hole * 2 <= tree_content.size() - 1) ); hole = child)
{
child = hole * 2;
if ((child != tree_content.size() -1) && (tree_content.at(child+1).key_dist < tree_content.at(child).key_dist))
{
++child;
}
if (tree_content.at(child).key_dist < temp.key_dist)
{
tree_content.at(hole) = tree_content.at(child);
}
else
{
break;
}
}
tree_content.at(hole) = temp;
return result;
}
class SamplePoint
{
private:
double sample_point[3]; //record the coordinate of the sample point
double surface_normal[3]; //record the surface normal at the sample point
list<Point> neighbour; //record the 55 nearest neighbour for the sample point
double fifth_nearest_dist; //record the fifth nearest distance for the sample point
SamplePoint(const double px, const double py, const double pz);
void set_surface_normal(const double nx, const double ny, const double nz);
double sample_point_x() const;
double sample_point_y() const;
double sample_point_z() const;
double surface_normal_x() const;
double surface_normal_y() const;
double surface_normal_z() const;
void set_fifth_nearest_dist(const double d);
//void add_neighbour(SamplePoint p);
void perturb_point();
friend class SamplePointSet;
//friend class Construct_coord_iterator;
friend class SamplePointKdTreeNode;
};
/**
class Construct_coord_iterator
{
public:
const double* operator()(const SamplePoint& sp) const
{
return static_cast<const double*>(sp.sample_point);
}
const double* operator()(const SamplePoint& sp, int) const
{
return static_cast<const double*>(sp.sample_point+3);
}
};
**/
//k-neighbour
typedef CGAL::Search_traits_3<Rep_class> Traits;
typedef CGAL::Orthogonal_k_neighbor_search<Traits> Neighbor_search;
typedef Neighbor_search::Tree Tree;
typedef Tree::iterator Tree_itr;
//typedef CGAL::Fuzzy_sphere<Traits> Fuzzy_sphere;
Tree sample_point_tree_row;
void SamplePoint::perturb_point()
{
Point p(sample_point[0], sample_point[1], sample_point[2]);
Neighbor_search search(sample_point_tree_row, p, 2);
Neighbor_search::iterator itr = search.begin();
++itr;
double dist = itr -> second;
double disturbance = dist * (0 + rand() * (0.5) / RAND_MAX);
//double disturbance = dist * ((rand() % 2 * 0.01) - (rand() % 2 * 0.01));
sample_point[0] += disturbance;
//disturbance = dist * (-0.1 + rand() * 0.1);
sample_point[1] += disturbance;
//disturbance = dist * (-0.1 + rand() * 0.1);
sample_point[2] += disturbance;
}
void SamplePoint::set_fifth_nearest_dist(const double d)
{
fifth_nearest_dist = d;
}
/**
void SamplePoint::add_neighbour(SamplePoint p)
{
neighbour.push_back(p);
}
**/
SamplePoint::SamplePoint(const double px, const double py, const double pz)
{
sample_point[0] = px;
sample_point[1] = py;
sample_point[2] = pz;
}
void SamplePoint::set_surface_normal(const double nx, const double ny, const double nz)
{
surface_normal[0] = nx;
surface_normal[1] = ny;
surface_normal[2] = nz;
}
double SamplePoint::sample_point_x() const
{
return sample_point[0];
}
double SamplePoint::sample_point_y() const
{
return sample_point[1];
}
double SamplePoint::sample_point_z() const
{
return sample_point[2];
}
double SamplePoint::surface_normal_x() const
{
return surface_normal[0];
}
double SamplePoint::surface_normal_y() const
{
return surface_normal[1];
}
double SamplePoint::surface_normal_z() const
{
return surface_normal[2];
}
class SelectedTriangleElement
{
private:
int facet_index;
vector<int> vertice_index;
SelectedTriangleElement(){}
friend class SelectedTriangleSet;
friend class SamplePointSet;
};
class SelectedTriangleElementConnection
{
private:
Point_tr vertice_1;
Point_tr vertice_2;
SelectedTriangleElementConnection(const Point_tr v1, const Point_tr v2);
friend class SelectedTriangleSet;
friend class SamplePointSet;
};
SelectedTriangleElementConnection::SelectedTriangleElementConnection(const Point_tr v1, const Point_tr v2)
{
vertice_1 = v1;
vertice_2 = v2;
}
class SelectedTriangleSet
{
private:
vector<SelectedTriangleElement> triangle_set;
vector<Point_tr> selected_vertice;
vector<SelectedTriangleElementConnection> connection_edge;
int num_of_facet;
void add_vertice_to_set(const Point_tr v1, const Point_tr v2, const Point_tr v3);
bool exist_vertice(const Point_tr query);
int get_vertice_index(const Point_tr v);
bool exist_connection_edge(const Point_tr query_1, const Point_tr query_2);
int get_num_of_vertices() const;
int get_num_of_edges() const;
int get_num_of_facets() const;
friend class SamplePointSet;
};
int SelectedTriangleSet::get_num_of_vertices() const
{
return selected_vertice.size();
}
int SelectedTriangleSet::get_num_of_edges() const
{
return (2 * (selected_vertice.size() - 3));
}
int SelectedTriangleSet::get_num_of_facets() const
{
return (2 * (selected_vertice.size() - 2));
}
int SelectedTriangleSet::get_vertice_index(const Point_tr v)
{
for(int i=0; i<selected_vertice.size(); ++i)
{
if ((selected_vertice.at(i).x() == v.x()) && (selected_vertice.at(i).y() == v.y()) && (selected_vertice.at(i).z() == v.z()))
{
return i;
}
}
return -1;
}
bool SelectedTriangleSet::exist_connection_edge(const Point_tr query_1, const Point_tr query_2)
{
vector<SelectedTriangleElementConnection>::iterator ec_itr;
bool b1;
bool b2;
bool b3;
bool b4;
bool b5;
bool b6;
for(ec_itr = connection_edge.begin(); ec_itr != connection_edge.end(); ++ec_itr)
{
b1 = (((ec_itr -> vertice_1).x() == query_1.x()) || ((ec_itr -> vertice_1).x() == query_2.x()));
b2 = (((ec_itr -> vertice_1).y() == query_1.y()) || ((ec_itr -> vertice_1).y() == query_2.y()));
b3 = (((ec_itr -> vertice_1).z() == query_1.z()) || ((ec_itr -> vertice_1).z() == query_2.z()));
b4 = (((ec_itr -> vertice_2).x() == query_1.x()) || ((ec_itr -> vertice_2).x() == query_2.x()));
b5 = (((ec_itr -> vertice_2).y() == query_1.y()) || ((ec_itr -> vertice_2).y() == query_2.y()));
b6 = (((ec_itr -> vertice_2).z() == query_1.z()) || ((ec_itr -> vertice_2).z() == query_2.z()));
if (b1 && b2 && b3 && b4 && b5 && b6)
{
return true;
}
}
return false;
}
bool SelectedTriangleSet::exist_vertice(const Point_tr query)
{
vector<Point_tr>::iterator v_itr;
for(v_itr = selected_vertice.begin(); v_itr != selected_vertice.end(); ++v_itr)
{
if (((v_itr -> x()) == query.x()) && ((v_itr -> y()) == query.y()) && ((v_itr -> z()) == query.z()))
{
return true;
}
}
return false;
}
void SelectedTriangleSet::add_vertice_to_set(const Point_tr v1, const Point_tr v2, const Point_tr v3)
{
SelectedTriangleElement e;
if (!(exist_vertice(v1)))
{
selected_vertice.push_back(v1);
}
if (!(exist_vertice(v2)))
{
selected_vertice.push_back(v2);
}
if (!(exist_vertice(v3)))
{
selected_vertice.push_back(v3);
}
e.vertice_index.push_back(get_vertice_index(v1));
e.vertice_index.push_back(get_vertice_index(v2));
e.vertice_index.push_back(get_vertice_index(v3));
triangle_set.push_back(e);
/**
if (!(exist_connection_edge(v1, v2)))
{
SelectedTriangleElementConnection ec(v1, v2);
connection_edge.push_back(ec);
}
if (!(exist_connection_edge(v1, v3)))
{
SelectedTriangleElementConnection ec(v1, v3);
connection_edge.push_back(ec);
}
if (!(exist_connection_edge(v2, v3)))
{
SelectedTriangleElementConnection ec(v2, v3);
connection_edge.push_back(ec);
}
**/
}
class SamplePointKdTreeNode
{
private:
SamplePoint * node;
SamplePointKdTreeNode * left;
SamplePointKdTreeNode * right;
SamplePointKdTreeNode() : node(NULL), left(NULL), right(NULL){}
void insert(SamplePointKdTreeNode * curr_node, SamplePoint * p, int dim=0);
SamplePoint * get_sample_point_by_coordinate(SamplePointKdTreeNode * curr_node, double x_1, double y_1, double z_i, int dim=0);
friend class SamplePointSet;
};
void SamplePointKdTreeNode::insert(SamplePointKdTreeNode * curr_node, SamplePoint * p, int dim)
{
SamplePointKdTreeNode * next_node = curr_node;
if (next_node -> node == NULL)
{
next_node -> node = p;
}
else
{
switch (dim)
{
case 0:
if (curr_node -> node -> sample_point_x() < p -> sample_point_x())
{
if (next_node -> left == NULL)
{
next_node -> left = new SamplePointKdTreeNode();
}
next_node = next_node -> left;
}
else
{
if (next_node -> right == NULL)
{
next_node -> right = new SamplePointKdTreeNode();
}
next_node = next_node -> right;
}
break;
case 1:
if (curr_node -> node -> sample_point_y() < p -> sample_point_y())
{
if (next_node -> left == NULL)
{
next_node -> left = new SamplePointKdTreeNode();
}
next_node = next_node -> left;
}
else
{
if (next_node -> right == NULL)
{
next_node -> right = new SamplePointKdTreeNode();
}
next_node = next_node -> right;
}
break;
case 2:
if (curr_node -> node -> sample_point_z() < p -> sample_point_z())
{
if (next_node -> left == NULL)
{
next_node -> left = new SamplePointKdTreeNode();
}
next_node = next_node -> left;
}
else
{
if (next_node -> right == NULL)
{
next_node -> right = new SamplePointKdTreeNode();
}
next_node = next_node -> right;
}
break;
}
insert(next_node, p, (dim + 1) % 3);
}
}
SamplePoint * SamplePointKdTreeNode::get_sample_point_by_coordinate(SamplePointKdTreeNode * curr_node, double x_i, double y_i, double z_i, int dim)
{
SamplePointKdTreeNode * next_node = curr_node;
if (next_node -> node == NULL)
{
cerr << "attempting to search for a sample point that does not exist" << endl;
system("pause");
exit(1);
}
else if ((next_node -> node -> sample_point_x() == x_i) &&
(next_node -> node -> sample_point_y() == y_i) &&
(next_node -> node -> sample_point_z() == z_i))
{
return next_node -> node;
}
else
{
switch (dim)
{
case 0:
if (curr_node -> node -> sample_point_x() < x_i)
{
next_node = next_node -> left;
}
else
{
next_node = next_node -> right;
}
break;
case 1:
if (curr_node -> node -> sample_point_y() < y_i)
{
next_node = next_node -> left;
}
else
{
next_node = next_node -> right;
}
break;
case 2:
if (curr_node -> node -> sample_point_z() < z_i)
{
next_node = next_node -> left;
}
else
{
next_node = next_node -> right;
}
break;
}
return get_sample_point_by_coordinate(next_node, x_i, y_i, z_i, (dim + 1) % 3);
}
}
class SamplePointSet
{
private:
list<Point_tr> sample_point_data;
list<SamplePoint> sample_point_list;
SelectedTriangleSet improved_triangulation;
Tree sample_point_tree_final;
Triangulations Tr;
SamplePointKdTreeNode * root;
void build_surface_normal();
double squared_distance(const Point p1, const Point p2) const;
void add_sample(const double px, const double py, const double pz);
void improve_triangulation();
void check_for_circumradius(list<finite_facets_itr>& temp);
void check_for_voronoi_edge(list<finite_facets_itr>& temp);
double max_of_three(const double d1, const double d2, const double d3) const;
void find_fifth_nearest_dist();
void search_nearest_neighbour(const int num_of_neighbour);
double implicit_function(list<SamplePoint * > pl, const Point query);
void create_triangulation_file(const char * const file) const;
void get_sample_point(const char * const file_source);
void compute_surface_normal(SamplePoint& sp, const double x_variance, const double y_variance, const double z_variance, const double x_y_covariance, const double y_z_covariance, const double x_z_covariance);
SamplePoint get_sample_point_from_point(Point p);
//SamplePoint * get_sample_point_by_coordinate(const double x_i, const double y_i, const double z_i);
//SamplePoint * get_sample_point_by_coordinate_operation(const double x_i, const double y_i, const double z_i, const int dim, Tree_itr ref);
public:
SamplePointSet();
void get_triangulation(const char * const file_source, const char * const file_result);
};
SamplePointSet::SamplePointSet()
{
root = new SamplePointKdTreeNode();
}
void SamplePointSet::create_triangulation_file(const char * const file_result) const
{
ofstream output_file;
output_file.open(file_result);
output_file << "OFF\n";
output_file << improved_triangulation.get_num_of_vertices() << ' ' << improved_triangulation.get_num_of_facets() << ' ' << improved_triangulation.get_num_of_edges() << '\n';
for(int i=0; i<improved_triangulation.selected_vertice.size(); ++i)
{
output_file << improved_triangulation.selected_vertice.at(i).x() << ' '
<< improved_triangulation.selected_vertice.at(i).y() << ' '
<< improved_triangulation.selected_vertice.at(i).z() << '\n';
}
for(int i = 0; i < improved_triangulation.triangle_set.size(); ++i)
{
output_file << "3 ";
for(int j = 0; j < improved_triangulation.triangle_set.at(i).vertice_index.size(); ++j)
{
output_file << improved_triangulation.triangle_set.at(i).vertice_index.at(j) << ' ';
}
output_file << '\n';
}
output_file.close();
}
void SamplePointSet::get_triangulation(const char * const file_source,const char * const file_result)
{
if( !mclInitializeApplication(NULL,0) )
{
fprintf(stderr, "Could not initialize the application.\n");
system("pause");
exit(1);
}
if (!find_surface_normalInitialize())
{
fprintf(stderr,"Could not initialize the library.\n");
system("pause");
exit(1);
}
cout << "getting sample point input ....." << endl;
get_sample_point(file_source);
list<SamplePoint>::iterator start = sample_point_list.begin();
list<SamplePoint>::iterator end = sample_point_list.end();
list<SamplePoint>::iterator curr;
Tr.insert(sample_point_data.begin(), sample_point_data.end());
sample_point_data.clear();
/**
cout << "perturbing sample points ..... " << endl;
for(curr = start; curr != end; ++curr)
{
curr -> perturb_point();
}
cout << "rebuilding search tree ....." << endl;
for(curr = start; curr != end; ++curr)
{
double tempX = curr -> sample_point_x();
double tempY = curr -> sample_point_y();
double tempZ = curr -> sample_point_z();
Point p(tempX, tempY, tempZ);
sample_point_tree_final.insert(p);
}
**/
cout << "initializing neighbour and distance ....." << endl;
search_nearest_neighbour(55);
find_fifth_nearest_dist();
//sample_point_tree_final.clear();
/**
for(curr = start; curr != end; ++curr)
{
cout << "here" << endl;
curr -> neighbour = search_nearest_neighbour(*curr, 55);
cout << "sdfsasdaa" << endl;
cout << ++i << endl;
curr -> set_fifth_nearest_dist(find_fifth_nearest_dist(*curr));
cout << "there" << endl;
}
**/
cout << "rebuilding search tree ....." << endl;
int i=0;
for(curr = start; curr != end; ++curr)
{ cout << ++i << "RST" << endl;
root -> insert(root, &(*curr));
}
/**
cout << "building row triangulation ....." << endl;
for(curr = start; curr != end; ++curr)
{
cout << "here" << endl;
double tempX = curr -> sample_point_x();
double tempY = curr -> sample_point_y();
double tempZ = curr -> sample_point_z();
Tr.insert(Point_tr(tempX, tempY, tempZ));
}
**/
cout << "computing surface normal ....." << endl;
build_surface_normal();
cout << "improving triangulation ..." << endl;
improve_triangulation();
cout << "creating off files ... " << endl;
create_triangulation_file(file_result);
cout << "Done" << endl;
find_surface_normalTerminate();
mclTerminateApplication();
}
void SamplePointSet::improve_triangulation()
{
list<finite_facets_itr> temp;
check_for_circumradius(temp);
check_for_voronoi_edge(temp);
/**
Triangle temp_triangle;
double temp_circumradius = 0;
double temp_max_inter_vertice_dist;
finite_facets_itr starting_facet = Tr.finite_facets_begin();
finite_facets_itr ending_facet = Tr.finite_facets_end();
finite_facets_itr current_facet = starting_facet;
int i=0;
for(current_facet = starting_facet; current_facet != ending_facet; ++current_facet)
{
cout << ++i << "CCR" << endl;
temp_triangle = Tr.triangle(*current_facet);
improved_triangulation.add_vertice_to_set(temp_triangle.vertex(0), temp_triangle.vertex(1), temp_triangle.vertex(2));
cout << improved_triangulation.selected_vertice.size() << "SSSSSSSSSSSSSSS" << endl;
}
**/
}
void SamplePointSet::check_for_circumradius(list<finite_facets_itr>& temp)
{
Triangle temp_triangle;
double temp_circumradius = 0;
double temp_max_inter_vertice_dist;
finite_facets_itr starting_facet = Tr.finite_facets_begin();
finite_facets_itr ending_facet = Tr.finite_facets_end();
finite_facets_itr current_facet = starting_facet;
double d1;
double d2;
double d3;
int i=0;
for(current_facet = starting_facet; current_facet != ending_facet; ++current_facet)
{
cout << ++i << "CCR" << endl;
temp_triangle = Tr.triangle(*current_facet);
/**
temp_max_inter_vertice_dist = max_of_three(find_fifth_nearest_dist(SamplePoint(temp_triangle.vertex(0).x(), temp_triangle.vertex(0).y(), temp_triangle.vertex(0).z())) * 3,
find_fifth_nearest_dist(SamplePoint(temp_triangle.vertex(1).x(), temp_triangle.vertex(0).y(), temp_triangle.vertex(0).z())) * 3,
find_fifth_nearest_dist(SamplePoint(temp_triangle.vertex(2).x(), temp_triangle.vertex(0).y(), temp_triangle.vertex(0).z())) * 3);
**/
/**
Fuzzy_sphere sp1(SamplePoint(temp_triangle.vertex(0).x(), temp_triangle.vertex(0).y(), temp_triangle.vertex(0).z()), 0, 0);
Fuzzy_sphere sp2(SamplePoint(temp_triangle.vertex(1).x(), temp_triangle.vertex(1).y(), temp_triangle.vertex(1).z()), 0, 0);
Fuzzy_sphere sp3(SamplePoint(temp_triangle.vertex(2).x(), temp_triangle.vertex(2).y(), temp_triangle.vertex(2).z()), 0, 0);
Tree::iterator temp_vertice1 = final_tree.search(Tree::iterator temp_vertice1, sp1);
std::ostream_iterator<SamplePoint> temp_vertice2 = final_tree.search(std::ostream_iterator<SamplePoint>(std::cout, "\n"), sp2);
std::ostream_iterator<SamplePoint> temp_vertice3 = final_tree.search(std::ostream_iterator<SamplePoint>(std::cout, "\n"), sp3);
temp_max_inter_vertice_dist = max_of_three(temp_vertice_1 -> fifth_nearest_dist,
temp_vertice_2 -> fifth_nearest_dist,
temp_vertice_3 -> fifth_nearest_dist);
**/
d1 = root -> get_sample_point_by_coordinate(root,
temp_triangle.vertex(0).x(),
temp_triangle.vertex(0).y(),
temp_triangle.vertex(0).z()) -> fifth_nearest_dist;
d2 = root -> get_sample_point_by_coordinate(root,
temp_triangle.vertex(1).x(),
temp_triangle.vertex(1).y(),
temp_triangle.vertex(1).z()) -> fifth_nearest_dist;
d3 = root -> get_sample_point_by_coordinate(root,
temp_triangle.vertex(2).x(),
temp_triangle.vertex(2).y(),
temp_triangle.vertex(2).z()) -> fifth_nearest_dist;
d1 *= 3;
d2 *= 3;
d3 *= 3;
temp_max_inter_vertice_dist = max_of_three(d1, d2, d3);
temp_circumradius = (squared_radius(temp_triangle.vertex(0), temp_triangle.vertex(1), temp_triangle.vertex(2)));
cout << "d1: " << d1 << endl;
cout << "d2: " << d2 << endl;
cout << "d3: " << d3 << endl;
cout << "temp_circumradius: " << temp_circumradius << endl;
cout << "temp_max_inter_vertice_dist: " << temp_max_inter_vertice_dist << endl;
cout << "***************" << temp.size() << endl;
if (temp_circumradius <= temp_max_inter_vertice_dist)
{
cout << "triangle filtered" << endl;
temp.push_back(current_facet);
/**
improved_triangulation.add_vertice_to_set(temp_triangle.vertex(0), temp_triangle.vertex(1), temp_triangle.vertex(2), facet_index);
++facet_index;
**/
//improved_triangulation.add_vertice_to_set(temp_triangle.vertex(0), temp_triangle.vertex(1), temp_triangle.vertex(2));
}
}
}
void SamplePointSet::check_for_voronoi_edge(list<finite_facets_itr>& temp)
{
cout << temp.size() << endl;
cout << "vornoi"<< endl;
system("pause");
Triangle temp_triangle;
Segment temp_dual_voronoi_edge_of_triangle;
CGAL::Object obj;
list<Point>::iterator pl_itr;
list<finite_facets_itr>::iterator starting_facet = temp.begin();
list<finite_facets_itr>::iterator ending_facet = temp.end();
list<finite_facets_itr>::iterator current_facet;
int i=0;
for(current_facet = starting_facet; current_facet != ending_facet; ++current_facet)
{
cout << ++i << "CFVE" << endl;
list<Point> temp_neighbour_list_1;
list<Point> temp_neighbour_list_2;
list<Point> temp_neighbour_list_3;
list<SamplePoint *> point_set;
temp_triangle = Tr.triangle(*(*current_facet));
/**
temp_neighbour_list_1 = (search_nearest_neighbour(SamplePoint(temp_triangle.vertex(0).x(), temp_triangle.vertex(0).y(), temp_triangle.vertex(0).z()), 15));
temp_neighbour_list_1.push_back(SamplePoint(temp_triangle.vertex(0).x(), temp_triangle.vertex(0).y(), temp_triangle.vertex(0).z()));
temp_neighbour_list_2 = (search_nearest_neighbour(SamplePoint(temp_triangle.vertex(1).x(), temp_triangle.vertex(1).y(), temp_triangle.vertex(1).z()), 15));
temp_neighbour_list_2.push_back(SamplePoint(temp_triangle.vertex(1).x(), temp_triangle.vertex(1).y(), temp_triangle.vertex(1).z()));
temp_neighbour_list_3 = (search_nearest_neighbour(SamplePoint(temp_triangle.vertex(2).x(), temp_triangle.vertex(2).y(), temp_triangle.vertex(2).z()), 15));
temp_neighbour_list_3.push_back(SamplePoint(temp_triangle.vertex(2).x(), temp_triangle.vertex(2).y(), temp_triangle.vertex(2).z()));
**/
/**
temp_neighbour_list_1 = final_tree.search(SamplePoint(temp_triangle.vertex(0).x(), temp_triangle.vertex(0).y(), temp_triangle.vertex(0).z()), 0);
//temp_neighbour_list_1.push_back(SamplePoint(temp_triangle.vertex(0).x(), temp_triangle.vertex(0).y(), temp_triangle.vertex(0).z()));
temp_neighbour_list_2 = final_tree.search(SamplePoint(temp_triangle.vertex(1).x(), temp_triangle.vertex(1).y(), temp_triangle.vertex(1).z()), 0);
//temp_neighbour_list_2.push_back(SamplePoint(temp_triangle.vertex(1).x(), temp_triangle.vertex(1).y(), temp_triangle.vertex(1).z()));
temp_neighbour_list_3 = final_tree.search(SamplePoint(temp_triangle.vertex(2).x(), temp_triangle.vertex(2).y(), temp_triangle.vertex(2).z()), 0);
//temp_neighbour_list_3.push_back(SamplePoint(temp_triangle.vertex(2).x(), temp_triangle.vertex(2).y(), temp_triangle.vertex(2).z()));
**/
temp_neighbour_list_1 = root -> get_sample_point_by_coordinate(root,
temp_triangle.vertex(0).x(),
temp_triangle.vertex(0).y(),
temp_triangle.vertex(0).z()) -> neighbour;
temp_neighbour_list_2 = root -> get_sample_point_by_coordinate(root,
temp_triangle.vertex(1).x(),
temp_triangle.vertex(1).y(),
temp_triangle.vertex(1).z()) -> neighbour;
temp_neighbour_list_3 = root -> get_sample_point_by_coordinate(root,
temp_triangle.vertex(2).x(),
temp_triangle.vertex(2).y(),
temp_triangle.vertex(2).z()) -> neighbour;
double tempX;
double tempY;
double tempZ;
for(pl_itr = temp_neighbour_list_1.begin(); pl_itr != temp_neighbour_list_1.end(); ++pl_itr)
{
tempX = pl_itr -> x();
tempY = pl_itr -> y();
tempZ = pl_itr -> z();
point_set.push_back( root -> get_sample_point_by_coordinate(root, tempX, tempY, tempZ));
//point_set.push_back(*pl_itr);
}
for(pl_itr = temp_neighbour_list_2.begin(); pl_itr != temp_neighbour_list_2.end(); ++pl_itr)
{
tempX = pl_itr -> x();
tempY = pl_itr -> y();
tempZ = pl_itr -> z();
point_set.push_back( root -> get_sample_point_by_coordinate(root, tempX, tempY, tempZ));
//point_set.push_back(*pl_itr);
}
for(pl_itr = temp_neighbour_list_3.begin(); pl_itr != temp_neighbour_list_3.end(); ++pl_itr)
{
tempX = pl_itr -> x();
tempY = pl_itr -> y();
tempZ = pl_itr -> z();
point_set.push_back( root -> get_sample_point_by_coordinate(root, tempX, tempY, tempZ));
//point_set.push_back(*pl_itr);
}
obj = Tr.dual(*(*current_facet));
double temp_start;
double temp_end;
if (assign(temp_dual_voronoi_edge_of_triangle, obj))
{
temp_start = implicit_function(point_set, Point(temp_dual_voronoi_edge_of_triangle.start().x(), temp_dual_voronoi_edge_of_triangle.start().y(), temp_dual_voronoi_edge_of_triangle.start().z()));
temp_end = implicit_function(point_set, Point(temp_dual_voronoi_edge_of_triangle.end().x(), temp_dual_voronoi_edge_of_triangle.end().y(), temp_dual_voronoi_edge_of_triangle.end().z()));
if (temp_start * temp_end < 0)
{
improved_triangulation.add_vertice_to_set(temp_triangle.vertex(0), temp_triangle.vertex(1), temp_triangle.vertex(2));
}
}
}
}
double SamplePointSet::implicit_function(list<SamplePoint *> pl, const Point query)
{
cout << "start IF" << endl;
double result = 0;
double result_temp_1 = 0;
double result_temp_2 = 0;
double weight_temp_1 = 0;
double weight_temp_2 = 0;
double weight_power_1 = 0;
double weight_power_2 = 0;
double fifth_nearest_dist_of_query= 0;
double total = 0;
/**
Neighbor_search search(sample_point_tree_final, query, 5);
Neighbor_search::iterator itr;
int i=0;
for(itr = search.begin(); itr != search.end(); ++itr)
{
++i;
if (i == 5)
{
break;
}
}
fifth_nearest_dist_of_query = itr -> second;
**/
//fifth_nearest_dist_of_query = find_fifth_nearest_dist(query);
list<SamplePoint* >::iterator p_itr;
//list<SamplePoint* > sp;
//list<SamplePoint>::iterator sp_itr;
/**
for(p_itr = pl.begin(); p_itr != pl.end(); ++p_itr)
{
sp.push_back(get_sample_point_from_point(*p_itr));
}
**/
for(p_itr = pl.begin(); p_itr != pl.end(); ++p_itr)
{
weight_power_1 = pow(( (*p_itr) -> sample_point_x() - query.x()), 2);
weight_power_1 += pow(( (*p_itr) -> sample_point_y() - query.y()), 2);
weight_power_1 += pow(( (*p_itr) -> sample_point_z() - query.z()), 2);
weight_power_1 = sqrt(weight_power_1);
weight_power_1 /= fifth_nearest_dist_of_query;
weight_power_1 = pow(weight_power_1, 2);
weight_power_1 *= -1;
weight_temp_1 += exp(weight_power_1);
}
for(p_itr = pl.begin(); p_itr != pl.end(); ++p_itr)
{
weight_power_2 = pow(( (*p_itr) -> sample_point_x() - query.x()), 2);
weight_power_2 += pow(( (*p_itr) -> sample_point_y() - query.y()), 2);
weight_power_2 += pow(( (*p_itr) -> sample_point_z() - query.z()), 2);
weight_power_2 = sqrt(weight_power_2);
weight_power_2 /= fifth_nearest_dist_of_query;
weight_power_2 = pow(weight_power_1, 2);
weight_power_2 *= -1;
weight_temp_2 = exp(weight_power_2);
result_temp_1 = weight_temp_2 / weight_temp_1;
result_temp_2 = ( (*p_itr) -> surface_normal_x()) * (query.x() - ( (*p_itr) -> sample_point_x()));
result_temp_2 += ( (*p_itr) -> surface_normal_y()) * (query.y() - ( (*p_itr) -> sample_point_y()));
result_temp_2 += ( (*p_itr) -> surface_normal_z()) * (query.z() - ( (*p_itr) -> sample_point_z()));
result += result_temp_1 * result_temp_2;
}
cout << "end_IF" << endl;
return result;
}
double SamplePointSet::max_of_three(const double d1, const double d2, const double d3) const
{
if ((d1 >= d2) && (d1 >= d3))
{
return d1;
}
else if ((d2 >= d1) && (d2 >= d3))
{
return d2;
}
else
{
return d3;
}
}
void SamplePointSet::get_sample_point(const char * const file_source)
{
ifstream inFile; //stream for reading sample data points from a file
inFile.open(file_source); //open the file storing the input sample points
if (!inFile) { //return erorr message and terminate the program iff the file containing the sample point data cannot be opened succssfully.
cerr << "Unable to open file containing the sample point data";
system("pause");
exit(1);
}
double tempI; //variable for holding the data currently read from the input file
double tempX = 0; //variable for string the x-copordinate of a sample point, which is read from the input file
double tempY = 0; //variable for string the y-copordinate of a sample point, which is read from the input file
double tempZ = 0; //variable for string the z-copordinate of a sample point, which is read from the input file
int tempCounter = 0; /**variable for recording wehther the program is current reading an x-coordinate, y-coordinate or z-coordinate from a file
If it has the value of 0, then it indicates that the porgramme is reading an x-coordinate
If it has the value of 1, then it indicates that the porgramme is reading an y-coordinate
If it has the value of 2, then it indicates that the porgramme is reading an z-coordinate
**/
int i=0;
while (!inFile.eof()) { //read the sample data points
inFile >> tempI;
if (tempCounter == 0)
{
tempX = tempI;
tempCounter = 1;
}
else if (tempCounter == 1)
{
tempY = tempI;
tempCounter = 2;
}
else if (tempCounter == 2)
{
cout << ++i << endl;
tempZ = tempI;
add_sample(tempX, tempY, tempZ);
tempCounter = 0;
}
}
inFile.close();
}
void SamplePointSet::add_sample(const double px, const double py, const double pz)
{
SamplePoint p(px, py, pz); //initalize a SamplePoint object for the inout point
sample_point_list.push_back(p); //add the SamplePoint object representing the sample point into the list of SamplePoint
Point_tr pi(px, py, pz);
sample_point_data.push_back(pi);
//Tr.insert(Point_tr(px, py, pz)); //add the sample point into the triangulation
sample_point_tree_final.insert(Point(px, py, pz));
}
double SamplePointSet::squared_distance(const Point p1, const Point p2) const
{
double result = 0;
result += p1.x() * p2.x();
result += p1.y() * p2.y();
result += p1.z() * p2.z();
return result;
}
void SamplePointSet::build_surface_normal()
{
cout << "BSN start" << endl;
vector<double> temp_neighbour_dist;
double temp_fifth_nearest_dist = 0;
list<Point> co_input; //input for covriance matrix() calculation
Point co_mean; //mean vector for covarience matrix() calculation
double x_total = 0; //mean vector calculation
double y_total = 0; //mean vector calculation
double z_total = 0; //mean vector calculation
double x_variance = 0;
double y_variance = 0;
double z_variance = 0;
double x_y_covariance = 0;
double y_z_covariance = 0;
double x_z_covariance = 0;
list<SamplePoint>::iterator sample_point_set_itr1;
list<SamplePoint *>::iterator sample_point_set_itr2;
list<Point>::iterator co_itr;
list<Point> neighbour_temp;
int i=0;
for(sample_point_set_itr1 = sample_point_list.begin(); sample_point_set_itr1 != sample_point_list.end(); ++sample_point_set_itr1)
{
cout << ++i << "CSN " << endl;
temp_fifth_nearest_dist = 0;
list<Point> co_input; //input for covriance matrix() calculation
x_total = 0; //mean vector calculation
y_total = 0; //mean vector calculation
z_total = 0; //mean vector calculation
x_variance = 0;
y_variance = 0;
z_variance = 0;
x_y_covariance = 0;
y_z_covariance = 0;
x_z_covariance = 0;
/**
temp_fifth_nearest_dist = find_fifth_nearest_dist(SamplePoint(sample_point_set_itr1 -> sample_point_x(),
sample_point_set_itr1 -> sample_point_y(),
sample_point_set_itr1 -> sample_point_z()));
neighbour = search_nearest_neighbour(*sample_point_set_itr1, 8);
**/
/**
neighbour = root -> get_sample_point_by_coordinate(root,
sample_point_set_itr1 -> sample_point_x(),
sample_point_set_itr1 -> sample_point_y(),
sample_point_set_itr1 -> sample_point_z()) -> neighbour;
**/
temp_fifth_nearest_dist = sample_point_set_itr1 -> fifth_nearest_dist;
neighbour_temp = sample_point_set_itr1 -> neighbour;
list<SamplePoint *> neighbour;
list<Point>::iterator itr;
for (itr = neighbour_temp.begin(); itr != neighbour_temp.end(); ++itr)
{
neighbour.push_back(root -> get_sample_point_by_coordinate(root, itr -> x(),itr -> y(), itr -> z()));
}
for(sample_point_set_itr2 = neighbour.begin(); sample_point_set_itr2 != neighbour.end(); ++sample_point_set_itr2)
{
double dist = pow((sample_point_set_itr1 -> sample_point_x() - (*sample_point_set_itr2) -> sample_point_x()), 2);
dist += pow((sample_point_set_itr1 -> sample_point_y() - (*sample_point_set_itr2) -> sample_point_y()), 2);
dist += pow((sample_point_set_itr1 -> sample_point_z() - (*sample_point_set_itr2) -> sample_point_z()), 2);
dist = sqrt(dist);
if((dist <= temp_fifth_nearest_dist * 3) && (dist != 0))
{
double diffX = ((*sample_point_set_itr2) -> sample_point_x()) - (sample_point_set_itr1 -> sample_point_x());
double diffY = ((*sample_point_set_itr2) -> sample_point_y()) - (sample_point_set_itr1 -> sample_point_y());
double diffZ = ((*sample_point_set_itr2) -> sample_point_z()) - (sample_point_set_itr1 -> sample_point_z());
co_input.push_back(Point(diffX, diffY, diffZ));
}
}
for(co_itr = co_input.begin(); co_itr != co_input.end(); ++co_itr)
{
x_total += co_itr -> x();
y_total += co_itr -> y();
z_total += co_itr -> z();
}
co_mean = Point((x_total/co_input.size()), (y_total/co_input.size()), (z_total/co_input.size()));
for(co_itr = co_input.begin(); co_itr != co_input.end(); ++co_itr)
{
x_variance += ((co_itr -> x()) - (co_mean.x())) * ((co_itr -> x()) - (co_mean.x()));
y_variance += ((co_itr -> y()) - (co_mean.y())) * ((co_itr -> y()) - (co_mean.y()));
z_variance += ((co_itr -> z()) - (co_mean.z())) * ((co_itr -> z()) - (co_mean.z()));
x_y_covariance += ((co_itr -> x()) - (co_mean.x())) * ((co_itr -> y()) - (co_mean.y()));
y_z_covariance += ((co_itr -> y()) - (co_mean.y())) * ((co_itr -> z()) - (co_mean.z()));
x_z_covariance += ((co_itr -> x()) - (co_mean.x())) * ((co_itr -> z()) - (co_mean.z()));
}
x_variance /= (co_input.size() -1);
y_variance /= (co_input.size() -1);
z_variance /= (co_input.size() -1);
x_y_covariance /= (co_input.size() -1);
y_z_covariance /= (co_input.size() -1);
x_z_covariance /= (co_input.size() -1);
sample_point_set_itr1 -> set_surface_normal(0,0,0);
compute_surface_normal(*(sample_point_set_itr1), x_variance, y_variance, z_variance, x_y_covariance, y_z_covariance, x_z_covariance);
}
cout << "BSN end" << endl;
}
void SamplePointSet::compute_surface_normal(SamplePoint& sp, const double x_variance, const double y_variance, const double z_variance, const double x_y_covariance, const double y_z_covariance, const double x_z_covariance)
{
mxArray * input[9];
mxArray * output = NULL;
double * result;
double normal_component_1;
double normal_component_2;
double normal_component_3;
input[0] = mxCreateDoubleScalar(x_variance);
input[1] = mxCreateDoubleScalar(x_y_covariance);
input[2] = mxCreateDoubleScalar(x_z_covariance);
input[3] = mxCreateDoubleScalar(x_y_covariance);
input[4] = mxCreateDoubleScalar(y_variance);
input[5] = mxCreateDoubleScalar(y_z_covariance);
input[6] = mxCreateDoubleScalar(x_z_covariance);
input[7] = mxCreateDoubleScalar(y_z_covariance);
input[8] = mxCreateDoubleScalar(z_variance);
mlfFind_surface_normal(1,&output,input[0], input[1], input[2], input[3], input[4], input[5], input[6], input[7], input[8]);
result = mxGetPr(output);
normal_component_1 = *(result);
normal_component_2 = *(result + 1);
normal_component_3 = *(result + 2);
sp.set_surface_normal(normal_component_1, normal_component_2, normal_component_3);
mxDestroyArray(input[0]);
mxDestroyArray(input[1]);
mxDestroyArray(input[2]);
mxDestroyArray(input[3]);
mxDestroyArray(input[4]);
mxDestroyArray(input[5]);
mxDestroyArray(input[6]);
mxDestroyArray(input[7]);
mxDestroyArray(input[8]);
mxDestroyArray(output);
}
void SamplePointSet::search_nearest_neighbour(const int num_of_neighbour)
{
cout << "SNN start" << endl;
list<SamplePoint>::iterator start = sample_point_list.begin();
list<SamplePoint>::iterator end = sample_point_list.end();
list<SamplePoint>::iterator curr;
Neighbor_search * search;
Point * query;
Neighbor_search::iterator itr;
int i=0;
for(curr = start; curr != end; ++curr)
{
cout << ++i << endl;
query = new Point(curr -> sample_point_x(), curr -> sample_point_y(), curr -> sample_point_z());
search = new Neighbor_search(sample_point_tree_final, *query, num_of_neighbour);
for (itr = search -> begin(); itr != search -> end(); ++itr)
{
(curr -> neighbour).push_back(itr -> first);
}
delete query;
delete search;
query = NULL;
search = NULL;
}
/**
Neighbor_search search(sample_point_tree_final, query_point, num_of_neighbour);
list<Point> result_neighbour;
Neighbor_search::iterator itr = search.begin();
int i=0;
for(i=0, itr = search.begin(); itr != search.end(); ++itr)
{
if (i == num_of_neighbour)
{
break;
}
result_neighbour.push_back(itr->first);
++i;
}
**/
cout << "SNN end" << endl;
//return result_neighbour;
}
void SamplePointSet::find_fifth_nearest_dist()
{
list<SamplePoint>::iterator start = sample_point_list.begin();
list<SamplePoint>::iterator end = sample_point_list.end();
list<SamplePoint>::iterator curr;
Neighbor_search * search;
Point * query;
Neighbor_search::iterator itr;
int i=0;
for(curr = start; curr != end; ++curr)
{
query = new Point(curr -> sample_point_x(), curr -> sample_point_y(), curr -> sample_point_z());
search = new Neighbor_search(sample_point_tree_final, *query, 6);
i=0;
for (itr = search -> begin(); itr != search -> end(); ++itr)
{
++i;
if (i == 5)
{
break;
}
}
(curr -> fifth_nearest_dist) = itr -> second;
delete query;
delete search;
query = NULL;
search = NULL;
}
/**
Neighbor_search search(sample_point_tree_final, p, 5);
Neighbor_search::iterator itr = search.begin();
int i=0;
for(itr = search.begin(); itr != search.end(); ++itr)
{
++i;
if(i == 5)
{
break;
}
}
return itr -> second;
**/
}
/********************************************************SamplePoint.h******************************************************/
int main()
{
SamplePointSet sample;
sample.get_triangulation("C:\\bunny.txt", "C:\\result3.off");
system("pause");
return 0;
}
- [cgal-discuss] advice about CGAL Delaunay Triangulation, Kwok Jasper, 02/02/2009
- Re: [cgal-discuss] advice about CGAL Delaunay Triangulation, Andreas Fabri, 02/02/2009
- RE: [cgal-discuss] advice about CGAL Delaunay Triangulation, Kwok Jasper, 02/02/2009
- Re: [cgal-discuss] advice about CGAL Delaunay Triangulation, Andreas Fabri, 02/02/2009
- RE: [cgal-discuss] advice about CGAL Delaunay Triangulation, Kwok Jasper, 02/02/2009
- RE: [cgal-discuss] advice about CGAL Delaunay Triangulation, Kwok Jasper, 02/02/2009
- Re: [cgal-discuss] advice about CGAL Delaunay Triangulation, Andreas Fabri, 02/02/2009
- RE: [cgal-discuss] advice about CGAL Delaunay Triangulation, Kwok Jasper, 02/03/2009
- RE: [cgal-discuss] advice about CGAL Delaunay Triangulation, Kwok Jasper, 02/02/2009
- Re: [cgal-discuss] advice about CGAL Delaunay Triangulation, Andreas Fabri, 02/02/2009
- RE: [cgal-discuss] advice about CGAL Delaunay Triangulation, Kwok Jasper, 02/02/2009
- Re: [cgal-discuss] advice about CGAL Delaunay Triangulation, Andreas Fabri, 02/02/2009
Archive powered by MHonArc 2.6.16.