Polygon2d.centroid

AUTHORS

@@ -10,4 +10,4 @@ Jakub Hampl <[email protected]> (convex hull, point in polygon)
 Dave Cameron <[email protected]> (centroids)
 Lalit Umbarkar <[email protected]> (bounding box expand/offset)
 Sebastian Kazenbroot-Guppy <[email protected]> (circle-bbox intersections)
-Andrey Kuzmin <[email protected]> (arc midpoint)
+Andrey Kuzmin <[email protected]> (arc midpoint, polygon centroid)

src/Polygon2d.elm

@@ -10,7 +10,7 @@
 module Polygon2d exposing
     ( Polygon2d
     , singleLoop, withHoles, convexHull
-    , outerLoop, innerLoops, vertices, edges, perimeter, area, boundingBox
+    , outerLoop, innerLoops, vertices, edges, perimeter, area, centroid, boundingBox
     , contains
     , scaleAbout, rotateAround, translateBy, translateIn, mirrorAcross
     , at, at_
@@ -37,7 +37,7 @@ holes. This module contains a variety of polygon-related functionality, such as
 
 # Properties
 
-@docs outerLoop, innerLoops, vertices, edges, perimeter, area, boundingBox
+@docs outerLoop, innerLoops, vertices, edges, perimeter, area, centroid, boundingBox
 
 
 # Queries
@@ -377,6 +377,199 @@ area polygon =
             (Quantity.sum (List.map counterclockwiseArea (innerLoops polygon)))
 
 
+{-| Get the centroid of a polygon. Returns `Nothing`
+if the polygon has no vertices or empty area.
+-}
+centroid : Polygon2d units coordinates -> Maybe (Point2d units coordinates)
+centroid polygon =
+    case outerLoop polygon of
+        first :: _ :: _ ->
+            let
+                offset =
+                    Point2d.unwrap first
+            in
+            centroidHelp
+                offset.x
+                offset.y
+                first
+                (outerLoop polygon)
+                (innerLoops polygon)
+                0
+                0
+                0
+
+        _ ->
+            Nothing
+
+
+centroidHelp :
+    Float
+    -> Float
+    -> Point2d units coordinates
+    -> List (Point2d units coordinates)
+    -> List (List (Point2d units coordinates))
+    -> Float
+    -> Float
+    -> Float
+    -> Maybe (Point2d units coordinates)
+centroidHelp x0 y0 firstPoint currentLoop remainingLoops xSum ySum areaSum =
+    case currentLoop of
+        [] ->
+            case remainingLoops of
+                loop :: newRemainingLoops ->
+                    case loop of
+                        first :: _ :: _ ->
+                            -- enqueue a new loop
+                            centroidHelp
+                                x0
+                                y0
+                                first
+                                loop
+                                newRemainingLoops
+                                xSum
+                                ySum
+                                areaSum
+
+                        _ ->
+                            -- skip a loop with < 2 points
+                            centroidHelp
+                                x0
+                                y0
+                                firstPoint
+                                []
+                                newRemainingLoops
+                                xSum
+                                ySum
+                                areaSum
+
+                [] ->
+                    if areaSum > 0 then
+                        Just
+                            (Point2d.unsafe
+                                { x = xSum / (areaSum * 3) + x0
+                                , y = ySum / (areaSum * 3) + y0
+                                }
+                            )
+
+                    else
+                        Nothing
+
+        point1 :: currentLoopRest ->
+            case currentLoopRest of
+                point2 :: _ ->
+                    let
+                        p1 =
+                            Point2d.unwrap point1
+
+                        p2 =
+                            Point2d.unwrap point2
+
+                        p1x =
+                            p1.x - x0
+
+                        p1y =
+                            p1.y - y0
+
+                        p2x =
+                            p2.x - x0
+
+                        p2y =
+                            p2.y - y0
+
+                        a =
+                            p1x * p2y - p2x * p1y
+
+                        newXSum =
+                            xSum + (p1x + p2x) * a
+
+                        newYSum =
+                            ySum + (p1y + p2y) * a
+
+                        newAreaSum =
+                            areaSum + a
+                    in
+                    centroidHelp
+                        x0
+                        y0
+                        firstPoint
+                        currentLoopRest
+                        remainingLoops
+                        newXSum
+                        newYSum
+                        newAreaSum
+
+                [] ->
+                    let
+                        p1 =
+                            Point2d.unwrap point1
+
+                        p2 =
+                            Point2d.unwrap firstPoint
+
+                        p1x =
+                            p1.x - x0
+
+                        p1y =
+                            p1.y - y0
+
+                        p2x =
+                            p2.x - x0
+
+                        p2y =
+                            p2.y - y0
+
+                        a =
+                            p1x * p2y - p2x * p1y
+
+                        newXSum =
+                            xSum + (p1x + p2x) * a
+
+                        newYSum =
+                            ySum + (p1y + p2y) * a
+
+                        newAreaSum =
+                            areaSum + a
+                    in
+                    case remainingLoops of
+                        loop :: newRemainingLoops ->
+                            case loop of
+                                first :: _ :: _ ->
+                                    -- enqueue a new loop
+                                    centroidHelp
+                                        x0
+                                        y0
+                                        first
+                                        loop
+                                        newRemainingLoops
+                                        newXSum
+                                        newYSum
+                                        newAreaSum
+
+                                _ ->
+                                    -- skip a loop with < 2 points
+                                    centroidHelp
+                                        x0
+                                        y0
+                                        firstPoint
+                                        []
+                                        newRemainingLoops
+                                        newXSum
+                                        newYSum
+                                        newAreaSum
+
+                        [] ->
+                            if newAreaSum > 0 then
+                                Just
+                                    (Point2d.unsafe
+                                        { x = newXSum / (newAreaSum * 3) + x0
+                                        , y = newYSum / (newAreaSum * 3) + y0
+                                        }
+                                    )
+
+                            else
+                                Nothing
+
+
 {-| Scale a polygon about a given center point by a given scale. If the given
 scale is negative, the order of the polygon's vertices will be reversed so that
 the resulting polygon still has its outer vertices in counterclockwise order and

tests/Tests/Polygon2d.elm

@@ -2,10 +2,14 @@ module Tests.Polygon2d exposing
     ( containsTest
     , convexHullContainsAllPoints
     , convexHullIsConvex
+    , rectangleCentroidIsInTheCenter
+    , rotatingAroundCentroidKeepsCentroid
     , triangulationHasCorrectArea
     , triangulationHasCorrectNumberOfTriangles
+    , triangulationHasCorrectWeightedCentroid
     )
 
+import Area
 import Expect
 import Fuzz
 import Geometry.Expect as Expect
@@ -15,6 +19,7 @@ import LineSegment2d
 import Point2d
 import Polygon2d
 import Quantity exposing (zero)
+import Rectangle2d
 import Test exposing (Test)
 import Triangle2d
 import TriangularMesh
@@ -189,3 +194,80 @@ triangulationHasCorrectNumberOfTriangles =
                 |> List.length
                 |> Expect.equal expectedNumberOfTriangles
         )
+
+
+triangulationHasCorrectWeightedCentroid : Test
+triangulationHasCorrectWeightedCentroid =
+    Test.fuzz Fuzz.polygon2d
+        "The centroid of the polygon before triangulation is the same as weighted centroid of all the resulting triangles"
+        (\polygon ->
+            let
+                triangles =
+                    Polygon2d.triangulate polygon
+                        |> TriangularMesh.faceVertices
+                        |> List.map Triangle2d.fromVertices
+
+                centroidsAsVectors =
+                    triangles
+                        |> List.map (Triangle2d.centroid >> Vector2d.from Point2d.origin)
+
+                areasInSquareMeters =
+                    triangles
+                        |> List.map (Triangle2d.area >> Area.inSquareMeters)
+
+                weightedCentroid =
+                    List.map2 Vector2d.scaleBy areasInSquareMeters centroidsAsVectors
+                        |> List.foldl Vector2d.plus Vector2d.zero
+                        |> Vector2d.scaleBy (1 / List.sum areasInSquareMeters)
+                        |> (\vector -> Point2d.translateBy vector Point2d.origin)
+            in
+            case Polygon2d.centroid polygon of
+                Just centroid ->
+                    Expect.point2d weightedCentroid centroid
+
+                Nothing ->
+                    Expect.fail "Original polygon needs centroid"
+        )
+
+
+rotatingAroundCentroidKeepsCentroid : Test
+rotatingAroundCentroidKeepsCentroid =
+    Test.fuzz2 Fuzz.polygon2d
+        Fuzz.angle
+        "Rotating a polygon around its centroid keeps the centroid point"
+        (\polygon angle ->
+            case Polygon2d.centroid polygon of
+                Just centroid ->
+                    case
+                        polygon
+                            |> Polygon2d.rotateAround centroid angle
+                            |> Polygon2d.centroid
+                    of
+                        Just rotatedCentroid ->
+                            Expect.point2d centroid rotatedCentroid
+
+                        Nothing ->
+                            Expect.fail "Rotated polygon needs centroid"
+
+                Nothing ->
+                    Expect.fail "Original polygon needs centroid"
+        )
+
+
+rectangleCentroidIsInTheCenter : Test
+rectangleCentroidIsInTheCenter =
+    Test.fuzz Fuzz.rectangle2d
+        "The centroid of rectangle is in the center point"
+        (\rectangle ->
+            case
+                rectangle
+                    |> Rectangle2d.vertices
+                    |> Polygon2d.singleLoop
+                    |> Polygon2d.centroid
+            of
+                Just centroid ->
+                    Expect.point2d centroid (Rectangle2d.centerPoint rectangle)
+
+                Nothing ->
+                    Expect.fail "Polygon needs centroid"
+        )