Skip to content

Implemeted SQL function center() that calculates the gravity center(centroid) from array of points #33

New issue

Have a question about this project? Sign up for a free GitHub account to open an issue and contact its maintainers and the community.

Sign up for GitHub

By clicking “Sign up for GitHub”, you agree to our terms of service and privacy statement. We’ll occasionally send you account related emails.

Already on GitHub? Sign in to your account

Closed
wants to merge 1 commit into from
Closed

Implemeted SQL function center() that calculates the gravity center(centroid) from array of points #33

Show file tree
Hide file tree
Changes from all commits
Commits
Show all changes
1 commit
Select commit
File filter

Filter by extension

Filter by extension
Conversations
Failed to load comments.
Loading
Jump to
Jump to file
Failed to load files.
Loading
Diff view
Diff view
3 changes: 2 additions & 1 deletion Makefile
View file
Open in desktop
Original file line number Diff line number Diff line change
Expand Up @@ -31,7 +31,8 @@ DATA_built = $(RELEASE_SQL) \

DOCS = README.pg_sphere COPYRIGHT.pg_sphere
REGRESS = init tables points euler circle line ellipse poly path box index \
contains_ops contains_ops_compat bounding_box_gist gnomo epochprop
contains_ops contains_ops_compat bounding_box_gist gnomo epochprop \
spherepoint_array_center

ifneq ($(USE_HEALPIX),0)
REGRESS += healpix moc mocautocast
Expand Down
2 changes: 1 addition & 1 deletion doc/functions.sgm
View file
Open in desktop
Original file line number Diff line number Diff line change
Expand Up @@ -143,7 +143,7 @@
</sect3>

</sect2>

<sect2 id="funcs.strans">
<title>
<type>strans</type> functions
Expand Down
73 changes: 71 additions & 2 deletions doc/operators.sgm
View file
Open in desktop
Original file line number Diff line number Diff line change
Expand Up @@ -393,8 +393,10 @@
The center operator <literal>@@</literal> is a non-boolean
unary operator returning the center of an object. In the
current implementation of <application>pgSphere</application>,
only centers of circles and ellipses are supported. Instead
of using the operator, you can use the function
only centers of circles, ellipses, points, arrays of points are
</para>
<para>
supported. Instead of using the operator, you can use the function
<literal>center(object)</literal>.
</para>
<example>
Expand All @@ -403,6 +405,73 @@
<![CDATA[sql> SELECT @@ scircle '<(0d,20d),30d>';]]>
</programlisting>
</example>
<example>
<title>Center array of points</title>
<programlisting>
<![CDATA[sql> SELECT @@ ARRAY[
spoint(40.7128, -74.0060),
spoint(34.0522, -118.2437),
spoint(37.7749, -122.4194)
] AS center;]]>
<![CDATA[ center ]]>
<![CDATA[----------------------------------------]]>
<![CDATA[ (3.04366980631979 , 0.858938068921891)]]>
<![CDATA[(1 row)]]>
</programlisting>
</example>

<example>
<title>Center of a spoint</title>
<para>
For the completeness of the interface, the function returns the point itself
</para>
<programlisting>
<![CDATA[sql> select @@ spoint(0,0) AS spoint;]]>
<![CDATA[ spoint ]]>
<![CDATA[---------]]>
<![CDATA[ (0 , 0) ]]>
<![CDATA[(1 row)]]>
</programlisting>
</example>

<example>
<title>Basic example array of points</title>
<programlisting>
<![CDATA[sql> SELECT center(ARRAY[
spoint(40.7128, -74.0060),
spoint(34.0522, -118.2437),
spoint(37.7749, -122.4194)
]);]]>
<![CDATA[ center ]]>
<![CDATA[----------------------------------------]]>
<![CDATA[ (3.04366980631979 , 0.858938068921891)]]>
<![CDATA[ (1 row) ]]>
</programlisting>
</example>

<example>
<title>Behavior with empty array of points</title>
<programlisting>
<![CDATA[sql> SELECT center('{}');]]>
<![CDATA[ center ]]>
<![CDATA[ ----------------------------- ]]>

<![CDATA[ (1 row) ]]>
</programlisting>
</example>
<example>
<title>Example with opposite array points(poles)</title>
<programlisting>
<![CDATA[sql> SELECT center(ARRAY[
spoint(0, 10),
spoint(0, -10),
]);]]>
<![CDATA[ center ]]>
<![CDATA[----------------------------------------]]>
<![CDATA[ (0, 0) ]]>
<![CDATA[ (1 row) ]]>
</programlisting>
</example>
</sect2>

