From acdf2a8b372aec1da09370fca77ff7dccac7646d Mon Sep 17 00:00:00 2001
From: Teodor Sigaev <teodor@sigaev.ru>
Date: Wed, 30 Mar 2016 18:42:36 +0300
Subject: [PATCH] Introduce SP-GiST operator class over box.

Patch implements quad-tree over boxes, naive approach of 2D quad tree will not
work for any non-point objects because splitting space on node is not
efficient. The idea of pathc is treating 2D boxes as 4D points, so,
object will not overlap (in 4D space).

The performance tests reveal that this technique especially beneficial
with too much overlapping objects, so called "spaghetti data".

Author: Alexander Lebedev with editorization by Emre Hasegeli and me
---
 doc/src/sgml/spgist.sgml                 |  18 +
 src/backend/utils/adt/Makefile           |   4 +-
 src/backend/utils/adt/geo_spgist.c       | 699 +++++++++++++++++++++++
 src/include/catalog/pg_amop.h            |  16 +
 src/include/catalog/pg_amproc.h          |   5 +
 src/include/catalog/pg_opclass.h         |   1 +
 src/include/catalog/pg_opfamily.h        |   1 +
 src/include/catalog/pg_proc.h            |  11 +
 src/include/utils/geo_decls.h            |   6 +
 src/test/regress/expected/box.out        | 240 ++++++++
 src/test/regress/expected/opr_sanity.out |   6 +-
 src/test/regress/sql/box.sql             |  62 ++
 12 files changed, 1066 insertions(+), 3 deletions(-)
 create mode 100644 src/backend/utils/adt/geo_spgist.c

diff --git a/doc/src/sgml/spgist.sgml b/doc/src/sgml/spgist.sgml
index 0d44d064aa..24c11bae71 100644
--- a/doc/src/sgml/spgist.sgml
+++ b/doc/src/sgml/spgist.sgml
@@ -112,6 +112,24 @@
        <literal>@&gt;</>
       </entry>
      </row>
+     <row>
+      <entry><literal>box_ops</></entry>
+      <entry>box</entry>
+      <entry>
+       <literal>&lt;&lt;</>
+       <literal>&amp;&lt;</>
+       <literal>&amp;&amp;</>
+       <literal>&amp;&gt;</>
+       <literal>&gt;&gt;</>
+       <literal>~=</>
+       <literal>@&gt;</>
+       <literal>&lt;@</>
+       <literal>&amp;&lt;|</>
+       <literal>&lt;&lt;|</>
+       <literal>|&gt;&gt;</literal>
+       <literal>|&amp;&gt;</>
+      </entry>
+     </row>
      <row>
       <entry><literal>text_ops</></entry>
       <entry><type>text</></entry>
diff --git a/src/backend/utils/adt/Makefile b/src/backend/utils/adt/Makefile
index 2cb7baba42..2b4ebc72b4 100644
--- a/src/backend/utils/adt/Makefile
+++ b/src/backend/utils/adt/Makefile
@@ -14,8 +14,8 @@ OBJS = acl.o arrayfuncs.o array_expanded.o array_selfuncs.o \
 	bool.o cash.o char.o date.o datetime.o datum.o dbsize.o domains.o \
 	encode.o enum.o expandeddatum.o \
 	float.o format_type.o formatting.o genfile.o \
-	geo_ops.o geo_selfuncs.o inet_cidr_ntop.o inet_net_pton.o int.o \
-	int8.o json.o jsonb.o jsonb_gin.o jsonb_op.o jsonb_util.o \
+	geo_ops.o geo_selfuncs.o geo_spgist.o inet_cidr_ntop.o inet_net_pton.o \
+	int.o int8.o json.o jsonb.o jsonb_gin.o jsonb_op.o jsonb_util.o \
 	jsonfuncs.o like.o lockfuncs.o mac.o misc.o nabstime.o name.o \
 	network.o network_gist.o network_selfuncs.o \
 	numeric.o numutils.o oid.o oracle_compat.o \
