builtins: Implement ST_GeneratePoints function

Release note (sql change): Implement geo builtin ST_GeneratePoints

docs/generated/sql/functions.md

@@ -1686,6 +1686,10 @@ Bottom Left.</p>
 <p>Smaller densify_frac gives a more accurate Fréchet distance. However, the computation time and memory usage increases with the square of the number of subsegments.</p>
 <p>This function utilizes the GEOS module.</p>
 </span></td></tr>
+<tr><td><a name="st_generatepoints"></a><code>st_generatepoints(geometry: geometry, npoints: int4) &rarr; geometry</code></td><td><span class="funcdesc"><p>Generates pseudo-random points until the requested number are found within the input area. Uses system time as a seed.</p>
+</span></td></tr>
+<tr><td><a name="st_generatepoints"></a><code>st_generatepoints(geometry: geometry, npoints: int4, seed: int4) &rarr; geometry</code></td><td><span class="funcdesc"><p>Generates pseudo-random points until the requested number are found within the input area.</p>
+</span></td></tr>
 <tr><td><a name="st_geogfromewkb"></a><code>st_geogfromewkb(val: <a href="bytes.html">bytes</a>) &rarr; geography</code></td><td><span class="funcdesc"><p>Returns the Geography from an EWKB representation.</p>
 </span></td></tr>
 <tr><td><a name="st_geogfromewkt"></a><code>st_geogfromewkt(val: <a href="string.html">string</a>) &rarr; geography</code></td><td><span class="funcdesc"><p>Returns the Geography from an EWKT representation.</p>

pkg/geo/geomfn/generate_points.go

