diff --git a/doc/09_voronoi.md b/doc/09_voronoi.md
new file mode 100644
index 0000000..c3ec501
--- /dev/null
+++ b/doc/09_voronoi.md
@@ -0,0 +1,35 @@
+## Voronoi
+
+Function to construct the [Voronoi Diagram](https://en.wikipedia.org/wiki/Voronoi_diagram) from a dataset of scatter points, clipped to the significant area
+
+PostGIS wil include this in future versions ([doc for dev branch](http://postgis.net/docs/manual-dev/ST_Voronoi.html)) and will perform faster for sure, but in the meantime...
+
+
+### CDB_Voronoi (geom geometry[], buffer numeric DEFAULT 0.5, tolerance numeric DEFAULT 1e-9)
+
+#### Arguments
+
+| Name | Type | Description |
+|------|------|-------------|
+| geom   | geometry[]  | Array of points's geometries |
+| buffer | numeric   | enlargment ratio for the envelope area used for the restraints|
+| tolerance   | numeric |  Delaunay tolerance, optional |
+
+### Returns
+
+| Column Name | Type | Description |
+|-------------|------|-------------|
+| geom  | geometry collection | Collection of polygons of the Voronoi cells|
+
+
+#### Example Usage
+
+```sql
+WITH a AS (
+    SELECT
+    ARRAY[ST_GeomFromText('POINT(2.1744 41.403)'),ST_GeomFromText('POINT(2.1228 41.380)'),ST_GeomFromText('POINT(2.1511 41.374)'),ST_GeomFromText('POINT(2.1528 41.413)'),ST_GeomFromText('POINT(2.165 41.391)'),ST_GeomFromText('POINT(2.1498 41.371)'),ST_GeomFromText('POINT(2.1533 41.368)'),ST_GeomFromText('POINT(2.131386 41.41399)')] AS geomin
+)
+SELECT
+    CDB_voronoi(geomin, 0.2, 1e-9) as result
+FROM a;
+```
diff --git a/src/pg/sql/09_voronoi.sql b/src/pg/sql/09_voronoi.sql
new file mode 100644
index 0000000..2bb268d
--- /dev/null
+++ b/src/pg/sql/09_voronoi.sql
@@ -0,0 +1,215 @@
+-- =============================================================================================
+--
+-- CDB_Voronoi
+--
+-- =============================================================================================
+CREATE OR REPLACE FUNCTION CDB_voronoi(
+    IN geomin geometry[],
+    IN buffer numeric DEFAULT 0.5,
+    IN tolerance numeric DEFAULT 1e-9
+    )
+RETURNS geometry  AS $$
+DECLARE
+    geomout geometry;
+BEGIN
+    WITH
+    convexhull_1 as (
+        SELECT
+            ST_ConvexHull(ST_Collect(geomin)) as g,
+            buffer * |/ st_area(ST_ConvexHull(ST_Collect(geomin))) as r
+    ),
+    clipper as(
+        SELECT
+            st_buffer(ST_MinimumBoundingCircle(a.g), buffer*a.r)  as g
+        FROM convexhull_1 a
+    ),
+    env0 as (
+        SELECT
+            (st_dumppoints(st_expand(a.g, buffer*a.r))).geom as e
+        FROM  convexhull_1 a
+    ),
+    env as (
+        SELECT
+            array_agg(env0.e) as e
+        FROM  env0
+    ),
+    sample AS (
+        SELECT
+            ST_Collect(geomin || env.e) as geom
+        FROM env
+    ),
+    convexhull as (
+        SELECT
+            ST_ConvexHull(ST_Collect(geomin)) as cg
+    ),
+    tin as (
+        SELECT
+            ST_Dump(ST_DelaunayTriangles(geom, tolerance, 0)) as gd
+        FROM
+            sample
+    ),
+    tin_polygons as (
+        SELECT
+            (gd).Path as id,
+            (gd).Geom as pg,
+            ST_Centroid(ST_MinimumBoundingCircle((gd).Geom, 180)) as ct
+        FROM tin
+    ),
+    tin_lines as (
+        SELECT
+            id,
+            ST_ExteriorRing(pg) as lg
+        FROM tin_polygons
+    ),
+    tin_nodes as (
+        SELECT
+            id,
+            ST_PointN(lg,1) p1,
+            ST_PointN(lg,2) p2,
+            ST_PointN(lg,3) p3
+        FROM tin_lines
+    ),
+    tin_edges AS (
+        SELECT
+            p.id,
+            UNNEST(ARRAY[
+            ST_MakeLine(n.p1,n.p2) ,
+            ST_MakeLine(n.p2,n.p3) ,
+            ST_MakeLine(n.p3,n.p1)]) as Edge,
+            ST_Force2D(cdb_crankshaft._Find_Circle(n.p1,n.p2,n.p3)) as ct,
+            CASE WHEN st_distance(p.ct, ST_ExteriorRing(p.pg)) < tolerance THEN
+                TRUE
+            ELSE  FALSE END AS ctx,
+            p.pg,
+            ST_within(p.ct, convexhull.cg) as ctin
+        FROM
+            tin_polygons p,
+            tin_nodes n,
+            convexhull
+        WHERE p.id = n.id
+    ),
+    voro_nodes as (
+        SELECT
+            CASE WHEN x.ctx = TRUE THEN
+                ST_Centroid(x.edge)
+            ELSE
+                x.ct
+            END as xct,
+            CASE WHEN y.id is null THEN
+                CASE WHEN x.ctin = TRUE THEN
+                    ST_SetSRID(ST_MakePoint(
+                        ST_X(x.ct) + ((ST_X(ST_Centroid(x.edge)) - ST_X(x.ct)) * (1+buffer)),
+                        ST_Y(x.ct) + ((ST_Y(ST_Centroid(x.edge)) - ST_Y(x.ct)) * (1+buffer))
+                    ), ST_SRID(x.ct))
+                END
+            ELSE
+                y.ct
+            END as yct
+        FROM
+            tin_edges x
+        LEFT OUTER JOIN
+            tin_edges y
+        ON x.id <> y.id AND ST_Equals(x.edge, y.edge)
+    ),
+    voro_edges as(
+        SELECT
+            ST_LineMerge(ST_Collect(ST_MakeLine(xct, yct))) as v
+        FROM
+            voro_nodes
+    ),
+    voro_cells as(
+        SELECT
+            ST_Polygonize(
+                ST_Node(
+                    ST_LineMerge(
+                        ST_Union(v, ST_ExteriorRing(
+                            ST_Convexhull(v)
+                            )
+                        )
+                    )
+                )
+            ) as g
+        FROM
+            voro_edges
+    ),
+    voro_set as(
+        SELECT
+        (st_dump(v.g)).geom as g
+        FROM voro_cells v
+    )
+    SELECT
+        st_collect(ST_intersection(c.g, v.g))
+       -- ST_intersection(c.g, v.g)
+    INTO geomout
+    FROM
+        voro_set v,
+        clipper c;
+    RETURN geomout;
+END;
+$$ language plpgsql IMMUTABLE;
+
+/** ----------------------------------------------------------------------------------------
+  * @function   : FindCircle
+  * @precis     : Function that determines if three points form a circle. If so a table containing
+  *               centre and radius is returned. If not, a null table is returned.
+  * @version    : 1.0.1
+  * @param      : p_pt1        : First point in curve
+  * @param      : p_pt2        : Second point in curve
+  * @param      : p_pt3        : Third point in curve
+  * @return     : geometry     : In which X,Y ordinates are the centre X, Y and the Z being the radius of found circle
+  *                              or NULL if three points do not form a circle.
+  * @history    : Simon Greener - Feb 2012 - Original coding.
+  *               Rafa de la Torre - Aug 2016 - Small fix for type checking
+  * @copyright  : Simon Greener @ 2012
+  *               Licensed under a Creative Commons Attribution-Share Alike 2.5 Australia License. (http://creativecommons.org/licenses/by-sa/2.5/au/)
+**/
+CREATE OR REPLACE FUNCTION _Find_Circle(
+    IN p_pt1 geometry,
+    IN p_pt2 geometry,
+    IN p_pt3 geometry)
+  RETURNS geometry AS
+$BODY$
+DECLARE
+   v_Centre geometry;
+   v_radius NUMERIC;
+   v_CX     NUMERIC;
+   v_CY     NUMERIC;
+   v_dA     NUMERIC;
+   v_dB     NUMERIC;
+   v_dC     NUMERIC;
+   v_dD     NUMERIC;
+   v_dE     NUMERIC;
+   v_dF     NUMERIC;
+   v_dG     NUMERIC;
+BEGIN
+   IF ( p_pt1 IS NULL OR p_pt2 IS NULL OR p_pt3 IS NULL ) THEN
+      RAISE EXCEPTION 'All supplied points must be not null.';
+      RETURN NULL;
+   END IF;
+   IF ( ST_GeometryType(p_pt1) <> 'ST_Point' OR
+        ST_GeometryType(p_pt2) <> 'ST_Point' OR
+        ST_GeometryType(p_pt3) <> 'ST_Point' ) THEN
+      RAISE EXCEPTION 'All supplied geometries must be points.';
+      RETURN NULL;
+   END IF;
+   v_dA := ST_X(p_pt2) - ST_X(p_pt1);
+   v_dB := ST_Y(p_pt2) - ST_Y(p_pt1);
+   v_dC := ST_X(p_pt3) - ST_X(p_pt1);
+   v_dD := ST_Y(p_pt3) - ST_Y(p_pt1);
+   v_dE := v_dA * (ST_X(p_pt1) + ST_X(p_pt2)) + v_dB * (ST_Y(p_pt1) + ST_Y(p_pt2));
+   v_dF := v_dC * (ST_X(p_pt1) + ST_X(p_pt3)) + v_dD * (ST_Y(p_pt1) + ST_Y(p_pt3));
+   v_dG := 2.0  * (v_dA * (ST_Y(p_pt3) - ST_Y(p_pt2)) - v_dB * (ST_X(p_pt3) - ST_X(p_pt2)));
+   -- If v_dG is zero then the three points are collinear and no finite-radius
+   -- circle through them exists.
+   IF ( v_dG = 0 ) THEN
+      RETURN NULL;
+   ELSE
+      v_CX := (v_dD * v_dE - v_dB * v_dF) / v_dG;
+      v_CY := (v_dA * v_dF - v_dC * v_dE) / v_dG;
+      v_Radius := SQRT(POWER(ST_X(p_pt1) - v_CX,2) + POWER(ST_Y(p_pt1) - v_CY,2) );
+   END IF;
+   RETURN ST_SetSRID(ST_MakePoint(v_CX, v_CY, v_radius),ST_Srid(p_pt1));
+END;
+$BODY$
+  LANGUAGE plpgsql VOLATILE STRICT;
+
diff --git a/src/pg/test/expected/09_voronoi_test.out b/src/pg/test/expected/09_voronoi_test.out
new file mode 100644
index 0000000..02a0400
--- /dev/null
+++ b/src/pg/test/expected/09_voronoi_test.out
@@ -0,0 +1,7 @@
+SET client_min_messages TO WARNING;
+\set ECHO none
+       avg_area       
+----------------------
+ 0.000139602712846995
+(1 row)
+
diff --git a/src/pg/test/sql/09_voronoi_test.sql b/src/pg/test/sql/09_voronoi_test.sql
new file mode 100644
index 0000000..a827350
--- /dev/null
+++ b/src/pg/test/sql/09_voronoi_test.sql
@@ -0,0 +1,15 @@
+SET client_min_messages TO WARNING;
+\set ECHO none
+
+WITH a AS (
+    SELECT
+    ARRAY[ST_GeomFromText('POINT(2.1744 41.403)'),ST_GeomFromText('POINT(2.1228 41.380)'),ST_GeomFromText('POINT(2.1511 41.374)'),ST_GeomFromText('POINT(2.1528 41.413)'),ST_GeomFromText('POINT(2.165 41.391)'),ST_GeomFromText('POINT(2.1498 41.371)'),ST_GeomFromText('POINT(2.1533 41.368)'),ST_GeomFromText('POINT(2.131386 41.41399)')] AS geomin
+),
+b as(
+    SELECT
+        (st_dump(cdb_crankshaft.CDB_voronoi(geomin, 0.2, 1e-9))).geom as result
+    FROM a
+)
+SELECT
+    avg(st_area(result)) as avg_area
+FROM b;