diff --git a/doc/reference_processing.xml b/doc/reference_processing.xml
index 92acc3a95..af17d0db4 100644
--- a/doc/reference_processing.xml
+++ b/doc/reference_processing.xml
@@ -898,12 +898,13 @@ Return a <ulink
 url="http://en.wikipedia.org/wiki/Delaunay_triangulation">Delaunay
 triangulation</ulink> around the vertices of the input geometry.
 Output is a COLLECTION of polygons (for flags=0) or a MULTILINESTRING
-(for flags=1).  The tolerance, if any, is used to snap input vertices
+(for flags=1) or TIN (for flags=2).  The tolerance, if any, is used to snap input vertices
 togheter.
 			</para>
 
 			<para>Availability: 2.1.0 - requires GEOS &gt;= 3.4.0.</para>
 			<para>&Z_support;</para>
+			<para>&T_support;</para>
 
 		</refsection>
 			  <refsection>
diff --git a/liblwgeom/lwgeom_geos.c b/liblwgeom/lwgeom_geos.c
index 64906cfde..50bdba9c5 100644
--- a/liblwgeom/lwgeom_geos.c
+++ b/liblwgeom/lwgeom_geos.c
@@ -17,6 +17,8 @@
 
 #include <stdlib.h>
 
+LWTIN *lwtin_from_geos(const GEOSGeometry *geom, int want3d);
+
 #undef LWGEOM_PROFILE_BUILDAREA
 
 #define LWGEOM_GEOS_ERRMSG_MAXSIZE 256
@@ -1298,9 +1300,10 @@ lwgeom_offsetcurve(const LWLINE *lwline, double size, int quadsegs, int joinStyl
 #endif /* POSTGIS_GEOS_VERSION < 32 */
 }
 
