Delaunay Triangulation Function for TSQL and SQL Server Spatial

Previously I have written an article on creating a function that implements a Delaunay triangulation for PostGIS.

This article is about creating TSQL procedures that:

• Create a Delaunay triangulation (mesh) from 3D geometry data;
• Computes the total slope (3D) area of a Delaunay triangulation (mesh).

Create Delaunay Triangulation

The procedure that follows requires R Services for SQL Server, and Rolf Turner’s DELDIR R package, to be installed.

CREATE PROCEDURE dbo.STDelaunay( @p_geom   geometry,
                                 @p_facets geometry out)
AS
BEGIN
  Declare
    @vSrid   int,
    @polygon nvarchar(max),
    @vFacets nvarchar(max);

  if ( @p_geom is null or @p_geom.STGeometryType() = 'Point' )
    return;

  set @polygon = @p_geom.AsTextZM();
  set @vSrid = @p_geom.STSrid;

  execute sp_execute_external_script 
    @language = N'R',
    @script = N'
     library(deldir)
     # select the point x/y/z coordinates into a data frame
     points <- InputDataSet
     # calculate an approprate buffer distance (~10%):
     buffer_distance = (
         (
             abs(max(points$x) - min(points$x)) +
             abs(max(points$y) - min(points$y))
         ) / 2
     ) * (0.10)
  
     # the following use of deldir uses high precision and digits to prevent
     # slivers between the output triangles, and uses a relatively large bounding
     # box with four dummy points included to ensure that points in the
     # peripheral areas of the dataset are appropriately enveloped by their
     # corresponding triangles:
     #    points$x,
     #    points$y,
     #    points$z,
     voro = deldir(
         list( x=points$x, y=points$y, z=points$z),
         digits=22,
         frac=0.00000000000000000000000001,
         list(ndx=2,ndy=2),
         rw=c(
             min(points$x) - abs(min(points$x) - max(points$x)),
             max(points$x) + abs(min(points$x) - max(points$x)),
             min(points$y) - abs(min(points$y) - max(points$y)),
             max(points$y) + abs(min(points$y) - max(points$y))
         )
     )
     # Do the triangulation
     #
     triangles = triang.list(voro)
  
     # We will output all the triangles as a single Geometry Collection
     #
     geomCollection = "GEOMETRYCOLLECTION("
  
     # Now create the output
     # 1. Construct each triangular facet as a polygon WKT 
     # 2. Append to the GeometryCollection 
     #
     for (i in 1:length(triangles)) {
         triangle = triangles[[i]]
         wktpoly = "POLYGON(("
         for (j in 1:length(triangle$x)) {
              wktpoly = sprintf("%s%.8f %.8f %.2f,",wktpoly,triangle$x[[j]],triangle$y[[j]],triangle$z[[j]])
         }
         wktpoly = sprintf("%s%.8f %.8f %.2f))",wktpoly,triangle$x[[1]],triangle$y[[1]],triangle$z[[1]])
         if ( i < length(triangles) ) {
           geomCollection = sprintf("%s%s,",geomCollection,wktpoly)
         } else {
           geomCollection = sprintf("%s%s)",geomCollection,wktpoly)
         }
     }',
    @input_data_1 = N'SELECT distinct 
                             A.X as x,
                             A.Y as y,
                             A.Z as z
                        FROM [dbo].[STDumpPoints](geometry::STGeomFromText(@wkt,@srid)) as a',
    @params = N'@wkt varchar(max), @srid int, @geomCollection nvarchar(max) OUTPUT',
    @wkt = @polygon,
    @srid = @vSrid,
    @geomCollection = @vFacets OUTPUT;
  SET @p_facets = geometry::STGeomFromText(@vFacets,@p_geom.STSrid);
END;
GO

Example Delaunay Triangulation of Points

An example of calling this procedure is as follows:

