Skip to content

Instantly share code, notes, and snippets.

@pranavkantgaur
Last active July 4, 2016 12:20
Show Gist options
  • Select an option

  • Save pranavkantgaur/bef140830b89942e6bb7 to your computer and use it in GitHub Desktop.

Select an option

Save pranavkantgaur/bef140830b89942e6bb7 to your computer and use it in GitHub Desktop.
This code reads input from a PLY file, computes its Delaunay triangulation, imports it into Linear Cell Complex and then writes it back to a PLY file.
#include <iostream>
#include <vector>
#include <math.h>
#include <map>
#include <unordered_set>
#include <CGAL/Exact_predicates_inexact_constructions_kernel.h>
#include <CGAL/Triangulation_vertex_base_with_info_3.h>
#include <CGAL/Delaunay_triangulation_3.h>
#include <CGAL/Linear_cell_complex.h>
#include <CGAL/Linear_cell_complex_constructors.h>
#include "rply.h"
#define Pi 22.0/7.0
#define INVALID_VALUE -1.0f // used in context of distances
using namespace std;
using namespace CGAL;
typedef Exact_predicates_inexact_constructions_kernel K;
typedef Linear_cell_complex_traits<3, K> Traits;
typedef Linear_cell_complex<3, 3, Traits> LCC;
typedef LCC::Dart_handle DartHandle;
struct MyDartInfo
{
template<class Refs>
struct Dart_wrapper
{
typedef CGAL::Dart<3, Refs > Dart;
typedef Cell_attribute_with_point<Refs, DartHandle > VertexAttribute;
typedef cpp11::tuple<VertexAttribute> Attributes;
};
};
struct MyIntInfo
{
template<class Refs>
struct Dart_wrapper
{
typedef CGAL::Dart<3, Refs > Dart;
typedef Cell_attribute_with_point<Refs, size_t > VertexAttribute;
typedef cpp11::tuple<VertexAttribute> Attributes;
};
};
typedef Linear_cell_complex<3, 3, Traits, MyDartInfo> LCCWithDartInfo;
typedef Linear_cell_complex<3, 3, Traits, MyIntInfo> LCCWithIntInfo;
typedef LCCWithDartInfo::Dart_handle DartHandleWithDartInfo;
typedef Triangulation_vertex_base_with_info_3<DartHandle, K> Vb;
typedef Triangulation_data_structure_3<Vb> Tds;
typedef Delaunay_triangulation_3<K, Tds, Fast_location> Delaunay;
typedef Delaunay::Vertex_handle VertexHandle;
typedef Delaunay::Cell_handle CellHandle;
typedef Point_3<K> CGALPoint;
/*! \class Triangle
\brief Represents triangle
Triangle class stores indices of vertices.
*/
class Triangle
{
public:
size_t pointIds[3]; /*< Indices of points */
};
class DegenerateVertexSet
{
public:
LCC::Dart_handle vertHandles[5];
};
class CDTGenerator
{
public:
void generate();
LCC plc; /*!< Input piecewise linear cell complex representing the input */
Delaunay DT; /*!< Intermidiate structure used for storing Delaunay tetrahedralization */
LCC cdtMesh; /*!< Output mesh */
void markInfiniteVertexDart(LCCWithIntInfo::Dart_handle, LCCWithIntInfo&, int);
bool isInfinite(LCCWithIntInfo::Dart_handle, LCCWithIntInfo&, int, size_t);
void computeDelaunayTetrahedralization();
void writePLYOutput(LCCWithIntInfo::Dart_handle, LCCWithIntInfo&, string);
void readPLCInput();
bool areGeometricallySameSegments(DartHandle, DartHandle, LCC);
bool areGeometricallySameSegmentsWithDartInfo(LCCWithDartInfo::Dart_handle, LCCWithDartInfo::Dart_handle, LCCWithDartInfo);
void sew2CellsFromEdge(LCC &);
};
#include "cdtGen.h"
static size_t pointId = 0;
static float tempPoint[3];
static size_t dimensionId = 0;
vector<CGALPoint> plcVertexVector;
vector<Triangle> plcFaceVector;
/*! \fn static int vertex_cb(p_ply_argument argument)
\brief Callback for reading vertex from PLY file
\param [in] argument represents PLY file
*/
static int vertex_cb(p_ply_argument argument)
{
long eol;
ply_get_argument_user_data(argument, NULL, &eol);
tempPoint[dimensionId++] = ply_get_argument_value(argument);
if (eol)
{
plcVertexVector.push_back(CGALPoint(tempPoint[0],tempPoint[1],tempPoint[2]));
dimensionId = 0;
}
return 1;
}
/*! \fn static int face_cb(p_ply_argument argument)
\brief Reads face information from PLY file
\param [in] argument Represents PLY file
*/
static int face_cb(p_ply_argument argument)
{
long length, value_index;
static Triangle tempFace;
ply_get_argument_property(argument, NULL, &length, &value_index);
switch (value_index)
{
case 0:
case 1:
tempFace.pointIds[pointId++] = ply_get_argument_value(argument);
break;
case 2:
tempFace.pointIds[pointId] = ply_get_argument_value(argument);
pointId = 0;
plcFaceVector.push_back(tempFace);
break;
default:
break;
}
return 1;
}
/*! \fn bool CDTGenerator::areGeometricallySameSegments(DartHandle d1, DartHandle d2)
\brief Tests whether segments represented by d1 and d2 are _geometrically_ same.
\param [in] d1 DartHandle for first segment
\param [in] d2 DartHandle for second segment
*/
bool CDTGenerator::areGeometricallySameSegments(DartHandle d1, DartHandle d2, LCC lcc)
{
if (lcc.point(d1) == lcc.point(lcc.beta(d2, 1)))
if (lcc.point(lcc.beta(d1, 1)) == lcc.point(d2))
return true;
return false;
}
/*! \fn void CDTGenerator::sew2CellsFromEdge(LCC &lcc)
* \brief Sews 2-cells sharing common edge.
* \param [in, out] lcc Linear cell complex containing 2-cells.
*/
void CDTGenerator::sew2CellsFromEdge(LCC &lcc)
{
vector<pair<DartHandle, DartHandle> > twoCellsToBeSewed;
int sewedMark = plc.get_new_mark();
if (sewedMark == -1)
{
cout << "\nNo free mark available!!";
exit(0);
}
// sew facets sharing edge
for (LCC::Dart_range::iterator segIter1 = lcc.darts().begin(), segIterEnd1 = lcc.darts().end(); segIter1 != segIterEnd1; segIter1++)
{
if (!lcc.is_marked(segIter1, sewedMark)) // not sewed till now
{
for (LCC::Dart_range::iterator segIter2 = lcc.darts().begin(), segIterEnd2 = lcc.darts().end(); segIter2 != segIterEnd2; segIter2++)
{
if (!lcc.is_marked(segIter2, sewedMark) && lcc.is_sewable<2>(segIter1, segIter2))
{
if (areGeometricallySameSegments(segIter1, segIter2, lcc)) // checks the geometry of segments
{
lcc.mark(segIter1, sewedMark);
lcc.mark(segIter2, sewedMark);
twoCellsToBeSewed.push_back(pair<DartHandle, DartHandle>(segIter1, segIter2));
break;
}
}
else
continue;
}
}
else
continue;
}
// sew the faces sharing an edge
size_t k = 0;
for (vector<pair<DartHandle, DartHandle> >::iterator dIter = twoCellsToBeSewed.begin(), dIterEnd = twoCellsToBeSewed.end(); dIter != dIterEnd; dIter++)
if (lcc.is_sewable<2>(dIter->first, dIter->second))
{
lcc.sew<2>(dIter->first, dIter->second);
k++;
}
}
/*! \fn void CDTGenerator::readPLCInput()
\brief Reads PLC from input(only PLY supported currently) file.
Reads vertex and face information and initializes corresponding linear cell complex _plc_.
*/
void CDTGenerator::readPLCInput()
{
// read PLY file(assumed to contain the PLC)
string fileName;
cout << "\nPlease enter input filename:\t";
cin >> fileName;
p_ply inputPLY = ply_open(fileName.c_str(), NULL, 0, NULL);
if (!inputPLY) exit(0);
if (!ply_read_header(inputPLY)) exit(0);
ply_set_read_cb(inputPLY, "vertex", "x", vertex_cb, NULL, 0);
ply_set_read_cb(inputPLY, "vertex", "y", vertex_cb, NULL, 0);
ply_set_read_cb(inputPLY, "vertex", "z", vertex_cb, NULL, 1);
ply_set_read_cb(inputPLY, "face", "vertex_indices", face_cb, NULL, 0);
dimensionId = 0;
if (!ply_read(inputPLY))
{
cout << "Cannot read the PLY file :(";
exit(0);
}
ply_close(inputPLY);
// Initialize PLC
size_t vertexIds[3];
CGALPoint trianglePoints[3];
for (size_t n = 0, m = plcFaceVector.size(); n < m; n++)
{
for (size_t k = 0; k < 3; k++)
vertexIds[k] = plcFaceVector[n].pointIds[k];
for (size_t i = 0; i < 3; i++)
trianglePoints[i] = CGALPoint(plcVertexVector[vertexIds[i]].x(), plcVertexVector[vertexIds[i]].y(), plcVertexVector[vertexIds[i]].z());
plc.make_triangle(trianglePoints[0], trianglePoints[1], trianglePoints[2]);
}
// sew facets sharing edge
sew2CellsFromEdge(plc);
}
/*! void CDTGenerator::markInfinteVertexDart(LCC::Dart_handle d)
* \brief Marks all darts defining the infinite vertex in LCC.
* \param [in] d Handle to the dart representing infinite vertex.
* \param [in] lcc Linear cell complex to be marked.
* \param [in] infiniteVertexMark Mark representing a dart defining the infinite vertex.
*/
void CDTGenerator::markInfiniteVertexDart(LCCWithIntInfo::Dart_handle d, LCCWithIntInfo &lcc, int infiniteVertexMark)
{
if (!lcc.is_marked(d, infiniteVertexMark))
{
lcc.mark(d, infiniteVertexMark);
d = lcc.beta(d, 2, 1);
markInfiniteVertexDart(d, lcc, infiniteVertexMark);
d = lcc.beta(d, 3, 1);
markInfiniteVertexDart(d, lcc, infiniteVertexMark);
}
return;
}
/** \fn isInfinite(LCC::Dart_handle adart, LCC lcc, size_t cell_dimension)
* \brief Tests whether the given i-cell is infinite.
* \param [in] adart a dart handle to the i-cell.
* \param [in] lcc Linear cell complex to be checked.
* \param [in] infiniteVertexMark mark representing the infinite vertex.
* \param [in] cell_dimension dimension of i-cell.
* \return True if input vertex or face is infinite
**/
bool CDTGenerator::isInfinite(LCCWithIntInfo::Dart_handle adart, LCCWithIntInfo& lcc, int infiniteVertexMark, size_t cell_dimension)
{
bool isInfinite = false;
if (cell_dimension == 0)
{
// TODO: Add code for traversing all darts associated with input 0-cell
if (lcc.is_marked(adart, infiniteVertexMark))
return true;
}
if (cell_dimension == 2)
{
for (LCCWithIntInfo::One_dart_per_incident_cell_range<0, 2>::iterator pIter = lcc.one_dart_per_incident_cell<0, 2>(adart).begin(), pIterEnd = lcc.one_dart_per_incident_cell<0, 2>(adart).end(); pIter != pIterEnd; pIter++)
{
if (lcc.is_marked(pIter, infiniteVertexMark))
{
isInfinite = true;
break;
}
}
}
return isInfinite;
}
/*! \fn void CDTGenerator::writePLYOutput(LCC::Dart_handle dartToInfiniteVertex, LCC &lcc, string fileName)
* \brief Writes mesh represented as LCC to PLY file
* \param [in] dartToInfiniteVertex Dart handle to the infinite vertex.
* \param [in] lcc Linear cell complex.
* \param [in] fileName Name of output PLY file.
*/
void CDTGenerator::writePLYOutput(LCCWithIntInfo::Dart_handle dartToInfiniteVertex, LCCWithIntInfo &lcc, string fileName)
{
p_ply lccOutputPLY;
if ((lccOutputPLY = ply_create(fileName.c_str(), PLY_ASCII, NULL, 0, NULL)) == NULL)
{
cout << "\nCannot open file for writing!!";
exit(0);
}
// mark all darts defining infinite vertex
int infiniteVertexMark = lcc.get_new_mark();
if (infiniteVertexMark == -1)
exit(0);
// Marking all darts associated with infinite vertices
markInfiniteVertexDart(dartToInfiniteVertex, lcc, infiniteVertexMark);
// count number of vertices and faces in LCC
size_t nVertices = 0, nFaces = 0;
for (LCCWithIntInfo::One_dart_per_cell_range<0>::iterator pointCountIter = lcc.one_dart_per_cell<0>().begin(), pointCountIterEnd = lcc.one_dart_per_cell<0>().end(); pointCountIter != pointCountIterEnd; pointCountIter++)
if (!isInfinite(pointCountIter, lcc, infiniteVertexMark, 0))
nVertices++;
for (LCCWithIntInfo::One_dart_per_cell_range<2>::iterator faceCountIter = lcc.one_dart_per_cell<2>().begin(), faceCountIterEnd = lcc.one_dart_per_cell<2>().end(); faceCountIter != faceCountIterEnd; faceCountIter++)
if(!isInfinite(faceCountIter, lcc, infiniteVertexMark, 2))
nFaces++;
ply_add_element(lccOutputPLY, "vertex", nVertices);
ply_add_scalar_property(lccOutputPLY, "x", PLY_FLOAT);
ply_add_scalar_property(lccOutputPLY, "y", PLY_FLOAT);
ply_add_scalar_property(lccOutputPLY, "z", PLY_FLOAT);
ply_add_element(lccOutputPLY, "face", nFaces);
ply_add_list_property(lccOutputPLY, "vertex_indices", PLY_UCHAR, PLY_INT32);
if (!ply_write_header(lccOutputPLY))
{
cout << "Header cannot be written!!";
exit(0);
}
cout << "#### Writing vertices..." << endl;
// write vertices
size_t pointId = 0;
for (LCCWithIntInfo::One_dart_per_cell_range<0>::iterator pointIter = lcc.one_dart_per_cell<0>().begin(), pointIterEnd = lcc.one_dart_per_cell<0>().end(); pointIter != pointIterEnd; pointIter++)
{
if (!isInfinite(pointIter, lcc, infiniteVertexMark, 0))
{
CGALPoint pt = lcc.point(pointIter);
ply_write(lccOutputPLY, pt.x());
ply_write(lccOutputPLY, pt.y());
ply_write(lccOutputPLY, pt.z());
lcc.info<0>(pointIter) = pointId++;
}
}
cout << "#### Writing polygons..." << endl;
// write polygons
size_t nFiniteFaces = 0, nInfiniteFaces = 0;
for (LCCWithIntInfo::One_dart_per_cell_range<2>::iterator faceIter = lcc.one_dart_per_cell<2>().begin(), faceIterEnd = lcc.one_dart_per_cell<2>().end(); faceIter != faceIterEnd; faceIter++)
{
if (!isInfinite(faceIter, lcc, infiniteVertexMark, 2))
{
ply_write(lccOutputPLY, 3);
for (LCCWithIntInfo::One_dart_per_incident_cell_range<0, 2>::iterator pointInFaceIter = lcc.one_dart_per_incident_cell<0, 2>(faceIter).begin(), pointInFaceIterEnd = lcc.one_dart_per_incident_cell<0, 2>(faceIter).end(); pointInFaceIter != pointInFaceIterEnd; pointInFaceIter++)
ply_write(lccOutputPLY, lcc.info<0>(pointInFaceIter));
nFiniteFaces++;
}
else
nInfiniteFaces++;
}
cout << "Output file written successfully!" << endl;
lcc.free_mark(infiniteVertexMark);
ply_close(lccOutputPLY);
}
/*! \fn void CDTGenerator::computeDelaunayTetrahedralization()
\brief Computes Delaunay tetrahedralization.
*/
void CDTGenerator::computeDelaunayTetrahedralization()
{
vector <pair <CGALPoint, DartHandle> > lccVertexVector;
for (LCC::One_dart_per_cell_range<0>::iterator pIter = plc.one_dart_per_cell<0>().begin(), pIterEnd = plc.one_dart_per_cell<0>().end(); pIter != pIterEnd; pIter++)
lccVertexVector.push_back(pair<CGALPoint, DartHandle>(plc.point(pIter), pIter));
DT.insert(lccVertexVector.begin(), lccVertexVector.end());
// cout << "\nDelaunay tetrahedralization computed!!";
LCCWithIntInfo DTLCC;
string fileName("delaunay.ply");
LCCWithIntInfo::Dart_handle dartToInfiniteVertex = import_from_triangulation_3(DTLCC, DT);
cout << "Writing Delaunay triangulation to output file..." << endl;
writePLYOutput(dartToInfiniteVertex, DTLCC, fileName);
}
int main()
{
CDTGenerator cdt;
cdt.readPLCInput();
cdt.computeDelaunayTetrahedralization();
return 0;
}
@Ntinmoradiya
Copy link

Hi I am unable to compile this code...
could you please help me..

Sign up for free to join this conversation on GitHub. Already have an account? Sign in to comment