Subject: CGAL users discussion list
List archive
- From: Jonas Schild <>
- To:
- Subject: Re: [cgal-discuss] Normal generation out of 3d alpha shapes
- Date: Wed, 21 May 2008 18:34:24 +0200
Hi again,
with all your help I could gladly solve the issue. I first tried Santosh's code which worked. Though this was quite inconvenient, by writing and reading streams, it brought me to looking how the alpha shape is put into the ostream in alpha_shape_3.h. There I could find the same solution Monique and Mariette kindly proposed. So thank you all!
For anyone interested here is the code:
int number_of_vertices = 0;
std::vector<CgalPoint3> vertices_list;
std::map<CgalVertexHandle,int> vertices_map;
CgalVerticesIterator vit = mAlphaShape->alpha_shape_vertices_begin();
for( ; vit != mAlphaShape->alpha_shape_vertices_end(); ++vit)
{
vertices_list.push_back((*vit)->point());
vertices_map[*vit] = number_of_vertices++;
}
CgalCellHandle c;
CgalFacetsIterator fit = mAlphaShape->alpha_shape_facets_begin();
std::vector<osg::Vec3> triangle(3);
CgalPoint3 v1, v2, v3;
int i, i0, i1, i2, j;
for( ; fit != mAlphaShape->alpha_shape_facets_end(); ++fit)
{
c = fit->first;
i = fit->second;
// the following ensures that regular facets are output
// in ccw order
if ((mAlphaShape->classify(*fit) == CgalAlphaShape3::REGULAR) &&
(mAlphaShape->classify(c) == CgalAlphaShape3::INTERIOR))
{
c = c->neighbor(i);
i = c->index(fit->first);
}
i0 = CgalAlphaShape3::Triangulation_utils_3::vertex_triple_index(i,0);
i1 = CgalAlphaShape3::Triangulation_utils_3::vertex_triple_index(i,1);
i2 = CgalAlphaShape3::Triangulation_utils_3::vertex_triple_index(i,2);
v1 = vertices_list[vertices_map[c->vertex(i0)]];
v2 = vertices_list[vertices_map[c->vertex(i1)]];
v3 = vertices_list[vertices_map[c->vertex(i2)]];
triangle[0] = osg::Vec3(v1[0], v1[1], v1[2]);
triangle[1] = osg::Vec3(v2[0], v2[1], v2[2]);
triangle[2] = osg::Vec3(v3[0], v3[1], v3[2]);
osg::Vec3 a = triangle[1] - triangle[0];
osg::Vec3 b = triangle[2] - triangle[1];
osg::Vec3 n = a ^ b;
n.normalize();
}
Best regards
Jonas
Santosh Tiwari wrote:
Attached with this email is the complete code.
I have modified the code and removed unnecessary stuff.
I have not attempted to compile it after the modification. But, you can
use the code as a blueprint to get the correct normals.
Hope that works.
Santosh Tiwari wrote:
This is how I got the correct normals for all the boundary facets.Hi Santosh,
I tried using your code, but failed with correctly setting up the Facet
and Vertex data types, so I was not able to create the lists correctly.
I then only used your normals to my vertices but they don't seem to be
correct. My vertices seem to be in a different order. Could you help me
out with your typedefs?
Thanks for kindly providing your code!
Jonas
The relevant code.--
Alpha_shape_3 *as;
list<Facet> facets;
list<Vertex> vertices;
ostringstream facetStream;
as->get_alpha_shape_facets(back_inserter(facets),
Alpha_shape_3::REGULAR);
as->get_alpha_shape_vertices(back_inserter(vertices),
Alpha_shape_3::REGULAR);
cout<<"\n Number of boundary facets = "<<facets.size();
cout<<"\n Number of vertices = "<<vertices.size();
facetStream<<vertices.size()<<endl;
facetStream<<facets.size()<<endl;
facetStream<<*as;
//You can now read the facetStream. Here is some code to read it.
string facetString(facetstream.str());
istringstream in;
in.str(facetString);
in>>numVertices;
in>>numFacets;
vertices = new double*[numVertices];
for (i=0; i<numVertices; i++) {
vertices[i] = new double[3];
in>>x>>y>>z;
vertices[i][0] = x;
vertices[i][1] = y;
vertices[i][2] = z;
}
facets = new unsigned long int*[numFacets];
for (i=0; i<numFacets; i++) {
facets[i] = new unsigned long int[3];
in>>v1>>v2>>v3;
facets[i][0] = v1;
facets[i][1] = v2;
facets[i][2] = v3;
}
for(i=0; i<numFacets; i++) {
x1 = vertices[facets[i][0]][0];
y1 = vertices[facets[i][0]][1];
z1 = vertices[facets[i][0]][2];
x2 = vertices[facets[i][1]][0];
y2 = vertices[facets[i][1]][1];
z2 = vertices[facets[i][1]][2];
x3 = vertices[facets[i][2]][0];
y3 = vertices[facets[i][2]][1];
z3 = vertices[facets[i][2]][2];
}
//You now have all the facets. You can write them to a STL file.
//You can compute the normals for every facet now.
//Here is some code to compute the normals.
a1 = x3 - x2;
a2 = y3 - y2;
a3 = z3 - z2;
b1 = x2 - x1;
b2 = y2 - y1;
b3 = z2 - z1;
n1 = a2*b3 - a3*b2;
n2 = a1*b3 - a3*b1;
n3 = a1*b2 - a2*b1;
res = sqrt(n1*n1 + n2*n2 + n3*n3);
n1 /= res;
n2 /= res;
n3 /= res;
I personally do not prefer accessing the internal data structure of a
library even if it permits that. It is very difficult to debug unless
you
know the data structure really well.
__
Santosh Tiwari
EIB-326 Clemson University
http://www.clemson.edu/~stiwari/
-----Original Message-----
From:
[mailto:]
Sent: Sunday, May 18, 2008 10:54 AM
To:
Subject: Re: [cgal-discuss] Normal generation out of 3d alpha shapes
Andreas Fabri wrote:
Jonas Schild wrote:http://www.cgal.org/Manual/3.3/doc_html/cgal_manual/TriangulationDS_3_ref/Cl
Hi all,Hi Jonas,
I have the following problem:
I'm constructing a 3d alpha shape from a set of points using CGAL.
After creation of the alpha shape, I'm iterating over all facets to
extract the triangles. But:
The orientation of the normal (either by using cross product or using
CGAL's ortho_vector(point, point, point) of which results are the
same) is partially wrong. I've read that the orientation of the
triangle is based on orientation of the cell (i.e. pointing into the
cell). However, I'm not able to retrieve the normals pointing outside
the alpha shape. Part of the problem is, that for certain alpha
values, the constructed shape is no solid object. Is it possible to
iterate only over alphas where the shape is solid, not necessarily
convex?
Furthermore, I don't understand how ccw(int) on the triangulation
structure is supposed to work regarding finding the correct
orientation thus the normals are pointing outside the volume.
Here comes a partial answer.
With ccw(int) you refer probably to the base class of all 3D
traingulations:
ass_Triangulation_utils_3.html#Cross_link_anchor_994
As the manual says ccw(int) only makes sense when it is 2D.really? I would say it is more like a matter of taste...
What you should use is the undocumented funcyion:
static int Triangulation_utils_3<..>::vertex_triple_index(const int i,
const int j)
{
// indexes of the jth vertex of the facet of a cell
// opposite to vertx i
CGAL_triangulation_precondition( ( i >= 0 && i < 4 ) &&
( j >= 0 && j < 3 ) );
return tab_vertex_triple_index[i][j];
}
The fact that it is undocumented is a bug.
Monique
You are currently subscribed to cgal-discuss.
To unsubscribe or access the archives, go to
https://lists-sop.inria.fr/wws/info/cgal-discuss
------------------------------------------------------------------------
#include <CGAL/Exact_predicates_inexact_constructions_kernel.h>
#include <CGAL/Delaunay_triangulation_3.h>
#include <CGAL/Triangulation_hierarchy_3.h>
#include <CGAL/Alpha_shape_3.h>
#include <CGAL/IO/io.h>
#include <iostream>
#include <sstream>
#include <string>
#include <fstream>
#include <list>
using namespace std;
struct K : CGAL::Exact_predicates_inexact_constructions_kernel {};
typedef CGAL::Alpha_shape_vertex_base_3<K>
Vb;
typedef CGAL::Triangulation_hierarchy_vertex_base_3<Vb> Vbh;
typedef CGAL::Alpha_shape_cell_base_3<K> Fb;
typedef CGAL::Triangulation_data_structure_3<Vbh,Fb> Tds;
typedef CGAL::Delaunay_triangulation_3<K,Tds>
Delaunay;
typedef CGAL::Triangulation_hierarchy_3<Delaunay>
Delaunay_hierarchy;
typedef CGAL::Alpha_shape_3<Delaunay_hierarchy>
Alpha_shape_3;
typedef K::Point_3
Point;
typedef Alpha_shape_3::Alpha_iterator
Alpha_iterator;
typedef Alpha_shape_3::NT
NT;
typedef Alpha_shape_3::Facet
Facet;
typedef Alpha_shape_3::Vertex_handle
Vertex;
void writeData(string& facetString);
int main(int argc, const char* const argv[]) {
if(argc!=2) {
cerr<<"\n Usage: "<<argv[0]<<" point_cloud.dat"<<endl;
exit(1);
}
int mode;
mode = static_cast<int>(atoi(argv[1]));
// Read the points and compute delaunay triangulation
Delaunay_hierarchy dt;
ifstream is(argv[1]);
if(!is) {
cerr<<"\n Could not open file "<<argv[1]<<" for reading, exiting
\n";
exit(1);
}
register unsigned long int n;
is>>n;
Point p;
cout<<"\n "<<n<<" points in the file";
for(register unsigned long int i=0; i<n; i++) {
is>>p;
dt.insert(p);
if(i%10000==0) {
cout<<"\n "<<i<<" points read and added to delaunay
triangulation";
}
}
cout<<"\n "<<n<<" points read and added to delaunay triangulation";
cout<<"\n Delaunay computed";
// compute alpha shape
Alpha_shape_3 *as = new Alpha_shape_3(dt);
cout<<"\n Alpha shape computed";
// find optimal alpha values
Alpha_shape_3::NT alpha_solid = as->find_alpha_solid();
Alpha_iterator opt = as->find_optimal_alpha(1);
cout<<"\n Smallest alpha value to get a solid through data points is "
<<alpha_solid;
cout<<"\n Optimal alpha value to get one connected component is
"<<*opt;
as->set_alpha(*opt);
assert(as->number_of_solid_components() == 1);
//write the facets
list<Facet> facets;
list<Vertex> vertices;
as->get_alpha_shape_facets(back_inserter(facets),
Alpha_shape_3::REGULAR);
as->get_alpha_shape_vertices(back_inserter(vertices),
Alpha_shape_3::REGULAR);
cout<<"\n Number of boundary facets = "<<facets.size();
cout<<"\n Number of vertices = "<<vertices.size();
ostringstream facetStream;
facetStream<<vertices.size()<<endl;
facetStream<<facets.size()<<endl;
facetStream<<*as;
writeData(facetStream.str());
cout<<endl;
return 0;
}
void writeData(string& facetString) {
istringstream in;
double** vertices;
unsigned long int** facets;
register unsigned long int i;
double x, y, z;
unsigned long int v1, v2, v3;
double n1, n2, n3, x1, y1, z1, x2, y2, z2, x3, y3, z3;
double a1, a2, a3, b1, b2, b3, res;
register unsigned long int numVertices;
register unsigned long int numFacets;
in.str(facetString);
in>>numVertices;
in>>numFacets;
vertices = new double*[numVertices];
for (i=0; i<numVertices; i++) {
vertices[i] = new double[3];
in>>x>>y>>z;
vertices[i][0] = x;
vertices[i][1] = y;
vertices[i][2] = z;
//cout<<vertices[i][0]<<", "<<vertices[i][1]<<",
"<<vertices[i][2]<<endl;
}
facets = new unsigned long int*[numFacets];
for (i=0; i<numFacets; i++) {
facets[i] = new unsigned long int[3];
in>>v1>>v2>>v3;
facets[i][0] = v1;
facets[i][1] = v2;
facets[i][2] = v3;
}
for(i=0; i<numFacets; i++) {
x1 = vertices[facets[i][0]][0];
y1 = vertices[facets[i][0]][1];
z1 = vertices[facets[i][0]][2];
x2 = vertices[facets[i][1]][0];
y2 = vertices[facets[i][1]][1];
z2 = vertices[facets[i][1]][2];
x3 = vertices[facets[i][2]][0];
y3 = vertices[facets[i][2]][1];
z3 = vertices[facets[i][2]][2];
a1 = x2 - x1;
a2 = y2 - y1;
a3 = z2 - z1;
b1 = x3 - x2;
b2 = y3 - y2;
b3 = z3 - z2;
n1 = a2*b3 - a3*b2;
n2 = a3*b1 - a1*b3;
n3 = a1*b2 - a2*b1;
res = sqrt(n1*n1 + n2*n2 + n3*n3);
n1 = n1/res;
n2 = n2/res;
n3 = n3/res;
//Do whatever you want with facets and normals.
}
return;
}
- Normal generation out of 3d alpha shapes, Jonas Schild, 05/16/2008
- Re: [cgal-discuss] Normal generation out of 3d alpha shapes, Andreas Fabri, 05/16/2008
- Re: [cgal-discuss] Normal generation out of 3d alpha shapes, Monique . Teillaud, 05/18/2008
- RE: [cgal-discuss] Normal generation out of 3d alpha shapes, Santosh Tiwari, 05/18/2008
- Re: [cgal-discuss] Normal generation out of 3d alpha shapes, Jonas Schild, 05/20/2008
- Re: [cgal-discuss] Normal generation out of 3d alpha shapes, Santosh Tiwari, 05/20/2008
- Re: [cgal-discuss] Normal generation out of 3d alpha shapes, Jonas Schild, 05/21/2008
- Re: [cgal-discuss] Normal generation out of 3d alpha shapes, Santosh Tiwari, 05/20/2008
- Re: [cgal-discuss] Normal generation out of 3d alpha shapes, Jonas Schild, 05/20/2008
- RE: [cgal-discuss] Normal generation out of 3d alpha shapes, Santosh Tiwari, 05/18/2008
- Re: [cgal-discuss] Normal generation out of 3d alpha shapes, Monique . Teillaud, 05/18/2008
- Re: [cgal-discuss] Normal generation out of 3d alpha shapes, Monique . Teillaud, 05/18/2008
- Re: [cgal-discuss] Normal generation out of 3d alpha shapes, Jonas Schild, 05/20/2008
- Re: [cgal-discuss] Normal generation out of 3d alpha shapes, Mariette Yvinec, 05/20/2008
- Re: [cgal-discuss] Normal generation out of 3d alpha shapes, Andreas Fabri, 05/16/2008
Archive powered by MHonArc 2.6.16.