[geometry] Spherical Barycenters

classic Classic list List threaded Threaded
4 messages Options
Reply | Threaded
Open this post in threaded view
|

[geometry] Spherical Barycenters

Matt Juntunen
Hi all,

I'm attempting to wrap up the refactoring work I've been doing with GEOMETRY-32 but I've hit a snag with the code for computing barycenters in 2D spherical space. The results I'm getting seem mostly reasonable but they aren't as consistent as I'd like, which makes me wonder if I'm doing something wrong or if I've made an invalid assumption somewhere. The previous barycenter code was not really unit tested so I don't have much to compare to.

Here is what I have currently for computing the barycenter of convex regions with 3 or more sides:
- Split the region into a series of triangles using a triangle fan approach (eg, if the boundary points are p0, p1, p2, p3, split the region into the triangles p0, p1, p2 and p0, p2, p3)
- Compute the area of each triangle by calculating the interior angles A, B, C and using Girard's Theorem (area = A + B + C - pi)
- Compute the barycenter of each triangle using the formula for the center of mass found on the first page of the article here [1]. This is done by multiplying the pole of each side's great circle by its arc length and summing the result vectors. The sum vector is normalized to place it on the surface of the unit sphere. (This is the approach used in the previous code for the barycenter [2])
- Compute the overall barycenter for the region by scaling each triangle barycenter vector by the triangle area, adding the results, and normalizing.

As mentioned above, this seems to produce mostly reasonable results. The issue I'm having is with verifying their consistency. In my unit tests, I'm doing the following to check the barycenters:
- Compute a number of great circles that pass through the region barycenter
- Use each circle to split the area into two smaller convex regions - a plus and a minus region
- Compute the barycenters and areas for the plus and minus regions using the same logic as the larger region
- Combine the two barycenters by scaling their vectors by the region's area, summing, and normalizing
- Verify that the barycenter computed from the split parts matches the one computed for the region as a whole

This basically checks that the region is the sum of its parts in terms of barycenter. However, I'm finding that the barycenter computed from the split regions is always slightly off from the overall region barycenter, generally by a margin of between 0.001 to 0.01 radians in both azimuth and polar. So, I'm wondering if this is a math issue, a floating point issue, or something else.

Any insight would be greatly appreciated. My code can be found here [3] and here [4], although I haven't yet pushed the changes for the full algorithm I outlined above.

Regards,
Matt J


[1] https://asmedigitalcollection.asme.org/appliedmechanics/article/42/1/239/387501/The-Inertia-Tensor-for-a-Spherical-Triangle
[2] https://github.com/darkma773r/commons-geometry/blob/master/commons-geometry-spherical/src/main/java/org/apache/commons/geometry/spherical/twod/PropertiesComputer.java#L135
[3] https://github.com/darkma773r/commons-geometry/blob/geometry-32-working/commons-geometry-spherical/src/main/java/org/apache/commons/geometry/spherical/twod/ConvexArea2S.java
[4] https://github.com/darkma773r/commons-geometry/blob/geometry-32-working/commons-geometry-spherical/src/test/java/org/apache/commons/geometry/spherical/twod/ConvexArea2STest.java#L188
Reply | Threaded
Open this post in threaded view
|

Re: [geometry] Spherical Barycenters

Gilles Sadowski-2
Hello.

Le ven. 15 nov. 2019 à 15:34, Matt Juntunen
<[hidden email]> a écrit :
>
> Hi all,
>
> I'm attempting to wrap up the refactoring work I've been doing with GEOMETRY-32 but I've hit a snag with the code for computing barycenters in 2D spherical space. The results I'm getting seem mostly reasonable but they aren't as consistent as I'd like, which makes me wonder if I'm doing something wrong or if I've made an invalid assumption somewhere. The previous barycenter code was not really unit tested so I don't have much to compare to.

What does the CM code produce?

Gilles

