Merge 7d298d7 into b58339c

Author: mngr777

NEWS

@@ -77,6 +77,7 @@ PostGIS 3.2.0
   - #4986, GIST indexes on Postgres 14 are now created faster using Hilbert-sorting method.
            (Han Wang, Aliaksandr Kalenik, Darafei Praliaskouski, Giuseppe Broccolo)
   - #4949, Use proj_normalize_for_visualization to hand "axis swap" decisions (Paul Ramsey)
+  - ST_PixelAsCentroids, ST_PixelAsCentroid reimplemented on top of a C function (Sergei Shoulbakov)
 
  * New features*
   - #4923, topology.ValidateTopologyRelation (Sandro Santilli)

raster/rt_core/librtcore.h

@@ -1549,6 +1549,17 @@ rt_raster_compute_skewed_raster(
  */
 LWPOLY* rt_raster_pixel_as_polygon(rt_raster raster, int x, int y);
 
+/**
+ * Get a raster pixel centroid point.
+ *
+ * @param raster : the raster to get pixel from
+ * @param x : the column number
+ * @param y : the row number
+ *
+ * @return the pixel centroid point, or NULL on error.
+ */
+LWPOINT* rt_raster_pixel_as_centroid_point(rt_raster rast, int x, int y);
+
 /**
  * Get a raster as a surface (multipolygon).  If a band is specified,
  * those pixels with value (not NODATA) contribute to the area

raster/rt_core/rt_geometry.c

@@ -654,6 +654,44 @@ rt_raster_pixel_as_polygon(rt_raster rast, int x, int y)
     return poly;
 }
 
+/******************************************************************************
+* rt_raster_pixel_as_centroid_point()
+******************************************************************************/
+
+/**
+ * Get a raster pixel centroid point.
+ *
+ * @param raster : the raster to get pixel from
+ * @param x : the column number
+ * @param y : the row number
+ *
+ * @return the pixel centroid point, or NULL on error.
+ */
+LWPOINT*
+rt_raster_pixel_as_centroid_point(rt_raster rast, int x, int y)
+{
+    double scale_x, scale_y;
+    double skew_x, skew_y;
+    double ul_x, ul_y;
+    int32_t srid;
+    double center_x, center_y;
+    LWPOINT* point;
+
+    scale_x = rt_raster_get_x_scale(rast);
+    scale_y = rt_raster_get_y_scale(rast);
+    skew_x = rt_raster_get_x_skew(rast);
+    skew_y = rt_raster_get_y_skew(rast);
+    ul_x = rt_raster_get_x_offset(rast);
+    ul_y = rt_raster_get_y_offset(rast);
+    srid = rt_raster_get_srid(rast);
+
+    center_x = scale_x * x + skew_x * y + ul_x + (scale_x + skew_x) * 0.5;
+    center_y = scale_y * y + skew_y * x + ul_y + (scale_y + skew_y) * 0.5;
+    point = lwpoint_make2d(srid, center_x, center_y);
+
+    return point;
+}
+
 /******************************************************************************
 * rt_raster_get_envelope_geom()
 ******************************************************************************/