CGAL users discussion list

[Daniel Duque Campayo <>]

Hello everyone,

I am having some troubles with finding the Cartesian z values of an alpha 
shape given (x,y). I have developed the code below in which I

a) iterate over a list of alpha-facets
b) find the z value of the plane containing the facet:  -(d+a*x+b*x)/c 
c) check whether this point is contained in the facet, by a Tr.has_on(Pc). 
This Tr is a triangle build from the facet's vertices, and according to the 
manual, 
http://www.cgal.org/Manual/3.3/doc_html/cgal_manual/Kernel_23_ref/Class_Triangle_3.html

" t.has_on ( Point_3<Kernel> p)          A point is on t, if it is on a 
vertex, an 
edge or the face of t."

This does not work (unless the coordinates precisely correspond to an alpha 
vertex). In fact:

Point Pc=circumcenter(Tr);      assert(Tr.has_on(Pc));

throws an error, so I must have understood it wrong. Perhaps this is due to 
rounding? (I'm using exact predicates, inexact constructions.)

My question is then: any idea on how to efficiently compute if a (x,y) 
pair "projects" onto a 3d triangle? (or any idea in general to approach this 
efficiently). I know that, depending on the alpha value, there a facet might 
be found or not, but I can work on that later. This alpha shape is pretty 
dense anyway, that's not the problem for sure (as the circumcenter example 
makes clear).

Thanks very much. I copy relevant portions of my messy code below, attach the 
whole of it.

Best,

Daniel




//...
    list<Facet> al_fs;
    as.get_alpha_shape_facets(back_inserter(al_fs),
                                Alpha_shape_3::REGULAR);
//...

      double x= //...
      double y= //...
      double z= //... get from somewhere, it does not matter

      for(list<Facet>::iterator fit=al_fs.begin();
          fit!=al_fs.end();fit++) {

        Cell_handle c=fit->first;
        int i=fit->second;
        Vertex_handle V0=c->vertex((i+1)%4);
        Vertex_handle V1=c->vertex((i+2)%4);
        Vertex_handle V2=c->vertex((i+3)%4);
        Point P0=V0->point();
        Point P1=V1->point();
        Point P2=V2->point();
        Plane Pl(P0,P1,P2);

        assert(Pl.c());
        double zint= -(Pl.d()+Pl.a()*x+Pl.b()*y)/Pl.c();
        Triangle Tr(P0,P1,P2);

        Point inter(x,y,zint);

        if( Tr.has_on(inter) ) {
          cout << fabs(zint-z)) << endl;
          break;
        }
      }

-- 
Daniel Duque
http://rincon.uam.es/dir?cw=950067138671875
// Alpha shape calculation.
// z histogram of pivots, and list of them
// Connected histogram

#include <CGAL/Exact_predicates_inexact_constructions_kernel.h>
#include <CGAL/Delaunay_triangulation_3.h>
#include <CGAL/Alpha_shape_3.h>
//#include <CGAL/intersections.h>

#include <fstream>
#include <list>
#include <vector>
#include <algorithm> // for find
#include <cstdlib> // For atof
#include <cassert>

#include "histo.h"
#include "my_vertex_base.h"

#define NDEBUG


//typedef CGAL::Alpha_shape_vertex_base_3<K>          Vb;

struct K : CGAL::Exact_predicates_inexact_constructions_kernel {};
typedef my_vertex_base<K>                           Vb;
typedef CGAL::Alpha_shape_cell_base_3<K>            Fb;
typedef CGAL::Triangulation_data_structure_3<Vb,Fb> Tds;
typedef CGAL::Delaunay_triangulation_3<K,Tds>       Triangulation_3;
typedef CGAL::Alpha_shape_3<Triangulation_3>        Alpha_shape_3;

typedef Alpha_shape_3::Alpha_iterator       Alpha_iterator;
typedef Alpha_shape_3::Cell_handle          Cell_handle;
typedef Alpha_shape_3::Vertex_handle        Vertex_handle;
typedef Alpha_shape_3::Facet                Facet;
typedef Alpha_shape_3::Edge                 Edge;