>
> Here is what I have currently for computing the barycenter of convex regions with 3 or more sides:
> - Split the region into a series of triangles using a triangle fan approach (eg, if the boundary points are p0, p1, p2, p3, split the region into the triangles p0, p1, p2 and p0, p2, p3)
> - Compute the area of each triangle by calculating the interior angles A, B, C and using Girard's Theorem (area = A + B + C - pi)
> - Compute the barycenter of each triangle using the formula for the center of mass found on the first page of the article here [1]. This is done by multiplying the pole of each side's great circle by its arc length and summing the result vectors. The sum vector is normalized to place it on the surface of the unit sphere. (This is the approach used in the previous code for the barycenter [2])
> - Compute the overall barycenter for the region by scaling each triangle barycenter vector by the triangle area, adding the results, and normalizing.
>
> As mentioned above, this seems to produce mostly reasonable results. The issue I'm having is with verifying their consistency. In my unit tests, I'm doing the following to check the barycenters:
> - Compute a number of great circles that pass through the region barycenter
> - Use each circle to split the area into two smaller convex regions - a plus and a minus region
> - Compute the barycenters and areas for the plus and minus regions using the same logic as the larger region
> - Combine the two barycenters by scaling their vectors by the region's area, summing, and normalizing
> - Verify that the barycenter computed from the split parts matches the one computed for the region as a whole
>
> This basically checks that the region is the sum of its parts in terms of barycenter. However, I'm finding that the barycenter computed from the split regions is always slightly off from the overall region barycenter, generally by a margin of between 0.001 to 0.01 radians in both azimuth and polar. So, I'm wondering if this is a math issue, a floating point issue, or something else.
>
> Any insight would be greatly appreciated. My code can be found here [3] and here [4], although I haven't yet pushed the changes for the full algorithm I outlined above.
>
> Regards,
> Matt J
>
>
> [1] https://asmedigitalcollection.asme.org/appliedmechanics/article/42/1/239/387501/The-Inertia-Tensor-for-a-Spherical-Triangle
> [2] https://github.com/darkma773r/commons-geometry/blob/master/commons-geometry-spherical/src/main/java/org/apache/commons/geometry/spherical/twod/PropertiesComputer.java#L135
> [3] https://github.com/darkma773r/commons-geometry/blob/geometry-32-working/commons-geometry-spherical/src/main/java/org/apache/commons/geometry/spherical/twod/ConvexArea2S.java
> [4] https://github.com/darkma773r/commons-geometry/blob/geometry-32-working/commons-geometry-spherical/src/test/java/org/apache/commons/geometry/spherical/twod/ConvexArea2STest.java#L188

---------------------------------------------------------------------
To unsubscribe, e-mail: [hidden email]
For additional commands, e-mail: [hidden email]

Reply | Threaded
Open this post in threaded view
|

Re: [geometry] Spherical Barycenters

Matt Juntunen
Gilles,

The old commons-math code produces the same results. So, I'm not sure if this is a bug or a known limitation. There isn't any documentation in the old code in regard to the barycenter computation and there are only two trivial unit tests for it. I suppose if there any no further insights into this, I could just stick with the existing algorithm, use a large epsilon in the unit tests, and leave any fixes for future JIRA issues.

-Matt
________________________________
From: Gilles Sadowski <[hidden email]>
Sent: Friday, November 15, 2019 11:15 AM
To: Commons Developers List <[hidden email]>
Subject: Re: [geometry] Spherical Barycenters

Hello.

Le ven. 15 nov. 2019 à 15:34, Matt Juntunen
<[hidden email]> a écrit :
>
> Hi all,
>
> I'm attempting to wrap up the refactoring work I've been doing with GEOMETRY-32 but I've hit a snag with the code for computing barycenters in 2D spherical space. The results I'm getting seem mostly reasonable but they aren't as consistent as I'd like, which makes me wonder if I'm doing something wrong or if I've made an invalid assumption somewhere. The previous barycenter code was not really unit tested so I don't have much to compare to.

What does the CM code produce?

Gilles