diff --git a/src/backend/utils/adt/geo_spgist.c b/src/backend/utils/adt/geo_spgist.c
new file mode 100644
index 0000000000..7604f1e8ea
--- /dev/null
+++ b/src/backend/utils/adt/geo_spgist.c
@@ -0,0 +1,699 @@
+/*-------------------------------------------------------------------------
+ *
+ * geo_spgist.c
+ *	  SP-GiST implementation of 4-dimensional quad tree over boxes
+ *
+ * This module provides SP-GiST implementation for boxes using quad tree
+ * analogy in 4-dimensional space.  SP-GiST doesn't allow indexing of
+ * overlapping objects.  We are making 2D objects never-overlapping in
+ * 4D space.  This technique has some benefits compared to traditional
+ * R-Tree which is implemented as GiST.  The performance tests reveal
+ * that this technique especially beneficial with too much overlapping
+ * objects, so called "spaghetti data".
+ *
+ * Unlike the original quad tree, we are splitting the tree into 16
+ * quadrants in 4D space.  It is easier to imagine it as splitting space
+ * two times into 4:
+ *
+ *				|	   |
+ *				|	   |
+ *				| -----+-----
+ *				|	   |
+ *				|	   |
+ * -------------+-------------
+ *				|
+ *				|
+ *				|
+ *				|
+ *				|
+ *
+ * We are using box datatype as the prefix, but we are treating them
+ * as points in 4-dimensional space, because 2D boxes are not enough
+ * to represent the quadrant boundaries in 4D space.  They however are
+ * sufficient to point out the additional boundaries of the next
+ * quadrant.
+ *
+ * We are using traversal values provided by SP-GiST to calculate and
+ * to store the bounds of the quadrants, while traversing into the tree.
+ * Traversal value has all the boundaries in the 4D space, and is is
+ * capable of transferring the required boundaries to the following
+ * traversal values.  In conclusion, three things are necessary
+ * to calculate the next traversal value:
+ *
+ *	(1) the traversal value of the parent
+ *	(2) the quadrant of the current node
+ *	(3) the prefix of the current node
+ *
+ * If we visualize them on our simplified drawing (see the drawing above);
+ * transfered boundaries of (1) would be the outer axis, relevant part
+ * of (2) would be the up right part of the other axis, and (3) would be
+ * the inner axis.
+ *
+ * For example, consider the case of overlapping.  When recursion
+ * descends deeper and deeper down the tree, all quadrants in
+ * the current node will be checked for overlapping.  The boundaries
+ * will be re-calculated for all quadrants.  Overlap check answers
+ * the question: can any box from this quadrant overlap with the given
+ * box?  If yes, then this quadrant will be walked.  If no, then this
+ * quadrant will be skipped.
+ *
+ * This method provides restrictions for minimum and maximum values of
+ * every dimension of every corner of the box on every level of the tree
+ * except the root.  For the root node, we are setting the boundaries
+ * that we don't yet have as infinity.
+ *
+ * Portions Copyright (c) 1996-2016, PostgreSQL Global Development Group
+ * Portions Copyright (c) 1994, Regents of the University of California
+ *
+ * IDENTIFICATION
+ *			src/backend/utils/adt/geo_spgist.c
+ *
+ *-------------------------------------------------------------------------
+ */
+
+#include "postgres.h"
+
+#include "access/spgist.h"
+#include "access/stratnum.h"
+#include "catalog/pg_type.h"
+#include "utils/builtins.h"
+#include "utils/geo_decls.h"
+
+/*
+ * Comparator for qsort
+ *
+ * We don't need to use the floating point macros in here, because this
+ * is going only going to be used in a place to effect the performance
+ * of the index, not the correctness.
+ */
+static int
+compareDoubles(const void *a, const void *b)
+{
+	double		x = *(double *) a;
+	double		y = *(double *) b;
+
+	if (x == y)
+		return 0;
+	return (x > y) ? 1 : -1;
+}
+
+typedef struct
+{
+	double		low;
+	double		high;
+}	Range;
+
+typedef struct
+{
+	Range		left;
+	Range		right;
+}	RangeBox;
+
+typedef struct
+{
+	RangeBox	range_box_x;
+	RangeBox	range_box_y;
+}	RectBox;
+
+/*
+ * Calculate the quadrant
+ *
+ * The quadrant is 8 bit unsigned integer with 4 least bits in use.
+ * This function accepts BOXes as input.  They are not casted to
+ * RangeBoxes, yet.  All 4 bits are set by comparing a corner of the box.
+ * This makes 16 quadrants in total.
+ */
+static uint8
+getQuadrant(BOX *centroid, BOX *inBox)
+{
+	uint8		quadrant = 0;
+
+	if (inBox->low.x > centroid->low.x)
+		quadrant |= 0x8;
+
+	if (inBox->high.x > centroid->high.x)
+		quadrant |= 0x4;
+
+	if (inBox->low.y > centroid->low.y)
+		quadrant |= 0x2;
+
+	if (inBox->high.y > centroid->high.y)
+		quadrant |= 0x1;
+
+	return quadrant;
+}
+
+/*
+ * Get RangeBox using BOX
+ *
+ * We are turning the BOX to our structures to emphasize their function
+ * of representing points in 4D space.  It also is more convenient to
+ * access the values with this structure.
+ */
+static RangeBox *
+getRangeBox(BOX *box)
+{
+	RangeBox   *range_box = (RangeBox *) palloc(sizeof(RangeBox));
+
+	range_box->left.low = box->low.x;
+	range_box->left.high = box->high.x;
+
+	range_box->right.low = box->low.y;
+	range_box->right.high = box->high.y;
+
+	return range_box;
+}
+
+/*
+ * Initialize the traversal value
+ *
+ * In the beginning, we don't have any restrictions.  We have to
+ * initialize the struct to cover the whole 4D space.
+ */
+static RectBox *
+initRectBox()
+{
+	RectBox	   *rect_box = (RectBox *) palloc(sizeof(RectBox));
+	double		infinity = get_float8_infinity();
+
+	rect_box->range_box_x.left.low = -infinity;
+	rect_box->range_box_x.left.high = infinity;
+
+	rect_box->range_box_x.right.low = -infinity;
+	rect_box->range_box_x.right.high = infinity;
+
+	rect_box->range_box_y.left.low = -infinity;
+	rect_box->range_box_y.left.high = infinity;
+
+	rect_box->range_box_y.right.low = -infinity;
+	rect_box->range_box_y.right.high = infinity;
+
+	return rect_box;
+}
+
+/*
+ * Calculate the next traversal value
+ *
+ * All centroids are bounded by RectBox, but SP-GiST only keeps
+ * boxes.  When we are traversing the tree, we must calculate RectBox,
+ * using centroid and quadrant.
+ */
+static RectBox *
+nextRectBox(RectBox *rect_box, RangeBox *centroid, uint8 quadrant)
+{
+	RectBox	   *next_rect_box = (RectBox *) palloc(sizeof(RectBox));
+
+	memcpy(next_rect_box, rect_box, sizeof(RectBox));
+
+	if (quadrant & 0x8)
+		next_rect_box->range_box_x.left.low = centroid->left.low;
+	else
+		next_rect_box->range_box_x.left.high = centroid->left.low;
+
+	if (quadrant & 0x4)
+		next_rect_box->range_box_x.right.low = centroid->left.high;
+	else
+		next_rect_box->range_box_x.right.high = centroid->left.high;
+
+	if (quadrant & 0x2)
+		next_rect_box->range_box_y.left.low = centroid->right.low;
+	else
+		next_rect_box->range_box_y.left.high = centroid->right.low;
+
+	if (quadrant & 0x1)
+		next_rect_box->range_box_y.right.low = centroid->right.high;
+	else
+		next_rect_box->range_box_y.right.high = centroid->right.high;
+
+	return next_rect_box;
+}
+
+/* Can any range from range_box overlap with this argument? */
+static bool
+overlap2D(RangeBox *range_box, Range *query)
+{
+	return FPge(range_box->right.high, query->low) &&
+		   FPle(range_box->left.low, query->high);
+}
+
+/* Can any rectangle from rect_box overlap with this argument? */
+static bool
+overlap4D(RectBox *rect_box, RangeBox *query)
+{
+	return overlap2D(&rect_box->range_box_x, &query->left) &&
+		   overlap2D(&rect_box->range_box_y, &query->right);
+}
+
+/* Can any range from range_box contain this argument? */
+static bool
+contain2D(RangeBox *range_box, Range *query)
+{
+	return FPge(range_box->right.high, query->high) &&
+		   FPle(range_box->left.low, query->low);
+}
+
+/* Can any rectangle from rect_box contain this argument? */
+static bool
+contain4D(RectBox *rect_box, RangeBox * query)
+{
+	return contain2D(&rect_box->range_box_x, &query->left) &&
+		   contain2D(&rect_box->range_box_y, &query->right);
+}
+
+/* Can any range from range_box be contained by this argument? */
+static bool
+contained2D(RangeBox *range_box, Range *query)
+{
+	return FPle(range_box->left.low, query->high) &&
+		   FPge(range_box->left.high, query->low) &&
+		   FPle(range_box->right.low, query->high) &&
+		   FPge(range_box->right.high, query->low);
+}
+
+/* Can any rectangle from rect_box be contained by this argument? */
+static bool
+contained4D(RectBox *rect_box, RangeBox *query)
+{
+	return contained2D(&rect_box->range_box_x, &query->left) &&
+		   contained2D(&rect_box->range_box_y, &query->right);
+}
+
+/* Can any range from range_box to be lower than this argument? */
+static bool
+lower2D(RangeBox *range_box, Range *query)
+{
+	return FPlt(range_box->left.low, query->low) &&
+		   FPlt(range_box->right.low, query->low);
+}
+
+/* Can any range from range_box to be higher than this argument? */
+static bool
+higher2D(RangeBox *range_box, Range *query)
+{
+	return FPgt(range_box->left.high, query->high) &&
+		   FPgt(range_box->right.high, query->high);
+}
+
+/* Can any rectangle from rect_box be left of this argument? */
+static bool
+left4D(RectBox *rect_box, RangeBox *query)
+{
+	return lower2D(&rect_box->range_box_x, &query->left);
+}
+
+/* Can any rectangle from rect_box does not extend the right of this argument? */
+static bool
+overLeft4D(RectBox *rect_box, RangeBox *query)
+{
+	return lower2D(&rect_box->range_box_x, &query->right);
+}
+
+/* Can any rectangle from rect_box be right of this argument? */
+static bool
+right4D(RectBox *rect_box, RangeBox *query)
+{
+	return higher2D(&rect_box->range_box_x, &query->left);
+}
+
+/* Can any rectangle from rect_box does not extend the left of this argument? */
+static bool
+overRight4D(RectBox *rect_box, RangeBox *query)
+{
+	return higher2D(&rect_box->range_box_x, &query->right);
+}
+
+/* Can any rectangle from rect_box be below of this argument? */
+static bool
+below4D(RectBox *rect_box, RangeBox *query)
+{
+	return lower2D(&rect_box->range_box_y, &query->right);
+}
+
+/* Can any rectangle from rect_box does not extend above this argument? */
+static bool
+overBelow4D(RectBox *rect_box, RangeBox *query)
+{
+	return lower2D(&rect_box->range_box_y, &query->left);
+}
+
+/* Can any rectangle from rect_box be above of this argument? */
+static bool
+above4D(RectBox *rect_box, RangeBox *query)
+{
+	return higher2D(&rect_box->range_box_y, &query->right);
+}
+
+/* Can any rectangle from rect_box does not extend below of this argument? */
+static bool
+overAbove4D(RectBox *rect_box, RangeBox *query)
+{
+	return higher2D(&rect_box->range_box_y, &query->right);
+}
+
+/*
+ * SP-GiST config function
+ */
+Datum
+spg_box_quad_config(PG_FUNCTION_ARGS)
+{
+	spgConfigOut *cfg = (spgConfigOut *) PG_GETARG_POINTER(1);
+
+	cfg->prefixType = BOXOID;
+	cfg->labelType = VOIDOID;	/* We don't need node labels. */
+	cfg->canReturnData = true;
+	cfg->longValuesOK = false;
+
+	PG_RETURN_VOID();
+}
+
+/*
+ * SP-GiST choose function
+ */
+Datum
+spg_box_quad_choose(PG_FUNCTION_ARGS)
+{
+	spgChooseIn *in = (spgChooseIn *) PG_GETARG_POINTER(0);
+	spgChooseOut *out = (spgChooseOut *) PG_GETARG_POINTER(1);
+	BOX		   *centroid = DatumGetBoxP(in->prefixDatum),
+			   *box = DatumGetBoxP(in->datum);
+
+	out->resultType = spgMatchNode;
+	out->result.matchNode.restDatum = BoxPGetDatum(box);
+
+	/* nodeN will be set by core, when allTheSame. */
+	if (!in->allTheSame)
+		out->result.matchNode.nodeN = getQuadrant(centroid, box);
+
+	PG_RETURN_VOID();
+}
+
+/*
+ * SP-GiST pick-split function
+ *
+ * It splits a list of boxes into quadrants by choosing a central 4D
+ * point as the median of the coordinates of the boxes.
+ */
+Datum
+spg_box_quad_picksplit(PG_FUNCTION_ARGS)
+{
+	spgPickSplitIn	*in = (spgPickSplitIn *) PG_GETARG_POINTER(0);
+	spgPickSplitOut	*out = (spgPickSplitOut *) PG_GETARG_POINTER(1);
+	BOX		   *centroid;
+	int			median,
+				i;
+	double	   *lowXs = palloc(sizeof(double) * in->nTuples);
+	double	   *highXs = palloc(sizeof(double) * in->nTuples);
+	double	   *lowYs = palloc(sizeof(double) * in->nTuples);
+	double	   *highYs = palloc(sizeof(double) * in->nTuples);
+
+	/* Calculate median of all 4D coordinates */
+	for (i = 0; i < in->nTuples; i++)
+	{
+		BOX  *box = DatumGetBoxP(in->datums[i]);
+
+		lowXs[i] = box->low.x;
+		highXs[i] = box->high.x;
+		lowYs[i] = box->low.y;
+		highYs[i] = box->high.y;
+	}
+
+	qsort(lowXs, in->nTuples, sizeof(double), compareDoubles);
+	qsort(highXs, in->nTuples, sizeof(double), compareDoubles);
+	qsort(lowYs, in->nTuples, sizeof(double), compareDoubles);
+	qsort(highYs, in->nTuples, sizeof(double), compareDoubles);
+
+	median = in->nTuples / 2;
+
+	centroid = palloc(sizeof(BOX));
+
+	centroid->low.x = lowXs[median];
+	centroid->high.x = highXs[median];
+	centroid->low.y = lowYs[median];
+	centroid->high.y = highYs[median];
+
+	/* Fill the output */
+	out->hasPrefix = true;
+	out->prefixDatum = BoxPGetDatum(centroid);
+
+	out->nNodes = 16;
+	out->nodeLabels = NULL;		/* We don't need node labels. */
+
+	out->mapTuplesToNodes = palloc(sizeof(int) * in->nTuples);
+	out->leafTupleDatums = palloc(sizeof(Datum) * in->nTuples);
+
+	/*
+	 * Assign ranges to corresponding nodes according to quadrants
+	 * relative to the "centroid" range
+	 */
+	for (i = 0; i < in->nTuples; i++)
+	{
+		BOX  *box = DatumGetBoxP(in->datums[i]);
+		uint8 quadrant = getQuadrant(centroid, box);
+
+		out->leafTupleDatums[i] = BoxPGetDatum(box);
+		out->mapTuplesToNodes[i] = quadrant;
+	}
+
+	PG_RETURN_VOID();
+}
+
+/*
+ * SP-GiST inner consistent function
+ */
+Datum
+spg_box_quad_inner_consistent(PG_FUNCTION_ARGS)
+{
+	spgInnerConsistentIn *in = (spgInnerConsistentIn *) PG_GETARG_POINTER(0);
+	spgInnerConsistentOut *out = (spgInnerConsistentOut *) PG_GETARG_POINTER(1);
+	int				i;
+	MemoryContext	old_ctx;
+	RectBox		   *rect_box;
+	uint8			quadrant;
+	RangeBox	   *centroid,
+				  **queries;
+
+	if (in->allTheSame)
+	{
+		/* Report that all nodes should be visited */
+		out->nNodes = in->nNodes;
+		out->nodeNumbers = (int *) palloc(sizeof(int) * in->nNodes);
+		for (i = 0; i < in->nNodes; i++)
+			out->nodeNumbers[i] = i;
+
+		PG_RETURN_VOID();
+	}
+
+	/*
+	 * We are saving the traversal value or initialize it an unbounded
+	 * one, if we have just begun to walk the tree.
+	 */
+	if (in->traversalValue)
+		rect_box = in->traversalValue;
+	else
+		rect_box = initRectBox();
+
+	/*
+	 * We are casting the prefix and queries to RangeBoxes for ease of
+	 * the following operations.
+	 */
+	centroid = getRangeBox(DatumGetBoxP(in->prefixDatum));
+	queries = (RangeBox **) palloc(in->nkeys * sizeof(RangeBox *));
+	for (i = 0; i < in->nkeys; i++)
+		queries[i] = getRangeBox(DatumGetBoxP(in->scankeys[i].sk_argument));
+
+	/* Allocate enough memory for nodes */
+	out->nNodes = 0;
+	out->nodeNumbers = (int *) palloc(sizeof(int) * in->nNodes);
+	out->traversalValues = (void **) palloc(sizeof(void *) * in->nNodes);
+
+	/*
+	 * We switch memory context, because we want to allocate memory for
+	 * new traversal values (next_rect_box) and pass these pieces of
+	 * memory to further call of this function.
+	 */
+	old_ctx = MemoryContextSwitchTo(in->traversalMemoryContext);
+
+	for (quadrant = 0; quadrant < in->nNodes; quadrant++)
+	{
+		bool		flag;
+		RectBox	   *next_rect_box = nextRectBox(rect_box, centroid, quadrant);
+
+		for (i = 0; i < in->nkeys; i++)
+		{
+			StrategyNumber strategy = in->scankeys[i].sk_strategy;
+
+			switch (strategy)
+			{
+				case RTOverlapStrategyNumber:
+					flag = overlap4D(next_rect_box, queries[i]);
+					break;
+
+				case RTContainsStrategyNumber:
+					flag = contain4D(next_rect_box, queries[i]);
+					break;
+
+				case RTSameStrategyNumber:
+				case RTContainedByStrategyNumber:
+					flag = contained4D(next_rect_box, queries[i]);
+					break;
+
+				case RTLeftStrategyNumber:
+					flag = left4D(next_rect_box, queries[i]);
+					break;
+
+				case RTOverLeftStrategyNumber:
+					flag = overLeft4D(next_rect_box, queries[i]);
+					break;
+
+				case RTRightStrategyNumber:
+					flag = right4D(next_rect_box, queries[i]);
+					break;
+
+				case RTOverRightStrategyNumber:
+					flag = overRight4D(next_rect_box, queries[i]);
+					break;
+
+				case RTAboveStrategyNumber:
+					flag = above4D(next_rect_box, queries[i]);
+					break;
+
+				case RTOverAboveStrategyNumber:
+					flag = overAbove4D(next_rect_box, queries[i]);
+					break;
+
+				case RTBelowStrategyNumber:
+					flag = below4D(next_rect_box, queries[i]);
+					break;
+
+				case RTOverBelowStrategyNumber:
+					flag = overBelow4D(next_rect_box, queries[i]);
+					break;
+
+				default:
+					elog(ERROR, "unrecognized strategy: %d", strategy);
+			}
+
+			/* If any check is failed, we have found our answer. */
+			if (!flag)
+				break;
+		}
+
+		if (flag)
+		{
+			out->traversalValues[out->nNodes] = next_rect_box;
+			out->nodeNumbers[out->nNodes] = quadrant;
+			out->nNodes++;
+		}
+		else
+		{
+			/*
+			 * If this node is not selected, we don't need to keep
+			 * the next traversal value in the memory context.
+			 */
+			pfree(next_rect_box);
+		}
+	}
+
+	/* Switch back */
+	MemoryContextSwitchTo(old_ctx);
+
+	PG_RETURN_VOID();
+}
+
+/*
+ * SP-GiST inner consistent function
+ */
+Datum
+spg_box_quad_leaf_consistent(PG_FUNCTION_ARGS)
+{
+	spgLeafConsistentIn *in = (spgLeafConsistentIn *) PG_GETARG_POINTER(0);
+	spgLeafConsistentOut *out = (spgLeafConsistentOut *) PG_GETARG_POINTER(1);
+	Datum		leaf = in->leafDatum;
+	bool		flag;
+	int			i;
+
+	/* All tests are exact. */
+	out->recheck = false;
+
+	/* leafDatum is what it is... */
+	out->leafValue = in->leafDatum;
+
+	/* Perform the required comparison(s) */
+	for (i = 0; i < in->nkeys; i++)
+	{
+		StrategyNumber strategy = in->scankeys[i].sk_strategy;
+		Datum		query = in->scankeys[i].sk_argument;
+
+		switch (strategy)
+		{
+			case RTOverlapStrategyNumber:
+				flag = DatumGetBool(DirectFunctionCall2(box_overlap, leaf,
+														query));
+				break;
+
+			case RTContainsStrategyNumber:
+				flag = DatumGetBool(DirectFunctionCall2(box_contain, leaf,
+														query));
+				break;
+
+			case RTContainedByStrategyNumber:
+				flag = DatumGetBool(DirectFunctionCall2(box_contained, leaf,
+														query));
+				break;
+
+			case RTSameStrategyNumber:
+				flag = DatumGetBool(DirectFunctionCall2(box_same, leaf,
+														query));
+				break;
+
+			case RTLeftStrategyNumber:
+				flag = DatumGetBool(DirectFunctionCall2(box_left, leaf,
+														query));
+				break;
+
+			case RTOverLeftStrategyNumber:
+				flag = DatumGetBool(DirectFunctionCall2(box_overleft, leaf,
+														query));
+				break;
+
+			case RTRightStrategyNumber:
+				flag = DatumGetBool(DirectFunctionCall2(box_right, leaf,
+														query));
+				break;
+
+			case RTOverRightStrategyNumber:
+				flag = DatumGetBool(DirectFunctionCall2(box_overright, leaf,
+														query));
+				break;
+
+			case RTAboveStrategyNumber:
+				flag = DatumGetBool(DirectFunctionCall2(box_above, leaf,
+														query));
+				break;
+
+			case RTOverAboveStrategyNumber:
+				flag = DatumGetBool(DirectFunctionCall2(box_overabove, leaf,
+														query));
+				break;
+
+			case RTBelowStrategyNumber:
+				flag = DatumGetBool(DirectFunctionCall2(box_below, leaf,
+														query));
+				break;
+
+			case RTOverBelowStrategyNumber:
+				flag = DatumGetBool(DirectFunctionCall2(box_overbelow, leaf,
+														query));
+				break;
+
+			default:
+				elog(ERROR, "unrecognized strategy: %d", strategy);
+		}
+
+		/* If any check is failed, we have found our answer. */
+		if (!flag)
+			break;
+	}
+
+	PG_RETURN_BOOL(flag);
+}
diff --git a/src/include/catalog/pg_amop.h b/src/include/catalog/pg_amop.h
index c642c3927e..a15b0ec309 100644
--- a/src/include/catalog/pg_amop.h
+++ b/src/include/catalog/pg_amop.h
@@ -832,6 +832,22 @@ DATA(insert (	3474   3831 3831 8 s	3892 4000 0 ));
 DATA(insert (	3474   3831 2283 16 s	3889 4000 0 ));
 DATA(insert (	3474   3831 3831 18 s	3882 4000 0 ));
 