typedef K::Point_3                          Point;
typedef K::Line_3   Line;
typedef K::Plane_3  Plane;
typedef K::Triangle_3  Triangle;
typedef CGAL::Object      Object;

using std::list;
using std::vector;
using std::ifstream;
using std::ofstream;
using std::cout;
using std::endl;
using std::back_inserter;
//using std::assert;

typedef vector<double> vctr;
typedef list<Vertex_handle>::const_iterator Vl_it;

typedef std::pair<Vertex_handle, double> v_w_dist;

bool bydist (const v_w_dist& v1, const v_w_dist& v2) {
  return v1.second < v2.second;
}


//namespace params {
  double dtime;
  vctr LL(3,0);
  unsigned N=0;
//};

bool header(ifstream& hist);
bool every_time(ifstream& hist);
void every_part(ifstream& hist, vctr & , vctr & );
bool periodic(vctr&, const int);
bool crop(const Point& P);
int get_pivots(ifstream & pp, list<int>& L);

//template<class T> int norepeat_push(list<T>&,const T&);

const double scale=1.0;

int main(int argc, char * argv[]) {

  ifstream hist("HISTORY");

  assert(hist);

//  header(hist);

  histo Hz(10); // Histogram for z values
  histo profile(10); // Histogram for the intrinsic profile

  histo Hcoin1(10); // Histogram for coincidence
  histo Hcoin2(10); // Histogram for coincidences

  //histo len(20);

  int step=0;

  ofstream pivots("alphapivots.res");

  while (!hist.eof()) {

    step++;

    cout << step << endl;

//    if(step==20) break;

    if(! every_time(hist)) break;  // TODO: not quite right for dtime

//    list<Point> L;

    assert(N!=0);

    Triangulation_3 T; // (L.begin(),L.end());

    for(uint i=1;i<=N;i++) {

      vctr rr(3),vv(3);

      every_part(hist,rr,vv);  // get vx

      Vertex_handle vh = T.insert(Point(rr[0],rr[1],scale*rr[2]));
      vh->indx()=i;


//      Hz.insert(rr[2]+LL[2]/2.0);

//      Periodic b.c.:

      vctr rr0(rr);

      if(periodic(rr,0))  T.insert(Point(rr[0],rr[1],scale*rr[2]));

      rr=rr0;

      if(periodic(rr,1))  T.insert(Point(rr[0],rr[1],scale*rr[2]));

      rr=rr0;

      if(periodic(rr,3))  T.insert(Point(rr[0],rr[1],scale*rr[2]));

    }

//  Build alpha shape.-

    Triangulation_3 T2(T); // Keep copy

    Alpha_shape_3 as(T);

    T=T2;

    T2.clear(); // Thanks, T2

    Alpha_iterator aa = as.alpha_upper_bound(1.4);
    as.set_alpha(*aa);

    list<Vertex_handle> al_vs;
    as.get_alpha_shape_vertices(back_inserter(al_vs),
				Alpha_shape_3::REGULAR);

    list<Facet> al_fs;
    as.get_alpha_shape_facets(back_inserter(al_fs),
				Alpha_shape_3::REGULAR);

//    cout << al_fs.size() << endl;

    int size=0;

    for(Vl_it it=al_vs.begin();it!=al_vs.end(); it++) {
      int theindex=(*it)->indx();
      if( theindex ) size++;  // TODO: build list here
    }

    pivots << size << "     " << step << endl;

    for(Vl_it it=al_vs.begin();it!=al_vs.end(); it++) {
      int theindex=(*it)->indx();
      if( theindex )  {
	double z=((*it)->point()[2])/scale;
	Hz.insert(z+LL[2]/2.0);
	pivots << theindex << ' ' << z << endl;
      }
    }

    pivots << LL[0] << "   " << LL[1] << endl;

//  z intrinsic profile
    typedef Triangulation_3::Finite_vertices_iterator Fvi;

    for(Fvi vit=T.finite_vertices_begin();
        vit!=T.finite_vertices_end(); vit++) {

      Point P=vit->point();
      double x=P[0];
      double y=P[1];
      double z=P[2];
      if ((z>0) || (!(vit->indx()))) continue;  // TODO

      cout << z << endl;
//      Line ll(P,Point(x,y,1));

//      bool found=false;
      for(list<Facet>::iterator fit=al_fs.begin();
	  fit!=al_fs.end();fit++) {

	Cell_handle c=fit->first;
	int i=fit->second;
	Vertex_handle V0=c->vertex((i+1)%4);
	Vertex_handle V1=c->vertex((i+2)%4);
	Vertex_handle V2=c->vertex((i+3)%4);

        if( (!V0->indx() ) || (!V1->indx() ) || (!V2->indx() ) ) continue;
        Point P0=V0->point();
	if (P0[2]>0) continue; // TODO
        Point P1=V1->point();
        Point P2=V2->point();
	Plane Pl(P0,P1,P2);

	assert(Pl.c());
        double zint= -(Pl.d()+Pl.a()*x+Pl.b()*y)/Pl.c();
	Triangle Tr(P0,P1,P2);
	Point Pc=circumcenter(Tr);
	assert(Tr.has_on(Pc));

//	Object result=CGAL::intersection(ll,Pl);
	Point inter(x,y,zint);
//	bool bb=CGAL::assign(inter,result);
//	assert(bb);
	//    if(bb) cout << inter << endl;

	if( Tr.has_on(inter) ) {
	  //	  profile.insert(sqrt(squared_distance(P,inter)));
	  cout << z << " -> " << zint << endl;
	  profile.insert(fabs(zint-z));
	  //	  found=true;
	  break;
	}
	//	assert(found);
      }

    }

    continue;

// Network profiles.-
    list<Vertex_handle> v_dists;

    for(Vl_it vit=al_vs.begin();vit!=al_vs.end(); vit++) {
      Vertex_handle vv=T.nearest_vertex( (*vit)->point() );
      //      assert(T.is_vertex (vv));
      int theindex=(*vit)->indx();
      if(theindex) {
	vv->dist()=0;
	v_dists.push_back(vv);
      }
    }

    for(;;) {
      list<Vertex_handle> v2(v_dists); // copy
      v_dists.clear();

      int done=0;

      for(Vl_it vit=v2.begin();vit!=v2.end(); vit++) {

	Point P=(*vit)->point();
	double d0=(*vit)->dist();

	list<Vertex_handle> neigh;
	T.incident_vertices (*vit, back_inserter(neigh));

	for(Vl_it nit=neigh.begin();nit!=neigh.end(); nit++) {
	  if((*nit)==T.infinite_vertex()) continue;

	  double dd=d0+sqrt(CGAL::squared_distance(P,(*nit)->point()));

	  if(dd < (*nit)->dist() ) {
	    (*nit)->dist()=dd;
	    v_dists.push_back(*nit);
	    done++;
	  }

	}

      }

      //      cout << done << "  done \n";
      if(!done) break;

    }

    for(Fvi vit=T.finite_vertices_begin();
	vit!=T.finite_vertices_end(); vit++) {
      int theindex=vit->indx();
      //      double z=((*it)->point()[2])/scale;
      if(theindex) {
	double dd=vit->dist();
	profile.insert(dd);
      }
    }

  }

  hist.close();
  pivots.close();

  ofstream Hout("hist.dat");
  //  Hz.relativize();
  Hout << Hz;
  Hout.close();

  Hout.open("profile.dat");
  //  profile.relativize();
  Hout << profile;
  Hout.close();


  return 0;
}