declare @vGeom geometry = geometry::STGeomFromText(
'GEOMETRYCOLLECTION (
POINT (-90.2618913358086 33.2184179805337 403),
POINT (-90.2534717690058 33.2235960587115 332),
POINT (-90.2526790106307 33.2219983038141 393),
POINT (-90.2548360643114 33.2220152472292 349),
POINT (-90.2569931215599 33.2220321532301 375),
POINT (-90.2591702280608 33.2202360124111 329),
POINT (-90.2570132117573 33.2202191464654 242),
POINT (-90.2548561990173 33.2202022431079 311),
POINT (-90.2526991898448 33.2201853023387 296),
POINT (-90.2527193654874 33.2183723003489 375),
POINT (-90.2548763301557 33.2183892384722 331),
POINT (-90.2570332983911 33.2184061391864 246),
POINT (-90.2591902701899 33.2184230024914 292),
POINT (-90.2602038624534 33.2165408049195 228),
POINT (-90.2580468665807 33.2165321889039 295),
POINT (-90.2558898738654 33.2165235354584 305),
POINT (-90.2537328843112 33.2165148445831 396),
POINT (-90.2537432444696 33.2147017821533 329),
POINT (-90.255900189497 33.2147104705514 213),
POINT (-90.2580571376854 33.2147191215224 218),
POINT (-90.2602140890308 33.2147277350661 226),
POINT (-90.2591415442112 33.2133037851107 310),
POINT (-90.254879330931 33.2133111535788 228)
)',4326);
declare @vCollection geometry;
exec dbo.STDelaunay @p_geom = @vGeom, @p_facets = @vCollection output;
select @vCollection
union all
select @vGeom;
GO

Which looks like this:

Computing 3D Area

The triangulation can be extended to allow for the computation of the 3D area of each Delaunay triangle facet (polygon), and the total area (all facets) of the whole triangulation.

First we create a function to calculate the 3D area of a single triangular facet (polygon):

CREATE FUNCTION [dbo].[STTriangleArea3D] (
   @p_a geometry, 
   @p_b geometry, 
   @p_c geometry
)
returns float
as
/****f* STTriangleArea3D - Computes slope area of a 3D triangle
*******/
Begin
  Declare
    @ux     float,
    @uy     float,
    @uz     float,
    @vx     float,
    @vy     float,
    @vz     float,
    @crossx float,
    @crossy float,
    @crossz float,
    @absSq  float,
    @area3D float;

  IF ( @p_a is null or @p_b is null or @p_c is null ) 
    RETURN NULL;

  IF ( NOT ( @p_a.STGeometryType() = 'Point' and @p_b.STGeometryType() = 'Point' and @p_c.STGeometryType() = 'Point' ) )
    RETURN NULL;

    /**
     * Uses the formula 1/2 * | u x v | where u,v are the side vectors of the
     * triangle x is the vector cross-product
     */
    -- side vectors u and v
    set @ux = @p_b.STX - @p_a.STX;
    set @uy = @p_b.STY - @p_a.STY;
    set @uz = @p_b.Z   - @p_a.Z;

    set @vx = @p_c.STX - @p_a.STX;
    set @vy = @p_c.STY - @p_a.STY;
    set @vz = @p_c.Z   - @p_a.Z;

    -- cross-product = u x v
    set @crossx = @uy * @vz - @uz * @vy;
    set @crossy = @uz * @vx - @ux * @vz;
    set @crossz = @ux * @vy - @uy * @vx;

    -- tri area = 1/2 * | u x v |
    set @absSq  = @crossx * @crossx + @crossy * @crossy + @crossz * @crossz;
    set @area3D = sqrt(@absSq) / 2.0;

    return @area3D;
End;
GO

-- Example
SELECT [dbo].[STTriangleArea3D] (
             geometry::STGeomFromText('POINT(0 0 0)',0),
             geometry::STGeomFromText('POINT(1 0 1)',0),
             geometry::STGeomFromText('POINT(0.5 0.5 2)',0)) as Area3D;
GO

Area3D
0.82915619758885

Now create a simple polygon wrapper over this.

CREATE FUNCTION [dbo].[STTriangleArea3DByPolygon](
   @p_geom geometry
)
returns float
as
Begin
  if ( @p_geom is null or @p_geom.STGeometryType() <> 'Polygon' )
    return null;

  if ( @p_geom.STNumPoints() <> 4 )
    return null;

  return [dbo].[STTriangleArea3D]( @p_geom.STPointN(1),
                                   @p_geom.STPointN(2),
                                   @p_geom.STPointN(3) );
End;
GO

-- Example
select [dbo].[STTriangleArea3DByPolygon] (geometry::STGeomFromText('POLYGON((0 0 0, 1 0 1, 0.5 0.5 2, 0 0 0))',0)) as Area3d;
go

Area3D
0.82915619758885

Compute 3D Area of complete Delaunay Triangulation/Mesh

Finally, we can create a function that creates a Delaunay Mesh from an input geometry and computes the total 3D area.

CREATE PROCEDURE [dbo].[STDelaunayArea3D](@p_geometry geometry,
                                          @p_area3d   float out)