@@ -0,0 +1,228 @@
+// Copyright 2020 The Cockroach Authors.
+//
+// Use of this software is governed by the Business Source License
+// included in the file licenses/BSL.txt.
+//
+// As of the Change Date specified in that file, in accordance with
+// the Business Source License, use of this software will be governed
+// by the Apache License, Version 2.0, included in the file
+// licenses/APL.txt.
+
+package geomfn
+
+import (
+	"math"
+	"math/rand"
+
+	"github.com/cockroachdb/cockroach/pkg/geo"
+	"github.com/cockroachdb/cockroach/pkg/geo/geopb"
+	"github.com/cockroachdb/cockroach/pkg/geo/geos"
+	"github.com/cockroachdb/errors"
+	"github.com/twpayne/go-geom"
+)
+
+// GeneratePoints generates provided number of pseudo-random points for the input area.
+func GeneratePoints(g geo.Geometry, npoints int, seed int64) (geo.Geometry, error) {
+	if npoints < 0 {
+		return geo.Geometry{}, nil
+	}
+	if seed < 1 {
+		return geo.Geometry{}, errors.New("seed must be greater than zero")
+	}
+	pointsGeometry, err := geometryToPoints(g, npoints, seed)
+	if err != nil {
+		return geo.Geometry{}, errors.Newf("generating random points error: %v", err)
+	}
+	return pointsGeometry, nil
+}
+
+// geometryToPoints returns a MultiPoint geometry consisting of randomly generated points
+// that are covered by the geometry provided.
+// npoints is the number of points to return.
+// seed is the seed used by the random numbers generator.
+func geometryToPoints(g geo.Geometry, npoints int, seed int64) (geo.Geometry, error) {
+	var generatePointsFunction func(g geo.Geometry, npoints int, seed int64) (*geom.MultiPoint, error)
+	switch g.ShapeType() {
+	case geopb.ShapeType_Polygon:
+		generatePointsFunction = polygonGeometryToPoints
+	case geopb.ShapeType_MultiPolygon:
+		generatePointsFunction = multipolygonGeometryToPoints
+	default:
+		return geo.Geometry{}, errors.New("unsupported type")
+	}
+	// This is to be checked once we know Geometry type is supported,
+	// so that we can keep consistency with PostGIS implementation.
+	if npoints == 0 {
+		return geo.Geometry{}, nil
+	}
+	empty, err := IsEmpty(g)
+	if err != nil {
+		return geo.Geometry{}, errors.Newf("could not check if Geometry is empty: %v", err)
+	}
+	if empty {
+		return geo.Geometry{}, nil
+	}
+	mpt, err := generatePointsFunction(g, npoints, seed)
+	if err != nil {
+		return geo.Geometry{}, err
+	}
+	srid := g.SRID()
+	mpt.SetSRID(int(srid))
+	out, err := geo.MakeGeometryFromGeomT(mpt)
+	if err != nil {
+		return geo.Geometry{}, errors.Newf("could not transform geom.T into Geometry: %v", err)
+	}
+	return out, nil
+}
+
+func polygonGeometryToPoints(g geo.Geometry, npoints int, seed int64) (*geom.MultiPoint, error) {
+	area, err := Area(g)
+	if err != nil {
+		return nil, errors.Newf("could not calculate Polygon area: %v", err)
+	}
+	bbox := g.CartesianBoundingBox()
+	bboxWidth := bbox.HiX - bbox.LoX
+	bboxHeight := bbox.HiY - bbox.LoY
+	bboxArea := bboxHeight * bboxWidth
+	if area == 0.0 || bboxArea == 0.0 {
+		return nil, errors.New("zero area input Polygon")
+	}
+	// Gross up our test set a bit to increase odds of getting coverage in one pass.
+	sampleNpoints := float64(npoints) * bboxArea / area
+
+	sampleSqrt := math.Round(math.Sqrt(sampleNpoints))
+	if sampleSqrt == 0 {
+		sampleSqrt = 1
+	}
+	var sampleHeight, sampleWidth int
+	var sampleCellSize float64
+	// Calculate the grids we're going to randomize within.
+	if bboxWidth > bboxHeight {
+		sampleWidth = int(sampleSqrt)
+		sampleHeight = int(math.Ceil(sampleNpoints / float64(sampleWidth)))
+		sampleCellSize = bboxWidth / float64(sampleWidth)
+	} else {
+		sampleHeight = int(sampleSqrt)
+		sampleWidth = int(math.Ceil(sampleNpoints / float64(sampleHeight)))
+		sampleCellSize = bboxHeight / float64(sampleHeight)
+	}
+	// Prepare the polygon for fast true/false testing.
+	gprep, err := geos.PrepareGeometry(g.EWKB())
+	if err != nil {
+		return nil, errors.Newf("could not prepare Geometry: %v", err)
+	}
+	defer func() {
+		err = geos.PreparedGeomDestroy(gprep)
+	}()
+
+	// Initiate random number generator.
+	randomNumberGenerator := rand.New(rand.NewSource(seed))
+	// Generate a slice of points - for every cell on a grid store coordinates.
+	n := sampleHeight * sampleWidth
+	cells := make([]geom.Coord, n)
+	for i := 0; i < sampleWidth; i++ {
+		for j := 0; j < sampleHeight; j++ {
+			cells[i*sampleHeight+j] = geom.Coord{float64(i), float64(j)}
+		}
+	}
+	// Shuffle the points. Without shuffling, the generated point will
+	// always be adjacent to the previous one (in terms of grid cells).
+	if n > 1 {
+		randomNumberGenerator.Shuffle(n, func(i int, j int) {
+			temp := cells[j]
+			cells[j] = cells[i]
+			cells[i] = temp
+		})
+	}
+	results := geom.NewMultiPoint(geom.XY)
+	// Generate points and test them.
+	npointsGenerated := 0
+	iterations := 0
+	for npointsGenerated < npoints {
+		iterations++
+		for _, cell := range cells {
+			y := bbox.LoY + cell.X()*sampleCellSize
+			x := bbox.LoX + cell.Y()*sampleCellSize
+			x += randomNumberGenerator.Float64() * sampleCellSize
+			y += randomNumberGenerator.Float64() * sampleCellSize
+			if x > bbox.HiX || y > bbox.HiY {
+				continue
+			}
+			gseq, err := geos.CoordSequenceCreate(1, 2)
+			if err != nil {
+				return nil, errors.Newf("could not create coord sequence: %v", err)
+			}
+			gseq, err = geos.CoordSequenceSetXY(gseq, 0, x, y)
+			if err != nil {
+				return nil, errors.Newf("could not set X and Y for coord sequence: %v", err)
+			}
+			defer func() {
+				// err = geos.CoordSequenceDestroy(gseq)
+			}()
+			gpt, err := geos.GeometryPointCreate(gseq)
+			if err != nil {
+				return nil, errors.Newf("could not create Geometry Point: %v", err)
+			}
+			intersects, err := geos.PreparedIntersects(gprep, gpt)
+			if err != nil {
+				return nil, errors.Newf("could not check prepared intersection: %v", err)
+			}
+			if intersects {
+				npointsGenerated++
+				p := geom.NewPointFlat(geom.XY, []float64{x, y})
+				srid := g.SRID()
+				p.SetSRID(int(srid))
+				err = results.Push(p)
+				if err != nil {
+					return nil, errors.Newf("could not add point to the results: %v", err)
+				}
+				if npointsGenerated == npoints {
+					break
+				}
+			}
+		}
+		if npointsGenerated == npoints || iterations > 100 {
+			break
+		}
+	}
+	return results, nil
+}
+
+func multipolygonGeometryToPoints(
+	g geo.Geometry, npoints int, seed int64,
+) (*geom.MultiPoint, error) {
+	results := geom.NewMultiPoint(geom.XY)
+
+	area, err := Area(g)
+	if err != nil {
+		return nil, errors.Newf("could not calculate MultiPolygon area: %v", err)
+	}
+
+	gt, err := g.AsGeomT()
+	if err != nil {
+		return nil, errors.Newf("could not transform MultiPolygon into geom.T: %v", err)
+	}
+
+	gmp := gt.(*geom.MultiPolygon)
+	for i := 0; i < gmp.NumPolygons(); i++ {
+		poly := gmp.Polygon(i)
+		subarea := poly.Area()
+		subNpoints := int(math.Round(float64(npoints) * subarea / area))
+		if subNpoints > 0 {
+			g, err := geo.MakeGeometryFromGeomT(poly)
+			if err != nil {
+				return nil, errors.Newf("could not transform geom.T into Geometry: %v", err)
+			}
+			subMpt, err := polygonGeometryToPoints(g, subNpoints, seed)
+			if err != nil {
+				return nil, errors.Newf("error generating points for Polygon: %v", err)
+			}
+			for j := 0; j < subMpt.NumPoints(); j++ {
+				if err := results.Push(subMpt.Point(j)); err != nil {
+					return nil, errors.Newf("could not push point to the results: %v", err)
+				}
+			}
+		}
+	}
+	return results, nil
+}