bool crop(const Point& P) {

//   double dx=fabs(x)-pp.LL[0]/2.0;
//   double dy=fabs(y)-pp.LL[1]/2.0;
//   return (dx>0.0)||(dy>0.0);

  double x=P[0];
  double y=P[1];

  return (x> LL[0]/2.0)||
         (x<-LL[0]/2.0)||
         (y> LL[1]/2.0)||
         (y<-LL[1]/2.0);
}




bool header(ifstream& hist)
{

//  using namespace params;  // for N


  {
    std::string line;

    getline(hist,line); 
  }

  {
    int tr;  // for "trash"
    hist >> tr;

    assert(tr==0);
    hist >> tr;
  }

  hist >> N;
  N/=3;

  return true;
}


bool every_time(ifstream& hist)
{

//  using namespace params;  // for dtime, Lx, N

  // fix dtime stuff:  static int prev_time;
  static long prev_step=0;

  {
    std::string word;
    hist >> word;  // "timestep"

    if(hist.eof()) return false;
  }

  int ts;
  hist >> ts;

# ifdef DEBUG
  cout << "timestep " << ts << endl;
# endif

  hist >> N;
//  N/=3;
  {
    int tr;
    hist >> tr;
    hist >> tr;
  }

  double dt;

  hist >> dt;

  dtime=dt*(ts-prev_step);

  prev_step=ts;

  {
    hist >> LL[0];

    double tr;

    for(int i=0;i<3;i++) hist >> tr;

    hist >> LL[1];

    for(int i=0;i<3;i++) hist >> tr;

    hist >> LL[2];

  }


  return true;
}

