CGAL users discussion list

[Farzaneh <>]

Hi,

In my program I have a hierarchical delaunay triangulation (representing a
terrian) and I have some query points. I should find the face that this
query point locates and compute the distance of the point from that face.

I locate the point with:
Face_handle     fh1 = dt.locate(pl[i].p);

But my programs gets runtime error when I want to access the vertices of
this face:
Vertex_handle v1=Face->vertex(1);
(In compute "computeError" function that which I wrote for computing the
distance of the point from the face)

does any one has any idea whats the reason and how I can resolve it?
Maybe the reason is that the query point is on an edge or vertex! If so, how
could I handle this cases?
Is the problem with my triangulation or what?
(In addition I appreciate it if some one could tell me a way for computing
the distance of the point from the face, I'm going to compute it by
accessing the vertices of the face and ..., but I don't know how to access
the vertices coordinates yet . Is there an easier way to do this?)

Thanks in advance,

Here is my code:

#include <CGAL/Polytope_distance_d.h>
#include <CGAL/Polytope_distance_d_traits_3.h>
#include <CGAL/Euclidean_distance.h>
#include <CGAL/Triangulation_vertex_base_with_info_2.h>
#include <CGAL/Exact_predicates_inexact_constructions_kernel.h>
#include <CGAL/Triangulation_2.h>
#include <CGAL/Triangulation_euclidean_traits_xy_3.h>
#include <CGAL/Delaunay_triangulation_2.h>
#include <CGAL/Triangulation_hierarchy_2.h>
#include <vector>
#include <fstream>
#include <CGAL/Exact_predicates_inexact_constructions_kernel.h>
#include <CGAL/Delaunay_triangulation_2.h>
#include <CGAL/Triangulation_hierarchy_2.h>
#include <CGAL/point_generators_2.h>
#include <CGAL/algorithm.h>
#include <cassert>

#ifdef CGAL_USE_GMP
#include <CGAL/Gmpzf.h>
typedef CGAL::Gmpzf ET;
#else
#include <CGAL/MP_Float.h>
typedef CGAL::MP_Float ET;
#endif

using namespace std;
using namespace CGAL;

typedef CGAL::Exact_predicates_inexact_constructions_kernel K;
typedef CGAL::Triangulation_vertex_base_with_info_2<int,K>      Vbb;
typedef CGAL::Triangulation_hierarchy_vertex_base_2<Vbb> Vb;
typedef CGAL::Triangulation_face_base_2<K>               Fb;
typedef CGAL::Triangulation_data_structure_2<Vb,Fb>      Tds;
typedef CGAL::Delaunay_triangulation_2<K,Tds>            Dt;
typedef CGAL::Triangulation_hierarchy_2<Dt>              Triangulation;
typedef Triangulation::Point                             Point;
typedef Triangulation::Vertex_handle                                    
Vertex_handle;
typedef Triangulation::Finite_faces_iterator                    
Finite_faces_iterator;
typedef Triangulation::Face_handle                                            
  Face_handle;
typedef Triangulation::Vertex_handle                                          
  Vertex_handle;
typedef K::Point_3                                                            
                  Point_3;
typedef CGAL::Polytope_distance_d_traits_3<K, ET, double> Traits;
typedef CGAL::Polytope_distance_d<Traits>                 Polytope_distance;

struct PointData{
        Point p;
        double z;
        double er;
        bool in;
};

double computeError(Point query,double z,Face_handle Face){
        double x=query.x();
        double y=query.y();
        K::Point_3 p(x,y,z);
        K::Point_3 P[1] = { p};
        
        K::Point_3 Q[3] = { p,p,p};

        Vertex_handle v1=Face->vertex(1);
        Vertex_handle v2=Face->vertex(2);
        Vertex_handle v3=Face->vertex(0);
  Polytope_distance pd(P, P+1, Q, Q+3);
  assert (pd.is_valid());

  // get squared distance 
  double dis=0;
  dis=sqrt(CGAL::to_double (pd.squared_distance_numerator())
/CGAL::to_double (pd.squared_distance_denominator()));


        return dis;
}
int main()
{
        std::ifstream in("terrain.txt");
        std::istream_iterator<Point> begin(in);
        std::istream_iterator<Point> end;
        Triangulation dt;
        double x,y,z;
        x=y=z=0;
        double maxX,minX,maxY,minY;
        double curMaxError=100;//this is computed
        double maxEr=0;//this must be input
        double maxZ=1000;
        vector<PointData> pl;

        bool first=true;

        while(in>>x>>y>>z){
                Point tempP(x,y);
                PointData pd;
                pd.er=0;
                pd.in=false;
                pd.z=z;
                if(first){
                        maxX=minX=x;
                        maxY=minY=y;
                        first=false;
                }
                pd.p=tempP;
                pl.push_back(pd);
        }
        //adding the 4 corner
        //*******************
        //for now
        Point t00(minX,minY);
        Point t01(minX,maxY);
        
        Point t10(maxX,minY);

        Point t11(maxX,maxY);
        dt.insert(t00);
        dt.insert(t01);
        dt.insert(t10);
        dt.insert(t11);
        //------------------

        while(curMaxError>maxEr){
                int candidIndex=0;
                double tempMax=0;
                int i;
                for(i=0;i<pl.size();i++){
                        if(pl[i].in)continue;
                        Face_handle     fh1 = dt.locate(pl[i].p);
                        Vertex_handle h1=fh1->vertex(0);
                        pl[i].er=computeError(pl[i].p,pl[i].z,fh1);
                        if(pl[i].er>tempMax){
                                tempMax=pl[i].er;
                                candidIndex=i;
                        }
                }
                dt.insert(pl[candidIndex].p);
                pl[candidIndex].in=true;
                curMaxError=tempMax;
        }

        cout << "# of all points: "<<pl.size() << endl;
        cout << "# of simplification points: "<<dt.number_of_vertices() << 
endl;
        cout << "Error: "<< curMaxError<< endl;


        return 0;
}

-- 
View this message in context: 
http://cgal-discuss.949826.n4.nabble.com/face-handle-error-tp3319970p3319970.html
Sent from the cgal-discuss mailing list archive at Nabble.com.