+/*
+ * SP-GiST box_ops
+ */
+DATA(insert (	5000	603  603  1 s	 493	4000 0 ));
+DATA(insert (	5000	603  603  2 s	 494	4000 0 ));
+DATA(insert (	5000	603  603  3 s	 500	4000 0 ));
+DATA(insert (	5000	603  603  4 s	 495	4000 0 ));
+DATA(insert (	5000	603  603  5 s	 496	4000 0 ));
+DATA(insert (	5000	603  603  6 s	 499	4000 0 ));
+DATA(insert (	5000	603  603  7 s	 498	4000 0 ));
+DATA(insert (	5000	603  603  8 s	 497	4000 0 ));
+DATA(insert (	5000	603  603  9 s	2571	4000 0 ));
+DATA(insert (	5000	603  603 10 s	2570	4000 0 ));
+DATA(insert (	5000	603  603 11 s	2573	4000 0 ));
+DATA(insert (	5000	603  603 12 s	2572	4000 0 ));
+
 /*
  * GiST inet_ops
  */
diff --git a/src/include/catalog/pg_amproc.h b/src/include/catalog/pg_amproc.h
index f0ae008704..00320b4c33 100644
--- a/src/include/catalog/pg_amproc.h
+++ b/src/include/catalog/pg_amproc.h
@@ -443,6 +443,11 @@ DATA(insert (	4017   25 25 2 4028 ));
 DATA(insert (	4017   25 25 3 4029 ));
 DATA(insert (	4017   25 25 4 4030 ));
 DATA(insert (	4017   25 25 5 4031 ));