void every_part(ifstream& hist, vctr & r, vctr & v)
{

 const int nm=1; // Particles per molecule; 3 for SPC water

 for(int k=0;k<nm;k++) {
  std::string word;
  hist >> word;

  int tr;
  double d1,d2;

  hist >> tr;
  hist >> d1;
  hist >> d2;

  if(k==0) {
    hist >> r[0] >> r[1] >> r[2] ;
    double f[3];
    hist >> f[0] >> f[1] >> f[2] ; 
    hist >> f[0] >> f[1] >> f[2] ; }
  else {
    double f[3];
    hist >> f[0] >> f[1] >> f[2] ; }
  }

  return;
}

void normalize( vctr & vv, double delta ) {
  double cc=0;

  for(uint i=0; i < vv.size() ; i++) cc+=vv[i];

  for(uint i=0; i < vv.size() ; i++) vv[i]*= (delta/cc);

  return;
}



// Apply p.b.c.:

bool periodic(vctr &rr, const int i) {

  const double edg=4.0;

  int j=i;

  if(i==3) j=0;

  double& x=rr[j];
  double L = LL[j];
  double L2= LL[j]/2.0;

  bool appl=false;

  if     (L2-x<edg) {x-=L; appl=true;}
  else if(x+L2<edg) {x+=L; appl=true;}
  
  if((i!=3)||(!appl)) return appl;

  // Special: _both_ directions

  //  if(!appl) return;

  double& y=rr[1];
  L = LL[1];
  L2= LL[1]/2.0;

  appl=false;

  if     (L2-y<edg) {y-=L; appl=true;}
  else if(y+L2<edg) {y+=L; appl=true;}

  return appl;

}


template<class T>
int norepeat_push(list<T>& listV,const T& V1) {

  typename list<T>::const_iterator li=listV.begin();

  for ( ; li !=listV.end(); ++li)
    if(V1 == *li) return 0;

  listV.push_back(V1);

  return 1;
}

#include <CGAL/Compact_container.h>
#include <CGAL/Triangulation_vertex_base_3.h>
#include <CGAL/Alpha_shape_cell_base_3.h>

const double Large=1.0e10;
const double large=1.0e9;