AS
BEGIN
  Declare
    @vInterior  bit,
    @vFacetArea float,
    @vGeomId    integer,
    @vNumGeoms  integer,
    @vGeomN     geometry,
    @vGeom      geometry;

  -- Currently we support only MultiPoint and Polygon objects (or containing GeometryCollection).
  --
  if ( @p_geometry is null 
       OR
       @p_geometry.STGeometryType() NOT IN ('MultiPoint','Polygon','GeometryCollection' ) )
    return;

  -- Set up Interior Check
  SET @vInterior = case when @p_geometry.STGeometryType() = 'Polygon' then 1 else 0 end;

  exec [dbo].[STDelaunay] @p_geom = @p_geometry, @p_facets = @vGeom output;

  set @vFacetArea = 0.0;
  set @vNumGeoms = @vGeom.STNumGeometries();
  set @vGeomId   = 1;
  while (@vGeomID <= @vNumGeoms)
  Begin
    SET @vGeomN = @vGeom.STGeometryN(@vGeomId);
    -- Check it overlaps input geometry (if polygon)
    -- Touches means it is inside a interiorRing or hole
    -- PRINT cast(@vInterior as varchar(1)) + ' ---> '+CAST(ISNULL(dbo.STDetermine(@vGeomN,@p_geometry,1),-1) as varchar(20));
    if ( ( @vInterior = 1 
           AND 
           @vGeomN.STDisjoint(@p_geometry)=0
           AND
           @vGeomN.STTouches(@p_geometry)=0 ) 
        OR @vInterior=0 )
      SET @vFacetArea += [dbo].[STTriangleArea3DByPolygon](@vGeomN)
    SET @vGeomID += 1;
  End;
  SET @p_area3d = @vFacetArea;
END;
GO

Here are some examples:

declare @vArea3D float;
declare @vGeom geometry = geometry::STGeomFromText (
'MULTIPOINT(
(-90.2618913358086 33.2184179805337 403),
(-90.2534717690058 33.2235960587115 332),
(-90.2526790106307 33.2219983038141 393),
(-90.2548360643114 33.2220152472292 349),
(-90.2569931215599 33.2220321532301 375),
(-90.2591702280608 33.2202360124111 329),
(-90.2570132117573 33.2202191464654 242),
(-90.2548561990173 33.2202022431079 311),
(-90.2526991898448 33.2201853023387 296),
(-90.2527193654874 33.2183723003489 375),
(-90.2548763301557 33.2183892384722 331),
(-90.2570332983911 33.2184061391864 246),
(-90.2591902701899 33.2184230024914 292),
(-90.2602038624534 33.2165408049195 228),
(-90.2580468665807 33.2165321889039 295),
(-90.2558898738654 33.2165235354584 305),
(-90.2537328843112 33.2165148445831 396),
(-90.2537432444696 33.2147017821533 329),
(-90.255900189497 33.2147104705514 213),
(-90.2580571376854 33.2147191215224 218),
(-90.2602140890308 33.2147277350661 226),
(-90.2591415442112 33.2133037851107 310),
(-90.254879330931 33.2133111535788 228)
)',4326);
exec [dbo].[STDelaunayArea3D] @p_geometry = @vGeom, @p_area3d = @vArea3D out;
SELECT @vArea3D as Area3D;
go

Area3D
3.21589444680464

declare @vArea3D float;
declare @vGeom geometry = geometry::STGeomFromText('GEOMETRYCOLLECTION(POINT(5 5 5),POLYGON((0 0 0,10 0 0,10 10 0,0 10 0,0 0 0)))',0);
exec [dbo].[STDelaunayArea3D] @p_geometry = @vGeom, @p_area3d = @vArea3D out;
select @vGeom.STGeometryN(2).STArea() as Area2D, @vArea3D as Area3D;
go

Area2D Area3D
100    141.42135623731

declare @vArea3D float;
declare @vGeom geometry = geometry::STGeomFromText('POLYGON((0 0 0,10 0 0,10 10 0,0 10 0,0 0 0),(1 1 1,9 1 1,9 9 1,1 9 1,1 1 1))',0);
exec [dbo].[STDelaunayArea3D] @p_geometry = @vGeom, @p_area3d = @vArea3D out;
select @vGeom.STArea() as Area2D, @vArea3D as Area3D;
go

Area2D Area3D
36     50.9116882454314

I hope this is of interest to someone out there.