+DATA(insert (	5000   603 603 1 5012 ));
+DATA(insert (	5000   603 603 2 5013 ));
+DATA(insert (	5000   603 603 3 5014 ));
+DATA(insert (	5000   603 603 4 5015 ));
+DATA(insert (	5000   603 603 5 5016 ));
 
 /* BRIN opclasses */
 /* minmax bytea */
diff --git a/src/include/catalog/pg_opclass.h b/src/include/catalog/pg_opclass.h
index a9446b7ec9..b564046dea 100644
--- a/src/include/catalog/pg_opclass.h
+++ b/src/include/catalog/pg_opclass.h
@@ -228,6 +228,7 @@ DATA(insert (	403		range_ops			PGNSP PGUID 3901  3831 t 0 ));
 DATA(insert (	405		range_ops			PGNSP PGUID 3903  3831 t 0 ));
 DATA(insert (	783		range_ops			PGNSP PGUID 3919  3831 t 0 ));
 DATA(insert (	4000	range_ops			PGNSP PGUID 3474  3831 t 0 ));
+DATA(insert (	4000	box_ops			    PGNSP PGUID 5000  603  t 0 ));
 DATA(insert (	4000	quad_point_ops		PGNSP PGUID 4015  600 t 0 ));
 DATA(insert (	4000	kd_point_ops		PGNSP PGUID 4016  600 f 0 ));
 DATA(insert (	4000	text_ops			PGNSP PGUID 4017  25 t 0 ));
