  We present a numerically robust algorithm for computing the constrained
Delaunay tetrahedrization (CDT) of a piecewise-linear complex, which has a 100%
success rate on the 4408 valid models in the Thingi10k dataset. We build on the
underlying theory of the well-known TetGen software, but use a floating-point
implementation based on indirect geometric predicates to implicitly represent
Steiner points: this new approach dramatically simplifies the implementation,
removing the need for ad-hoc tolerances in geometric operations. Our approach
leads to a robust and parameter-free implementation, with an empirically
manageable number of added Steiner points. Furthermore, our algorithm addresses
a major gap in TetGen's theory which may lead to algorithmic failure on valid
models, even when assuming perfect precision in the calculations. Our output
tetrahedrization conforms with the input geometry without approximations. We
can further round our output to floating-point coordinates for downstream
applications, which almost always results in valid floating-point meshes unless
the input triangulation is very close to being degenerate.
