CGAL users discussion list

Hi,
I tried to follow the suggestion of Laurent:
<Here you do not test the equality between opposite_vertex_index and V[vit]
<(anyway, I just said they can only be equal by chance). Instead, you do
<something similar to:
< opposite_vertex_index = V[vit];
 <if(opposite_vertex_index != 0) { ....}
I replaced in the attached code the condition:  if(opposite_vertex_index=
V[vit])  with
if(vit->point()==cell->vertex(opposite_vertex_index)->point()). With this
condition for each vertex in the iteration of Finite_vertices_iterator I
want select the intersected facets,in the list of intersected primitives,
which have this vertex as opposite vertex.
 (I recall that I combined a 3D mesh generation procedure with an AABB
tree. My query object
is a line with fixed direction and crossing a vertex of the mesh (I
consider parallel rays crossing each vertex for different intersection
queries) and I use a kind of AABB_primitive, such as  Id is
Tr::Facet_finite_iterator)
I would expect to find just two cells incident to the vertex and
containing the intersected facets and sort along the direction of the
intersection line query the two intersection points lying in the facets
opposite to the vertex.
This should allow me to define an intersection point "before" and an
intersection point "after" the vertex of interest for each step.

Conversely I detect almost everywhere just one intersection and not 2 as
expected.Two intersection are very seldom detected.

1) So my criterion is not sufficient to choose the two opposite facet? Or
I make other mistakes in applying it?

2) Moreover, if I print the vertexes when the condition
if(vit->point()==cell->vertex(opposite_vertex_index)->point())
is fulfilled I retrieve the point coordinates but also an integer. Why its
value is always 0?
Below the code.....

Thanks for your attention!


        for( Finite_vertices_iterator vit =
c3t3.triangulation().finite_vertices_begin();
                vit != c3t3.triangulation().finite_vertices_end();
                ++vit)
        {
                Point c = vit->point();
                Line line_vertex(c,dir);
                std::cout << 
tree.number_of_intersected_primitives(line_vertex)
                << " intersected primitives"<<std::endl;
                // computes all intersections with line query
                std::list<Object_and_primitive_id> intersections;
                tree.all_intersections(line_vertex, 
std::back_inserter(intersections));
std::list<Object_and_primitive_id>::iterator it_intersections=
intersections.begin();
                std::list<Primitive_id> primitives;
                std::vector<std::pair<Point,Id> > myvector;
                std::vector<std::pair<Point,Id> >::iterator it;
       tree.all_intersected_primitives(line_vertex,
std::back_inserter(primitives));
                for(; it_intersections != intersections.end() ; 
++it_intersections)

                {
                        // get intersection object
                        Object_and_primitive_id op = *it_intersections;
                        CGAL::Object object = op.first;
                        Id id=op.second;


                        Point point;
                        Segment segment;
                        if(CGAL::assign(point,object))
                        {
                                it=myvector.begin();
                                myvector.insert(it, make_pair(point,id));
                        }
                        else if(CGAL::assign(segment,object))
                                std::cout << "intersection object is a 
segment" <<std::endl;
                }
                std::sort(myvector.begin(), myvector.end(), sort_pred());

for ( it=myvector.begin() ; it != myvector.end(); it++ ){
                        C3t3::Cell_handle cell    = (*it).second->first;
                        int opposite_vertex_index = (*it).second->second;
const Vertex_handle& vopp=cell->vertex(opposite_vertex_index);
//if(opposite_vertex_index==V[vit]) {OLD TEST
if(vit->point()==cell->vertex(opposite_vertex_index)->point()) {
std::cout<<vit->point()<< " "<<" "<<" "<<"
"<<cell->vertex(opposite_vertex_index)->point()<<std::endl;
if(c3t3.is_in_complex(cell)==true) {
                                std::cout << " intersection/coordinates"<<"
"<<std::setprecision(15)<<(*it).first << " => " << std::endl;
for(int i = 0; i < 4; i++){
const Vertex_handle& vh = (*((*it).second)).first->vertex(i);
                                                                std::cout<<" 
vertex index:"<<" " <<V[vh] << " "<<std::endl;
                                        std::cout <<" vertex coordinates:"<<" 
"<<vh->point().x()<<"
"<<vh->point().y()<<" "<<vh->point().z()<<std::endl;
}
}
}
}








On lundi 23 mai 2011 16:01:09 

 wrote:

1) Is my criterion sufficient to choose the two opposite facet? Or I have
to select all the cells incident to the vertex before setting the
condition on the "opposite vertex_index"? (In this case I have an error in
the first argument of the function
c3t3.triangulation().incident_cells(vh,std::back_inserter(cells)), where I
set std::vector<Cell_handle> cells. I see that vh=Vertex_handle but I
don't know how set vh to the vertex pointed by the finite_vertex iterator
I'm using).

I do not know. I had a look at your code but I do not know how you fill the
vector myvector. Maybe the error is here. Actually I found other errors (see
below)...

2) How  can I compare the two intersection points lying on these facets in
order to establish the point on the line "before" and "after" the vertex
under study? Is a sort function a good choice or I have to evaluate the
sign of the scalar product between the direction of the oriented line and
a vector connecting the two points?