diff --git a/src/include/catalog/pg_opfamily.h b/src/include/catalog/pg_opfamily.h
index fa444469cd..1499a502f3 100644
--- a/src/include/catalog/pg_opfamily.h
+++ b/src/include/catalog/pg_opfamily.h
@@ -182,5 +182,6 @@ DATA(insert OID = 4081 (	3580	uuid_minmax_ops			PGNSP PGUID ));
 DATA(insert OID = 4103 (	3580	range_inclusion_ops		PGNSP PGUID ));
 DATA(insert OID = 4082 (	3580	pg_lsn_minmax_ops		PGNSP PGUID ));
 DATA(insert OID = 4104 (	3580	box_inclusion_ops		PGNSP PGUID ));
+DATA(insert OID = 5000 (	4000	box_ops		PGNSP PGUID ));
 
 #endif   /* PG_OPFAMILY_H */
diff --git a/src/include/catalog/pg_proc.h b/src/include/catalog/pg_proc.h
index 7619c40eb3..c86b920cb6 100644
--- a/src/include/catalog/pg_proc.h
+++ b/src/include/catalog/pg_proc.h
@@ -5097,6 +5097,17 @@ DESCR("SP-GiST support for quad tree over range");
 DATA(insert OID = 3473 (  spg_range_quad_leaf_consistent	PGNSP PGUID 12 1 0 0 0 f f f f t f i s 2 0 16 "2281 2281" _null_ _null_ _null_ _null_  _null_ spg_range_quad_leaf_consistent _null_ _null_ _null_ ));
 DESCR("SP-GiST support for quad tree over range");
 