pkg/geo/geomfn/generate_points_test.go

@@ -0,0 +1,121 @@
+// Copyright 2020 The Cockroach Authors.
+//
+// Use of this software is governed by the Business Source License
+// included in the file licenses/BSL.txt.
+//
+// As of the Change Date specified in that file, in accordance with
+// the Business Source License, use of this software will be governed
+// by the Apache License, Version 2.0, included in the file
+// licenses/APL.txt.
+
+package geomfn
+
+import (
+	"testing"
+
+	"github.com/cockroachdb/cockroach/pkg/geo"
+	"github.com/stretchr/testify/require"
+)
+
+func TestGeneratePoints(t *testing.T) {
+	type args struct {
+		g       geo.Geometry
+		npoints int
+		seed    int64
+	}
+	tests := []struct {
+		name    string
+		args    args
+		want    geo.Geometry
+		wantErr bool
+	}{
+		{
+			"number of points to generate less than minimum",
+			args{geo.MustParseGeometry("POLYGON((0 0, 1 0, 1 1.1, 0 0))"), -1, 1},
+			geo.Geometry{},
+			false,
+		},
+		{
+			"supported geometry, zero points to generate",
+			args{geo.MustParseGeometry("POLYGON((0 0, 1 0, 1 1.1, 0 0))"), 0, 1},
+			geo.Geometry{},
+			false,
+		},
+		{
+			"unsupported geometry, zero points to generate",
+			args{geo.MustParseGeometry("LINESTRING(50 50,150 150,150 50)"), 0, 1},
+			geo.Geometry{},
+			true,
+		},
+		{
+			"provided seed less than minimum",
+			args{geo.MustParseGeometry("POLYGON((0 0, 1 0, 1 1.1, 0 0))"), 1, -1},
+			geo.Geometry{},
+			true,
+		},
+		{
+			"unsupported geometry type: LineString",
+			args{geo.MustParseGeometry("LINESTRING(50 50,150 150,150 50)"), 4, 1},
+			geo.Geometry{},
+			true,
+		},
+		{
+			"Polygon with zero area",
+			args{geo.MustParseGeometry("POLYGON((0 0, 1 1, 1 1, 0 0))"), 4, 1},
+			geo.Geometry{},
+			true,
+		},
+		{
+			"empty geometry",
+			args{geo.MustParseGeometry("POLYGON EMPTY"), 4, 1},
+			geo.Geometry{},
+			false,
+		},
+		{
+			"Polygon - square",
+			args{geo.MustParseGeometry("POLYGON((1 1,1 2,2 2,2 1,1 1))"), 4, 2},
+			geo.MustParseGeometry("MULTIPOINT(1.55948870460174560115 1.30713530828859747501, 1.4638274793965568854 1.71250027479276534237, 1.6030771445518701146 1.60597601333107697918, 1.06344848690006488212 1.30988904426268293335)"),
+			false,
+		},
+		{
+			"Polygon, width > height",
+			args{geo.MustParseGeometry("POLYGON((0 0,1 2,3 2,2 0,0 0))"), 7, 14},
+			geo.MustParseGeometry("MULTIPOINT(0.3830566789833752539 0.23117069178437729682, 2.51179632398947294547 1.69060434239713908156, 1.94222431495883451902 0.56512577777958861169, 1.87408814538545565043 1.21241169406013726828, 2.29499038937696608897 1.19218670122114289711, 1.93144882885377144888 0.79976266805657403314, 1.33675111888047548625 1.70583597131752906506)"),
+			false,
+		},
+		{
+			"Polygon, width < height",
+			args{geo.MustParseGeometry("POLYGON((0 0,2 5,2.5 4,3 5,3 1,0 0))"), 5, 1996},
+			geo.MustParseGeometry("MULTIPOINT(0.69465779472271460548 0.92001319545446269554, 2.4417593921811042712 3.71642371685872197062, 2.79787890688424933927 3.8425013135166361522, 1.05776032659919683176 1.77173131482243051416, 1.79695770199420046254 2.42853164217679839965)"),
+			false,
+		},
+		{
+			"Polygon with a hole",
+			args{geo.MustParseGeometry("POLYGON((-1 -1, -1 1, 0 2, 1 1, 1 -1, 0 -2, -1 -1),(-0.5 -0.5, -0.5 0.5, 0.5 0.5, 0.5 -0.5, 0 -0.5, -0.5 -0.5))"), 7, 42},
+			geo.MustParseGeometry("MULTIPOINT(0.62392209943416365725 1.14891034830739080519, 0.5066157824317096825 -1.39898866703452817717, -0.71853651145541486134 1.15646037366762399756, -0.03771709502060871522 -0.80984418343170694321, 0.56879738017966463559 -0.74897977355716416348, 0.05783156505961972726 1.37709385168907050279, 0.71334734167548885519 -0.61452515240937377605)"),
+			false,
+		},
+		{
+			"MultiPolygon",
+			args{geo.MustParseGeometry("MULTIPOLYGON(((0 0,4 0,4 4,0 4,0 0),(1 1,2 1,2 2,1 2,1 1)), ((0 0,-4 0,-4 -4,0 -4,0 0),(1 1,2 1,2 2,1 2,1 1)))"), 9, 15},
+			geo.MustParseGeometry("MULTIPOINT(3.22867259458672961614 3.29323429726158023456, 1.70359136534834165744 0.36049598512632602398, 2.71072342994228554502 1.20150277742601008235, 1.9059667543503933107 3.67598273180139756278, 2.54757976155147858321 3.35228130507990673692, -0.1495732283473727442 -3.65380934845110161291, -2.3391608436532123072 -0.37011413317119368216, -3.6902632348743011903 -3.45769500629712034367, -1.72346872742639400933 -1.83187553147969461875, -0.04275659844954238231 -2.64469612218461502806)"),
+			false,
+		},
+		{
+			"Polygon with a specified SRID",
+			args{geo.MustParseGeometry("SRID=4326;POLYGON((0 0,2 5,2.5 4,3 5,3 1,0 0))"), 5, 1996},
+			geo.MustParseGeometry("SRID=4326;MULTIPOINT(0.69465779472271460548 0.92001319545446269554, 2.4417593921811042712 3.71642371685872197062, 2.79787890688424933927 3.8425013135166361522, 1.05776032659919683176 1.77173131482243051416, 1.79695770199420046254 2.42853164217679839965)"),
+			false,
+		},
+	}
+	for _, tt := range tests {
+		t.Run(tt.name, func(t *testing.T) {
+			got, err := GeneratePoints(tt.args.g, tt.args.npoints, tt.args.seed)
+			if (err != nil) != tt.wantErr {
+				t.Errorf("GeneratePoints() error = %v, wantErr %v", err, tt.wantErr)
+				return
+			}
+			require.Equal(t, tt.want, got)
+		})
+	}
+}