<sect2 id="op.direction">
Expand Down
14 changes: 7 additions & 7 deletions expected/init_test.out.in
View file
Open in desktop
Original file line number Diff line number Diff line change
Expand Up @@ -24,12 +24,12 @@ psql:pg_sphere.test.sql:159: NOTICE: argument type spath is only a shell
psql:pg_sphere.test.sql:178: NOTICE: type "sbox" is not yet defined
DETAIL: Creating a shell type definition.
psql:pg_sphere.test.sql:185: NOTICE: argument type sbox is only a shell
psql:pg_sphere.test.sql:8540: NOTICE: type "spherekey" is not yet defined
psql:pg_sphere.test.sql:8570: NOTICE: type "spherekey" is not yet defined
DETAIL: Creating a shell type definition.
psql:pg_sphere.test.sql:8547: NOTICE: argument type spherekey is only a shell
psql:pg_sphere.test.sql:8561: NOTICE: type "pointkey" is not yet defined
psql:pg_sphere.test.sql:8577: NOTICE: argument type spherekey is only a shell
psql:pg_sphere.test.sql:8591: NOTICE: type "pointkey" is not yet defined
DETAIL: Creating a shell type definition.
psql:pg_sphere.test.sql:8568: NOTICE: argument type pointkey is only a shell
psql:pg_sphere.test.sql:8574: NOTICE: argument type pointkey is only a shell
psql:pg_sphere.test.sql:8580: NOTICE: argument type pointkey is only a shell
psql:pg_sphere.test.sql:8586: NOTICE: argument type pointkey is only a shell
psql:pg_sphere.test.sql:8598: NOTICE: argument type pointkey is only a shell
psql:pg_sphere.test.sql:8604: NOTICE: argument type pointkey is only a shell
psql:pg_sphere.test.sql:8610: NOTICE: argument type pointkey is only a shell
psql:pg_sphere.test.sql:8616: NOTICE: argument type pointkey is only a shell
4 changes: 2 additions & 2 deletions expected/init_test_healpix.out.in
View file
Open in desktop
Original file line number Diff line number Diff line change
@@ -1,2 +1,2 @@
psql:pg_sphere.test.sql:9153: NOTICE: return type smoc is only a shell
psql:pg_sphere.test.sql:9159: NOTICE: argument type smoc is only a shell
psql:pg_sphere.test.sql:9183: NOTICE: return type smoc is only a shell
psql:pg_sphere.test.sql:9189: NOTICE: argument type smoc is only a shell
7 changes: 7 additions & 0 deletions expected/points.out
View file
Open in desktop
Original file line number Diff line number Diff line change
Expand Up @@ -648,3 +648,10 @@ SELECT '( 0h 2m 30s , -90d 0m 0s)'::spoint<->'( 12h 2m 30s , -90d 0m 0s)'::spoin
0
(1 row)

-- Center operator -------------------
select @@ spoint(0,0) AS spoint;
spoint
---------
(0 , 0)
(1 row)

73 changes: 73 additions & 0 deletions expected/spherepoint_array_center.out
View file
Open in desktop
Original file line number Diff line number Diff line change
@@ -0,0 +1,73 @@
SELECT center(ARRAY[
spoint(40.7128, -74.0060),
spoint(34.0522, -118.2437),
spoint(37.7749, -122.4194)
]);
center
----------------------------------------
(3.04366980631979 , 0.858938068921891)
(1 row)

SELECT center('{}'::SPoint[]);
center
--------

(1 row)

CREATE FUNCTION spoint_from_xyz(FLOAT8, FLOAT8, FLOAT8)
RETURNS spoint
AS 'pg_sphere', 'spoint_from_xyz'
LANGUAGE 'c'
IMMUTABLE STRICT PARALLEL SAFE;
SELECT spoint_from_xyz(1, 0, 0);
spoint_from_xyz
-----------------
(0 , 0)
(1 row)

SELECT spoint_from_xyz(0, 0, 0);
spoint_from_xyz
-----------------
(0 , 0)
(1 row)

SELECT center(ARRAY[
spoint_from_xyz(1, 0, 0),
spoint_from_xyz(-1, 0, 0)
]);
center
-----------------------
(1.5707963267949 , 0)
(1 row)