+DATA(insert OID = 5012 (  spg_box_quad_config PGNSP PGUID 12 1 0 0 0 f f f f t f i s 2 0 2278 "2281 2281" _null_ _null_ _null_ _null_  _null_ spg_box_quad_config _null_ _null_ _null_ ));
+DESCR("SP-GiST support for quad tree over box");
+DATA(insert OID = 5013 (  spg_box_quad_choose PGNSP PGUID 12 1 0 0 0 f f f f t f i s 2 0 2278 "2281 2281" _null_ _null_ _null_ _null_  _null_ spg_box_quad_choose _null_ _null_ _null_ ));
+DESCR("SP-GiST support for quad tree over box");
+DATA(insert OID = 5014 (  spg_box_quad_picksplit	PGNSP PGUID 12 1 0 0 0 f f f f t f i s 2 0 2278 "2281 2281" _null_ _null_ _null_ _null_  _null_ spg_box_quad_picksplit _null_ _null_ _null_ ));
+DESCR("SP-GiST support for quad tree over box");
+DATA(insert OID = 5015 (  spg_box_quad_inner_consistent	PGNSP PGUID 12 1 0 0 0 f f f f t f i s 2 0 2278 "2281 2281" _null_ _null_ _null_ _null_  _null_ spg_box_quad_inner_consistent _null_ _null_ _null_ ));
+DESCR("SP-GiST support for quad tree over box");
+DATA(insert OID = 5016 (  spg_box_quad_leaf_consistent	PGNSP PGUID 12 1 0 0 0 f f f f t f i s 2 0 16 "2281 2281" _null_ _null_ _null_ _null_  _null_ spg_box_quad_leaf_consistent _null_ _null_ _null_ ));
+DESCR("SP-GiST support for quad tree over box");
+
 /* replication slots */
 DATA(insert OID = 3779 (  pg_create_physical_replication_slot PGNSP PGUID 12 1 0 0 0 f f f f t f v u 2 0 2249 "19 16" "{19,16,19,3220}" "{i,i,o,o}" "{slot_name,immediately_reserve,slot_name,xlog_position}" _null_ _null_ pg_create_physical_replication_slot _null_ _null_ _null_ ));
 DESCR("create a physical replication slot");
diff --git a/src/include/utils/geo_decls.h b/src/include/utils/geo_decls.h
index 9d8d660d15..acf320207c 100644
--- a/src/include/utils/geo_decls.h
+++ b/src/include/utils/geo_decls.h
@@ -426,6 +426,12 @@ extern Datum gist_point_consistent(PG_FUNCTION_ARGS);
 extern Datum gist_point_distance(PG_FUNCTION_ARGS);
 extern Datum gist_point_fetch(PG_FUNCTION_ARGS);
 
+/* utils/adt/geo_spgist.c */
+Datum spg_box_quad_config(PG_FUNCTION_ARGS);
+Datum spg_box_quad_choose(PG_FUNCTION_ARGS);
+Datum spg_box_quad_picksplit(PG_FUNCTION_ARGS);
+Datum spg_box_quad_inner_consistent(PG_FUNCTION_ARGS);
+Datum spg_box_quad_leaf_consistent(PG_FUNCTION_ARGS);
 
 /* geo_selfuncs.c */
 extern Datum areasel(PG_FUNCTION_ARGS);
diff --git a/src/test/regress/expected/box.out b/src/test/regress/expected/box.out
index fc3154014b..300190f762 100644
--- a/src/test/regress/expected/box.out
+++ b/src/test/regress/expected/box.out
@@ -215,3 +215,243 @@ SELECT '' AS four, height(f1), width(f1) FROM BOX_TBL;
       |      0 |     0
 (4 rows)
 
