From: Bborie Park <bkpark at ucdavis.edu>
Date: Thu, 22 Sep 2011 15:07:25 +0000 (+0000)
Subject: Further tweaks to rt_raster_gdal_rasterize to correctly handle auto-computed extents... 
X-Git-Tag: 2.0.0alpha1~970
X-Git-Url: https://granicus.if.org/sourcecode?a=commitdiff_plain;h=c5d7903487d5559315bf709d70f620b045804882;p=postgis

Further tweaks to rt_raster_gdal_rasterize to correctly handle auto-computed extents of multipoints.  Addition of ST_Intersects for two rasters.  Refactored the one raster and one geometry version of ST_Intersects.

Associated ticket is #1176


git-svn-id: http://svn.osgeo.org/postgis/trunk@7884 b70326c6-7e19-0410-871a-916f4a2858ee
---

diff --git a/raster/rt_core/rt_api.c b/raster/rt_core/rt_api.c
index c1608f9f6..0935fe6d8 100644
--- a/raster/rt_core/rt_api.c
+++ b/raster/rt_core/rt_api.c
@@ -768,29 +768,6 @@ rt_pixtype_name(rt_pixtype pixtype) {
 
 /*- rt_band ----------------------------------------------------------*/
 
-struct rt_extband_t {
-    uint8_t bandNum;
-    char* path; /* externally owned ? */
-};
-
-struct rt_band_t {
-    rt_pixtype pixtype;
-    int32_t offline;
-    uint16_t width;
-    uint16_t height;
-    int32_t hasnodata; /* a flag indicating if this band contains nodata values */
-    int32_t isnodata;   /* a flag indicating if this band is filled only with
-                           nodata values */
-    double nodataval; /* int will be converted ... */
-    int32_t ownsData; /* XXX mloskot: its behaviour needs to be documented */
-
-    union {
-        void* mem; /* actual data, externally owned */
-        struct rt_extband_t offline;
-    } data;
-
-};
-
 rt_band
 rt_band_new_inline(uint16_t width, uint16_t height,
         rt_pixtype pixtype, uint32_t hasnodata, double nodataval,
@@ -1051,7 +1028,9 @@ rt_band_get_isnodata_flag(rt_band band) {
 int
 rt_band_set_nodata(rt_band band, double val) {
     rt_pixtype pixtype = PT_END;
-    //double oldnodataval = band->nodataval;
+    /*
+		double oldnodataval = band->nodataval;
+		*/
 
     int32_t checkvalint = 0;
     uint32_t checkvaluint = 0;
@@ -1146,7 +1125,7 @@ rt_band_set_nodata(rt_band band, double val) {
     RASTER_DEBUGF(3, "rt_band_set_nodata: band->nodataval = %f", band->nodataval);
 
 
-    // the nodata value was just set, so this band has NODATA
+    /* the nodata value was just set, so this band has NODATA */
     rt_band_set_hasnodata_flag(band, 1);
 
 #ifdef POSTGIS_RASTER_WARN_ON_TRUNCATION
@@ -4010,26 +3989,112 @@ rt_raster_generate_new_band(rt_raster raster, rt_pixtype pixtype,
     return index;
 }
 
-
+/**
+ * Convert an xr, yr raster point to an xw, yw point on map
+ *
+ * @param raster : the raster to get info from
+ * @param xr : the pixel's column
+ * @param yr : the pixel's row
+ * @param xw : output parameter, X ordinate of the geographical point
+ * @param yw : output parameter, Y ordinate of the geographical point
+ */
 void
 rt_raster_cell_to_geopoint(rt_raster raster,
-        double x, double y,
-        double* x1, double* y1) {
+	double xr, double yr,
+	double* xw, double* yw
+) {
+	assert(NULL != raster);
+	assert(NULL != xw);
+	assert(NULL != yw);
 
+	/* Six parameters affine transformation */
+	*xw = raster->scaleX * xr + raster->skewX * yr + raster->ipX;
+	*yw = raster->scaleY * yr + raster->skewY * xr + raster->ipY;
 
-    assert(NULL != raster);
-    assert(NULL != x1);
-    assert(NULL != y1);
+	RASTER_DEBUGF(3, "rt_raster_cell_to_geopoint(%g, %g) = %g, %g",
+		xr, yr, *xw, *yw);
 
-    /* Six parameters affine transformation */
-    *x1 = raster->scaleX * x + raster->skewX * y + raster->ipX;
-    *y1 = raster->scaleY * y + raster->skewY * x + raster->ipY;
+	/*
+	{
+		double gt[] = {0.0, 1.0, 0.0, 0.0, 0.0, 1.0};
+		double x1;
+		double y1;
+
+		gt[0] = rt_raster_get_x_offset(raster);
+		gt[1] = rt_raster_get_x_scale(raster);
+		gt[2] = rt_raster_get_x_skew(raster);
+		gt[3] = rt_raster_get_y_offset(raster);
+		gt[4] = rt_raster_get_y_skew(raster);
+		gt[5] = rt_raster_get_y_scale(raster);
+
+		GDALApplyGeoTransform(gt, xr, yr, &x1, &y1);
+		RASTER_DEBUGF(4, "GDALApplyGeoTransform for (%f, %f) = (%f, %f)",
+			xr, yr, x1, y1);
+	}
+	*/
+}
 
-    RASTER_DEBUGF(3, "rt_raster_cell_to_geopoint(%g,%g)", x, y);
-    RASTER_DEBUGF(3, " ipx/y:%g/%g", raster->ipX, raster->ipY);
-    RASTER_DEBUGF(3, "cell_to_geopoint: ipX:%g, ipY:%g, %g,%g -> %g,%g",
-            raster->ipX, raster->ipY, x, y, *x1, *y1);
+/**
+ * Convert an xw,yw map point to a xr,yr raster point
+ *
+ * @param raster : the raster to get info from
+ * @param xw : X ordinate of the geographical point
+ * @param yw : Y ordinate of the geographical point
+ * @param xr : output parameter, the pixel's column
+ * @param yr : output parameter, the pixel's row
+ *
+ * @return if zero, error occurred in function
+ */
+int
+rt_raster_geopoint_to_cell(rt_raster raster,
+	double xw, double yw,
+	double *xr, double *yr
+) {
+	double d;
+
+	assert(NULL != raster);
+	assert(NULL != xr);
+	assert(NULL != yr);
+
+	/* Six parameters affine transformation */
+	d = raster->skewX * raster->skewY - raster->scaleX * raster->scaleY;
+	if (FLT_EQ(d, 0)) {
+		rterror("rt_raster_geopoint_to_cell: Attempting to compute coordinates on a raster with scale equal to 0");
+		return 0;
+	}
+
+	*xr = (raster->skewX * (yw - raster->ipY) - raster->scaleY * (xw - raster->ipX)) / d;
+	*yr = (raster->scaleX * (yw - raster->ipY) + raster->skewY * (raster->ipX - xw)) / (-1 * d);
+
+	*xr = floor(*xr);
+	*yr = floor(*yr);
+
+	RASTER_DEBUGF(3, "rt_raster_geopoint_to_cell(%g, %g) = %g, %g",
+		xw, yw, *xr, *yr);
+
+	/*
+	{
+		double gt[] = {0.0, 1.0, 0.0, 0.0, 0.0, 1.0};
+		double igt[6];
+		double x1;
+		double y1;
+
+		gt[0] = rt_raster_get_x_offset(raster);
+		gt[1] = rt_raster_get_x_scale(raster);
+		gt[2] = rt_raster_get_x_skew(raster);
+		gt[3] = rt_raster_get_y_offset(raster);
+		gt[4] = rt_raster_get_y_skew(raster);
+		gt[5] = rt_raster_get_y_scale(raster);
 
+		GDALInvGeoTransform(gt, igt);
+
+		GDALApplyGeoTransform(igt, xw, yw, &x1, &y1);
+		RASTER_DEBUGF(4, "GDALApplyGeoTransform for (%f, %f) = (%f, %f)",
+			xw, yw, x1, y1);
+	}
+	*/
+
+	return 1;
 }
 
 rt_geomval
@@ -4270,7 +4335,7 @@ rt_raster_dump_as_wktpolygons(rt_raster raster, int nband, int * pnElements)
          * postgresql memory context to work with, and the memory created
          * for pszSrcText is created outside this context.
          **/
-         //rtdealloc(pszSrcText);
+        /*rtdealloc(pszSrcText);*/
         free(pszSrcText);
         pszSrcText = NULL;
 
@@ -4883,7 +4948,7 @@ rt_band_from_wkb(uint16_t width, uint16_t height,
 
     /* And we should check if the band is a nodata band */
     /* TODO: No!! This is too slow */
-    //rt_band_check_is_nodata(band);
+    /*rt_band_check_is_nodata(band);*/
 
     return band;
 }
@@ -5172,7 +5237,7 @@ rt_raster_to_wkb(rt_raster raster, uint32_t *wkbsize) {
         if (band->isnodata) *ptr |= BANDTYPE_FLAG_ISNODATA;
         ptr += 1;
 
-#if 0 // no padding required for WKB
+#if 0 /* no padding required for WKB */
         /* Add padding (if needed) */
         if (pixbytes > 1) {
             memset(ptr, '\0', pixbytes - 1);
@@ -5239,7 +5304,7 @@ rt_raster_to_wkb(rt_raster raster, uint32_t *wkbsize) {
                 return 0;
         }
 
-#if 0 // no padding for WKB
+#if 0 /* no padding for WKB */
         /* Consistency checking (ptr is pixbytes-aligned) */
         assert(!((uint64_t) ptr % pixbytes));
 #endif
@@ -5262,7 +5327,7 @@ rt_raster_to_wkb(rt_raster raster, uint32_t *wkbsize) {
             }
         }
 
-#if 0 // no padding for WKB
+#if 0 /* no padding for WKB */
         /* Pad up to 8-bytes boundary */
         while ((uint64_t) ptr % 8) {
             *ptr = 0;
@@ -5361,7 +5426,7 @@ rt_raster_serialized_size(rt_raster raster) {
         /* Align size to 8-bytes boundary (trailing padding) */
         /* XXX jorgearevalo: bug here. If the size is actually 8-bytes aligned,
          this line will add 8 bytes trailing padding, and it's not necessary */
-        //size += 8 - (size % 8);
+        /*size += 8 - (size % 8);*/
         if (size % 8)
             size += 8 - (size % 8);
 
@@ -7324,6 +7389,7 @@ rt_raster rt_raster_gdal_warp(
  * @param grid_yw : the Y value of point on grid to align raster to
  * @param skew_x : the X skew of the raster
  * @param skew_y : the Y skew of the raster
+ * @param options : array of options.  only option is "ALL_TOUCHED"
  *
  * @return the raster of the provided geometry
  */
@@ -7337,7 +7403,8 @@ rt_raster_gdal_rasterize(const unsigned char *wkb,
 	double *scale_x, double *scale_y,
 	double *ul_xw, double *ul_yw,
 	double *grid_xw, double *grid_yw,
-	double *skew_x, double *skew_y
+	double *skew_x, double *skew_y,
+	char **options
 ) {
 	rt_raster rast;
 	int i = 0;
@@ -7362,6 +7429,7 @@ rt_raster_gdal_rasterize(const unsigned char *wkb,
 	OGRSpatialReferenceH src_sr = NULL;
 	OGRGeometryH src_geom;
 	OGREnvelope src_env;
+	OGRwkbGeometryType wkbtype = wkbUnknown;
 
 	int ul_user = 0;
 	double djunk = 0;
@@ -7499,10 +7567,17 @@ rt_raster_gdal_rasterize(const unsigned char *wkb,
 	RASTER_DEBUGF(3, "scale (x, y) = %f, %f", _scale_x, _scale_y);
 	RASTER_DEBUGF(3, "dim (x, y) = %d, %d", _width, _height);
 
-	if (
-		(wkbFlatten(OGR_G_GetGeometryType(src_geom)) == wkbPoint) &&
-		FLT_NEQ(_scale_x, 0) &&
-		FLT_NEQ(_scale_y, 0)
+	/*
+	 	if geometry is a point or point-like and bounds not set, increase
+		extent by half-pixel to avoid missing points on border
+	*/
+	wkbtype = wkbFlatten(OGR_G_GetGeometryType(src_geom));
+	if ((
+			(wkbtype == wkbPoint) ||
+			(wkbtype == wkbMultiPoint)
+		) &&
+		FLT_EQ(_width, 0) &&
+		FLT_EQ(_height, 0)
 	) {
 		src_env.MinX -= (_scale_x / 2.);
 		src_env.MaxX += (_scale_x / 2.);
@@ -7795,7 +7870,7 @@ rt_raster_gdal_rasterize(const unsigned char *wkb,
 		1, &src_geom,
 		NULL, NULL,
 		_value,
-		NULL,
+		options,
 		NULL, NULL
 	);
 	rtdealloc(band_list);
@@ -7851,3 +7926,657 @@ rt_raster_gdal_rasterize(const unsigned char *wkb,
 
 	return rast;
 }
+
+static
+int rt_raster_intersects_algorithm(
+	rt_raster rast1, rt_raster rast2,
+	rt_band band1, rt_band band2,
+	int hasnodata1, int hasnodata2,
+	double nodata1, double nodata2
+) {
+	int i;
+	int byHeight = 1;
+	uint32_t dimValue;
+
+	uint32_t row;
+	uint32_t rowoffset;
+	uint32_t col;
+	uint32_t coloffset;
+
+	enum line_points {X1, Y1, X2, Y2};
+	enum point {pX, pY};
+	double line1[4];
+	double line2[4];
+	double P[2];
+	double Qw[2];
+	double Qr[2];
+	double d;
+	double val1;
+	int noval1;
+	double val2;
+	int noval2;
+	uint32_t adjacent[8] = {0};
+
+	double xscale;
+	double yscale;
+
+	uint16_t width1;
+	uint16_t height1;
+	uint16_t width2;
+	uint16_t height2;
+
+	width1 = rt_raster_get_width(rast1);
+	height1 = rt_raster_get_height(rast1);
+	width2 = rt_raster_get_width(rast2);
+	height2 = rt_raster_get_height(rast2);
+
+	/* sampling scale */
+	xscale = fmin(rt_raster_get_x_scale(rast1), rt_raster_get_x_scale(rast2)) / 10.;
+	yscale = fmin(rt_raster_get_y_scale(rast1), rt_raster_get_y_scale(rast2)) / 10.;
+
+	/* see if skew made rast2's rows are parallel to rast1's cols */
+	rt_raster_cell_to_geopoint(
+		rast1,
+		0, 0,
+		&(line1[X1]), &(line1[Y1])
+	);
+
+	rt_raster_cell_to_geopoint(
+		rast1,
+		0, height1,
+		&(line1[X2]), &(line1[Y2])
+	);
+
+	rt_raster_cell_to_geopoint(
+		rast2,
+		0, 0,
+		&(line2[X1]), &(line2[Y1])
+	);
+
+	rt_raster_cell_to_geopoint(
+		rast2,
+		width2, 0,
+		&(line2[X2]), &(line2[Y2])
+	);
+
+	/* parallel vertically */
+	if (FLT_EQ(line1[X2] - line1[X1], 0.) && FLT_EQ(line2[X2] - line2[X1], 0.))
+		byHeight = 0;
+	/* parallel */
+	else if (FLT_EQ(((line1[Y2] - line1[Y1]) / (line1[X2] - line1[X1])), ((line2[Y2] - line2[Y1]) / (line2[X2] - line2[X1]))))
+		byHeight = 0;
+
+	if (byHeight)
+		dimValue = height2;
+	else
+		dimValue = width2;
+	RASTER_DEBUGF(4, "byHeight: %d, dimValue: %d", byHeight, dimValue);
+
+	/* 3 x 3 search */
+	for (coloffset = 0; coloffset < 3; coloffset++) {
+		for (rowoffset = 0; rowoffset < 3; rowoffset++) {
+			/* smaller raster */
+			for (col = coloffset; col <= width1; col += 3) {
+
+				rt_raster_cell_to_geopoint(
+					rast1,
+					col, 0,
+					&(line1[X1]), &(line1[Y1])
+				);
+
+				rt_raster_cell_to_geopoint(
+					rast1,
+					col, height1,
+					&(line1[X2]), &(line1[Y2])
+				);
+
+				/* larger raster */
+				for (row = rowoffset; row <= dimValue; row += 3) {
+
+					if (byHeight) {
+						rt_raster_cell_to_geopoint(
+							rast2,
+							0, row,
+							&(line2[X1]), &(line2[Y1])
+						);
+
+						rt_raster_cell_to_geopoint(
+							rast2,
+							width2, row,
+							&(line2[X2]), &(line2[Y2])
+						);
+					}
+					else {
+						rt_raster_cell_to_geopoint(
+							rast2,
+							row, 0,
+							&(line2[X1]), &(line2[Y1])
+						);
+
+						rt_raster_cell_to_geopoint(
+							rast2,
+							row, height2,
+							&(line2[X2]), &(line2[Y2])
+						);
+					}
+
+					RASTER_DEBUGF(4, "(col, row) = (%d, %d)", col, row);
+					RASTER_DEBUGF(4, "line1(x1, y1, x2, y2) = (%f, %f, %f, %f)",
+						line1[X1], line1[Y1], line1[X2], line1[Y2]);
+					RASTER_DEBUGF(4, "line2(x1, y1, x2, y2) = (%f, %f, %f, %f)",
+						line2[X1], line2[Y1], line2[X2], line2[Y2]);
+
+					/* intersection */
+					/* http://en.wikipedia.org/wiki/Line-line_intersection */
+					d = ((line1[X1] - line1[X2]) * (line2[Y1] - line2[Y2])) - ((line1[Y1] - line1[Y2]) * (line2[X1] - line2[X2]));
+					/* no intersection */
+					if (FLT_EQ(d, 0.)) {
+						continue;
+					}
+
+					P[pX] = (((line1[X1] * line1[Y2]) - (line1[Y1] * line1[X2])) * (line2[X1] - line2[X2])) -
+						((line1[X1] - line1[X2]) * ((line2[X1] * line2[Y2]) - (line2[Y1] * line2[X2])));
+					P[pX] = P[pX] / d;
+
+					P[pY] = (((line1[X1] * line1[Y2]) - (line1[Y1] * line1[X2])) * (line2[Y1] - line2[Y2])) -
+						((line1[Y1] - line1[Y2]) * ((line2[X1] * line2[Y2]) - (line2[Y1] * line2[X2])));
+					P[pY] = P[pY] / d;
+
+					RASTER_DEBUGF(4, "P(x, y) = (%f, %f)", P[pX], P[pY]);
+
+					/* intersection within bounds */
+					if ((
+							(FLT_EQ(P[pX], line1[X1]) || FLT_EQ(P[pX], line1[X2])) ||
+								(P[pX] > fmin(line1[X1], line1[X2]) && (P[pX] < fmax(line1[X1], line1[X2])))
+						) && (
+							(FLT_EQ(P[pY], line1[Y1]) || FLT_EQ(P[pY], line1[Y2])) ||
+								(P[pY] > fmin(line1[Y1], line1[Y2]) && (P[pY] < fmax(line1[Y1], line1[Y2])))
+						) && (
+							(FLT_EQ(P[pX], line2[X1]) || FLT_EQ(P[pX], line2[X2])) ||
+								(P[pX] > fmin(line2[X1], line2[X2]) && (P[pX] < fmax(line2[X1], line2[X2])))
+						) && (
+							(FLT_EQ(P[pY], line2[Y1]) || FLT_EQ(P[pY], line2[Y2])) ||
+								(P[pY] > fmin(line2[Y1], line2[Y2]) && (P[pY] < fmax(line2[Y1], line2[Y2])))
+					)) {
+						RASTER_DEBUG(4, "within bounds");
+
+						for (i = 0; i < 8; i++) adjacent[i] = 0;
+
+						/* test points around intersection */
+						for (i = 0; i < 8; i++) {
+							switch (i) {
+								case 7:
+									Qw[pX] = P[pX] - xscale;
+									Qw[pY] = P[pY] + yscale;
+									break;
+								/* 270 degrees = 09:00 */
+								case 6:
+									Qw[pX] = P[pX] - xscale;
+									Qw[pY] = P[pY];
+									break;
+								case 5:
+									Qw[pX] = P[pX] - xscale;
+									Qw[pY] = P[pY] - yscale;
+									break;
+								/* 180 degrees = 06:00 */
+								case 4:
+									Qw[pX] = P[pX];
+									Qw[pY] = P[pY] - yscale;
+									break;
+								case 3:
+									Qw[pX] = P[pX] + xscale;
+									Qw[pY] = P[pY] - yscale;
+									break;
+								/* 90 degrees = 03:00 */
+								case 2:
+									Qw[pX] = P[pX] + xscale;
+									Qw[pY] = P[pY];
+									break;
+								/* 45 degrees */
+								case 1:
+									Qw[pX] = P[pX] + xscale;
+									Qw[pY] = P[pY] + yscale;
+									break;
+								/* 0 degrees = 00:00 */
+								case 0:
+									Qw[pX] = P[pX];
+									Qw[pY] = P[pY] + yscale;
+									break;
+							}
+
+							/* unable to convert point to cell */
+							noval1 = 0;
+							if (!rt_raster_geopoint_to_cell(
+								rast1,
+								Qw[pX], Qw[pY],
+								&(Qr[pX]), &(Qr[pY])
+							)) {
+								noval1 = 1;
+							}
+							/* cell is outside bounds of grid */
+							else if (
+								(Qr[pX] < 0 || Qr[pX] > width1 || FLT_EQ(Qr[pX], width1)) ||
+								(Qr[pY] < 0 || Qr[pY] > height1 || FLT_EQ(Qr[pY], height1))
+							) {
+								noval1 = 1;
+							}
+							else if (hasnodata1 == FALSE)
+								val1 = 1;
+							/* unable to get value at cell */
+							else if (rt_band_get_pixel(band1, Qr[pX], Qr[pY], &val1) < 0)
+								noval1 = 1;
+
+							/* unable to convert point to cell */
+							noval2 = 0;
+							if (!rt_raster_geopoint_to_cell(
+								rast2,
+								Qw[pX], Qw[pY],
+								&(Qr[pX]), &(Qr[pY])
+							)) {
+								noval2 = 1;
+							}
+							/* cell is outside bounds of grid */
+							else if (
+								(Qr[pX] < 0 || Qr[pX] > width2 || FLT_EQ(Qr[pX], width2)) ||
+								(Qr[pY] < 0 || Qr[pY] > height2 || FLT_EQ(Qr[pY], height2))
+							) {
+								noval2 = 1;
+							}
+							else if (hasnodata2 == FALSE)
+								val2 = 1;
+							/* unable to get value at cell */
+							else if (rt_band_get_pixel(band2, Qr[pX], Qr[pY], &val2) < 0)
+								noval2 = 1;
+
+							if (!noval1) {
+								RASTER_DEBUGF(4, "val1 = %f", val1);
+							}
+							if (!noval2) {
+								RASTER_DEBUGF(4, "val2 = %f", val2);
+							}
+
+							/* pixels touch */
+							if (!noval1 && (
+								(hasnodata1 == FALSE) || (
+									(hasnodata1 != FALSE) &&
+									FLT_NEQ(val1, nodata1)
+								)
+							)) {
+								adjacent[i]++;
+							}
+							if (!noval2 && (
+								(hasnodata2 == FALSE) || (
+									(hasnodata2 != FALSE) &&
+									FLT_NEQ(val2, nodata2)
+								)
+							)) {
+								adjacent[i] += 3;
+							}
+
+							/* two pixel values not present */
+							if (noval1 || noval2) {
+								RASTER_DEBUGF(4, "noval1 = %d, noval2 = %d", noval1, noval2);
+								continue;
+							}
+
+							/* pixels valid, so intersect */
+							if ((
+									(hasnodata1 == FALSE) || (
+										(hasnodata1 != FALSE) &&
+										FLT_NEQ(val1, nodata1)
+									)
+								) && (
+									(hasnodata2 == FALSE) || (
+										(hasnodata2 != FALSE) &&
+										FLT_NEQ(val2, nodata2)
+									)
+							)) {
+								RASTER_DEBUG(3, "The two rasters do intersect");
+
+								return 1;
+							}
+						}
+
+						/* pixels touch */
+						for (i = 0; i < 4; i++) {
+							RASTER_DEBUGF(4, "adjacent[%d] = %d, adjacent[%d] = %d"
+								, i, adjacent[i], i + 4, adjacent[i + 4]);
+							if (adjacent[i] == 0) continue;
+
+							if (adjacent[i] + adjacent[i + 4] == 4) {
+								RASTER_DEBUG(3, "The two rasters touch");
+
+								return 1;
+							}
+						}
+					}
+					else {
+						RASTER_DEBUG(4, "outside of bounds");
+					}
+				}
+			}
+		}
+	}
+
+	return 0;
+}
+
+/**
+ * Return zero if error occurred in function.
+ * Parameter intersects returns non-zero if two rasters intersect
+ *
+ * @param rast1 : the first raster whose band will be tested
+ * @param nband1 : the 0-based band of raster rast1 to use
+ *   if value is less than zero, bands are ignored.
+ *   if nband1 gte zero, nband2 must be gte zero
+ * @param rast2 : the second raster whose band will be tested
+ * @param nband2 : the 0-based band of raster rast2 to use
+ *   if value is less than zero, bands are ignored
+ *   if nband2 gte zero, nband1 must be gte zero
+ * @param intersects : non-zero value if the two rasters' bands intersects
+ *
+ * @return if zero, an error occurred in function
+ */
+int
+rt_raster_intersects(
+	rt_raster rast1, int nband1,
+	rt_raster rast2, int nband2,
+	int *intersects
+) {
+	int i;
+	int j;
+	int within = 0;
+
+	LWPOLY *hull[2] = {NULL};
+	GEOSGeometry *ghull[2] = {NULL};
+
+	uint16_t width1;
+	uint16_t height1;
+	uint16_t width2;
+	uint16_t height2;
+	double area1;
+	double area2;
+	double pixarea1;
+	double pixarea2;
+	rt_raster rastS = NULL;
+	rt_raster rastL = NULL;
+	uint16_t *widthS = NULL;
+	uint16_t *heightS = NULL;
+	uint16_t *widthL = NULL;
+	uint16_t *heightL = NULL;
+	int nbandS;
+	int nbandL;
+	rt_band bandS = NULL;
+	rt_band bandL = NULL;
+	int hasnodataS = FALSE;
+	int hasnodataL = FALSE;
+	double nodataS = 0;
+	double nodataL = 0;
+
+	uint32_t row;
+	uint32_t rowoffset;
+	uint32_t col;
+	uint32_t coloffset;
+
+	enum line_points {X1, Y1, X2, Y2};
+	enum point {pX, pY};
+	double lineS[4];
+	double Qr[2];
+	double valS;
+	double valL;
+
+	RASTER_DEBUG(3, "Starting");
+
+	assert(NULL != rast1);
+	assert(NULL != rast2);
+
+	if (nband1 < 0 && nband2 < 0) {
+		nband1 = -1;
+		nband2 = -1;
+	}
+	else {
+		assert(nband1 >= 0 && nband1 < rt_raster_get_num_bands(rast1));
+		assert(nband2 >= 0 && nband2 < rt_raster_get_num_bands(rast2));
+	}
+
+	/* same srid */
+	if (rt_raster_get_srid(rast1) != rt_raster_get_srid(rast2)) {
+		rterror("rt_raster_intersects: The two rasters provided have different SRIDs");
+		*intersects = 0;
+		return 0;
+	}
+
+	/* raster extents need to intersect */
+	do {
+		int rtn;
+
+		initGEOS(lwnotice, lwgeom_geos_error);
+
+		rtn = 1;
+		for (i = 0; i < 2; i++) {
+			hull[i] = rt_raster_get_convex_hull(i < 1 ? rast1 : rast2);
+			if (NULL == hull[i]) {
+				for (j = 0; j < i; j++) {
+					GEOSGeom_destroy(ghull[j]);
+					lwpoly_free(hull[j]);
+				}
+				rtn = 0;
+				break;
+			}
+			ghull[i] = (GEOSGeometry *) LWGEOM2GEOS((LWGEOM *) hull[i]);
+			if (NULL == ghull[i]) {
+				for (j = 0; j < i; j++) {
+					GEOSGeom_destroy(ghull[j]);
+					lwpoly_free(hull[j]);
+				}
+				lwpoly_free(hull[i]);
+				rtn = 0;
+				break;
+			}
+		}
+		if (!rtn) break;
+
+		/* test to see if raster within the other */
+		within = 0;
+		if (GEOSWithin(ghull[0], ghull[1]) == 1)
+			within = -1;
+		else if (GEOSWithin(ghull[1], ghull[0]) == 1)
+			within = 1;
+
+		if (within != 0)
+			rtn = 1;
+		else
+			rtn = GEOSIntersects(ghull[0], ghull[1]);
+
+		for (i = 0; i < 2; i++) {
+			GEOSGeom_destroy(ghull[i]);
+			lwpoly_free(hull[i]);
+		}
+
+		if (rtn != 2) {
+			RASTER_DEBUGF(4, "convex hulls of rasters do %sintersect", rtn != 1 ? "NOT " : "");
+			if (rtn != 1) {
+				*intersects = 0;
+				return 1;
+			}
+			/* band isn't specified */
+			else if (nband1 < 0) {
+				*intersects = 1;
+				return 1;
+			}
+		}
+	}
+	while (0);
+
+	/* smaller raster by area or width */
+	width1 = rt_raster_get_width(rast1);
+	height1 = rt_raster_get_height(rast1);
+	width2 = rt_raster_get_width(rast2);
+	height2 = rt_raster_get_height(rast2);
+	pixarea1 = fabs(rt_raster_get_x_scale(rast1) * rt_raster_get_y_scale(rast1));
+	pixarea2 = fabs(rt_raster_get_x_scale(rast2) * rt_raster_get_y_scale(rast2));
+	area1 = fabs(width1 * height1 * pixarea1);
+	area2 = fabs(width2 * height2 * pixarea2);
+	RASTER_DEBUGF(4, "pixarea1, pixarea2, area1, area2 = %f, %f, %f, %f",
+		pixarea1, pixarea2, area1, area2);
+	if (
+		(within <= 0) ||
+		(area1 < area2) ||
+		FLT_EQ(area1, area2) ||
+		(area1 < pixarea2) || /* area of rast1 smaller than pixel area of rast2 */
+		FLT_EQ(area1, pixarea2)
+	) {
+		rastS = rast1;
+		nbandS = nband1;
+		widthS = &width1;
+		heightS = &height1;
+
+		rastL = rast2;
+		nbandL = nband2;
+		widthL = &width2;
+		heightL = &height2;
+	}
+	else {
+		rastS = rast2;
+		nbandS = nband2;
+		widthS = &width2;
+		heightS = &height2;
+
+		rastL = rast1;
+		nbandL = nband1;
+		widthL = &width1;
+		heightL = &height1;
+	}
+
+	/* no band to use, set band to zero */
+	if (nband1 < 0) {
+		nbandS = 0;
+		nbandL = 0;
+	}
+
+	RASTER_DEBUGF(4, "rast1 @ %p", rast1);
+	RASTER_DEBUGF(4, "rast2 @ %p", rast2);
+	RASTER_DEBUGF(4, "rastS @ %p", rastS);
+	RASTER_DEBUGF(4, "rastL @ %p", rastL);
+
+	/* load band of smaller raster */
+	bandS = rt_raster_get_band(rastS, nbandS);
+	if (NULL == bandS) {
+		rterror("rt_raster_intersects: Unable to get band %d of the first raster\n", nbandS);
+		*intersects = 0;
+		return 0;
+	}
+	if (bandS->offline) {
+		rterror("rt_raster_intersects not implemented yet for OFFDB bands");
+		*intersects = 0;
+		return 0;
+	}
+
+	hasnodataS = rt_band_get_hasnodata_flag(bandS);
+	if (hasnodataS != FALSE)
+		nodataS = rt_band_get_nodata(bandS);
+
+	/* load band of larger raster */
+	bandL = rt_raster_get_band(rastL, nbandL);
+	if (NULL == bandL) {
+		rterror("rt_raster_intersects: Unable to get band %d of the first raster\n", nbandL);
+		*intersects = 0;
+		return 0;
+	}
+	if (bandL->offline) {
+		rterror("rt_raster_intersects not implemented yet for OFFDB bands");
+		*intersects = 0;
+		return 0;
+	}
+
+	hasnodataL = rt_band_get_hasnodata_flag(bandL);
+	if (hasnodataL != FALSE)
+		nodataL = rt_band_get_nodata(bandL);
+
+	/* no band to use, ignore nodata */
+	if (nband1 < 0) {
+		hasnodataS = FALSE;
+		hasnodataL = FALSE;
+	}
+
+	/* special case where a raster can fit inside another raster's pixel */
+	if (within != 0 && ((pixarea1 > area2) || (pixarea2 > area1))) {
+		RASTER_DEBUG(4, "Using special case of raster fitting into another raster's pixel");
+		/* 3 x 3 search */
+		for (coloffset = 0; coloffset < 3; coloffset++) {
+			for (rowoffset = 0; rowoffset < 3; rowoffset++) {
+				for (col = coloffset; col < *widthS; col += 3) {
+					for (row = rowoffset; row < *heightS; row += 3) {
+						if (hasnodataS == FALSE)
+							valS = 1;
+						else if (rt_band_get_pixel(bandS, col, row, &valS) < 0)
+							continue;
+
+						if ((hasnodataS == FALSE) || (
+							(hasnodataS != FALSE) &&
+							FLT_NEQ(valS, nodataS)
+						)) {
+							rt_raster_cell_to_geopoint(
+								rastS,
+								col, row,
+								&(lineS[X1]), &(lineS[Y1])
+							);
+
+							if (!rt_raster_geopoint_to_cell(
+								rastL,
+								lineS[X1], lineS[Y1],
+								&(Qr[pX]), &(Qr[pY])
+							)) {
+								continue;
+							}
+
+							if (
+								(Qr[pX] < 0 || Qr[pX] > *widthL || FLT_EQ(Qr[pX], *widthL)) ||
+								(Qr[pY] < 0 || Qr[pY] > *heightL || FLT_EQ(Qr[pY], *heightL))
+							) {
+								continue;
+							}
+
+							if (hasnodataS == FALSE)
+								valL = 1;
+							else if (rt_band_get_pixel(bandL, Qr[pX], Qr[pY], &valL) < 0)
+								continue;
+
+							if ((hasnodataL == FALSE) || (
+								(hasnodataL != FALSE) &&
+								FLT_NEQ(valL, nodataL)
+							)) {
+								RASTER_DEBUG(3, "The two rasters do intersect");
+								*intersects = 1;
+								return 1;
+							}
+						}
+					}
+				}
+			}
+		}
+	}
+
+	*intersects = rt_raster_intersects_algorithm(
+		rastS, rastL,
+		bandS, bandL,
+		hasnodataS, hasnodataL,
+		nodataS, nodataL
+	);
+
+	if (*intersects) return 1;
+
+	*intersects = rt_raster_intersects_algorithm(
+		rastL, rastS,
+		bandL, bandS,
+		hasnodataL, hasnodataS,
+		nodataL, nodataS
+	);
+
+	if (*intersects) return 1;
+
+	RASTER_DEBUG(3, "The two rasters do not intersect");
+
+	*intersects = 0;
+	return 1;
+}
diff --git a/raster/rt_core/rt_api.h b/raster/rt_core/rt_api.h
index 7d60c43ed..4c8e67139 100644
--- a/raster/rt_core/rt_api.h
+++ b/raster/rt_core/rt_api.h
@@ -90,10 +90,21 @@
 #endif
 #endif
 
+/* Microsoft does not support fmax or fmin in math.h */
+#if defined(WIN32)
+#if !defined(HAVE_MAX)
+#define fmax max
+#endif
+#if !defined(HAVE_MIN)
+#define fmin min
+#endif
+#endif
+
 #include <stdlib.h> /* For size_t, srand and rand */
 #include <stdint.h> /* For C99 int types */
 
 #include "liblwgeom.h"
+#include "lwgeom_geos.h"
 
 #include "gdal_alg.h"
 #include "gdal_frmts.h"
@@ -137,170 +148,6 @@ typedef struct rt_valuecount_t* rt_valuecount;
 typedef struct rt_gdaldriver_t* rt_gdaldriver;
 typedef struct rt_reclassexpr_t* rt_reclassexpr;
 
-
-/**
- * Struct definitions
- *
- * These structs are defined here as they are needed elsewhere
- * including rt_pg/rt_pg.c and reduce duplicative declarations
- *
- */
-struct rt_raster_serialized_t {
-    /*---[ 8 byte boundary ]---{ */
-    uint32_t size; /* required by postgresql: 4 bytes */
-    uint16_t version; /* format version (this is version 0): 2 bytes */
-    uint16_t numBands; /* Number of bands: 2 bytes */
-
-    /* }---[ 8 byte boundary ]---{ */
-    double scaleX; /* pixel width: 8 bytes */
-
-    /* }---[ 8 byte boundary ]---{ */
-    double scaleY; /* pixel height: 8 bytes */
-
-    /* }---[ 8 byte boundary ]---{ */
-    double ipX; /* insertion point X: 8 bytes */
-
-    /* }---[ 8 byte boundary ]---{ */
-    double ipY; /* insertion point Y: 8 bytes */
-
-    /* }---[ 8 byte boundary ]---{ */
-    double skewX; /* skew about the X axis: 8 bytes */
-
-    /* }---[ 8 byte boundary ]---{ */
-    double skewY; /* skew about the Y axis: 8 bytes */
-
-    /* }---[ 8 byte boundary ]--- */
-    int32_t srid; /* Spatial reference id: 4 bytes */
-    uint16_t width; /* pixel columns: 2 bytes */
-    uint16_t height; /* pixel rows: 2 bytes */
-};
-
-/* NOTE: the initial part of this structure matches the layout
- *       of data in the serialized form version 0, starting
- *       from the numBands element
- */
-struct rt_raster_t {
-    uint32_t size;
-    uint16_t version;
-
-    /* Number of bands, all share the same dimension
-     * and georeference */
-    uint16_t numBands;
-
-    /* Georeference (in projection units) */
-    double scaleX; /* pixel width */
-    double scaleY; /* pixel height */
-    double ipX; /* geo x ordinate of the corner of upper-left pixel */
-    double ipY; /* geo y ordinate of the corner of bottom-right pixel */
-    double skewX; /* skew about the X axis*/
-    double skewY; /* skew about the Y axis */
-
-    int32_t srid; /* spatial reference id */
-    uint16_t width; /* pixel columns - max 65535 */
-    uint16_t height; /* pixel rows - max 65535 */
-    rt_band *bands; /* actual bands */
-
-};
-
-/* WKT string representing each polygon in WKT format acompagned by its
-correspoding value */
-struct rt_geomval_t {
-    int srid;
-    double val;
-    char * geom;
-};
-
-/* summary stats of specified band */
-struct rt_bandstats_t {
-	double sample;
-	uint32_t count;
-
-	double min;
-	double max;
-	double sum;
-	double mean;
-	double stddev;
-
-	double *values;
-	int sorted; /* flag indicating that values is sorted ascending by value */
-};
-
-/* histogram bin(s) of specified band */
-struct rt_histogram_t {
-	uint32_t count;
-	double percent;
-
-	double min;
-	double max;
-
-	int inc_min;
-	int inc_max;
-};
-
-/* quantile(s) of the specified band */
-struct rt_quantile_t {
-	double quantile;
-	double value;
-};
-
-/* listed-list structures for rt_band_get_quantiles_stream */
-struct quantile_llist {
-	uint8_t algeq; /* AL-GEQ (1) or AL-GT (0) */
-	double quantile;
-	uint64_t tau; /* position in sequence */
-
-	struct quantile_llist_element *head; /* H index 0 */
-	struct quantile_llist_element *tail; /* H index last */
-	uint32_t count; /* # of elements in H */
-
-	/* faster access to elements at specific intervals */
-	struct quantile_llist_index *index;
-	uint32_t index_max; /* max # of elements in index */
-
-	uint64_t sum1; /* N1H */
-	uint64_t sum2; /* N2H */
-};
-
-struct quantile_llist_element {
-	double value;
-	uint32_t count;
-
-	struct quantile_llist_element *prev;
-	struct quantile_llist_element *next;
-};
-
-struct quantile_llist_index {
-	struct quantile_llist_element *element;
-	uint32_t index;
-};
-
-/* number of times a value occurs */
-struct rt_valuecount_t {
-	double value;
-	uint32_t count;
-	double percent;
-};
-
-/* reclassification expression */
-struct rt_reclassexpr_t {
-	struct rt_reclassrange {
-		double min;
-		double max;
-		int inc_min; /* include min */
-		int inc_max; /* include max */
-		int exc_min; /* exceed min */
-		int exc_max; /* exceed max */
-	} src, dst;
-};
-
-/* gdal driver information */
-struct rt_gdaldriver_t {
-    int idx;
-    char *short_name;
-    char *long_name;
-		char *create_options;
-};
-
 /**
 * Global functions for memory/logging handlers.
 */
@@ -676,6 +523,7 @@ rt_histogram rt_band_get_histogram(rt_bandstats stats,
 rt_quantile rt_band_get_quantiles(rt_bandstats stats,
 	double *quantiles, int quantiles_count, uint32_t *rtn_count);
 
+struct quantile_llist;
 int quantile_llist_destroy(struct quantile_llist **list,
 	uint32_t list_count);
 
@@ -966,17 +814,32 @@ void rt_raster_set_srid(rt_raster raster, int32_t srid);
 int32_t rt_raster_get_srid(rt_raster raster);
 
 /**
- * Convert an x,y raster point to an x1,y1 point on map
+ * Convert an xr, yr raster point to an xw, yw point on map
  *
  * @param raster : the raster to get info from
- * @param x : the pixel's column
- * @param y : the pixel's row
- * @param x1 : output parameter, X ordinate of the geographical point
- * @param y1 : output parameter, Y ordinate of the geographical point
+ * @param xr : the pixel's column
+ * @param yr : the pixel's row
+ * @param xw : output parameter, X ordinate of the geographical point
+ * @param yw : output parameter, Y ordinate of the geographical point
  */
 void rt_raster_cell_to_geopoint(rt_raster raster,
-                                double x, double y,
-                                double* x1, double* y1);
+	double xr, double yr,
+	double* xw, double* yw);
+
+/**
+ * Convert an xw, yw map point to a xr, yr raster point
+ *
+ * @param raster : the raster to get info from
+ * @param xw : X ordinate of the geographical point
+ * @param yw : Y ordinate of the geographical point
+ * @param xr : output parameter, the pixel's column
+ * @param yr : output parameter, the pixel's row
+ *
+ * @return if zero, error occurred in function
+ */
+int rt_raster_geopoint_to_cell(rt_raster raster,
+	double xw, double yw,
+	double *xr, double *yr);
 
 /**
  * Get raster's polygon convex hull.
@@ -1184,6 +1047,7 @@ rt_raster rt_raster_gdal_warp(rt_raster raster, const char *src_srs,
  * @param grid_yw : the Y value of point on grid to align raster to
  * @param skew_x : the X skew of the raster
  * @param skew_y : the Y skew of the raster
+ * @param options : array of options.  only option is "ALL_TOUCHED"
  *
  * @return the raster of the provided geometry
  */
@@ -1196,7 +1060,31 @@ rt_raster rt_raster_gdal_rasterize(const unsigned char *wkb,
 	double *scale_x, double *scale_y,
 	double *ul_xw, double *ul_yw,
 	double *grid_xw, double *grid_yw,
-	double *skew_x, double *skew_y);
+	double *skew_x, double *skew_y,
+	char **options
+);
+
+/**
+ * Return zero if error occurred in function.
+ * Parameter intersects returns non-zero if two rasters intersect
+ *
+ * @param rast1 : the first raster whose band will be tested
+ * @param nband1 : the 0-based band of raster rast1 to use
+ *   if value is less than zero, bands are ignored.
+ *   if nband1 gte zero, nband2 must be gte zero
+ * @param rast2 : the second raster whose band will be tested
+ * @param nband2 : the 0-based band of raster rast2 to use
+ *   if value is less than zero, bands are ignored
+ *   if nband2 gte zero, nband1 must be gte zero
+ * @param intersects : non-zero value if the two rasters' bands intersects
+ *
+ * @return if zero, an error occurred in function
+ */
+int rt_raster_intersects(
+	rt_raster rast1, int nband1,
+	rt_raster rast2, int nband2,
+	int *intersects
+);
 
 /*- utilities -------------------------------------------------------*/
 
@@ -1279,4 +1167,190 @@ rt_util_pixtype_to_gdal_datatype(rt_pixtype pt);
 */
 #define ROUND(x, y) (((x > 0.0) ? floor((x * pow(10, y) + 0.5)) : ceil((x * pow(10, y) - 0.5))) / pow(10, y));
 
+/**
+ * Struct definitions
+ *
+ * These structs are defined here as they are needed elsewhere
+ * including rt_pg/rt_pg.c and reduce duplicative declarations
+ *
+ */
+struct rt_raster_serialized_t {
+    /*---[ 8 byte boundary ]---{ */
+    uint32_t size; /* required by postgresql: 4 bytes */
+    uint16_t version; /* format version (this is version 0): 2 bytes */
+    uint16_t numBands; /* Number of bands: 2 bytes */
+
+    /* }---[ 8 byte boundary ]---{ */
+    double scaleX; /* pixel width: 8 bytes */
+
+    /* }---[ 8 byte boundary ]---{ */
+    double scaleY; /* pixel height: 8 bytes */
+
+    /* }---[ 8 byte boundary ]---{ */
+    double ipX; /* insertion point X: 8 bytes */
+
+    /* }---[ 8 byte boundary ]---{ */
+    double ipY; /* insertion point Y: 8 bytes */
+
+    /* }---[ 8 byte boundary ]---{ */
+    double skewX; /* skew about the X axis: 8 bytes */
+
+    /* }---[ 8 byte boundary ]---{ */
+    double skewY; /* skew about the Y axis: 8 bytes */
+
+    /* }---[ 8 byte boundary ]--- */
+    int32_t srid; /* Spatial reference id: 4 bytes */
+    uint16_t width; /* pixel columns: 2 bytes */
+    uint16_t height; /* pixel rows: 2 bytes */
+};
+
+/* NOTE: the initial part of this structure matches the layout
+ *       of data in the serialized form version 0, starting
+ *       from the numBands element
+ */
+struct rt_raster_t {
+    uint32_t size;
+    uint16_t version;
+
+    /* Number of bands, all share the same dimension
+     * and georeference */
+    uint16_t numBands;
+
+    /* Georeference (in projection units) */
+    double scaleX; /* pixel width */
+    double scaleY; /* pixel height */
+    double ipX; /* geo x ordinate of the corner of upper-left pixel */
+    double ipY; /* geo y ordinate of the corner of bottom-right pixel */
+    double skewX; /* skew about the X axis*/
+    double skewY; /* skew about the Y axis */
+
+    int32_t srid; /* spatial reference id */
+    uint16_t width; /* pixel columns - max 65535 */
+    uint16_t height; /* pixel rows - max 65535 */
+    rt_band *bands; /* actual bands */
+
+};
+
+struct rt_extband_t {
+    uint8_t bandNum;
+    char* path; /* externally owned ? */
+};
+
+struct rt_band_t {
+    rt_pixtype pixtype;
+    int32_t offline;
+    uint16_t width;
+    uint16_t height;
+    int32_t hasnodata; /* a flag indicating if this band contains nodata values */
+    int32_t isnodata;   /* a flag indicating if this band is filled only with
+                           nodata values */
+    double nodataval; /* int will be converted ... */
+    int32_t ownsData; /* XXX mloskot: its behaviour needs to be documented */
+
+    union {
+        void* mem; /* actual data, externally owned */
+        struct rt_extband_t offline;
+    } data;
+
+};
+
+/* WKT string representing each polygon in WKT format accompanied by its
+corresponding value */
+struct rt_geomval_t {
+    int srid;
+    double val;
+    char * geom;
+};
+
+/* summary stats of specified band */
+struct rt_bandstats_t {
+	double sample;
+	uint32_t count;
+
+	double min;
+	double max;
+	double sum;
+	double mean;
+	double stddev;
+
+	double *values;
+	int sorted; /* flag indicating that values is sorted ascending by value */
+};
+
+/* histogram bin(s) of specified band */
+struct rt_histogram_t {
+	uint32_t count;
+	double percent;
+
+	double min;
+	double max;
+
+	int inc_min;
+	int inc_max;
+};
+
+/* quantile(s) of the specified band */
+struct rt_quantile_t {
+	double quantile;
+	double value;
+};
+
+/* listed-list structures for rt_band_get_quantiles_stream */
+struct quantile_llist {
+	uint8_t algeq; /* AL-GEQ (1) or AL-GT (0) */
+	double quantile;
+	uint64_t tau; /* position in sequence */
+
+	struct quantile_llist_element *head; /* H index 0 */
+	struct quantile_llist_element *tail; /* H index last */
+	uint32_t count; /* # of elements in H */
+
+	/* faster access to elements at specific intervals */
+	struct quantile_llist_index *index;
+	uint32_t index_max; /* max # of elements in index */
+
+	uint64_t sum1; /* N1H */
+	uint64_t sum2; /* N2H */
+};
+
+struct quantile_llist_element {
+	double value;
+	uint32_t count;
+
+	struct quantile_llist_element *prev;
+	struct quantile_llist_element *next;
+};
+
+struct quantile_llist_index {
+	struct quantile_llist_element *element;
+	uint32_t index;
+};
+
+/* number of times a value occurs */
+struct rt_valuecount_t {
+	double value;
+	uint32_t count;
+	double percent;
+};
+
+/* reclassification expression */
+struct rt_reclassexpr_t {
+	struct rt_reclassrange {
+		double min;
+		double max;
+		int inc_min; /* include min */
+		int inc_max; /* include max */
+		int exc_min; /* exceed min */
+		int exc_max; /* exceed max */
+	} src, dst;
+};
+
+/* gdal driver information */
+struct rt_gdaldriver_t {
+	int idx;
+	char *short_name;
+	char *long_name;
+	char *create_options;
+};
+
 #endif /* RT_API_H_INCLUDED */
diff --git a/raster/rt_pg/rt_pg.c b/raster/rt_pg/rt_pg.c
index ab7973dde..cfb9d41de 100644
--- a/raster/rt_pg/rt_pg.c
+++ b/raster/rt_pg/rt_pg.c
@@ -236,6 +236,9 @@ Datum RASTER_metadata(PG_FUNCTION_ARGS);
 /* get raster band's meta data */
 Datum RASTER_bandmetadata(PG_FUNCTION_ARGS);
 
+/* determine if two rasters intersect */
+Datum RASTER_intersects(PG_FUNCTION_ARGS);
+
 /* Replace function taken from
  * http://ubuntuforums.org/showthread.php?s=aa6f015109fd7e4c7e30d2fd8b717497&t=141670&page=3
  */
@@ -6497,6 +6500,9 @@ Datum RASTER_asRaster(PG_FUNCTION_ARGS)
 	double *skew_x = NULL;
 	double *skew_y = NULL;
 
+	char **options = NULL;
+	int options_len = 0;
+
 	uint32_t num_bands = 0;
 
 	int srid = SRID_UNKNOWN;
@@ -6907,6 +6913,27 @@ Datum RASTER_asRaster(PG_FUNCTION_ARGS)
 	}
 	POSTGIS_RT_DEBUGF(3, "RASTER_asRaster: skew (x, y) = %f, %f", skew[0], skew[1]);
 
+	/* all touched */
+	if (!PG_ARGISNULL(14) && PG_GETARG_BOOL(14) == TRUE) {
+		if (options_len < 1) {
+			options_len = 1;
+			options = (char **) palloc(sizeof(char *) * options_len);
+		}
+		else {
+			options_len++;
+			options = (char **) repalloc(options, sizeof(char *) * options_len);
+		}
+
+		options[options_len - 1] = palloc(sizeof(char*) * (strlen("ALL_TOUCHED=TRUE") + 1));
+		options[options_len - 1] = "ALL_TOUCHED=TRUE";
+	}
+
+	if (options_len) {
+		options_len++;
+		options = (char **) repalloc(options, sizeof(char *) * options_len);
+		options[options_len - 1] = NULL;
+	}
+
 	/* get geometry's srid */
 #ifdef GSERIALIZED_ON
 	srid = gserialized_get_srid(pggeom);
@@ -6926,6 +6953,7 @@ Datum RASTER_asRaster(PG_FUNCTION_ARGS)
 				pfree(hasnodatas);
 				pfree(nodatavals);
 			}
+			if (options_len) pfree(options);
 
 			lwgeom_free(geom);
 			PG_FREE_IF_COPY(pggeom, 0);
@@ -6979,7 +7007,8 @@ Datum RASTER_asRaster(PG_FUNCTION_ARGS)
 		scale_x, scale_y,
 		ul_xw, ul_yw,
 		grid_xw, grid_yw,
-		skew_x, skew_y
+		skew_x, skew_y,
+		options
 	);
 
 	if (pixtypes_len) pfree(pixtypes);
@@ -6988,9 +7017,10 @@ Datum RASTER_asRaster(PG_FUNCTION_ARGS)
 		pfree(hasnodatas);
 		pfree(nodatavals);
 	}
+	if (options_len) pfree(options);
 
 	if (!rast) {
-		elog(ERROR, "RASTER_asRaster: Could not create rasterize geometry");
+		elog(ERROR, "RASTER_asRaster: Could not rasterize geometry");
 		PG_RETURN_NULL();
 	}
 
@@ -7446,6 +7476,151 @@ Datum RASTER_bandmetadata(PG_FUNCTION_ARGS)
 	PG_RETURN_DATUM(result);
 }
 
+/**
+ * See if two rasters intersect
+ */
+PG_FUNCTION_INFO_V1(RASTER_intersects);
+Datum RASTER_intersects(PG_FUNCTION_ARGS)
+{
+	const int set_count = 2;
+	rt_pgraster *pgrast;
+	rt_raster rast[2] = {NULL};
+	uint32_t bandindex[2] = {0};
+	uint32_t hasbandindex[2] = {0};
+
+	uint32_t i;
+	uint32_t j;
+	uint32_t k;
+	uint32_t numBands;
+	int rtn;
+	int intersects;
+
+	LWPOLY *hull[2] = {NULL};
+	GEOSGeometry *ghull[2] = {NULL};
+
+	for (i = 0, j = 0; i < set_count; i++) {
+		/* pgrast is null, return null */
+		if (PG_ARGISNULL(j)) {
+			for (k = 0; k < i; k++) rt_raster_destroy(rast[k]);
+			PG_RETURN_NULL();
+		}
+		pgrast = (rt_pgraster *) PG_DETOAST_DATUM(PG_GETARG_DATUM(j));
+		j++;
+
+		/* raster */
+		rast[i] = rt_raster_deserialize(pgrast, FALSE);
+		if (!rast[i]) {
+			elog(ERROR, "RASTER_intersects: Could not deserialize the %s raster", i < 1 ? "first" : "second");
+			for (k = 0; k < i; k++) rt_raster_destroy(rast[k]);
+			PG_RETURN_NULL();
+		}
+
+		/* numbands */
+		numBands = rt_raster_get_num_bands(rast[i]);
+		if (numBands < 1) {
+			elog(NOTICE, "The %s raster provided has no bands", i < 1 ? "first" : "second");
+			for (k = 0; k < i; k++) rt_raster_destroy(rast[k]);
+			PG_RETURN_NULL();
+		}
+
+		/* band index */
+		if (!PG_ARGISNULL(j)) {
+			bandindex[i] = PG_GETARG_INT32(j);
+			if (bandindex[i] < 1 || bandindex[i] > numBands) {
+				elog(NOTICE, "Invalid band index (must use 1-based) for the %s raster. Returning NULL", i < 1 ? "first" : "second");
+				for (k = 0; k < i; k++) rt_raster_destroy(rast[k]);
+				PG_RETURN_NULL();
+			}
+			hasbandindex[i] = 1;
+		}
+		else
+			hasbandindex[i] = 0;
+		POSTGIS_RT_DEBUGF(4, "hasbandindex[%d] = %d", i, hasbandindex[i]);
+		POSTGIS_RT_DEBUGF(4, "bandindex[%d] = %d", i, bandindex[i]);
+		j++;
+	}
+
+	/* hasbandindex must be balanced */
+	if (
+		(hasbandindex[0] && !hasbandindex[1]) ||
+		(!hasbandindex[0] && hasbandindex[1])
+	) {
+		elog(NOTICE, "Missing band index.  Band indices must be provided for both rasters if any one is provided");
+		for (k = 0; k < i; k++) rt_raster_destroy(rast[k]);
+		PG_RETURN_NULL();
+	}
+
+	/* SRID must match */
+	if (rt_raster_get_srid(rast[0]) != rt_raster_get_srid(rast[1])) {
+		elog(NOTICE, "The two rasters provided have different SRIDs");
+		for (k = 0; k < set_count; k++) rt_raster_destroy(rast[k]);
+		PG_RETURN_NULL();
+	}
+
+	/* raster extents need to intersect */
+	do {
+		initGEOS(lwnotice, lwgeom_geos_error);
+
+		rtn = 1;
+		for (i = 0; i < 2; i++) {
+			hull[i] = rt_raster_get_convex_hull(rast[i]);
+			if (NULL == hull[i]) {
+				for (j = 0; j < i; j++) {
+					GEOSGeom_destroy(ghull[j]);
+					lwpoly_free(hull[j]);
+				}
+				rtn = 0;
+				break;
+			}
+			ghull[i] = (GEOSGeometry *) LWGEOM2GEOS((LWGEOM *) hull[i]);
+			if (NULL == ghull[i]) {
+				for (j = 0; j < i; j++) {
+					GEOSGeom_destroy(ghull[j]);
+					lwpoly_free(hull[j]);
+				}
+				lwpoly_free(hull[i]);
+				rtn = 0;
+				break;
+			}
+		}
+		if (!rtn) break;
+
+		rtn = GEOSIntersects(ghull[0], ghull[1]);
+
+		for (i = 0; i < 2; i++) {
+			GEOSGeom_destroy(ghull[i]);
+			lwpoly_free(hull[i]);
+		}
+
+		if (rtn != 2) {
+			if (rtn != 1) {
+				for (k = 0; k < set_count; k++) rt_raster_destroy(rast[k]);
+				PG_RETURN_BOOL(0);
+			}
+			/* band isn't specified */
+			else if (!hasbandindex[0]) {
+				for (k = 0; k < set_count; k++) rt_raster_destroy(rast[k]);
+				PG_RETURN_BOOL(1);
+			}
+		}
+	}
+	while (0);
+
+	rtn = rt_raster_intersects(
+		rast[0], (hasbandindex[0] ? bandindex[0] - 1 : -1),
+		rast[1], (hasbandindex[1] ? bandindex[1] - 1 : -1),
+		&intersects
+	);
+	for (k = 0; k < set_count; k++) rt_raster_destroy(rast[k]);
+
+	if (!rtn) {
+		elog(ERROR, "RASTER_intersects: Unable to test for intersection on the two rasters");
+		PG_RETURN_NULL();
+	}
+
+	PG_RETURN_BOOL(intersects);
+}
+
 /* ---------------------------------------------------------------- */
 /*  Memory allocation / error reporting hooks                       */
 /* ---------------------------------------------------------------- */
diff --git a/raster/rt_pg/rtpostgis.sql.in.c b/raster/rt_pg/rtpostgis.sql.in.c
index 7875cea59..61a83fc32 100644
--- a/raster/rt_pg/rtpostgis.sql.in.c
+++ b/raster/rt_pg/rtpostgis.sql.in.c
@@ -1302,7 +1302,8 @@ CREATE OR REPLACE FUNCTION _st_asraster(
 	nodataval double precision[] DEFAULT ARRAY[0]::double precision[],
 	upperleftx double precision DEFAULT NULL, upperlefty double precision DEFAULT NULL,
 	gridx double precision DEFAULT NULL, gridy double precision DEFAULT NULL,
-	skewx double precision DEFAULT 0, skewy double precision DEFAULT 0
+	skewx double precision DEFAULT 0, skewy double precision DEFAULT 0,
+	touched boolean DEFAULT FALSE
 )
 	RETURNS raster
 	AS 'MODULE_PATHNAME', 'RASTER_asRaster'
@@ -1315,10 +1316,11 @@ CREATE OR REPLACE FUNCTION st_asraster(
 	pixeltype text[] DEFAULT ARRAY['8BUI']::text[],
 	value double precision[] DEFAULT ARRAY[1]::double precision[],
 	nodataval double precision[] DEFAULT ARRAY[0]::double precision[],
-	skewx double precision DEFAULT 0, skewy double precision DEFAULT 0
+	skewx double precision DEFAULT 0, skewy double precision DEFAULT 0,
+	touched boolean DEFAULT FALSE
 )
 	RETURNS raster
-	AS $$ SELECT _st_asraster($1, $2, $3, NULL, NULL, $6, $7, $8, NULL, NULL, $4, $5, $9, $10) $$
+	AS $$ SELECT _st_asraster($1, $2, $3, NULL, NULL, $6, $7, $8, NULL, NULL, $4, $5, $9, $10, $11) $$
 	LANGUAGE 'sql' STABLE;
 
 CREATE OR REPLACE FUNCTION st_asraster(
@@ -1328,10 +1330,11 @@ CREATE OR REPLACE FUNCTION st_asraster(
 	value double precision[] DEFAULT ARRAY[1]::double precision[],
 	nodataval double precision[] DEFAULT ARRAY[0]::double precision[],
 	upperleftx double precision DEFAULT NULL, upperlefty double precision DEFAULT NULL,
-	skewx double precision DEFAULT 0, skewy double precision DEFAULT 0
+	skewx double precision DEFAULT 0, skewy double precision DEFAULT 0,
+	touched boolean DEFAULT FALSE
 )
 	RETURNS raster
-	AS $$ SELECT _st_asraster($1, $2, $3, NULL, NULL, $4, $5, $6, $7, $8, NULL, NULL,	$9, $10) $$
+	AS $$ SELECT _st_asraster($1, $2, $3, NULL, NULL, $4, $5, $6, $7, $8, NULL, NULL,	$9, $10, $11) $$
 	LANGUAGE 'sql' STABLE;
 
 CREATE OR REPLACE FUNCTION st_asraster(
@@ -1341,10 +1344,11 @@ CREATE OR REPLACE FUNCTION st_asraster(
 	pixeltype text[] DEFAULT ARRAY['8BUI']::text[],
 	value double precision[] DEFAULT ARRAY[1]::double precision[],
 	nodataval double precision[] DEFAULT ARRAY[0]::double precision[],
-	skewx double precision DEFAULT 0, skewy double precision DEFAULT 0
+	skewx double precision DEFAULT 0, skewy double precision DEFAULT 0,
+	touched boolean DEFAULT FALSE
 )
 	RETURNS raster
-	AS $$ SELECT _st_asraster($1, NULL, NULL, $2, $3, $6, $7, $8, NULL, NULL, $4, $5, $9, $10) $$
+	AS $$ SELECT _st_asraster($1, NULL, NULL, $2, $3, $6, $7, $8, NULL, NULL, $4, $5, $9, $10, $11) $$
 	LANGUAGE 'sql' STABLE;
 
 CREATE OR REPLACE FUNCTION st_asraster(
@@ -1354,10 +1358,11 @@ CREATE OR REPLACE FUNCTION st_asraster(
 	value double precision[] DEFAULT ARRAY[1]::double precision[],
 	nodataval double precision[] DEFAULT ARRAY[0]::double precision[],
 	upperleftx double precision DEFAULT NULL, upperlefty double precision DEFAULT NULL,
-	skewx double precision DEFAULT 0, skewy double precision DEFAULT 0
+	skewx double precision DEFAULT 0, skewy double precision DEFAULT 0,
+	touched boolean DEFAULT FALSE
 )
 	RETURNS raster
-	AS $$ SELECT _st_asraster($1, NULL, NULL, $2, $3, $4, $5, $6, $7, $8, NULL, NULL,	$9, $10) $$
+	AS $$ SELECT _st_asraster($1, NULL, NULL, $2, $3, $4, $5, $6, $7, $8, NULL, NULL,	$9, $10, $11) $$
 	LANGUAGE 'sql' STABLE;
 
 CREATE OR REPLACE FUNCTION st_asraster(
@@ -1367,10 +1372,11 @@ CREATE OR REPLACE FUNCTION st_asraster(
 	pixeltype text,
 	value double precision DEFAULT 1,
 	nodataval double precision DEFAULT 0,
-	skewx double precision DEFAULT 0, skewy double precision DEFAULT 0
+	skewx double precision DEFAULT 0, skewy double precision DEFAULT 0,
+	touched boolean DEFAULT FALSE
 )
 	RETURNS raster
-	AS $$ SELECT _st_asraster($1, $2, $3, NULL, NULL, ARRAY[$6]::text[], ARRAY[$7]::double precision[], ARRAY[$8]::double precision[], NULL, NULL, $4, $5, $9, $10) $$
+	AS $$ SELECT _st_asraster($1, $2, $3, NULL, NULL, ARRAY[$6]::text[], ARRAY[$7]::double precision[], ARRAY[$8]::double precision[], NULL, NULL, $4, $5, $9, $10, $11) $$
 	LANGUAGE 'sql' STABLE;
 
 CREATE OR REPLACE FUNCTION st_asraster(
@@ -1380,10 +1386,11 @@ CREATE OR REPLACE FUNCTION st_asraster(
 	value double precision DEFAULT 1,
 	nodataval double precision DEFAULT 0,
 	upperleftx double precision DEFAULT NULL, upperlefty double precision DEFAULT NULL,
-	skewx double precision DEFAULT 0, skewy double precision DEFAULT 0
+	skewx double precision DEFAULT 0, skewy double precision DEFAULT 0,
+	touched boolean DEFAULT FALSE
 )
 	RETURNS raster
-	AS $$ SELECT _st_asraster($1, $2, $3, NULL, NULL, ARRAY[$4]::text[], ARRAY[$5]::double precision[], ARRAY[$6]::double precision[], $7, $8, NULL, NULL, $9, $10) $$
+	AS $$ SELECT _st_asraster($1, $2, $3, NULL, NULL, ARRAY[$4]::text[], ARRAY[$5]::double precision[], ARRAY[$6]::double precision[], $7, $8, NULL, NULL, $9, $10, $11) $$
 	LANGUAGE 'sql' STABLE;
 
 CREATE OR REPLACE FUNCTION st_asraster(
@@ -1393,10 +1400,11 @@ CREATE OR REPLACE FUNCTION st_asraster(
 	pixeltype text,
 	value double precision DEFAULT 1,
 	nodataval double precision DEFAULT 0,
-	skewx double precision DEFAULT 0, skewy double precision DEFAULT 0
+	skewx double precision DEFAULT 0, skewy double precision DEFAULT 0,
+	touched boolean DEFAULT FALSE
 )
 	RETURNS raster
-	AS $$ SELECT _st_asraster($1, NULL, NULL, $2, $3, ARRAY[$6]::text[], ARRAY[$7]::double precision[], ARRAY[$8]::double precision[], NULL, NULL, $4, $5, $9, $10) $$
+	AS $$ SELECT _st_asraster($1, NULL, NULL, $2, $3, ARRAY[$6]::text[], ARRAY[$7]::double precision[], ARRAY[$8]::double precision[], NULL, NULL, $4, $5, $9, $10, $11) $$
 	LANGUAGE 'sql' STABLE;
 
 CREATE OR REPLACE FUNCTION st_asraster(
@@ -1406,10 +1414,11 @@ CREATE OR REPLACE FUNCTION st_asraster(
 	value double precision DEFAULT 1,
 	nodataval double precision DEFAULT 0,
 	upperleftx double precision DEFAULT NULL, upperlefty double precision DEFAULT NULL,
-	skewx double precision DEFAULT 0, skewy double precision DEFAULT 0
+	skewx double precision DEFAULT 0, skewy double precision DEFAULT 0,
+	touched boolean DEFAULT FALSE
 )
 	RETURNS raster
-	AS $$ SELECT _st_asraster($1, NULL, NULL, $2, $3, ARRAY[$4]::text[], ARRAY[$5]::double precision[], ARRAY[$6]::double precision[], $7, $8, NULL, NULL,$9, $10) $$
+	AS $$ SELECT _st_asraster($1, NULL, NULL, $2, $3, ARRAY[$4]::text[], ARRAY[$5]::double precision[], ARRAY[$6]::double precision[], $7, $8, NULL, NULL,$9, $10, $11) $$
 	LANGUAGE 'sql' STABLE;
 
 CREATE OR REPLACE FUNCTION st_asraster(
@@ -1417,7 +1426,8 @@ CREATE OR REPLACE FUNCTION st_asraster(
 	ref raster,
 	pixeltype text[] DEFAULT ARRAY['8BUI']::text[],
 	value double precision[] DEFAULT ARRAY[1]::double precision[],
-	nodataval double precision[] DEFAULT ARRAY[0]::double precision[]
+	nodataval double precision[] DEFAULT ARRAY[0]::double precision[],
+	touched boolean DEFAULT FALSE
 )
 	RETURNS raster
 	AS $$
@@ -1445,7 +1455,7 @@ CREATE OR REPLACE FUNCTION st_asraster(
 			g := geom;
 		END IF;
 	
-		RETURN _st_asraster(g, scale_x, scale_y, NULL, NULL, $3, $4, $5, NULL, NULL, ul_x, ul_y, skew_x, skew_y);
+		RETURN _st_asraster(g, scale_x, scale_y, NULL, NULL, $3, $4, $5, NULL, NULL, ul_x, ul_y, skew_x, skew_y, $6);
 	END;
 	$$ LANGUAGE 'plpgsql' STABLE;
 
@@ -1454,10 +1464,11 @@ CREATE OR REPLACE FUNCTION st_asraster(
 	ref raster,
 	pixeltype text,
 	value double precision DEFAULT 1,
-	nodataval double precision DEFAULT 0
+	nodataval double precision DEFAULT 0,
+	touched boolean DEFAULT FALSE
 )
 	RETURNS raster
-	AS $$ SELECT st_asraster($1, $2, ARRAY[$3]::text[], ARRAY[$4]::double precision[], ARRAY[$5]::double precision[]) $$
+	AS $$ SELECT st_asraster($1, $2, ARRAY[$3]::text[], ARRAY[$4]::double precision[], ARRAY[$5]::double precision[], $6) $$
 	LANGUAGE 'sql' STABLE;
 
 -----------------------------------------------------------------------
@@ -2105,7 +2116,7 @@ CREATE OR REPLACE FUNCTION st_world2rastercoordx(rast raster, xw float8, yw floa
     LANGUAGE 'plpgsql' IMMUTABLE STRICT;
 
 ---------------------------------------------------------------------------------
--- ST_World2RasteCoordX(rast raster, xw float8)
+-- ST_World2RasterCoordX(rast raster, xw float8)
 -- Returns the column number of the pixels covering the provided world X coordinate
 -- for a non-rotated raster.
 -- This function works even if the world coordinate is outside the raster extent.
@@ -2508,183 +2519,200 @@ CREATE OPERATOR ~ (
     RESTRICT = contsel, JOIN = contjoinsel
     );
 
+-----------------------------------------------------------------------
+-- Raster/Raster Spatial Relationship
+-----------------------------------------------------------------------
+CREATE OR REPLACE FUNCTION _st_intersects(rastA raster, nbandA integer, rastB raster, nbandB integer)
+	RETURNS boolean
+	AS 'MODULE_PATHNAME', 'RASTER_intersects'
+	LANGUAGE 'C' IMMUTABLE;
+
+CREATE OR REPLACE FUNCTION st_intersects(rastA raster, nbandA integer, rastB raster, nbandB integer)
+	RETURNS boolean
+	AS $$ SELECT $1 && $3 AND _st_intersects($1, $2, $3, $4) $$
+	LANGUAGE 'SQL' IMMUTABLE;
+
+CREATE OR REPLACE FUNCTION st_intersects(rastA raster, rastB raster, nbandA integer DEFAULT NULL, nbandB integer DEFAULT NULL)
+	RETURNS boolean
+	AS $$ SELECT $1 && $2 AND _st_intersects($1, $3, $2, $4) $$
+	LANGUAGE 'SQL' IMMUTABLE;
+
 -----------------------------------------------------------------------
 -- Raster/Geometry Spatial Relationship
 -----------------------------------------------------------------------
 
+CREATE OR REPLACE FUNCTION _st_intersects(rast raster, geom geometry, nband integer DEFAULT NULL)
+	RETURNS boolean
+	AS $$
+	DECLARE
+		gr raster;
+		scale double precision;
+	BEGIN
+		IF ST_Intersects(geom, ST_ConvexHull(rast)) IS NOT TRUE THEN
+			RETURN FALSE;
+		ELSEIF nband IS NULL THEN
+			RETURN TRUE;
+		END IF;
+
+		-- scale is set to 1/100th of raster for granularity
+		SELECT CASE WHEN scaley < scalex THEN scaley ELSE scalex END / 100 INTO scale FROM ST_Metadata(rast);
+		gr := ST_AsRaster(geom, scale, scale);
+		IF gr IS NULL THEN
+			RAISE EXCEPTION 'Unable to convert geometry to a raster';
+			RETURN FALSE;
+		END IF;
+
+		RETURN ST_Intersects(rast, 1, gr, 1);
+	END;
+	$$ LANGUAGE 'plpgsql' IMMUTABLE;
+
+CREATE OR REPLACE FUNCTION st_intersects(rast raster, geom geometry, nband integer DEFAULT NULL)
+	RETURNS boolean
+	AS $$ SELECT $1 && $2 AND _st_intersects($1, $2, $3) $$
+	LANGUAGE 'SQL' IMMUTABLE;
+
+CREATE OR REPLACE FUNCTION st_intersects(rast raster, nband integer, geom geometry)
+	RETURNS boolean
+	AS $$ SELECT $1 && $3 AND _st_intersects($1, $3, $2) $$
+	LANGUAGE 'SQL' IMMUTABLE STRICT;
+
 -----------------------------------------------------------------------
--- _st_intersects(geomin geometry, rast raster, band integer, hasnodata boolean)
--- If hasnodata is true, check for the presence of withvalue pixels in the area
+-- _st_intersects(geom geometry, rast raster, nband integer)
+-- If band is provided, check for the presence of withvalue pixels in the area
 -- shared by the raster and the geometry. If only nodata value pixels are found, the
 -- geometry does not intersect with the raster.
 -----------------------------------------------------------------------
-CREATE OR REPLACE FUNCTION _st_intersects(geomin geometry, rast raster, band integer, hasnodata boolean)
-    RETURNS boolean AS
-    $$
-    DECLARE
-        nodata float8 := 0.0;
-        geomintersect geometry;
-        x1w double precision := 0.0;
-        x2w double precision := 0.0;
-        y1w double precision := 0.0;
-        y2w double precision := 0.0;
-        x1 integer := 0;
-        x2 integer := 0;
-        x3 integer := 0;
-        x4 integer := 0;
-        y1 integer := 0;
-        y2 integer := 0;
-        y3 integer := 0;
-        y4 integer := 0;
-        x integer := 0;
-        y integer := 0;
-        xinc integer := 0;
-        yinc integer := 0;
-        pixelval double precision;
-        bintersect boolean := FALSE;
-        gtype text;
-        scale float8;
-				w int;
-				h int;
-    BEGIN
+-- This function can not be STRICT
+CREATE OR REPLACE FUNCTION _st_intersects(geom geometry, rast raster, nband integer DEFAULT NULL)
+	RETURNS boolean AS $$
+	DECLARE
+		hasnodata boolean := TRUE;
+		nodata float8 := 0.0;
+		convexhull geometry;
+		geomintersect geometry;
+		x1w double precision := 0.0;
+		x2w double precision := 0.0;
+		y1w double precision := 0.0;
+		y2w double precision := 0.0;
+		x1 integer := 0;
+		x2 integer := 0;
+		x3 integer := 0;
+		x4 integer := 0;
+		y1 integer := 0;
+		y2 integer := 0;
+		y3 integer := 0;
+		y4 integer := 0;
+		x integer := 0;
+		y integer := 0;
+		xinc integer := 0;
+		yinc integer := 0;
+		pixelval double precision;
+		bintersect boolean := FALSE;
+		gtype text;
+		scale float8;
+		w int;
+		h int;
+	BEGIN
+		convexhull := ST_ConvexHull(rast);
+		IF nband IS NOT NULL THEN
+			SELECT hasnodatavalue INTO hasnodata FROM ST_BandMetaData(rast, nband);
+		END IF;
 
-        -- Get the intersection between with the geometry.
-        -- We will search for withvalue pixel only in this area.
-        geomintersect := st_intersection(geomin, st_convexhull(rast));
+		IF ST_Intersects(geom, convexhull) IS NOT TRUE THEN
+			RETURN FALSE;
+		ELSEIF nband IS NULL OR hasnodata IS FALSE THEN
+			RETURN TRUE;
+		END IF;
 
---RAISE NOTICE 'geomintersect1=%', astext(geomintersect);
+		-- Get the intersection between with the geometry.
+		-- We will search for withvalue pixel only in this area.
+		geomintersect := st_intersection(geom, convexhull);
 
-        -- If the intersection is empty, return false
-        IF st_isempty(geomintersect) THEN
-            RETURN FALSE;
-        END IF;
+--RAISE NOTICE 'geomintersect=%', st_astext(geomintersect);
 
-        -- If the band does not have a nodatavalue, there is no need to search for with value pixels
-        IF NOT hasnodata OR st_bandnodatavalue(rast, band) IS NULL THEN
-            RETURN TRUE;
-        END IF;
+		-- If the intersection is empty, return false
+		IF st_isempty(geomintersect) THEN
+			RETURN FALSE;
+		END IF;
 
-        -- We create a minimalistic buffer around the intersection in order to scan every pixels
-        -- that would touch the edge or intersect with the geometry
-        SELECT (scalex * skewy), width, height INTO scale, w, h FROM ST_Metadata(rast);
-        geomintersect := st_buffer(geomintersect, scale / 1000000);
+		-- We create a minimalistic buffer around the intersection in order to scan every pixels
+		-- that would touch the edge or intersect with the geometry
+		SELECT (scalex * skewy), width, height INTO scale, w, h FROM ST_Metadata(rast);
+		geomintersect := st_buffer(geomintersect, scale / 1000000);
 
 --RAISE NOTICE 'geomintersect2=%', astext(geomintersect);
 
-        -- Find the world coordinates of the bounding box of the intersecting area
-        x1w := st_xmin(geomintersect);
-        y1w := st_ymin(geomintersect);
-        x2w := st_xmax(geomintersect);
-        y2w := st_ymax(geomintersect);
-        nodata := st_bandnodatavalue(rast, band);
+		-- Find the world coordinates of the bounding box of the intersecting area
+		x1w := st_xmin(geomintersect);
+		y1w := st_ymin(geomintersect);
+		x2w := st_xmax(geomintersect);
+		y2w := st_ymax(geomintersect);
+		nodata := st_bandnodatavalue(rast, nband);
 
 --RAISE NOTICE 'x1w=%, y1w=%, x2w=%, y2w=%', x1w, y1w, x2w, y2w;
 
-        -- Convert world coordinates to raster coordinates
-        x1 := st_world2rastercoordx(rast, x1w, y1w);
-        y1 := st_world2rastercoordy(rast, x1w, y1w);
-        x2 := st_world2rastercoordx(rast, x2w, y1w);
-        y2 := st_world2rastercoordy(rast, x2w, y1w);
-        x3 := st_world2rastercoordx(rast, x1w, y2w);
-        y3 := st_world2rastercoordy(rast, x1w, y2w);
-        x4 := st_world2rastercoordx(rast, x2w, y2w);
-        y4 := st_world2rastercoordy(rast, x2w, y2w);
+		-- Convert world coordinates to raster coordinates
+		x1 := st_world2rastercoordx(rast, x1w, y1w);
+		y1 := st_world2rastercoordy(rast, x1w, y1w);
+		x2 := st_world2rastercoordx(rast, x2w, y1w);
+		y2 := st_world2rastercoordy(rast, x2w, y1w);
+		x3 := st_world2rastercoordx(rast, x1w, y2w);
+		y3 := st_world2rastercoordy(rast, x1w, y2w);
+		x4 := st_world2rastercoordx(rast, x2w, y2w);
+		y4 := st_world2rastercoordy(rast, x2w, y2w);
 
 --RAISE NOTICE 'x1=%, y1=%, x2=%, y2=%, x3=%, y3=%, x4=%, y4=%', x1, y1, x2, y2, x3, y3, x4, y4;
 
-        -- Order the raster coordinates for the upcoming FOR loop.
-        x1 := int4smaller(int4smaller(int4smaller(x1, x2), x3), x4);
-        y1 := int4smaller(int4smaller(int4smaller(y1, y2), y3), y4);
-        x2 := int4larger(int4larger(int4larger(x1, x2), x3), x4);
-        y2 := int4larger(int4larger(int4larger(y1, y2), y3), y4);
+		-- Order the raster coordinates for the upcoming FOR loop.
+		x1 := int4smaller(int4smaller(int4smaller(x1, x2), x3), x4);
+		y1 := int4smaller(int4smaller(int4smaller(y1, y2), y3), y4);
+		x2 := int4larger(int4larger(int4larger(x1, x2), x3), x4);
+		y2 := int4larger(int4larger(int4larger(y1, y2), y3), y4);
 
-        -- Make sure the range is not lower than 1.
-        -- This can happen when world coordinate are exactly on the left border
-        -- of the raster and that they do not span on more than one pixel.
-        x1 := int4smaller(int4larger(x1, 1), w);
-        y1 := int4smaller(int4larger(y1, 1), h);
+		-- Make sure the range is not lower than 1.
+		-- This can happen when world coordinate are exactly on the left border
+		-- of the raster and that they do not span on more than one pixel.
+		x1 := int4smaller(int4larger(x1, 1), w);
+		y1 := int4smaller(int4larger(y1, 1), h);
 
-        -- Also make sure the range does not exceed the width and height of the raster.
-        -- This can happen when world coordinate are exactly on the lower right border
-        -- of the raster.
-        x2 := int4smaller(x2, w);
-        y2 := int4smaller(y2, h);
+		-- Also make sure the range does not exceed the width and height of the raster.
+		-- This can happen when world coordinate are exactly on the lower right border
+		-- of the raster.
+		x2 := int4smaller(x2, w);
+		y2 := int4smaller(y2, h);
 
 --RAISE NOTICE 'x1=%, y1=%, x2=%, y2=%', x1, y1, x2, y2;
 
-        -- Search exhaustively for withvalue pixel on a moving 3x3 grid
-        -- (very often more efficient than searching on a mere 1x1 grid)
-        FOR xinc in 0..2 LOOP
-            FOR yinc in 0..2 LOOP
-                FOR x IN x1+xinc..x2 BY 3 LOOP
-                    FOR y IN y1+yinc..y2 BY 3 LOOP
-                        -- Check first if the pixel intersects with the geometry. Often many won't.
-                        bintersect := NOT st_isempty(st_intersection(st_pixelaspolygon(rast, band, x, y), geomin));
-
-                        IF bintersect THEN
-                            -- If the pixel really intersects, check its value. Return TRUE if with value.
-                            pixelval := st_value(rast, band, x, y);
-                            IF pixelval != nodata THEN
-                                RETURN TRUE;
-                            END IF;
-                        END IF;
-                    END LOOP;
-                END LOOP;
-            END LOOP;
-        END LOOP;
-        RETURN FALSE;
-    END;
-    $$
-    LANGUAGE 'plpgsql' IMMUTABLE STRICT;
-
-
--- This function can not be STRICT
-CREATE OR REPLACE FUNCTION st_intersects(geometry, raster, integer)
-    RETURNS boolean AS
-    $$ SELECT $1 && $2 AND _st_intersects($1, $2, $3, TRUE);
-    $$ LANGUAGE 'SQL' IMMUTABLE;
-
--- This function can not be STRICT
-CREATE OR REPLACE FUNCTION st_intersects(raster, integer, geometry)
-    RETURNS boolean AS
-    $$ SELECT st_intersects($3, $1, $2);
-    $$ LANGUAGE 'SQL' IMMUTABLE;
-
--- This function can not be STRICT
-CREATE OR REPLACE FUNCTION st_intersects(geometry, raster)
-    RETURNS boolean AS
-    $$ SELECT st_intersects($1, $2, 1);
-    $$ LANGUAGE 'SQL' IMMUTABLE;
-
--- This function can not be STRICT
-CREATE OR REPLACE FUNCTION st_intersects(raster, geometry)
-    RETURNS boolean AS
-    $$ SELECT st_intersects($2, $1, 1);
-    $$ LANGUAGE 'SQL' IMMUTABLE;
-
--- This function can not be STRICT
-CREATE OR REPLACE FUNCTION st_intersects(geometry, raster, integer, boolean)
-    RETURNS boolean AS
-    $$ SELECT $1 && $2 AND _st_intersects($1, $2, $3, $4);
-    $$ LANGUAGE 'SQL' IMMUTABLE;
-
--- This function can not be STRICT
-CREATE OR REPLACE FUNCTION st_intersects(raster, integer, boolean, geometry)
-    RETURNS boolean AS
-    $$ SELECT st_intersects($4, $1, $2, $3);
-    $$ LANGUAGE 'SQL' IMMUTABLE;
+		-- Search exhaustively for withvalue pixel on a moving 3x3 grid
+		-- (very often more efficient than searching on a mere 1x1 grid)
+		FOR xinc in 0..2 LOOP
+			FOR yinc in 0..2 LOOP
+				FOR x IN x1+xinc..x2 BY 3 LOOP
+					FOR y IN y1+yinc..y2 BY 3 LOOP
+						-- Check first if the pixel intersects with the geometry. Often many won't.
+						bintersect := NOT st_isempty(st_intersection(st_pixelaspolygon(rast, nband, x, y), geom));
+						
+						IF bintersect THEN
+							-- If the pixel really intersects, check its value. Return TRUE if with value.
+							pixelval := st_value(rast, nband, x, y);
+							IF pixelval != nodata THEN
+								RETURN TRUE;
+							END IF;
+						END IF;
+					END LOOP;
+				END LOOP;
+			END LOOP;
+		END LOOP;
 
--- This function can not be STRICT
-CREATE OR REPLACE FUNCTION st_intersects(geometry, raster, boolean)
-    RETURNS boolean AS
-    $$ SELECT st_intersects($1, $2, 1, $3);
-    $$ LANGUAGE 'SQL' IMMUTABLE;
+		RETURN FALSE;
+	END;
+	$$ LANGUAGE 'plpgsql' IMMUTABLE;
 
 -- This function can not be STRICT
-CREATE OR REPLACE FUNCTION st_intersects(raster, boolean, geometry)
-    RETURNS boolean AS
-    $$ SELECT st_intersects($3, $1, 1, $2);
-    $$ LANGUAGE 'SQL' IMMUTABLE;
+CREATE OR REPLACE FUNCTION st_intersects(geom geometry, rast raster, nband integer DEFAULT NULL)
+	RETURNS boolean AS
+	$$ SELECT $1 && $2 AND _st_intersects($1, $2, $3); $$
+	LANGUAGE 'SQL' IMMUTABLE;
 
 -----------------------------------------------------------------------
 -- _st_intersection(geom geometry, rast raster, band integer)
@@ -2695,58 +2723,50 @@ CREATE OR REPLACE FUNCTION st_intersects(raster, boolean, geometry)
 -- Raster nodata value areas are not vectorized and hence do not intersect
 -- with any geometries.
 -----------------------------------------------------------------------
-CREATE OR REPLACE FUNCTION ST_Intersection(geomin geometry, rast raster, band integer)
-  RETURNS SETOF geomval AS
-$BODY$
-    DECLARE
-        intersects boolean := FALSE;
-    BEGIN
-        intersects := ST_Intersects(geomin, rast, band);
-        IF intersects THEN
-            -- Return the intersections of the geometry with the vectorized parts of
-            -- the raster and the values associated with those parts, if really their
-            -- intersection is not empty.
-            RETURN QUERY
-                SELECT intgeom,
-                       val
-                FROM (SELECT ST_Intersection((gv).geom, geomin) AS intgeom,
-                             (gv).val
-                      FROM ST_DumpAsPolygons(rast, band) gv
-                      WHERE ST_Intersects((gv).geom, geomin)
-                     ) foo
-                WHERE NOT ST_IsEmpty(intgeom);
-        ELSE
-            -- If the geometry does not intersect with the raster, return an empty
-            -- geometry and a null value
-            RETURN QUERY
-                SELECT emptygeom,
-                       NULL::float8
-                FROM ST_GeomCollFromText('GEOMETRYCOLLECTION EMPTY', ST_SRID($1)) emptygeom;
-        END IF;
-    END;
-    $BODY$
-  LANGUAGE 'plpgsql' IMMUTABLE STRICT;
+CREATE OR REPLACE FUNCTION st_intersection(geomin geometry, rast raster, band integer DEFAULT 1)
+	RETURNS SETOF geomval AS $$
+	DECLARE
+		intersects boolean := FALSE;
+	BEGIN
+		intersects := ST_Intersects(geomin, rast, band);
+		IF intersects THEN
+			-- Return the intersections of the geometry with the vectorized parts of
+			-- the raster and the values associated with those parts, if really their
+			-- intersection is not empty.
+			RETURN QUERY
+				SELECT
+					intgeom,
+					val
+				FROM (
+					SELECT
+						ST_Intersection((gv).geom, geomin) AS intgeom,
+						(gv).val
+					FROM ST_DumpAsPolygons(rast, band) gv
+					WHERE ST_Intersects((gv).geom, geomin)
+				) foo
+				WHERE NOT ST_IsEmpty(intgeom);
+		ELSE
+			-- If the geometry does not intersect with the raster, return an empty
+			-- geometry and a null value
+			RETURN QUERY
+				SELECT
+					emptygeom,
+					NULL::float8
+				FROM ST_GeomCollFromText('GEOMETRYCOLLECTION EMPTY', ST_SRID($1)) emptygeom;
+		END IF;
+	END;
+	$$
+	LANGUAGE 'plpgsql' IMMUTABLE STRICT;
 
 CREATE OR REPLACE FUNCTION st_intersection(rast raster, geom geometry)
-    RETURNS SETOF geomval AS
-    $$
-        SELECT (gv).geom, (gv).val FROM st_intersection($2, $1, 1) gv;
-    $$
-    LANGUAGE SQL IMMUTABLE STRICT;
-
-CREATE OR REPLACE FUNCTION st_intersection(geom geometry, rast raster)
-    RETURNS SETOF geomval AS
-    $$
-        SELECT (gv).geom, (gv).val FROM st_intersection($1, $2, 1) gv;
-    $$
-    LANGUAGE SQL IMMUTABLE STRICT;
+	RETURNS SETOF geomval
+	AS $$ SELECT (gv).geom, (gv).val FROM st_intersection($2, $1, 1) gv; $$
+	LANGUAGE SQL IMMUTABLE STRICT;
 
 CREATE OR REPLACE FUNCTION st_intersection(rast raster, band integer, geom geometry)
-    RETURNS SETOF geomval AS
-    $$
-        SELECT (gv).geom, (gv).val FROM st_intersection($3, $1, $2) gv;
-    $$
-    LANGUAGE SQL IMMUTABLE STRICT;
+	RETURNS SETOF geomval
+	AS $$ SELECT (gv).geom, (gv).val FROM st_intersection($3, $1, $2) gv; $$
+	LANGUAGE SQL IMMUTABLE STRICT;
 
 ------------------------------------------------------------------------------
 -- RASTER_COLUMNS
diff --git a/raster/test/core/testapi.c b/raster/test/core/testapi.c
index dd7ae4abb..d27f815a7 100644
--- a/raster/test/core/testapi.c
+++ b/raster/test/core/testapi.c
@@ -1138,6 +1138,8 @@ static void testBandStats() {
 	qlls = NULL;
 	qlls_count = 0;
 
+	deepRelease(raster);
+
 	xmax = 100;
 	ymax = 100;
 	raster = rt_raster_new(xmax, ymax);
@@ -1579,7 +1581,8 @@ static void testGDALRasterize() {
 		&scale_x, &scale_y,
 		NULL, NULL,
 		NULL, NULL,
-		NULL, NULL
+		NULL, NULL,
+		NULL
 	);
 
 	free(wkb);
@@ -1594,6 +1597,469 @@ static void testGDALRasterize() {
 	deepRelease(raster);
 }
 
+static void testIntersects() {
+	rt_raster rast1;
+	rt_raster rast2;
+	rt_band band1;
+	rt_band band2;
+	double nodata;
+	int rtn;
+	int intersects;
+
+	/*
+		rast1
+
+		(-1, -1)
+						+-+-+
+						|1|1|
+						+-+-+
+						|1|1|
+						+-+-+
+								(1, 1)
+	*/
+	rast1 = rt_raster_new(2, 2);
+	assert(rast1);
+	rt_raster_set_offsets(rast1, -1, -1);
+
+	band1 = addBand(rast1, PT_8BUI, 1, 0);
+	CHECK(band1);
+	rt_band_set_nodata(band1, 0);
+	rtn = rt_band_set_pixel(band1, 0, 0, 1);
+	rtn = rt_band_set_pixel(band1, 0, 1, 1);
+	rtn = rt_band_set_pixel(band1, 1, 0, 1);
+	rtn = rt_band_set_pixel(band1, 1, 1, 1);
+
+	nodata = rt_band_get_nodata(band1);
+	CHECK_EQUALS(nodata, 0);
+
+	/*
+		rast2
+
+		(0, 0)
+						+-+-+
+						|1|1|
+						+-+-+
+						|1|1|
+						+-+-+
+								(2, 2)
+	*/
+	rast2 = rt_raster_new(2, 2);
+	assert(rast2);
+
+	band2 = addBand(rast2, PT_8BUI, 1, 0);
+	CHECK(band2);
+	rt_band_set_nodata(band2, 0);
+	rtn = rt_band_set_pixel(band2, 0, 0, 1);
+	rtn = rt_band_set_pixel(band2, 0, 1, 1);
+	rtn = rt_band_set_pixel(band2, 1, 0, 1);
+	rtn = rt_band_set_pixel(band2, 1, 1, 1);
+
+	nodata = rt_band_get_nodata(band2);
+	CHECK_EQUALS(nodata, 0);
+
+	rtn = rt_raster_intersects(
+		rast1, 0,
+		rast2, 0,
+		&intersects
+	);
+	CHECK((rtn != 0));
+	CHECK((intersects == 1));
+
+	rtn = rt_raster_intersects(
+		rast1, -1,
+		rast2, -1,
+		&intersects
+	);
+	CHECK((rtn != 0));
+	CHECK((intersects == 1));
+
+	/*
+		rast2
+
+		(0, 0)
+						+-+-+
+						|0|1|
+						+-+-+
+						|1|1|
+						+-+-+
+								(2, 2)
+	*/
+	rtn = rt_band_set_pixel(band2, 0, 0, 0);
+
+	rtn = rt_raster_intersects(
+		rast1, 0,
+		rast2, 0,
+		&intersects
+	);
+	CHECK((rtn != 0));
+	CHECK((intersects == 1));
+
+	/*
+		rast2
+
+		(0, 0)
+						+-+-+
+						|1|0|
+						+-+-+
+						|1|1|
+						+-+-+
+								(2, 2)
+	*/
+	rtn = rt_band_set_pixel(band2, 0, 0, 1);
+	rtn = rt_band_set_pixel(band2, 1, 0, 0);
+
+	rtn = rt_raster_intersects(
+		rast1, 0,
+		rast2, 0,
+		&intersects
+	);
+	CHECK((rtn != 0));
+	CHECK((intersects == 1));
+
+	/*
+		rast2
+
+		(0, 0)
+						+-+-+
+						|0|0|
+						+-+-+
+						|0|1|
+						+-+-+
+								(2, 2)
+	*/
+	rtn = rt_band_set_pixel(band2, 0, 0, 0);
+	rtn = rt_band_set_pixel(band2, 1, 0, 0);
+	rtn = rt_band_set_pixel(band2, 0, 1, 0);
+
+	rtn = rt_raster_intersects(
+		rast1, 0,
+		rast2, 0,
+		&intersects
+	);
+	CHECK((rtn != 0));
+	CHECK((intersects == 1));
+
+	/*
+		rast2
+
+		(0, 0)
+						+-+-+
+						|0|0|
+						+-+-+
+						|0|0|
+						+-+-+
+								(2, 2)
+	*/
+	rtn = rt_band_set_pixel(band2, 0, 0, 0);
+	rtn = rt_band_set_pixel(band2, 1, 0, 0);
+	rtn = rt_band_set_pixel(band2, 0, 1, 0);
+	rtn = rt_band_set_pixel(band2, 1, 1, 0);
+
+	rtn = rt_raster_intersects(
+		rast1, 0,
+		rast2, 0,
+		&intersects
+	);
+	CHECK((rtn != 0));
+	CHECK((intersects != 1));
+
+	/*
+		rast2
+
+		(2, 0)
+						+-+-+
+						|1|1|
+						+-+-+
+						|1|1|
+						+-+-+
+								(4, 2)
+	*/
+	rt_raster_set_offsets(rast2, 2, 0);
+
+	rtn = rt_band_set_pixel(band2, 0, 0, 1);
+	rtn = rt_band_set_pixel(band2, 1, 0, 1);
+	rtn = rt_band_set_pixel(band2, 0, 1, 1);
+	rtn = rt_band_set_pixel(band2, 1, 1, 1);
+
+	rtn = rt_raster_intersects(
+		rast1, 0,
+		rast2, 0,
+		&intersects
+	);
+	CHECK((rtn != 0));
+	CHECK((intersects != 1));
+
+	/*
+		rast2
+
+		(0.1, 0.1)
+						+-+-+
+						|1|1|
+						+-+-+
+						|1|1|
+						+-+-+
+								(0.9, 0.9)
+	*/
+	rt_raster_set_offsets(rast2, 0.1, 0.1);
+	rt_raster_set_scale(rast2, 0.4, 0.4);
+
+	rtn = rt_band_set_pixel(band2, 0, 0, 1);
+	rtn = rt_band_set_pixel(band2, 1, 0, 1);
+	rtn = rt_band_set_pixel(band2, 0, 1, 1);
+	rtn = rt_band_set_pixel(band2, 1, 1, 1);
+
+	rtn = rt_raster_intersects(
+		rast1, 0,
+		rast2, 0,
+		&intersects
+	);
+	CHECK((rtn != 0));
+	CHECK((intersects == 1));
+
+	/*
+		rast2
+
+		(-0.1, 0.1)
+						+-+-+
+						|1|1|
+						+-+-+
+						|1|1|
+						+-+-+
+								(0.9, 0.9)
+	*/
+	rt_raster_set_offsets(rast2, -0.1, 0.1);
+
+	rtn = rt_raster_intersects(
+		rast1, 0,
+		rast2, 0,
+		&intersects
+	);
+	CHECK((rtn != 0));
+	CHECK((intersects == 1));
+
+	deepRelease(rast2);
+
+	/*
+		rast2
+
+		(0, 0)
+						+-+-+-+
+						|1|1|1|
+						+-+-+-+
+						|1|1|1|
+						+-+-+-+
+						|1|1|1|
+						+-+-+-+
+									(3, 3)
+	*/
+	rast2 = rt_raster_new(3, 3);
+	assert(rast2);
+
+	band2 = addBand(rast2, PT_8BUI, 1, 0);
+	CHECK(band2);
+	rt_band_set_nodata(band2, 0);
+	rtn = rt_band_set_pixel(band2, 0, 0, 1);
+	rtn = rt_band_set_pixel(band2, 0, 1, 1);
+	rtn = rt_band_set_pixel(band2, 0, 2, 1);
+	rtn = rt_band_set_pixel(band2, 1, 0, 1);
+	rtn = rt_band_set_pixel(band2, 1, 1, 1);
+	rtn = rt_band_set_pixel(band2, 1, 2, 1);
+	rtn = rt_band_set_pixel(band2, 2, 0, 1);
+	rtn = rt_band_set_pixel(band2, 2, 1, 1);
+	rtn = rt_band_set_pixel(band2, 2, 2, 1);
+
+	nodata = rt_band_get_nodata(band2);
+	CHECK_EQUALS(nodata, 0);
+
+	rtn = rt_raster_intersects(
+		rast1, 0,
+		rast2, 0,
+		&intersects
+	);
+	CHECK((rtn != 0));
+	CHECK((intersects == 1));
+
+	/*
+		rast2
+
+		(-2, -2)
+						+-+-+-+
+						|1|1|1|
+						+-+-+-+
+						|1|1|1|
+						+-+-+-+
+						|1|1|1|
+						+-+-+-+
+									(1, 1)
+	*/
+	rt_raster_set_offsets(rast2, -2, -2);
+
+	rtn = rt_band_set_pixel(band2, 0, 0, 1);
+	rtn = rt_band_set_pixel(band2, 0, 1, 1);
+	rtn = rt_band_set_pixel(band2, 0, 2, 1);
+	rtn = rt_band_set_pixel(band2, 1, 0, 1);
+	rtn = rt_band_set_pixel(band2, 1, 1, 1);
+	rtn = rt_band_set_pixel(band2, 1, 2, 1);
+	rtn = rt_band_set_pixel(band2, 2, 0, 1);
+	rtn = rt_band_set_pixel(band2, 2, 1, 1);
+	rtn = rt_band_set_pixel(band2, 2, 2, 1);
+
+	rtn = rt_raster_intersects(
+		rast1, 0,
+		rast2, 0,
+		&intersects
+	);
+	CHECK((rtn != 0));
+	CHECK((intersects == 1));
+
+	/*
+		rast2
+
+		(-2, -2)
+						+-+-+-+
+						|0|1|1|
+						+-+-+-+
+						|1|0|1|
+						+-+-+-+
+						|1|1|0|
+						+-+-+-+
+									(1, 1)
+	*/
+	rtn = rt_band_set_pixel(band2, 0, 0, 0);
+	rtn = rt_band_set_pixel(band2, 0, 1, 1);
+	rtn = rt_band_set_pixel(band2, 0, 2, 1);
+	rtn = rt_band_set_pixel(band2, 1, 0, 1);
+	rtn = rt_band_set_pixel(band2, 1, 1, 0);
+	rtn = rt_band_set_pixel(band2, 1, 2, 1);
+	rtn = rt_band_set_pixel(band2, 2, 0, 1);
+	rtn = rt_band_set_pixel(band2, 2, 1, 1);
+	rtn = rt_band_set_pixel(band2, 2, 2, 0);
+
+	rtn = rt_raster_intersects(
+		rast1, 0,
+		rast2, 0,
+		&intersects
+	);
+	CHECK((rtn != 0));
+	CHECK((intersects == 1));
+
+	/*
+		rast2
+
+		(-2, -2)
+						+-+-+-+
+						|0|1|1|
+						+-+-+-+
+						|1|0|0|
+						+-+-+-+
+						|1|0|0|
+						+-+-+-+
+									(1, 1)
+	*/
+	rtn = rt_band_set_pixel(band2, 0, 0, 0);
+	rtn = rt_band_set_pixel(band2, 0, 1, 1);
+	rtn = rt_band_set_pixel(band2, 0, 2, 1);
+	rtn = rt_band_set_pixel(band2, 1, 0, 1);
+	rtn = rt_band_set_pixel(band2, 1, 1, 0);
+	rtn = rt_band_set_pixel(band2, 1, 2, 0);
+	rtn = rt_band_set_pixel(band2, 2, 0, 1);
+	rtn = rt_band_set_pixel(band2, 2, 1, 0);
+	rtn = rt_band_set_pixel(band2, 2, 2, 0);
+
+	rtn = rt_raster_intersects(
+		rast1, 0,
+		rast2, 0,
+		&intersects
+	);
+	CHECK((rtn != 0));
+	CHECK((intersects == 1));
+
+	/*
+		rast2
+
+		(-2, -2)
+						+-+-+-+
+						|0|1|0|
+						+-+-+-+
+						|1|0|0|
+						+-+-+-+
+						|0|0|0|
+						+-+-+-+
+									(1, 1)
+	*/
+	rtn = rt_band_set_pixel(band2, 0, 0, 0);
+	rtn = rt_band_set_pixel(band2, 0, 1, 1);
+	rtn = rt_band_set_pixel(band2, 0, 2, 0);
+	rtn = rt_band_set_pixel(band2, 1, 0, 1);
+	rtn = rt_band_set_pixel(band2, 1, 1, 0);
+	rtn = rt_band_set_pixel(band2, 1, 2, 0);
+	rtn = rt_band_set_pixel(band2, 2, 0, 0);
+	rtn = rt_band_set_pixel(band2, 2, 1, 0);
+	rtn = rt_band_set_pixel(band2, 2, 2, 0);
+
+	rtn = rt_raster_intersects(
+		rast1, 0,
+		rast2, 0,
+		&intersects
+	);
+	CHECK((rtn != 0));
+	CHECK((intersects == 1));
+
+	deepRelease(rast2);
+
+	/* skew tests */
+	/* rast2 (skewed by -0.5, 0.5) */
+	rast2 = rt_raster_new(3, 3);
+	assert(rast2);
+	rt_raster_set_skews(rast2, -0.5, 0.5);
+
+	band2 = addBand(rast2, PT_8BUI, 1, 0);
+	CHECK(band2);
+	rt_band_set_nodata(band2, 0);
+	rtn = rt_band_set_pixel(band2, 0, 0, 1);
+	rtn = rt_band_set_pixel(band2, 0, 1, 2);
+	rtn = rt_band_set_pixel(band2, 0, 2, 3);
+	rtn = rt_band_set_pixel(band2, 1, 0, 1);
+	rtn = rt_band_set_pixel(band2, 1, 1, 2);
+	rtn = rt_band_set_pixel(band2, 1, 2, 3);
+	rtn = rt_band_set_pixel(band2, 2, 0, 1);
+	rtn = rt_band_set_pixel(band2, 2, 1, 2);
+	rtn = rt_band_set_pixel(band2, 2, 2, 3);
+
+	rtn = rt_raster_intersects(
+		rast1, 0,
+		rast2, 0,
+		&intersects
+	);
+	CHECK((rtn != 0));
+	CHECK((intersects == 1));
+
+	/* rast2 (skewed by -1, 1) */
+	rt_raster_set_skews(rast2, -1, 1);
+
+	rtn = rt_raster_intersects(
+		rast1, 0,
+		rast2, 0,
+		&intersects
+	);
+	CHECK((rtn != 0));
+	CHECK((intersects == 1));
+
+	/* rast2 (skewed by 1, -1) */
+	rt_raster_set_skews(rast2, 1, -1);
+
+	rtn = rt_raster_intersects(
+		rast1, 0,
+		rast2, 0,
+		&intersects
+	);
+	CHECK((rtn != 0));
+	CHECK((intersects == 1));
+
+	deepRelease(rast2);
+	deepRelease(rast1);
+}
+
 int
 main()
 {
@@ -2020,6 +2486,10 @@ main()
 		testGDALRasterize();
 		printf("Successfully tested rt_raster_gdal_rasterize\n");
 
+		printf("Testing rt_raster_intersects\n");
+		testIntersects();
+		printf("Successfully tested rt_raster_intersects\n");
+
     deepRelease(raster);
 
     return EXIT_SUCCESS;
diff --git a/raster/test/regress/Makefile.in b/raster/test/regress/Makefile.in
index 6b55e3929..f088a2c2c 100644
--- a/raster/test/regress/Makefile.in
+++ b/raster/test/regress/Makefile.in
@@ -110,6 +110,7 @@ TEST_GIST = \
 	
 TEST_SREL = \
 	rt_spatial_relationship.sql \
+	rt_intersects.sql \
 	$(NULL)
 	
 TEST_BUGS = \
diff --git a/raster/test/regress/rt_intersects.sql b/raster/test/regress/rt_intersects.sql
new file mode 100644
index 000000000..63c6e8c67
--- /dev/null
+++ b/raster/test/regress/rt_intersects.sql
@@ -0,0 +1,437 @@
+DROP TABLE IF EXISTS raster_intersects_rast;
+DROP TABLE IF EXISTS raster_intersects_geom;
+CREATE TABLE raster_intersects_rast (
+	rid integer,
+	rast raster
+);
+CREATE TABLE raster_intersects_geom (
+	gid integer,
+	geom geometry
+);
+CREATE OR REPLACE FUNCTION make_test_raster(rid integer, width integer DEFAULT 2, height integer DEFAULT 2, ul_x double precision DEFAULT 0, ul_y double precision DEFAULT 0, skew_x double precision DEFAULT 0, skew_y double precision DEFAULT 0)
+	RETURNS void
+	AS $$
+	DECLARE
+		x int;
+		y int;
+		rast raster;
+	BEGIN
+		rast := ST_MakeEmptyRaster(width, height, ul_x, ul_y, 1, 1, skew_x, skew_y, -1);
+		rast := ST_AddBand(rast, 1, '8BUI', 1, 0);
+
+
+		INSERT INTO raster_intersects_rast VALUES (rid, rast);
+
+		RETURN;
+	END;
+	$$ LANGUAGE 'plpgsql';
+SELECT make_test_raster(0, 2, 2, -1, -1);
+SELECT make_test_raster(1, 2, 2);
+SELECT make_test_raster(2, 3, 3);
+DROP FUNCTION make_test_raster(integer, integer, integer, double precision, double precision, double precision, double precision);
+
+INSERT INTO raster_intersects_rast VALUES (10, (
+	SELECT
+		ST_SetValue(rast, 1, 1, 1, 0)
+	FROM raster_intersects_rast
+	WHERE rid = 1
+));
+INSERT INTO raster_intersects_rast VALUES (11, (
+	SELECT
+		ST_SetValue(rast, 1, 2, 1, 0)
+	FROM raster_intersects_rast
+	WHERE rid = 1
+));
+INSERT INTO raster_intersects_rast VALUES (12, (
+	SELECT
+		ST_SetValue(
+			ST_SetValue(
+				ST_SetValue(rast, 1, 1, 1, 0),
+				1, 2, 1, 0
+			),
+			1, 1, 2, 0
+		)
+	FROM raster_intersects_rast
+	WHERE rid = 1
+));
+INSERT INTO raster_intersects_rast VALUES (13, (
+	SELECT
+		ST_SetValue(
+			ST_SetValue(
+				ST_SetValue(
+					ST_SetValue(rast, 1, 1, 1, 0),
+					1, 2, 1, 0
+				),
+				1, 1, 2, 0
+			),
+			1, 2, 2, 0
+		)
+	FROM raster_intersects_rast
+	WHERE rid = 1
+));
+INSERT INTO raster_intersects_rast VALUES (14, (
+	SELECT
+		ST_SetUpperLeft(rast, 2, 0)
+	FROM raster_intersects_rast
+	WHERE rid = 1
+));
+INSERT INTO raster_intersects_rast VALUES (15, (
+	SELECT
+		ST_SetScale(
+			ST_SetUpperLeft(rast, 0.1, 0.1),
+			0.4, 0.4
+		)
+	FROM raster_intersects_rast
+	WHERE rid = 1
+));
+INSERT INTO raster_intersects_rast VALUES (16, (
+	SELECT
+		ST_SetScale(
+			ST_SetUpperLeft(rast, -0.1, 0.1),
+			0.4, 0.4
+		)
+	FROM raster_intersects_rast
+	WHERE rid = 1
+));
+
+INSERT INTO raster_intersects_rast VALUES (20, (
+	SELECT
+		ST_SetUpperLeft(rast, -2, -2)
+	FROM raster_intersects_rast
+	WHERE rid = 2
+));
+INSERT INTO raster_intersects_rast VALUES (21, (
+	SELECT
+		ST_SetValue(
+			ST_SetValue(
+				ST_SetValue(rast, 1, 1, 1, 0),
+				1, 2, 2, 0
+			),
+			1, 3, 3, 0
+		)
+	FROM raster_intersects_rast
+	WHERE rid = 20
+));
+INSERT INTO raster_intersects_rast VALUES (22, (
+	SELECT
+		ST_SetValue(
+			ST_SetValue(
+				rast, 1, 3, 2, 0
+			),
+			1, 2, 3, 0
+		)
+	FROM raster_intersects_rast
+	WHERE rid = 21
+));
+INSERT INTO raster_intersects_rast VALUES (23, (
+	SELECT
+		ST_SetValue(
+			ST_SetValue(
+				rast, 1, 3, 1, 0
+			),
+			1, 1, 3, 0
+		)
+	FROM raster_intersects_rast
+	WHERE rid = 22
+));
+
+INSERT INTO raster_intersects_rast VALUES (30, (
+	SELECT
+		ST_SetSkew(rast, -0.5, 0.5)
+	FROM raster_intersects_rast
+	WHERE rid = 2
+));
+INSERT INTO raster_intersects_rast VALUES (31, (
+	SELECT
+		ST_SetSkew(rast, -1, 1)
+	FROM raster_intersects_rast
+	WHERE rid = 2
+));
+INSERT INTO raster_intersects_rast VALUES (32, (
+	SELECT
+		ST_SetSkew(rast, 1, -1)
+	FROM raster_intersects_rast
+	WHERE rid = 2
+));
+
+SELECT
+	'1.1',
+	r1.rid,
+	r2.rid,
+	ST_Intersects(r1.rast, r2.rast)
+FROM raster_intersects_rast r1
+JOIN raster_intersects_rast r2
+	ON r1.rid != r2.rid
+WHERE r1.rid = 0;
+
+SELECT
+	'1.2',
+	r1.rid,
+	r2.rid,
+	ST_Intersects(r1.rast, 1, r2.rast, 1)
+FROM raster_intersects_rast r1
+JOIN raster_intersects_rast r2
+	ON r1.rid != r2.rid
+WHERE r1.rid = 0;
+
+-- point
+INSERT INTO raster_intersects_geom VALUES (
+	1, (
+		SELECT ST_MakePoint(0, 0)
+	)
+), (
+	2, (
+		SELECT ST_MakePoint(0.1, 0.1)
+	)
+), (
+	3, (
+		SELECT ST_MakePoint(-0.1, -0.1)
+	)
+), (
+	4, (
+		SELECT ST_MakePoint(-1, -1)
+	)
+), (
+	5, (
+		SELECT ST_MakePoint(-1.1, -1)
+	)
+), (
+	6, (
+		SELECT ST_MakePoint(-1, -1.1)
+	)
+), (
+	7, (
+		SELECT ST_MakePoint(-1.5, -1.5)
+	)
+), (
+	8, (
+		SELECT ST_MakePoint(3, 3)
+	)
+);
+
+-- multipoint
+INSERT INTO raster_intersects_geom VALUES (
+	11, (
+		SELECT ST_Collect(geom) FROM raster_intersects_geom WHERE gid BETWEEN 1 AND 10
+	)
+), (
+	12, (
+		SELECT ST_Collect(geom) FROM raster_intersects_geom WHERE gid BETWEEN 3 AND 10
+	)
+), (
+	13, (
+		SELECT ST_Collect(geom) FROM raster_intersects_geom WHERE gid BETWEEN 4 AND 10
+	)
+), (
+	14, (
+		SELECT ST_Collect(geom) FROM raster_intersects_geom WHERE gid BETWEEN 5 AND 10
+	)
+), (
+	15, (
+		SELECT ST_Collect(geom) FROM raster_intersects_geom WHERE gid BETWEEN 6 AND 10
+	)
+);
+
+-- linestring
+INSERT INTO raster_intersects_geom VALUES (
+	21, (
+		SELECT ST_MakeLine(ARRAY[
+			ST_MakePoint(1, 1),
+			ST_MakePoint(1, 0)
+		])
+	)
+), (
+	22, (
+		SELECT ST_MakeLine(ARRAY[
+			ST_MakePoint(-1, -1),
+			ST_MakePoint(1, 1),
+			ST_MakePoint(1, 0)
+		])
+	)
+), (
+	23, (
+		SELECT ST_MakeLine(ARRAY[
+			ST_MakePoint(-1, -1),
+			ST_MakePoint(-1, 1),
+			ST_MakePoint(1, 1),
+			ST_MakePoint(1, -1)
+		])
+	)
+), (
+	24, (
+		SELECT ST_MakeLine(ARRAY[
+			ST_MakePoint(-1.1, 1.1),
+			ST_MakePoint(1.1, 1.1),
+			ST_MakePoint(1.1, -1.1),
+			ST_MakePoint(-1.1, -1.1),
+			ST_MakePoint(-1.1, 1.1)
+		])
+	)
+), (
+	25, (
+		SELECT ST_MakeLine(ARRAY[
+			ST_MakePoint(-2, 1),
+			ST_MakePoint(1, 2),
+			ST_MakePoint(2, -1),
+			ST_MakePoint(-1, -2),
+			ST_MakePoint(-2, 1)
+		])
+	)
+), (
+	26, (
+		SELECT ST_MakeLine(ARRAY[
+			ST_MakePoint(-0.5, 0.5),
+			ST_MakePoint(0, 0.5),
+			ST_MakePoint(0, 0),
+			ST_MakePoint(0, -0.5),
+			ST_MakePoint(-0.5, 0.5)
+		])
+	)
+), (
+	27, (
+		SELECT ST_MakeLine(ARRAY[
+			ST_MakePoint(0.5, 0.5),
+			ST_MakePoint(1, 1),
+			ST_MakePoint(1, 0),
+			ST_MakePoint(0.5, 0.5)
+		])
+	)
+), (
+	28, (
+		SELECT ST_MakeLine(ARRAY[
+			ST_MakePoint(1, 1),
+			ST_MakePoint(0, 2),
+			ST_MakePoint(1, 2),
+			ST_MakePoint(1, 1)
+		])
+	)
+), (
+	29, (
+		SELECT ST_MakeLine(ARRAY[
+			ST_MakePoint(0, 2),
+			ST_MakePoint(1, 2),
+			ST_MakePoint(1, 4),
+			ST_MakePoint(0, 2)
+		])
+	)
+);
+
+-- polygon
+INSERT INTO raster_intersects_geom VALUES (
+	31, (
+		SELECT ST_MakePolygon(geom) FROM raster_intersects_geom WHERE gid = 24
+	)
+), (
+	32, (
+		SELECT ST_MakePolygon(geom) FROM raster_intersects_geom WHERE gid = 25
+	)
+), (
+	33, (
+		SELECT ST_MakePolygon(geom) FROM raster_intersects_geom WHERE gid = 26
+	)
+), (
+	34, (
+		SELECT ST_MakePolygon(geom) FROM raster_intersects_geom WHERE gid = 27
+	)
+), (
+	35, (
+		SELECT ST_MakePolygon(geom) FROM raster_intersects_geom WHERE gid = 28
+	)
+), (
+	36, (
+		SELECT ST_MakePolygon(geom) FROM raster_intersects_geom WHERE gid = 29
+	)
+);
+
+-- multipolygon
+INSERT INTO raster_intersects_geom VALUES (
+	41, (
+		SELECT ST_Collect(geom) FROM raster_intersects_geom WHERE gid BETWEEN 31 and 40
+	)
+), (
+	42, (
+		SELECT ST_Collect(geom) FROM raster_intersects_geom WHERE gid BETWEEN 32 and 40
+	)
+), (
+	43, (
+		SELECT ST_Collect(geom) FROM raster_intersects_geom WHERE gid BETWEEN 33 and 40
+	)
+), (
+	44, (
+		SELECT ST_Collect(geom) FROM raster_intersects_geom WHERE gid BETWEEN 34 and 40
+	)
+), (
+	45, (
+		SELECT ST_Collect(geom) FROM raster_intersects_geom WHERE gid BETWEEN 35 and 40
+	)
+), (
+	46, (
+		SELECT ST_Collect(geom) FROM raster_intersects_geom WHERE gid BETWEEN 36 and 40
+	)
+);
+
+SELECT
+	'2.1',
+	r1.rid,
+	g1.gid,
+	ST_GeometryType(g1.geom),
+	ST_Intersects(r1.rast, g1.geom)
+FROM raster_intersects_rast r1
+JOIN raster_intersects_geom g1
+	ON 1 = 1
+WHERE r1.rid = 0;
+
+SELECT
+	'2.2',
+	r1.rid,
+	g1.gid,
+	ST_GeometryType(g1.geom),
+	ST_Intersects(g1.geom, r1.rast)
+FROM raster_intersects_rast r1
+JOIN raster_intersects_geom g1
+	ON 1 = 1
+WHERE r1.rid = 0;
+
+SELECT
+	'2.3',
+	r1.rid,
+	g1.gid,
+	ST_GeometryType(g1.geom),
+	ST_Intersects(r1.rast, g1.geom)
+FROM raster_intersects_rast r1
+JOIN raster_intersects_geom g1
+	ON 1 = 1
+WHERE r1.rid = 2;
+
+SELECT
+	'2.4',
+	r1.rid,
+	g1.gid,
+	ST_GeometryType(g1.geom),
+	ST_Intersects(g1.geom, r1.rast)
+FROM raster_intersects_rast r1
+JOIN raster_intersects_geom g1
+	ON 1 = 1
+WHERE r1.rid = 2;
+
+SELECT
+	'2.5',
+	r1.rid,
+	g1.gid,
+	ST_GeometryType(g1.geom),
+	ST_Intersects(r1.rast, g1.geom, 1)
+FROM raster_intersects_rast r1
+JOIN raster_intersects_geom g1
+	ON 1 = 1
+WHERE r1.rid = 0;
+
+SELECT
+	'2.6',
+	r1.rid,
+	g1.gid,
+	ST_GeometryType(g1.geom),
+	ST_Intersects(r1.rast, g1.geom, 1)
+FROM raster_intersects_rast r1
+JOIN raster_intersects_geom g1
+	ON 1 = 1
+WHERE r1.rid = 2;
diff --git a/raster/test/regress/rt_intersects_expected b/raster/test/regress/rt_intersects_expected
new file mode 100644
index 000000000..e47990fc3
--- /dev/null
+++ b/raster/test/regress/rt_intersects_expected
@@ -0,0 +1,238 @@
+NOTICE:  table "raster_intersects_rast" does not exist, skipping
+NOTICE:  table "raster_intersects_geom" does not exist, skipping
+1.1|0|1|t
+1.1|0|2|t
+1.1|0|10|t
+1.1|0|11|t
+1.1|0|12|t
+1.1|0|13|t
+1.1|0|14|f
+1.1|0|15|t
+1.1|0|16|t
+1.1|0|20|t
+1.1|0|21|t
+1.1|0|22|t
+1.1|0|23|t
+1.1|0|30|t
+1.1|0|31|t
+1.1|0|32|t
+1.2|0|1|t
+1.2|0|2|t
+1.2|0|10|t
+1.2|0|11|t
+1.2|0|12|t
+1.2|0|13|f
+1.2|0|14|f
+1.2|0|15|t
+1.2|0|16|t
+1.2|0|20|t
+1.2|0|21|t
+1.2|0|22|t
+1.2|0|23|t
+1.2|0|30|t
+1.2|0|31|t
+1.2|0|32|t
+2.1|0|1|ST_Point|t
+2.1|0|2|ST_Point|t
+2.1|0|3|ST_Point|t
+2.1|0|4|ST_Point|t
+2.1|0|5|ST_Point|f
+2.1|0|6|ST_Point|f
+2.1|0|7|ST_Point|f
+2.1|0|8|ST_Point|f
+2.1|0|11|ST_MultiPoint|t
+2.1|0|12|ST_MultiPoint|t
+2.1|0|13|ST_MultiPoint|t
+2.1|0|14|ST_MultiPoint|f
+2.1|0|15|ST_MultiPoint|f
+2.1|0|21|ST_LineString|t
+2.1|0|22|ST_LineString|t
+2.1|0|23|ST_LineString|t
+2.1|0|24|ST_LineString|f
+2.1|0|25|ST_LineString|f
+2.1|0|26|ST_LineString|t
+2.1|0|27|ST_LineString|t
+2.1|0|28|ST_LineString|t
+2.1|0|29|ST_LineString|f
+2.1|0|31|ST_Polygon|t
+2.1|0|32|ST_Polygon|t
+2.1|0|33|ST_Polygon|t
+2.1|0|34|ST_Polygon|t
+2.1|0|35|ST_Polygon|t
+2.1|0|36|ST_Polygon|f
+2.1|0|41|ST_MultiPolygon|t
+2.1|0|42|ST_MultiPolygon|t
+2.1|0|43|ST_MultiPolygon|t
+2.1|0|44|ST_MultiPolygon|t
+2.1|0|45|ST_MultiPolygon|t
+2.1|0|46|ST_MultiPolygon|f
+2.2|0|1|ST_Point|t
+2.2|0|2|ST_Point|t
+2.2|0|3|ST_Point|t
+2.2|0|4|ST_Point|t
+2.2|0|5|ST_Point|f
+2.2|0|6|ST_Point|f
+2.2|0|7|ST_Point|f
+2.2|0|8|ST_Point|f
+2.2|0|11|ST_MultiPoint|t
+2.2|0|12|ST_MultiPoint|t
+2.2|0|13|ST_MultiPoint|t
+2.2|0|14|ST_MultiPoint|f
+2.2|0|15|ST_MultiPoint|f
+2.2|0|21|ST_LineString|t
+2.2|0|22|ST_LineString|t
+2.2|0|23|ST_LineString|t
+2.2|0|24|ST_LineString|f
+2.2|0|25|ST_LineString|f
+2.2|0|26|ST_LineString|t
+2.2|0|27|ST_LineString|t
+2.2|0|28|ST_LineString|t
+2.2|0|29|ST_LineString|f
+2.2|0|31|ST_Polygon|t
+2.2|0|32|ST_Polygon|t
+2.2|0|33|ST_Polygon|t
+2.2|0|34|ST_Polygon|t
+2.2|0|35|ST_Polygon|t
+2.2|0|36|ST_Polygon|f
+2.2|0|41|ST_MultiPolygon|t
+2.2|0|42|ST_MultiPolygon|t
+2.2|0|43|ST_MultiPolygon|t
+2.2|0|44|ST_MultiPolygon|t
+2.2|0|45|ST_MultiPolygon|t
+2.2|0|46|ST_MultiPolygon|f
+2.3|2|1|ST_Point|t
+2.3|2|2|ST_Point|t
+2.3|2|3|ST_Point|f
+2.3|2|4|ST_Point|f
+2.3|2|5|ST_Point|f
+2.3|2|6|ST_Point|f
+2.3|2|7|ST_Point|f
+2.3|2|8|ST_Point|t
+2.3|2|11|ST_MultiPoint|t
+2.3|2|12|ST_MultiPoint|t
+2.3|2|13|ST_MultiPoint|t
+2.3|2|14|ST_MultiPoint|t
+2.3|2|15|ST_MultiPoint|t
+2.3|2|21|ST_LineString|t
+2.3|2|22|ST_LineString|t
+2.3|2|23|ST_LineString|t
+2.3|2|24|ST_LineString|t
+2.3|2|25|ST_LineString|t
+2.3|2|26|ST_LineString|t
+2.3|2|27|ST_LineString|t
+2.3|2|28|ST_LineString|t
+2.3|2|29|ST_LineString|t
+2.3|2|31|ST_Polygon|t
+2.3|2|32|ST_Polygon|t
+2.3|2|33|ST_Polygon|t
+2.3|2|34|ST_Polygon|t
+2.3|2|35|ST_Polygon|t
+2.3|2|36|ST_Polygon|t
+2.3|2|41|ST_MultiPolygon|t
+2.3|2|42|ST_MultiPolygon|t
+2.3|2|43|ST_MultiPolygon|t
+2.3|2|44|ST_MultiPolygon|t
+2.3|2|45|ST_MultiPolygon|t
+2.3|2|46|ST_MultiPolygon|t
+2.4|2|1|ST_Point|t
+2.4|2|2|ST_Point|t
+2.4|2|3|ST_Point|f
+2.4|2|4|ST_Point|f
+2.4|2|5|ST_Point|f
+2.4|2|6|ST_Point|f
+2.4|2|7|ST_Point|f
+2.4|2|8|ST_Point|t
+2.4|2|11|ST_MultiPoint|t
+2.4|2|12|ST_MultiPoint|t
+2.4|2|13|ST_MultiPoint|t
+2.4|2|14|ST_MultiPoint|t
+2.4|2|15|ST_MultiPoint|t
+2.4|2|21|ST_LineString|t
+2.4|2|22|ST_LineString|t
+2.4|2|23|ST_LineString|t
+2.4|2|24|ST_LineString|t
+2.4|2|25|ST_LineString|t
+2.4|2|26|ST_LineString|t
+2.4|2|27|ST_LineString|t
+2.4|2|28|ST_LineString|t
+2.4|2|29|ST_LineString|t
+2.4|2|31|ST_Polygon|t
+2.4|2|32|ST_Polygon|t
+2.4|2|33|ST_Polygon|t
+2.4|2|34|ST_Polygon|t
+2.4|2|35|ST_Polygon|t
+2.4|2|36|ST_Polygon|t
+2.4|2|41|ST_MultiPolygon|t
+2.4|2|42|ST_MultiPolygon|t
+2.4|2|43|ST_MultiPolygon|t
+2.4|2|44|ST_MultiPolygon|t
+2.4|2|45|ST_MultiPolygon|t
+2.4|2|46|ST_MultiPolygon|t
+2.5|0|1|ST_Point|t
+2.5|0|2|ST_Point|t
+2.5|0|3|ST_Point|t
+2.5|0|4|ST_Point|t
+2.5|0|5|ST_Point|f
+2.5|0|6|ST_Point|f
+2.5|0|7|ST_Point|f
+2.5|0|8|ST_Point|f
+2.5|0|11|ST_MultiPoint|t
+2.5|0|12|ST_MultiPoint|t
+2.5|0|13|ST_MultiPoint|t
+2.5|0|14|ST_MultiPoint|f
+2.5|0|15|ST_MultiPoint|f
+2.5|0|21|ST_LineString|t
+2.5|0|22|ST_LineString|t
+2.5|0|23|ST_LineString|t
+2.5|0|24|ST_LineString|f
+2.5|0|25|ST_LineString|f
+2.5|0|26|ST_LineString|t
+2.5|0|27|ST_LineString|t
+2.5|0|28|ST_LineString|t
+2.5|0|29|ST_LineString|f
+2.5|0|31|ST_Polygon|t
+2.5|0|32|ST_Polygon|t
+2.5|0|33|ST_Polygon|t
+2.5|0|34|ST_Polygon|t
+2.5|0|35|ST_Polygon|f
+2.5|0|36|ST_Polygon|f
+2.5|0|41|ST_MultiPolygon|t
+2.5|0|42|ST_MultiPolygon|t
+2.5|0|43|ST_MultiPolygon|t
+2.5|0|44|ST_MultiPolygon|t
+2.5|0|45|ST_MultiPolygon|f
+2.5|0|46|ST_MultiPolygon|f
+2.6|2|1|ST_Point|t
+2.6|2|2|ST_Point|t
+2.6|2|3|ST_Point|f
+2.6|2|4|ST_Point|f
+2.6|2|5|ST_Point|f
+2.6|2|6|ST_Point|f
+2.6|2|7|ST_Point|f
+2.6|2|8|ST_Point|t
+2.6|2|11|ST_MultiPoint|t
+2.6|2|12|ST_MultiPoint|t
+2.6|2|13|ST_MultiPoint|t
+2.6|2|14|ST_MultiPoint|t
+2.6|2|15|ST_MultiPoint|t
+2.6|2|21|ST_LineString|t
+2.6|2|22|ST_LineString|t
+2.6|2|23|ST_LineString|t
+2.6|2|24|ST_LineString|t
+2.6|2|25|ST_LineString|t
+2.6|2|26|ST_LineString|t
+2.6|2|27|ST_LineString|t
+2.6|2|28|ST_LineString|t
+2.6|2|29|ST_LineString|t
+2.6|2|31|ST_Polygon|t
+2.6|2|32|ST_Polygon|t
+2.6|2|33|ST_Polygon|t
+2.6|2|34|ST_Polygon|t
+2.6|2|35|ST_Polygon|t
+2.6|2|36|ST_Polygon|t
+2.6|2|41|ST_MultiPolygon|t
+2.6|2|42|ST_MultiPolygon|t
+2.6|2|43|ST_MultiPolygon|t
+2.6|2|44|ST_MultiPolygon|t
+2.6|2|45|ST_MultiPolygon|t
+2.6|2|46|ST_MultiPolygon|t
diff --git a/raster/test/regress/rt_spatial_relationship.sql b/raster/test/regress/rt_spatial_relationship.sql
index addd95307..2205d8dab 100644
--- a/raster/test/regress/rt_spatial_relationship.sql
+++ b/raster/test/regress/rt_spatial_relationship.sql
@@ -273,19 +273,19 @@ VALUES ( 46, 2, st_setsrid(st_makeline(ARRAY[st_makepoint(
 -- Test 1 - st_intersect(geometry, raster) 
 -----------------------------------------------------------------------
 
-SELECT 'test 1.1', rid, geomid, st_intersects(geom, rast) AS intersect
+SELECT 'test 1.1', rid, geomid, ST_Intersects(geom, rast) AS intersect
 	FROM rt_spatial_relationship_test, rt_spatial_relationship_test_geom
 	WHERE forrast = rid ORDER BY rid, geomid;
 	
-SELECT 'test 1.2', rid, geomid, st_intersects(geom, rast) AS intersect
+SELECT 'test 1.2', rid, geomid, ST_Intersects(geom, rast) AS intersect
 	FROM rt_spatial_relationship_test, rt_spatial_relationship_test_geom
 	WHERE forrast = rid AND ST_Intersects(geom, rast, 1) ORDER BY rid, geomid;
 	
-SELECT 'test 1.3', rid, geomid, st_intersects(geom, st_setbandnodatavalue(rast, NULL))
+SELECT 'test 1.3', rid, geomid, ST_Intersects(geom, st_setbandnodatavalue(rast, NULL), 1)
 	FROM rt_spatial_relationship_test, rt_spatial_relationship_test_geom
 	WHERE forrast = rid ORDER BY rid, geomid;
 
-SELECT 'test 1.4', rid, geomid, st_intersects(geom, st_setbandnodatavalue(rast, NULL))
+SELECT 'test 1.4', rid, geomid, ST_Intersects(geom, st_setbandnodatavalue(rast, NULL))
 	FROM rt_spatial_relationship_test, rt_spatial_relationship_test_geom
 	WHERE forrast = rid AND ST_Intersects(geom, rast, 1) ORDER BY rid, geomid;
 
@@ -301,7 +301,7 @@ FROM (SELECT *, ST_Intersection(geom, rast) gv
 SELECT 'test 2.2', rid, geomid, ST_AsText((gv).geom), (gv).val
 FROM (SELECT *, ST_Intersection(geom, rast) gv
       FROM rt_spatial_relationship_test, rt_spatial_relationship_test_geom
-      WHERE forrast = rid AND ST_Intersects(geom, rast)
+      WHERE forrast = rid AND ST_Intersects(geom, rast, 1)
      ) foo;
 
 SELECT 'test 2.3', rid, geomid, ST_AsText((gv).geom), (gv).val
@@ -313,6 +313,6 @@ FROM (SELECT *, ST_Intersection(geom, st_setbandnodatavalue(rast, NULL)) gv
 SELECT 'test 2.4', rid, geomid, ST_AsText((gv).geom), (gv).val
 FROM (SELECT *, ST_Intersection(geom, st_setbandnodatavalue(rast, NULL)) gv
       FROM rt_spatial_relationship_test, rt_spatial_relationship_test_geom
-      WHERE forrast = rid AND ST_Intersects(geom, st_setbandnodatavalue(rast, NULL))
+      WHERE forrast = rid AND ST_Intersects(geom, st_setbandnodatavalue(rast, NULL), 1)
      ) foo;
 
diff --git a/raster/test/regress/rt_spatial_relationship_expected b/raster/test/regress/rt_spatial_relationship_expected
index 3dff47213..66e61d0f3 100644
--- a/raster/test/regress/rt_spatial_relationship_expected
+++ b/raster/test/regress/rt_spatial_relationship_expected
@@ -1,28 +1,28 @@
 test 1.1|1|1|t
 test 1.1|1|2|t
 test 1.1|1|3|f
-test 1.1|1|4|f
+test 1.1|1|4|t
 test 1.1|1|5|t
 test 1.1|1|6|t
 test 1.1|1|7|f
 test 1.1|1|11|f
 test 1.1|1|12|t
 test 1.1|1|13|t
-test 1.1|1|14|f
+test 1.1|1|14|t
 test 1.1|1|15|t
-test 1.1|1|16|f
+test 1.1|1|16|t
 test 1.1|1|17|t
 test 1.1|1|18|t
 test 1.1|2|31|t
 test 1.1|2|32|t
 test 1.1|2|33|f
-test 1.1|2|34|f
-test 1.1|2|35|f
-test 1.1|2|36|f
+test 1.1|2|34|t
+test 1.1|2|35|t
+test 1.1|2|36|t
 test 1.1|2|41|f
 test 1.1|2|42|t
 test 1.1|2|43|t
-test 1.1|2|44|f
+test 1.1|2|44|t
 test 1.1|2|45|f
 test 1.1|2|46|f
 test 1.2|1|1|t
@@ -58,7 +58,7 @@ test 1.3|2|32|t
 test 1.3|2|33|f
 test 1.3|2|34|t
 test 1.3|2|35|t
-test 1.3|2|36|f
+test 1.3|2|36|t
 test 1.3|2|41|f
 test 1.3|2|42|t
 test 1.3|2|43|t
@@ -153,7 +153,6 @@ test 2.3|2|32|POINT(-75.5532328537098 49.2824585505576)|1
 test 2.3|2|33|GEOMETRYCOLLECTION EMPTY|
 test 2.3|2|34|POINT(-75.5518328537098 49.2814585505576)|4
 test 2.3|2|35|POINT(-75.5523150760919 49.2818643977074)|4
-test 2.3|2|36|GEOMETRYCOLLECTION EMPTY|
 test 2.3|2|41|GEOMETRYCOLLECTION EMPTY|
 test 2.3|2|42|LINESTRING(-75.5533120424461 49.2823793618213,-75.5531972042327 49.2824942000347)|1
 test 2.3|2|43|LINESTRING(-75.5529884291587 49.2811479840607,-75.5522356128797 49.2815620330142)|3