/* A vertex with an index; modified from CGAL/Alpha_shape_vertex_base_3.h */
template <class Gt, class Vb = CGAL::Alpha_shape_vertex_base_3<Gt> >
  class my_vertex_base
  : public Vb {
 public:

 typedef typename Vb::Cell_handle    Cell_handle;

 template < typename TDS2 >
 struct Rebind_TDS {
   typedef typename Vb::template Rebind_TDS<TDS2>::Other   Vb2;
   typedef my_vertex_base<Gt, Vb2>              Other;
 };

 typedef typename Gt::Point_3 Point;
 typedef typename Gt::FT      NT;
 typedef CGAL::Alpha_status<NT>     Alpha_status;
 typedef CGAL::Compact_container<Alpha_status>   Alpha_status_container;
 typedef typename Alpha_status_container::const_iterator Alpha_status_const_iterator;
 typedef typename Alpha_status_container::iterator   Alpha_status_iterator;
 
 private:
 int index;
 double dd;

 public:

 my_vertex_base() : Vb(), index(0), dd(Large) {}

 my_vertex_base(const Point& p) : Vb(p), index(0), dd(Large) {}

 my_vertex_base(const Point& p, Cell_handle c) : Vb(p, c), index(0), dd(Large) {}

 int& indx() {return index;}
 double& dist() {return dd;}

};

// A class to store histograms. All implemented here.

#include<vector>
#include<iostream>

using std::vector;
using std::ostream;
using std::endl;

struct data {
  data(const double vv=0,const unsigned long int ii=0):
          value(vv), times(ii)  {}
  double value;
  unsigned long int times;
};



class histo {

public:
  histo(const int dd=20) : delta(dd), h() , hm() {} // Constructor

  uint size(bool sign=true) const {
    if(sign) return h.size();
    else return hm.size();} // The lengths

  uint thedelta(void) const {return delta;} // Bins per unit length

  void insert(const double pos, const double val=1) { // Insert new point

    int slab= int( pos * delta  ) ;

    vector<data> *hh= &h; // alias

    if (slab<0) hh= &hm;

    slab= abs( slab );

    int more=slab+1-(hh->size());

    if(more>0) hh->insert(hh->end(),more,data()); // grow h if needed

    ((*hh)[slab]).value+=val;
    ((*hh)[slab]).times++;
  }

  void relativize(void) {
    for(uint i= hm.size() ; i>1; i--) {
       if(hm[i].times==0) continue;
       hm[i].value/=(double)hm[i].times;
    }
    for(uint i=0; i<h.size(); i++) {
       if(h[i].times==0) continue;
       h[i].value/=(double)h[i].times;
    }
  }


  void printout(vector<double>& x,vector<double>& y) const {

    x.clear();
    y.clear();

    for(uint i= hm.size() ; i>1; i--) {
      double xx= - double(i-1)/double(delta)  ;
      x.push_back( xx );
      y.push_back( hm[i-1].value );
    }

    for(uint i=0; i<h.size(); i++) {
      double xx= double(i)/double(delta)  ;
      x.push_back( xx );
      y.push_back( h[i].value );
    }
  }

  void printout(vector<double>& y, const bool sign=true)  {

    vector<data> *hh= &h; // alias

    if (!sign) hh= &hm;

    int more=hh->size()-y.size();

    if( more>0 )
      y.insert(y.end(),more,0.0);

    for(uint i=0; i<hh->size(); i++)
      y[i]+=(*hh)[i].value;

    return;
  }

  friend ostream& operator<<(ostream& fs, const histo & hh);

private:
  uint delta; // points per unit length

  vector<data> h; // data
  vector<data> hm; // data for -ve x values

  //  vector<unsigned long int> hm;
};

// Output to fstream:


inline ostream& operator<<(ostream& fs, const histo & hh) {

  for(uint i=0;i<hh.h.size();i++) {
    double xx= double(i)/double(hh.delta);
    fs << xx << ' ' << hh.h[i].value << ' ' << hh.h[i].times << endl;
  }

  return fs;
}

//  histo& operator/=(const double div) {
//    for(uint i=0; i< h.size(); i++)  h[i]/=div;
//    for(uint i=0; i<hm.size(); i++) hm[i]/=div;
//    return this;
//  }