+--
+-- Test the SP-GiST index
+--
+CREATE TEMPORARY TABLE box_temp (f1 box);
+INSERT INTO box_temp
+	SELECT box(point(i, i), point(i * 2, i * 2))
+	FROM generate_series(1, 50) AS i;
+CREATE INDEX box_spgist ON box_temp USING spgist (f1);
+INSERT INTO box_temp
+	VALUES (NULL),
+		   ('(-0,0)(0,100)'),
+		   ('(-3,4.3333333333)(40,1)'),
+		   ('(0,100)(0,infinity)'),
+		   ('(-infinity,0)(0,infinity)'),
+		   ('(-infinity,-infinity)(infinity,infinity)');
+SET enable_seqscan = false;
+SELECT * FROM box_temp WHERE f1 << '(10,20),(30,40)';
+        f1        
+------------------
+ (2,2),(1,1)
+ (4,4),(2,2)
+ (6,6),(3,3)
+ (8,8),(4,4)
+ (-0,100),(0,0)
+ (0,inf),(0,100)
+ (0,inf),(-inf,0)
+(7 rows)
+
+EXPLAIN (COSTS OFF) SELECT * FROM box_temp WHERE f1 << '(10,20),(30,40)';
+                  QUERY PLAN                  
+----------------------------------------------
+ Index Only Scan using box_spgist on box_temp
+   Index Cond: (f1 << '(30,40),(10,20)'::box)
+(2 rows)
+
+SELECT * FROM box_temp WHERE f1 &< '(10,4.333334),(5,100)';
+        f1        
+------------------
+ (2,2),(1,1)
+ (4,4),(2,2)
+ (6,6),(3,3)
+ (8,8),(4,4)
+ (10,10),(5,5)
+ (-0,100),(0,0)
+ (0,inf),(0,100)
+ (0,inf),(-inf,0)
+(8 rows)
+
+EXPLAIN (COSTS OFF) SELECT * FROM box_temp WHERE f1 &< '(10,4.333334),(5,100)';
+                     QUERY PLAN                     
+----------------------------------------------------
+ Index Only Scan using box_spgist on box_temp
+   Index Cond: (f1 &< '(10,100),(5,4.333334)'::box)
+(2 rows)
+
+SELECT * FROM box_temp WHERE f1 && '(15,20),(25,30)';
+          f1           
+-----------------------
+ (20,20),(10,10)
+ (22,22),(11,11)
+ (24,24),(12,12)
+ (26,26),(13,13)
+ (28,28),(14,14)
+ (30,30),(15,15)
+ (32,32),(16,16)
+ (34,34),(17,17)
+ (36,36),(18,18)
+ (38,38),(19,19)
+ (40,40),(20,20)
+ (42,42),(21,21)
+ (44,44),(22,22)
+ (46,46),(23,23)
+ (48,48),(24,24)
+ (50,50),(25,25)
+ (inf,inf),(-inf,-inf)
+(17 rows)
+
+EXPLAIN (COSTS OFF) SELECT * FROM box_temp WHERE f1 && '(15,20),(25,30)';
+                  QUERY PLAN                  
+----------------------------------------------
+ Index Only Scan using box_spgist on box_temp
+   Index Cond: (f1 && '(25,30),(15,20)'::box)
+(2 rows)
+
+SELECT * FROM box_temp WHERE f1 &> '(40,30),(45,50)';
+        f1         
+-------------------
+ (80,80),(40,40)
+ (82,82),(41,41)
+ (84,84),(42,42)
+ (86,86),(43,43)
+ (88,88),(44,44)
+ (90,90),(45,45)
+ (92,92),(46,46)
+ (94,94),(47,47)
+ (96,96),(48,48)
+ (98,98),(49,49)
+ (100,100),(50,50)
+(11 rows)
+
+EXPLAIN (COSTS OFF) SELECT * FROM box_temp WHERE f1 &> '(40,30),(45,50)';
+                  QUERY PLAN                  
+----------------------------------------------
+ Index Only Scan using box_spgist on box_temp
+   Index Cond: (f1 &> '(45,50),(40,30)'::box)
+(2 rows)
+
+SELECT * FROM box_temp WHERE f1 >> '(30,40),(40,30)';
+        f1         
+-------------------
+ (82,82),(41,41)
+ (84,84),(42,42)
+ (86,86),(43,43)
+ (88,88),(44,44)
+ (90,90),(45,45)
+ (92,92),(46,46)
+ (94,94),(47,47)
+ (96,96),(48,48)
+ (98,98),(49,49)
+ (100,100),(50,50)
+(10 rows)
+
+EXPLAIN (COSTS OFF) SELECT * FROM box_temp WHERE f1 >> '(30,40),(40,30)';
+                  QUERY PLAN                  
+----------------------------------------------
+ Index Only Scan using box_spgist on box_temp
+   Index Cond: (f1 >> '(40,40),(30,30)'::box)
+(2 rows)
+
+SELECT * FROM box_temp WHERE f1 <<| '(10,4.33334),(5,100)';
+            f1            
+--------------------------
+ (2,2),(1,1)
+ (4,4),(2,2)
+ (40,4.3333333333),(-3,1)
+(3 rows)
+
+EXPLAIN (COSTS OFF) SELECT * FROM box_temp WHERE f1 <<| '(10,4.33334),(5,100)';
+                     QUERY PLAN                     
+----------------------------------------------------
+ Index Only Scan using box_spgist on box_temp
+   Index Cond: (f1 <<| '(10,100),(5,4.33334)'::box)
+(2 rows)
+
+SELECT * FROM box_temp WHERE f1 &<| '(10,4.3333334),(5,1)';
+            f1            
+--------------------------
+ (2,2),(1,1)
+ (4,4),(2,2)
+ (40,4.3333333333),(-3,1)
+(3 rows)
+
+EXPLAIN (COSTS OFF) SELECT * FROM box_temp WHERE f1 &<| '(10,4.3333334),(5,1)';
+                     QUERY PLAN                     
+----------------------------------------------------
+ Index Only Scan using box_spgist on box_temp
+   Index Cond: (f1 &<| '(10,4.3333334),(5,1)'::box)
+(2 rows)
+
+SELECT * FROM box_temp WHERE f1 |&> '(49.99,49.99),(49.99,49.99)';
+        f1         
+-------------------
+ (100,100),(50,50)
+ (0,inf),(0,100)
+(2 rows)
+
+EXPLAIN (COSTS OFF) SELECT * FROM box_temp WHERE f1 |&> '(49.99,49.99),(49.99,49.99)';
+                        QUERY PLAN                         
+-----------------------------------------------------------
+ Index Only Scan using box_spgist on box_temp
+   Index Cond: (f1 |&> '(49.99,49.99),(49.99,49.99)'::box)
+(2 rows)
+
+SELECT * FROM box_temp WHERE f1 |>> '(37,38),(39,40)';
+        f1         
+-------------------
+ (82,82),(41,41)
+ (84,84),(42,42)
+ (86,86),(43,43)
+ (88,88),(44,44)
+ (90,90),(45,45)
+ (92,92),(46,46)
+ (94,94),(47,47)
+ (96,96),(48,48)
+ (98,98),(49,49)
+ (100,100),(50,50)
+ (0,inf),(0,100)
+(11 rows)
+
+EXPLAIN (COSTS OFF) SELECT * FROM box_temp WHERE f1 |>> '(37,38),(39,40)';
+                  QUERY PLAN                   
+-----------------------------------------------
+ Index Only Scan using box_spgist on box_temp
+   Index Cond: (f1 |>> '(39,40),(37,38)'::box)
+(2 rows)
+
+SELECT * FROM box_temp WHERE f1 @> '(10,11),(15,16)';
+          f1           
+-----------------------
+ (16,16),(8,8)
+ (18,18),(9,9)
+ (20,20),(10,10)
+ (inf,inf),(-inf,-inf)
+(4 rows)
+
+EXPLAIN (COSTS OFF) SELECT * FROM box_temp WHERE f1 @> '(10,11),(15,15)';
+                  QUERY PLAN                  
+----------------------------------------------
+ Index Only Scan using box_spgist on box_temp
+   Index Cond: (f1 @> '(15,15),(10,11)'::box)
+(2 rows)
+
+SELECT * FROM box_temp WHERE f1 <@ '(10,15),(30,35)';
+       f1        
+-----------------
+ (30,30),(15,15)
+(1 row)
+
+EXPLAIN (COSTS OFF) SELECT * FROM box_temp WHERE f1 <@ '(10,15),(30,35)';
+                  QUERY PLAN                  
+----------------------------------------------
+ Index Only Scan using box_spgist on box_temp
+   Index Cond: (f1 <@ '(30,35),(10,15)'::box)
+(2 rows)
+
+SELECT * FROM box_temp WHERE f1 ~= '(20,20),(40,40)';
+       f1        
+-----------------
+ (40,40),(20,20)
+(1 row)
+
+EXPLAIN (COSTS OFF) SELECT * FROM box_temp WHERE f1 ~= '(20,20),(40,40)';
+                  QUERY PLAN                  
+----------------------------------------------
+ Index Only Scan using box_spgist on box_temp
+   Index Cond: (f1 ~= '(40,40),(20,20)'::box)
+(2 rows)
+
+RESET enable_seqscan;
+DROP INDEX box_spgist;
diff --git a/src/test/regress/expected/opr_sanity.out b/src/test/regress/expected/opr_sanity.out
index 7c09fa3d59..a8dc0c150f 100644
--- a/src/test/regress/expected/opr_sanity.out
+++ b/src/test/regress/expected/opr_sanity.out
@@ -1728,15 +1728,19 @@ ORDER BY 1, 2, 3;
        4000 |            6 | ~=
        4000 |            7 | @>
        4000 |            8 | <@