-LWGEOM*
-lwgeom_delaunay_triangulation(const LWGEOM *lwgeom_in, double tolerance, int edgeOnly)
-{
+/*
+ * output = 1 for edges, 2 for TIN, 0 for polygons
+ */
+LWGEOM* lwgeom_delaunay_triangulation(const LWGEOM *lwgeom_in, double tolerance, int output) {
 #if POSTGIS_GEOS_VERSION < 34
 	lwerror("lwgeom_delaunay_triangulation: GEOS 3.4 or higher required");
 	return NULL;
@@ -1308,6 +1311,11 @@ lwgeom_delaunay_triangulation(const LWGEOM *lwgeom_in, double tolerance, int edg
 	GEOSGeometry *g1, *g3;
 	LWGEOM *lwgeom_result;
 
+	if (output < 0 || output > 2) {
+		lwerror("lwgeom_delaunay_triangulation: invalid output type specified %d", output);
+		return NULL;
+	}
+
 	initGEOS(lwnotice, lwgeom_geos_error);
 
 	g1 = (GEOSGeometry *)LWGEOM2GEOS(lwgeom_in);
@@ -1317,7 +1325,8 @@ lwgeom_delaunay_triangulation(const LWGEOM *lwgeom_in, double tolerance, int edg
 		return NULL;
 	}
 
-	g3 = GEOSDelaunayTriangulation(g1, tolerance, edgeOnly);
+	/* if output != 1 we want polys */
+	g3 = GEOSDelaunayTriangulation(g1, tolerance, output == 1);
 
 	/* Don't need input geometry anymore */
 	GEOSGeom_destroy(g1);
@@ -1331,12 +1340,21 @@ lwgeom_delaunay_triangulation(const LWGEOM *lwgeom_in, double tolerance, int edg
 	/* LWDEBUGF(3, "result: %s", GEOSGeomToWKT(g3)); */
 
 	GEOSSetSRID(g3, lwgeom_get_srid(lwgeom_in));
-	lwgeom_result = GEOS2LWGEOM(g3, lwgeom_has_z(lwgeom_in));
+
+	if (output == 2) {
+		lwgeom_result = (LWGEOM *)lwtin_from_geos(g3, lwgeom_has_z(lwgeom_in));
+	} else {
+		lwgeom_result = GEOS2LWGEOM(g3, lwgeom_has_z(lwgeom_in));
+	}
+
 	GEOSGeom_destroy(g3);
 
-	if (lwgeom_result == NULL)
-	{
-		lwerror("lwgeom_delaunay_triangulation: GEOS2LWGEOM returned null");
+	if (lwgeom_result == NULL) {
+		if (output != 2) {
+			lwerror("lwgeom_delaunay_triangulation: GEOS2LWGEOM returned null");
+		} else {
+			lwerror("lwgeom_delaunay_triangulation: lwtin_from_geos returned null");
+		}
 		return NULL;
 	}
 
diff --git a/liblwgeom/lwtin.c b/liblwgeom/lwtin.c
index d7fef8524..daf001204 100644
--- a/liblwgeom/lwtin.c
+++ b/liblwgeom/lwtin.c
@@ -13,9 +13,12 @@
 #include <stdio.h>
 #include <stdlib.h>
 #include <string.h>
+
+#include "lwgeom_geos.h"
 #include "liblwgeom_internal.h"
 #include "lwgeom_log.h"
 
+LWTIN *lwtin_from_geos(const GEOSGeometry *geom, int want3d);
 
 
 LWTIN* lwtin_add_lwtriangle(LWTIN *mobj, const LWTRIANGLE *obj)
@@ -23,6 +26,68 @@ LWTIN* lwtin_add_lwtriangle(LWTIN *mobj, const LWTRIANGLE *obj)
 	return (LWTIN*)lwcollection_add_lwgeom((LWCOLLECTION*)mobj, (LWGEOM*)obj);
 }
 
+LWTIN *lwtin_from_geos(const GEOSGeometry *geom, int want3d) {
+	int type = GEOSGeomTypeId(geom);
+	int hasZ;
+	int SRID = GEOSGetSRID(geom);
+
+	/* GEOS's 0 is equivalent to our unknown as for SRID values */
+	if ( SRID == 0 ) SRID = SRID_UNKNOWN;
+
+	if ( want3d ) {
+		hasZ = GEOSHasZ(geom);
+		if ( ! hasZ ) {
+			LWDEBUG(3, "Geometry has no Z, won't provide one");
+			want3d = 0;
+		}
+	}
+
+	switch (type) {
+		LWTRIANGLE **geoms;
+		uint32_t i, ngeoms;
+	case GEOS_GEOMETRYCOLLECTION:
+		LWDEBUG(4, "lwgeom_from_geometry: it's a Collection or Multi");
+
+		ngeoms = GEOSGetNumGeometries(geom);
+		geoms = NULL;
+		if ( ngeoms ) {
+			geoms = lwalloc(ngeoms * sizeof *geoms);
+			if (!geoms) {
+				lwerror("lwtin_from_geos: can't allocate geoms");
+				return NULL;
+			}
+			for (i=0; i<ngeoms; i++) {
+				const GEOSGeometry *poly, *ring;
+				const GEOSCoordSequence *cs;
+				POINTARRAY *pa;
+
+				poly = GEOSGetGeometryN(geom, i);
+				ring = GEOSGetExteriorRing(poly);
+				cs = GEOSGeom_getCoordSeq(ring);
+				pa = ptarray_from_GEOSCoordSeq(cs, want3d);
+
+				geoms[i] = lwtriangle_construct(SRID, NULL, pa);
+			}
+		}
+		return (LWTIN *)lwcollection_construct(TINTYPE, SRID, NULL, ngeoms, (LWGEOM **)geoms);
+	case GEOS_POLYGON:
+	case GEOS_MULTIPOINT:
+	case GEOS_MULTILINESTRING:
+	case GEOS_MULTIPOLYGON:
+	case GEOS_LINESTRING:
+	case GEOS_LINEARRING:
+	case GEOS_POINT:
+		lwerror("lwtin_from_geos: invalid geometry type for tin: %d", type);
+		break;
+
+	default:
+		lwerror("GEOS2LWGEOM: unknown geometry type: %d", type);
+		return NULL;
+	}
+
+	/* shouldn't get here */
+	return NULL;
+}
 
 void lwtin_free(LWTIN *tin)
 {
diff --git a/regress/delaunaytriangles.sql b/regress/delaunaytriangles.sql
index f87cdc66c..3f98221cd 100644
--- a/regress/delaunaytriangles.sql
+++ b/regress/delaunaytriangles.sql
@@ -3,3 +3,6 @@ SELECT 1,  ST_AsText(ST_DelaunayTriangles('MULTIPOINT(5 5, 6 0, 7 9)'::geometry)
 SELECT 2,  ST_AsText(ST_DelaunayTriangles('MULTIPOINT(5 5, 6 0, 7 9, 8 9)'::geometry));
 SELECT 3,  ST_AsText(ST_DelaunayTriangles('MULTIPOINT(5 5, 6 0, 7 9, 8 9)'::geometry, 2));
 SELECT 4,  ST_AsText(ST_DelaunayTriangles('MULTIPOINT(5 5, 6 0, 7 9, 8 9)'::geometry, 2, 1));
+SELECT 5,  ST_AsText(ST_DelaunayTriangles('MULTIPOINT(5 5, 6 0, 7 9)'::geometry, 0.0, 2));
+SELECT 6,  ST_AsText(ST_DelaunayTriangles('MULTIPOINT(5 5, 6 0, 7 9, 8 9)'::geometry, 0.0, 2));
+SELECT 7,  ST_AsText(ST_DelaunayTriangles('MULTIPOINT(5 5, 6 0, 7 9, 8 9)'::geometry, 2.0, 2));
diff --git a/regress/delaunaytriangles_expected b/regress/delaunaytriangles_expected
index 78e5c6d31..d69923b55 100644
--- a/regress/delaunaytriangles_expected
+++ b/regress/delaunaytriangles_expected
@@ -2,3 +2,6 @@
 2|GEOMETRYCOLLECTION(POLYGON((5 5,6 0,8 9,5 5)),POLYGON((5 5,8 9,7 9,5 5)))
 3|GEOMETRYCOLLECTION(POLYGON((5 5,6 0,7 9,5 5)))
 4|MULTILINESTRING((5 5,7 9),(5 5,6 0),(6 0,7 9))
+5|TIN(((5 5,6 0,7 9,5 5)))
+6|TIN(((5 5,6 0,8 9,5 5)),((5 5,8 9,7 9,5 5)))
+7|TIN(((5 5,6 0,7 9,5 5)))