SELECT center(NULL::SPoint[]);
center
--------

(1 row)

SELECT @@ ARRAY[
spoint(40.7128, -74.0060),
spoint(34.0522, -118.2437),
spoint(37.7749, -122.4194)
] AS center;
center
----------------------------------------
(3.04366980631979 , 0.858938068921891)
(1 row)

SELECT @@ ARRAY[]::spoint[] AS center;
center
--------

(1 row)

SELECT @@ ARRAY[
spoint_from_xyz(1, 0, 0),
spoint_from_xyz(-1, 0, 0)
] AS center;
center
-----------------------
(1.5707963267949 , 0)
(1 row)

30 changes: 30 additions & 0 deletions pgs_point.sql.in
View file
Open in desktop
Original file line number Diff line number Diff line change
Expand Up @@ -167,3 +167,33 @@ CREATE OPERATOR <-> (
COMMENT ON OPERATOR <-> (spoint, spoint) IS
'distance between spherical points';

CREATE FUNCTION center(dots_vector SPoint[])
Copy link
Contributor

Choose a reason for hiding this comment

The reason will be displayed to describe this comment to others. Learn more.

I thought we had settled on naming this function centroid? Why was center chosen?

I see it referred to as "centroid" at these URLs:

https://stackoverflow.com/a/38201499/5153779
https://gis.stackexchange.com/questions/43505/calculating-a-spherical-polygon-centroid

However, "center" is used here:
https://github.com/chrisveness/geodesy/blob/8f4ef33d/latlon-nvector-spherical.js#L783

So I'm not sure which term is correct.

Does your implementation give the same results as the Python and JavaScript implementations at https://stackoverflow.com/questions/19897187#answer-38201499 and https://github.com/chrisveness/geodesy/blob/8f4ef33d/latlon-nvector-spherical.js#L783 ? If not, why not?

Copy link
Author

Choose a reason for hiding this comment

The reason will be displayed to describe this comment to others. Learn more.

@vitcpp and I discussed the name "centroid" earlier, but in the end we came to the conclusion that the best name is the center for the sql function to match the interface with the operator @@

Copy link
Contributor

Choose a reason for hiding this comment

The reason will be displayed to describe this comment to others. Learn more.

@esabol The main reason to choose 'center' function because such function already exists. It returns the center of circle or ellipse and it is intended to be expanded to other object types. In the doc it is described as center of an object. It was decided to overload this function and allow to pass array of points and the point itself (for interface completeness). It is why we decided to put on display PR with 'center' function. Furthermore, we may utilize operator @@. My opinion is that center and centroid are interchangeable in most cases. So, our idea is to reuse center and to avoid introducing new terms.

Copy link
Contributor

Choose a reason for hiding this comment

The reason will be displayed to describe this comment to others. Learn more.

OK, that makes sense (the center name) then, but I'd still like an answer to this:

Does your implementation give the same results as the Python and JavaScript implementations at https://stackoverflow.com/questions/19897187#answer-38201499 and https://github.com/chrisveness/geodesy/blob/8f4ef33d/latlon-nvector-spherical.js#L783 ? If not, why not?

Copy link
Author

Choose a reason for hiding this comment

The reason will be displayed to describe this comment to others. Learn more.

@esabol But this is function, which

Calculates the centre of a spherical polygon where the sides of the polygon are great circle arcs joining the vertices.
In my task, I have an arbitrary set of points. That is why I have now tested and on any two input points their function gives 0, because the polygon starts with three points

Copy link
Author

Choose a reason for hiding this comment

The reason will be displayed to describe this comment to others. Learn more.

@esabol Sorry... Give me time, then I'll try to figure out why the answers are completely different

Copy link
Contributor

@vitcpp vitcpp Jul 28, 2023 •
edited
Loading

Choose a reason for hiding this comment

The reason will be displayed to describe this comment to others. Learn more.

Dear Colleagues, it was intended that the proposed center function calculates the mass center of a set of points (particles) that do not form a figure. Such function may be useful in some astronomical tasks when we have to deal with stellar clusters or with other groups of astronomical objects which mass center on the sphere we'd like to find. When we want to find the center of a figure we definitely should use the Stokes theorem as it was implemented in the referenced JS library (center(spoly) should be implemented that way). I guess it is why @stepan-neretin7 sees some differences when trying to compare the proposed function against the function to calculate polygon center from the referenced JS library. Adding some more points on a polygon edge doesn't change the shape of the polygon. But adding more points change the mass center of the set of particles (points). There is another question - is it reasonable to implement such function? I think yes. Probably, the name center is best suited to calculate the center of a figure than to calculate the mass center of a set of points (may be come back to centroid - ?). I'm not sure about having such similar function in the JS library. If yes, it is reasonable to compare our function against its counterpart.

Copy link
Contributor

Choose a reason for hiding this comment

The reason will be displayed to describe this comment to others. Learn more.

OK with me.

Copy link
Author

Choose a reason for hiding this comment

The reason will be displayed to describe this comment to others. Learn more.

@vitcpp @esabol
Dear colleagues, I want to revive the discussion of my pull-request. I have two questions that I would like to discuss? First; After all, we'd better leave such an implementation (the center of an arbitrary set of points, not a polygon.. If we add more points on the side of an imaginary polygon from points, then my implementation will shift to this side) or we will think in the direction of that implementation for polygons through the Stokes theorem
Second: When two points are at the poles: point(0,-pi()/2), point(0,pi()/2), then the middle point will be ~epsilon and it is incorrect to translate this into lat, lng, because such a point is not projected onto the sphere. What should I give out in this case? NULL? But let's say the case with the points 0,0 and 0, pi
The center here is any point on a large circle perpendicular to the line connecting the points
What is better to return here? NAN? NULL?

Copy link
Contributor

Choose a reason for hiding this comment

The reason will be displayed to describe this comment to others. Learn more.

Deal All, the final result of this operation is to find centroid of a set of points, not the center of a polygon. We do not find the center of a figure like polygon in this task. The averaged vector is pretty fine solution. There are two special cases. Case 1: we find the center of two opposite points. In this case, the equidistant place is a circle, not a point. The geometric
center in the Cartesian coordinate system is the center of the sphere. Case 2: when we have, for example, two points on poles and 4 points on equator with longitudes 0, 90, 180, 270 degrees. In this case the geometric center in 3D is (0,0,0) but there is no geometric center on the sphere. Case 3: we have all points on the same big circle. In this case, there are two points are equidistant. I'm not sure about other special cases. I guess, there are no other special cases.

This function may be useful to find the approximate center of a group of astronomical objects that are placed in a local area on the sphere: the center of a star cluster, the center of a selected group of objects, etc. I think the defined special cases are impractical; I'm not sure in what tasks it can happen.

I propose this function to return NULL if 3D geometric center is (0,0,0) +- EPSILON. It means, the algorithm can't find centroid of the point set. I think, NULL is better because one can use ifull or coalesce to define some result in such cases.

RETURNS spoint
AS 'MODULE_PATHNAME', 'center'
LANGUAGE 'c'
IMMUTABLE STRICT PARALLEL SAFE;

CREATE OPERATOR @@ (
RIGHTARG = SPoint[],
PROCEDURE = center
);

COMMENT ON OPERATOR @@ (NONE , SPoint[]) IS
'Return gravity center from array of spherical points';

CREATE FUNCTION center(spoint)
RETURNS spoint
AS 'MODULE_PATHNAME', 'spherepoint_center'
LANGUAGE 'c'
IMMUTABLE STRICT PARALLEL SAFE;

COMMENT ON FUNCTION center(spoint) IS
'For the correctness of the interface, the function returns the spoint itself';

CREATE OPERATOR @@ (
RIGHTARG = SPoint,
PROCEDURE = center
);

COMMENT ON OPERATOR @@ (NONE , SPoint) IS
'For the correctness of the interface, the function returns the spoint itself';
2 changes: 2 additions & 0 deletions sql/points.sql
View file
Open in desktop
Original file line number Diff line number Diff line change
Expand Up @@ -234,3 +234,5 @@ SELECT '( 0h 2m 30s , 90d 0m 0s)'::spoint<->'( 12h 2m 30s , 90d 0m 0s)'::spoint;

SELECT '( 0h 2m 30s , -90d 0m 0s)'::spoint<->'( 12h 2m 30s , -90d 0m 0s)'::spoint;

-- Center operator -------------------
select @@ spoint(0,0) AS spoint;
37 changes: 37 additions & 0 deletions sql/spherepoint_array_center.sql
View file
Open in desktop
Original file line number Diff line number Diff line change
@@ -0,0 +1,37 @@
SELECT center(ARRAY[
spoint(40.7128, -74.0060),
spoint(34.0522, -118.2437),
spoint(37.7749, -122.4194)
]);

SELECT center('{}'::SPoint[]);

CREATE FUNCTION spoint_from_xyz(FLOAT8, FLOAT8, FLOAT8)
RETURNS spoint
AS 'pg_sphere', 'spoint_from_xyz'
LANGUAGE 'c'
IMMUTABLE STRICT PARALLEL SAFE;

SELECT spoint_from_xyz(1, 0, 0);

SELECT spoint_from_xyz(0, 0, 0);

SELECT center(ARRAY[
spoint_from_xyz(1, 0, 0),
spoint_from_xyz(-1, 0, 0)
]);

SELECT center(NULL::SPoint[]);

SELECT @@ ARRAY[
spoint(40.7128, -74.0060),
spoint(34.0522, -118.2437),
spoint(37.7749, -122.4194)
] AS center;

SELECT @@ ARRAY[]::spoint[] AS center;

SELECT @@ ARRAY[
spoint_from_xyz(1, 0, 0),
spoint_from_xyz(-1, 0, 0)
] AS center;
67 changes: 67 additions & 0 deletions src/point.c
View file
Open in desktop
Original file line number Diff line number Diff line change
Expand Up @@ -5,13 +5,17 @@
PG_FUNCTION_INFO_V1(spherepoint_in);
PG_FUNCTION_INFO_V1(spherepoint_from_long_lat);
PG_FUNCTION_INFO_V1(spherepoint_distance);
PG_FUNCTION_INFO_V1(spherepoint_center);
PG_FUNCTION_INFO_V1(spherepoint_long);
PG_FUNCTION_INFO_V1(spherepoint_lat);
PG_FUNCTION_INFO_V1(spherepoint_x);
PG_FUNCTION_INFO_V1(spherepoint_y);
PG_FUNCTION_INFO_V1(spherepoint_z);
PG_FUNCTION_INFO_V1(spherepoint_xyz);
PG_FUNCTION_INFO_V1(spherepoint_equal);
PG_FUNCTION_INFO_V1(spoint_from_xyz);
PG_FUNCTION_INFO_V1(center);


bool
spoint_eq(const SPoint *p1, const SPoint *p2)
Expand Down Expand Up @@ -179,6 +183,13 @@ spoint_dist(const SPoint *p1, const SPoint *p2)
}
}