>
> Here is what I have currently for computing the barycenter of convex regions with 3 or more sides:
> - Split the region into a series of triangles using a triangle fan approach (eg, if the boundary points are p0, p1, p2, p3, split the region into the triangles p0, p1, p2 and p0, p2, p3)
> - Compute the area of each triangle by calculating the interior angles A, B, C and using Girard's Theorem (area = A + B + C - pi)
> - Compute the barycenter of each triangle using the formula for the center of mass found on the first page of the article here [1]. This is done by multiplying the pole of each side's great circle by its arc length and summing the result vectors. The sum vector is normalized to place it on the surface of the unit sphere. (This is the approach used in the previous code for the barycenter [2])
> - Compute the overall barycenter for the region by scaling each triangle barycenter vector by the triangle area, adding the results, and normalizing.
>
> As mentioned above, this seems to produce mostly reasonable results. The issue I'm having is with verifying their consistency. In my unit tests, I'm doing the following to check the barycenters:
> - Compute a number of great circles that pass through the region barycenter
> - Use each circle to split the area into two smaller convex regions - a plus and a minus region
> - Compute the barycenters and areas for the plus and minus regions using the same logic as the larger region
> - Combine the two barycenters by scaling their vectors by the region's area, summing, and normalizing
> - Verify that the barycenter computed from the split parts matches the one computed for the region as a whole
>
> This basically checks that the region is the sum of its parts in terms of barycenter. However, I'm finding that the barycenter computed from the split regions is always slightly off from the overall region barycenter, generally by a margin of between 0.001 to 0.01 radians in both azimuth and polar. So, I'm wondering if this is a math issue, a floating point issue, or something else.
>
> Any insight would be greatly appreciated. My code can be found here [3] and here [4], although I haven't yet pushed the changes for the full algorithm I outlined above.
>
> Regards,
> Matt J
>
>
> [1] https://asmedigitalcollection.asme.org/appliedmechanics/article/42/1/239/387501/The-Inertia-Tensor-for-a-Spherical-Triangle
> [2] https://github.com/darkma773r/commons-geometry/blob/master/commons-geometry-spherical/src/main/java/org/apache/commons/geometry/spherical/twod/PropertiesComputer.java#L135
> [3] https://github.com/darkma773r/commons-geometry/blob/geometry-32-working/commons-geometry-spherical/src/main/java/org/apache/commons/geometry/spherical/twod/ConvexArea2S.java
> [4] https://github.com/darkma773r/commons-geometry/blob/geometry-32-working/commons-geometry-spherical/src/test/java/org/apache/commons/geometry/spherical/twod/ConvexArea2STest.java#L188

---------------------------------------------------------------------
To unsubscribe, e-mail: [hidden email]
For additional commands, e-mail: [hidden email]

Reply | Threaded
Open this post in threaded view
|

Re: [geometry] Spherical Barycenters

Gilles Sadowski-2
Hello.

Le sam. 16 nov. 2019 à 05:40, Matt Juntunen
<[hidden email]> a écrit :
>
> Gilles,
>
> The old commons-math code produces the same results.

:-(

> So, I'm not sure if this is a bug or a known limitation. There isn't any documentation in the old code in regard to the barycenter computation and there are only two trivial unit tests for it.

I had a look, but it's not even a clear to me what trivial cases it is testing.
It would help (me, at least) if there were separate methods that indeed
spell out some trivial cases (where the expected answer should hopefully
be obtained with ultimate accuracy, whatever the implementation).

> I suppose if there any no further insights into this, I could just stick with the existing algorithm, use a large epsilon in the unit tests, and leave any fixes for future JIRA issues.

Yes, but please:
* add code comment in the unit tests (as to the unexpected large tolerance
required for them to pass),
* add Javadoc in the "main" code (signalling that the unit tests point to a
potential bug),
* file a JIRA report for this issue (with a link to the original CM code, in
the "MATH_3_X" branch)

When your work is merged in "master", we should not forget to add a link
to the code in the [Geometry] repository.


Thanks a lot for reporting the problem,
Gilles

>> [...]

---------------------------------------------------------------------
To unsubscribe, e-mail: [hidden email]
For additional commands, e-mail: [hidden email]