The scalar product is the best way to sort point along a direction.

3) Moreover I'm setting the condition c3t3.is_in_complex(cell) to be
assured of getting facets belonging the 3D complex. Does it correspond to
get facets both inside the volume and on the surface or surface facets are
excluded?

"setting"? Did you mean "testing"? c3t3.is_in_complex(cell) returns true if
the cell is in the volume (3D complex). The facets of a cell that is in the
volume can either be internal to the volume or be on the surface (2D
complex).

I'm sending just the part where my criterion is used. I hope it can be
understood:

I found errors... See below.

....DEFINE CLASS AABB PRIMITIVES (Tr::Finite_facets_iterator Id)
....DEFINE TYPES

int main()
{....MESH A SPHERE AND SET THE TREE
Point a(0.0, 0.0, 0.0);
Point b(0.0, 1.0, 0.0);
Ray ray_query(a,b);
Direction dir(ray_query);
std::map<Vertex_handle, int> V;
int inum = 0;
for( Finite_vertices_iterator vit =
c3t3.triangulation().finite_vertices_begin();
       vit != c3t3.triangulation().finite_vertices_end();
       ++vit)
{
       V[vit] = inum++;
}

Here you create an indexation of finite vertices of the triangulation, in the
order of iteration using finite_vertices_(begin|end)()... The numeration is
only for you. Internally the triangulation data structure does not keep such
an indexation (just handles or iterators).

for( Finite_vertices_iterator vit =
c3t3.triangulation().finite_vertices_begin();
       vit != c3t3.triangulation().finite_vertices_end();
       ++vit)
{
      Point c = vit->point();
       Line line_vertex(c,dir);
       ......COMPUTE INTERSECTIONS AND STORE A COUPLE
std::vector<std::pair<Point,Id> >myvector

       for ( it=myvector.begin() ; it != myvector.end(); it++ ){
               C3t3::Cell_handle cell    = (*it).second->first;
               int opposite_vertex_index = (*it).second->second;

Here opposite_vertex_index is an integer i, where:
 0 <= i < 4
that corresponds to the number of the facet in the cell. It means that the
vertex that is opposite to the facet (*it) is cell-
vertex(opposite_vertex_index). Note that this is very different to
V[opposite_vertex_index]!!

          if(opposite_vertex_index= V[vit]) {

Here you do not test the equality between opposite_vertex_index and V[vit]
(anyway, I just said they can only be equal by chance). Instead, you do
something similar to:
 opposite_vertex_index = V[vit];
 if(opposite_vertex_index != 0) { ....}

According to the indexation you filled in the map V, the condition will only
be false it vit==finite_vertices_begin(). That may explains why you got so
many vertices...

-- 
Laurent Rineau, PhD
R&D Engineer at GeometryFactory           http://www.geometryfactory.com/
Release Manager of the CGAL Project       http://www.cgal.org/

-- 
You are currently subscribed to cgal-discuss.
To unsubscribe or access the archives, go to
https://lists-sop.inria.fr/wws/info/cgal-discuss