Datum
spherepoint_center(PG_FUNCTION_ARGS)
{
SPoint *p = (SPoint *) PG_GETARG_POINTER(0);
PG_RETURN_POINTER(p);
}

Datum
spherepoint_distance(PG_FUNCTION_ARGS)
{
Expand Down Expand Up @@ -263,3 +274,59 @@ spherepoint_equal(PG_FUNCTION_ARGS)

PG_RETURN_BOOL(spoint_eq(p1, p2));
}

Datum spoint_from_xyz(PG_FUNCTION_ARGS)
{
Vector3D point_coords;
SPoint *p = (SPoint *) palloc(sizeof(SPoint));

point_coords.x = PG_GETARG_FLOAT8(0);
point_coords.y = PG_GETARG_FLOAT8(1);
point_coords.z = PG_GETARG_FLOAT8(2);
vector3d_spoint(p, &point_coords);

if (p == NULL)
{
PG_RETURN_NULL();
}

spoint_check(p);
PG_RETURN_POINTER(p);
}

Datum center(PG_FUNCTION_ARGS)
{
int num_elements, i;
SPoint * p = (SPoint *) palloc(sizeof(SPoint));
SPoint * array_data;
Vector3D v;
Vector3D point_coords = {0,0,0};
ArrayType *dots_vector;

dots_vector = PG_GETARG_ARRAYTYPE_P(0);
num_elements = ArrayGetNItems(ARR_NDIM(dots_vector), ARR_DIMS(dots_vector));
if (num_elements <= 0)
{
PG_RETURN_NULL();
}

array_data = (SPoint *) ARR_DATA_PTR(dots_vector);

for (i = 0; i < num_elements; i++)
{
spoint_vector3d(&v, &array_data[i]);
point_coords.x += v.x;
point_coords.y += v.y;
point_coords.z += v.z;
}

point_coords.x /= num_elements;
point_coords.y /= num_elements;
point_coords.z /= num_elements;

vector3d_spoint(p, &point_coords);

spoint_check(p);

PG_RETURN_POINTER(p);
}
Loading