+       4000 |            9 | &<|
+       4000 |           10 | <<|
        4000 |           10 | <^
        4000 |           11 | <
        4000 |           11 | >^
+       4000 |           11 | |>>
        4000 |           12 | <=
+       4000 |           12 | |&>
        4000 |           14 | >=
        4000 |           15 | >
        4000 |           16 | @>
        4000 |           18 | =
-(108 rows)
+(112 rows)
 
 -- Check that all opclass search operators have selectivity estimators.
 -- This is not absolutely required, but it seems a reasonable thing
diff --git a/src/test/regress/sql/box.sql b/src/test/regress/sql/box.sql
index 234c2f28a3..1c9e52e0a3 100644
--- a/src/test/regress/sql/box.sql
+++ b/src/test/regress/sql/box.sql
@@ -117,3 +117,65 @@ SELECT '' AS one, b1.*, b2.*
    WHERE b1.f1 @> b2.f1 and not b1.f1 ~= b2.f1;
 
 SELECT '' AS four, height(f1), width(f1) FROM BOX_TBL;
+
+--
+-- Test the SP-GiST index
+--
+
+CREATE TEMPORARY TABLE box_temp (f1 box);
+
+INSERT INTO box_temp
+	SELECT box(point(i, i), point(i * 2, i * 2))
+	FROM generate_series(1, 50) AS i;
+
+CREATE INDEX box_spgist ON box_temp USING spgist (f1);
+
+INSERT INTO box_temp
+	VALUES (NULL),
+		   ('(-0,0)(0,100)'),
+		   ('(-3,4.3333333333)(40,1)'),
+		   ('(0,100)(0,infinity)'),
+		   ('(-infinity,0)(0,infinity)'),
+		   ('(-infinity,-infinity)(infinity,infinity)');
+
+SET enable_seqscan = false;
+
+SELECT * FROM box_temp WHERE f1 << '(10,20),(30,40)';
+EXPLAIN (COSTS OFF) SELECT * FROM box_temp WHERE f1 << '(10,20),(30,40)';
+
+SELECT * FROM box_temp WHERE f1 &< '(10,4.333334),(5,100)';
+EXPLAIN (COSTS OFF) SELECT * FROM box_temp WHERE f1 &< '(10,4.333334),(5,100)';
+
+SELECT * FROM box_temp WHERE f1 && '(15,20),(25,30)';
+EXPLAIN (COSTS OFF) SELECT * FROM box_temp WHERE f1 && '(15,20),(25,30)';
+
+SELECT * FROM box_temp WHERE f1 &> '(40,30),(45,50)';
+EXPLAIN (COSTS OFF) SELECT * FROM box_temp WHERE f1 &> '(40,30),(45,50)';
+
+SELECT * FROM box_temp WHERE f1 >> '(30,40),(40,30)';
+EXPLAIN (COSTS OFF) SELECT * FROM box_temp WHERE f1 >> '(30,40),(40,30)';
+
+SELECT * FROM box_temp WHERE f1 <<| '(10,4.33334),(5,100)';
+EXPLAIN (COSTS OFF) SELECT * FROM box_temp WHERE f1 <<| '(10,4.33334),(5,100)';
+
+SELECT * FROM box_temp WHERE f1 &<| '(10,4.3333334),(5,1)';
+EXPLAIN (COSTS OFF) SELECT * FROM box_temp WHERE f1 &<| '(10,4.3333334),(5,1)';
+
+SELECT * FROM box_temp WHERE f1 |&> '(49.99,49.99),(49.99,49.99)';
+EXPLAIN (COSTS OFF) SELECT * FROM box_temp WHERE f1 |&> '(49.99,49.99),(49.99,49.99)';
+
+SELECT * FROM box_temp WHERE f1 |>> '(37,38),(39,40)';
+EXPLAIN (COSTS OFF) SELECT * FROM box_temp WHERE f1 |>> '(37,38),(39,40)';
+
+SELECT * FROM box_temp WHERE f1 @> '(10,11),(15,16)';
+EXPLAIN (COSTS OFF) SELECT * FROM box_temp WHERE f1 @> '(10,11),(15,15)';
+
+SELECT * FROM box_temp WHERE f1 <@ '(10,15),(30,35)';
+EXPLAIN (COSTS OFF) SELECT * FROM box_temp WHERE f1 <@ '(10,15),(30,35)';
+
+SELECT * FROM box_temp WHERE f1 ~= '(20,20),(40,40)';
+EXPLAIN (COSTS OFF) SELECT * FROM box_temp WHERE f1 ~= '(20,20),(40,40)';
+
+RESET enable_seqscan;
+
+DROP INDEX box_spgist;
-- 
2.40.0