Almost-Delaunay Simplices in CGAL : Detailed Code Description / Developer's Manual
Deepak Bandyopadhyay

%%% Here is the description of the core Almost-Delaunay code for CGAL, along with 
%%% some MATLAB utilities that interface with it, that are described more fully in
%%% file README.ADMATLAB.txt.   

%%% Parameters and calling ADCGAL:
ADCGAL <cutoff> <prune> <file prefix>
cutoff = maximum limit on perturbation (in the same units as the coordinates, i.e. 
         Angstrom for proteins).  This is the parameter epsilon in my definition 
	 of AD(epsilon).  To try it out, a typical default value is 1.0 Angstrom,
         though values from 0.01 to 2.0 Angstrom are useful for capturing different 
	 ranges of perturbation in proteins.
prune = maximum edge length allowed in a simplex.  This gets rid of long edges on the
        convex hull and within interior pockets and cavities.  A typical default
	value in proteins is 10 (Angstrom), though values from 6.0 to 15.0 Angstrom
	may be useful.  "SHORT" edges, triangles and tetrahedra are those with all 
	edges below the prune in length.
file prefix = 1xyz if your 3D coordinates are stored as rows of three coordinates 
              in a file "1xyz.out".  This parameter can include a full path if the 
	      .out file lies elsewhere.  If this parameter is missing, input comes
	      from the standard input (hence you can use AD as an input pipe, though
	      minor changes are required to use it as an output pipe).  Output is 
	      written to files as follows:
	      1xyz.del: short Delaunay edges, written as rows of SORTED point index 
			pairs starting at 0; example row "0 4" 
	      1xyz.AD: short AD edges with threshold>0, written as rows of sorted edge 
		       index pairs and a threshold. Example row: 5 6 0.7773 
	      The files below are generated only by the version of ADCGAL that calculates
	      simplices, and not by the edge version (the symbol ADSIMP_CGAL has to be
	      defined during compilation.  
	      1xyz.del3: Delaunay triangles with short edges, rows of 3 sorted indices
	      1xyz.AD3: AD triangles, rows of 3 sorted indices and a threshold
	      1xyz.del4: Delaunay tetrahedra, 4 indices
	      1xyz.AD4: AD tetrahedra, 4 indices and a threshold

   
%%% Code internals, description of Class hierarchy: 

The almost-Delaunay codebase can be roughly divided into four parts: 

* Core modules: 
	- classes ADmainCGAL, ADedgeCGAL, ADsimpCGAL 

* Support modules: 
	  - classes ADPointNeighbors, ADconvex, MyDelaunay  
	  - global functions: ADfunc, MySTLExtensions 
	  - global funcTORs:  ADutil  

* Library code: 
	  - CGAL-3.0 from www.cgal.org (I include patch files to replace static  
	    filters, intersection tests and more: see Installation manual on how to 
	    copy them).  
	  - Boost (graph, pairs, smart pointer in Linux, etc.):  from www.boost.org 

* Open source third party classes integrated in the project:  
	  - auto_ref (fast smart pointer by Itay Maman, available at
	    www.codeproject.com/cpp/auto_ref.asp) 
	  - sorted_vector (by Martin Holzherr, available at 
	    www.codeproject.com/vcpp/stl/sorted_vector.asp)
 	
Apart from this there is the useful MATLAB interface to this code, see
README.ADMATLAB.txt.

--------------------------------------------------------------------------------------------

Here is the current set of classes/files in the ADCGAL implementaion,
along  with the meaning and role of each module:

ADmainCGAL.cpp : instantiates and runs the other classes (there could
be  one of these running per point set if we are processing many in
parallel.   Right now I would just run many instances of ADCGAL
together)

ADedgeCGAL.cpp : ADedge objects are created in a loop (once for each
potential short AD edge), and they perform all computations for an AD
edge, and store the data that needs to be used by  ADsimpCGAL in an
EdgeInfo object before being destroyed.  Thus if we are doing parallel
computation, several  ADedges could be computed in parallel,  but the
output must be merged and only the valid edges  (with
threshold<cutoff) allowed to continue to the AD simplex stage.  It can
also be tested and compiled/runs independent of ADsimpCGAL, by not
defining the symbol ADSIMP_CGAL and compiling  it with ADmainCGAL,
ADconvex, ADPointNeighbors, and the necessary include files.

ADsimpCGAL.cpp: One such object is created to perform all the
computations for an AD simplex,  with the result for the AD edge whose
candidate centers are being checked, passed in.   There is one call to
a method compute_3() to compute AD triangles, followed immediately by
one call to compute_4() for AD tetrahedra, and then the output from
both methods (triples/quadruples  and their potential lowest
thresholds) are accumulated, and then the ADsimpCGAL object dies.
Finally the accumulated lists of AD triples/quads are sorted by index
and then by threshold, and only the minimum thrsehold is kept for the 
output stage.  This is since the threshold could be computed for any of the
non-Delaunay edges in the AD triangle, and we need only the lowest. 
Parallelizing this code would need it to be run on all edges that pass
cutoff, and then the accumulated triples and quads would be min-filtered
as described above.

ADfunc.cpp: AD epsilon compuation for different cases, and other such 
functions used by several objects, that I've chosen to keep at global 
(namespace) scope rather than put in any one object.

MyDelaunay.h: This class stores a "pruned" DT and caches and returns lists of short
edges, triangles and tetrahedra on demand.  CGAL stores the triangulation
data structure using Triangulation_hierarchy_3<Delaunay_triangulation_3> where
cells have vertices 0..3, and neighbor cells 0..3 opposite to corresponding 
vertices.  Edges are represented as Triple(Cell,index,index), and Faces as 
pair(Cell,index). 

Turns out that pruning the triangulation by deleting cells with long edges 
(which would lead to a pruned-Delaunay class that biologists working with CGAL 
could use) led to all kinds of weird crashes since algorithms for DT edge and 
face traversal were not written to handle random missing tetrahedra 
(or worse, disconnected components).  So I changed the code to introduce a 
Boolean (proxy class for boolean) info in each cell, and I report each edge 
or facet traversed only if one of the cells it touches is still live.  It won't 
work to test the cell included in the edge/facet representation, since the scheme used
to avoid edge/facet repetitions during traversal, which compares cell
pointer values, may end up reporting an edge or facet in a dead cell.  
There has to be a more efficient way to do this, but it was the most
obvious extension of what I'd already written that could work, and it's 
only done once at the beginning.

ADPointNeighbors.{cpp,h}.  This class returns sets of points near one or two
vertices, pairs of points near 2 vertices and each other, and in general
manages proximity information.  It uses adjacency list graphs (from Boost
Graph Library) of short and non-Delaunay short edges that I create
initially, and remove edges from these graphs if their thresolds are found
to be over the cutoff.  It could be modified to use kd-tree's spherical
range query eventually, though I thought the graph representation would be
simpler and allow persistent storage and removal of edges on demand.

ADutil.h: Mostly when I try to think how something that I did in MATLAB
should be done (fast) in C++, I come up with templated code, usually a
function object (sometimes an algorithm) that I then place here, since
templated code needs to be in the header.  There's stuff for indexing a
sequence with a range of indices to produce a subsequence, finding the
minimum over a range of indices, modeling a constraint equation and what
its test on points and segments, lifting, and verifying that a displaced
point was within the computed hull (which actually taught me that cases
occur where "center of bounding box" was not in the convex hull, so I
replaced it by "centroid" which always is).

ADconvex.{cpp,h}: This is a class, that computes the convex hull in two ways, depending 
on the value of the constant ADCGAL_DELAUNAY_HULL in ADcommon.h (1 for using
Delaunay insertion as described in the CGAL paper, and 0 for using QuickHull as 
implemented by CGAL's Convex_hull_3.  We could go back to this class once a
good static-filtered version of Convex_hull_3 becomes available.

MySTLExtensions.h: input/output operators and iterators for pair<T,T>, if 
we know how to handle T (useful to read/write various edge lists).  Various other 
small extensions.

SortedTuples.h: implementations of sorted_pair, sorted_triple and 
sorted_quad that get sorted on creation.

ADCGAL.h: Header file that defines all the main classes (ADmain, ADsimp, and ADedge)
and the data types of packets they exchange with one another, in one place.  

ADcommon.h: Include file with short names for kernel and data types defined.

ADdefs.h: Include file with many of the parameter settings to change the behavior
	  of AD.  Eventually, all of these should be out of this file and in the
	  metadata (i.e. read in from a config file at runtime).  Here are the 
	  switches in this file:  
	  ADCGAL_VERBOSE: Defines if overall the program should spew out trace messages as it 
	  runs (#defined 1), or not (0).  This global definition can be overridden by #define 
	  statements in individual .cpp files, though that is a gruesome hack.   
	  BOOST_NO_INTRINSIC_WCHAR_T had to be defined to make the Boost Graph Library compile
	  on Windows/MSVC7.0.  Perhaps it is not needed now, but it does not hurt.
	  AD_COMPARE_EPSILON, AD_ZERO and AD_INFINITY are appropriate small constants for fuzzy floating point
	  comparisons and checks for zero and infinity to succeed. 
	  ADCGAL_DELAUNAY_HULL if 1 leads to the Delaunay insertion convex hull implementation that is
	  stable but slower (needs Static_filters), and if 0 leads to the Convex_Hull_3 implementation that is
	  faster (can work with Cartesian<double>) but crashes a lot because of numerical issues.
  
External Libraries used: CGAL (tested with v3.0 downloaded Oct 2003), Boost Graph Library (v.2.2 downloaded Nov 2003)
External modules used: 
sorted_vector.hpp 
	      by Martin Holzherr, available from www.codeproject.com/vcpp/stl/sorted_vector.asp 
	      Used to implement sets in a fast way, as well as keep other integer lists and tuple lists sorted.  
	      There are two variants, I use the one that does not allow repeats.    
auto_ref.{cpp,hpp}      
	      by Itay Maman, available from www.codeproject.com/cpp/auto_ref.asp
	      Used to keep smart pointers to large objects (eg. ADedge,EdgeInfo,ADsimp)
	      within an STL container, rather than the object itself.  Smart pointers would
	      makes sure they get duplicated/destructed when needed by the internal
	      container mechanisms, whereas raw pointers would lead to a memory leak.
	      See "Effective STL" by Scott Meyers for more information on this.