Skip to content

Commit 348e7e1

Browse files
brunobeltranbrunobeltran
committed
add function to get center of mass of path
1 parent 8d1e757 commit 348e7e1

File tree

3 files changed

+253
-2
lines changed

3 files changed

+253
-2
lines changed

‎lib/matplotlib/bezier.py

+76
Original file line numberDiff line numberDiff line change
@@ -297,6 +297,73 @@ def arc_length(self, rtol=None, atol=None):
297297
return curr_est
298298
return curr_est
299299

300+
def arc_center_of_mass(self):
301+
r"""
302+
Center of mass of the (even-odd-rendered) area swept out by the ray
303+
from the origin to the path.
304+
305+
Summing this vector for each segment along a closed path will produce
306+
that area's center of mass.
307+
308+
Returns
309+
-------
310+
r_cm : (2,) np.array<float>
311+
the "arc's center of mass"
312+
313+
Notes
314+
-----
315+
A simple analytical form can be derived for general Bezier curves.
316+
Suppose the curve was closed, so :math:`B(0) = B(1)`. Call the area
317+
enclosed by :math:`B(t)` :math:`B_\text{int}`. The center of mass of
318+
:math:`B_\text{int}` is defined by the expected value of the position
319+
vector `\vec{r}`
320+
321+
.. math::
322+
323+
\vec{R}_\text{cm} = \int_{B_\text{int}} \vec{r} \left( \frac{1}{
324+
\int_{B_\text{int}}} d\vec{r} \right) d\vec{r}
325+
326+
where :math:`(1/\text{Area}(B_\text{int})` can be interpreted as a
327+
probability density.
328+
329+
In order to compute this integral, we choose two functions
330+
:math:`F_0(x,y) = [x^2/2, 0]` and :math:`F_1(x,y) = [0, y^2/2]` such
331+
that :math:`[\div \cdot F_0, \div \cdot F_1] = \vec{r}`. Then, applying
332+
the divergence integral (componentwise), we get that
333+
334+
.. math::
335+
\vec{R}_\text{cm} &= \oint_{B(t)} F \cdot \vec{n} dt \\
336+
&= \int_0^1 \left[ \begin{array}{1}
337+
B^{(0)}(t) \frac{dB^{(1)}(t)}{dt} \\
338+
- B^{(1)}(t) \frac{dB^{(0)}(t)}{dt} \end{array} \right] dt
339+
340+
After expanding in Berstein polynomials and moving the integral inside
341+
all the sums, we get that
342+
343+
.. math::
344+
\vec{R}_\text{cm} = \frac{1}{6} \sum_{i,j=0}^n\sum_{k=0}^{n-1}
345+
\frac{{n \choose i}{n \choose j}{{n-1} \choose k}}
346+
{{3n - 1} \choose {i + j + k}}
347+
\left(\begin{array}{1}
348+
P^{(0)}_i P^{(0)}_j (P^{(1)}_{k+1} - P^{(1)}_k)
349+
- P^{(1)}_i P^{(1)}_j (P^{(0)}_{k+1} - P^{(0)}_k)
350+
\right) \end{array}
351+
352+
where :math:`P_i = [P^{(0)}_i, P^{(1)}_i]` is the :math:`i`'th control
353+
point of the curve and :math:`n` is the degree of the curve.
354+
"""
355+
n = self.degree
356+
r_cm = np.zeros(2)
357+
P = self.control_points
358+
dP = np.diff(P, axis=0)
359+
Pn = np.array([[1, -1]])*dP[:, ::-1] # n = [y, -x]
360+
for i in range(n + 1):
361+
for j in range(n + 1):
362+
for k in range(n):
363+
r_cm += _comb(n, i) * _comb(n, j) * _comb(n - 1, k) \
364+
* P[i]*P[j]*Pn[k] / _comb(3*n - 1, i + j + k)
365+
return r_cm/6
366+
300367
def arc_area(self):
301368
r"""
302369
(Signed) area swept out by ray from origin to curve.
@@ -385,6 +452,15 @@ def arc_area(self):
385452
* (P[j, 0]*dP[k, 1] - P[j, 1]*dP[k, 0])
386453
return (1/4)*area
387454

455+
def center_of_mass(self):
456+
"""Return the center of mass of the curve (not the filled curve!)
457+
458+
Notes
459+
-----
460+
Computed as the mean of the control points.
461+
"""
462+
return np.mean(self._cpoints, axis=0)
463+
388464
@classmethod
389465
def differentiate(cls, B):
390466
"""Return the derivative of a BezierSegment, itself a BezierSegment"""

‎lib/matplotlib/path.py

+139-2
Original file line numberDiff line numberDiff line change
@@ -694,10 +694,147 @@ def signed_area(self, **kwargs):
694694
# add final implied CLOSEPOLY, if necessary
695695
if start_point is not None \
696696
and not np.all(np.isclose(start_point, prev_point)):
697-
B = BezierSegment(np.array([prev_point, start_point]))
698-
area += B.arc_area()
697+
Bclose = BezierSegment(np.array([prev_point, start_point]))
698+
area += Bclose.arc_area()
699699
return area
700700

701+
def center_of_mass(self, dimension=None, **kwargs):
702+
r"""
703+
Center of mass of the path, assuming constant density.
704+
705+
The center of mass is defined to be the expected value of a vector
706+
located uniformly within either the filled area of the path
707+
(:code:`dimension=2`) or the along path's edge (:code:`dimension=1`) or
708+
along isolated points of the path (:code:`dimension=0`). Notice in
709+
particular that for this definition, if the filled area is used, then
710+
any 0- or 1-dimensional components of the path will not contribute to
711+
the center of mass. Similarly, for if *dimension* is 1, then isolated
712+
points in the path (i.e. "0-dimensional" strokes made up of only
713+
:code:`Path.MOVETO`'s) will not contribute to the center of mass.
714+
715+
For the 2d case, the center of mass is computed using the same
716+
filling strategy as `signed_area`. So, if a path is self-intersecting,
717+
the drawing rule "even-odd" is used and only the filled area is
718+
counted, and all sub paths are treated as if they had been closed. That
719+
is, if there is a MOVETO without a preceding CLOSEPOLY, one is added.
720+
721+
For the 1d measure, the curve is averaged as-is (the implied CLOSEPOLY
722+
is not added).
723+
724+
For the 0d measure, any non-isolated points are ignored.
725+
726+
Parameters
727+
----------
728+
dimension : 2, 1, or 0 (optional)
729+
Whether to compute the center of mass by taking the expected value
730+
of a position uniformly distributed within the filled path
731+
(2D-measure), the path's edge (1D-measure), or between the
732+
discrete, isolated points of the path (0D-measure), respectively.
733+
By default, the intended dimension of the path is inferred by
734+
checking first if `Path.signed_area` is non-zero (implying a
735+
*dimension* of 2), then if the `Path.arc_length` is non-zero
736+
(implying a *dimension* of 1), and finally falling back to the
737+
counting measure (*dimension* of 0).
738+
kwargs : Dict[str, object]
739+
Passed thru to `Path.cleaned` via `Path.iter_bezier`.
740+
741+
Returns
742+
-------
743+
r_cm : (2,) np.array<float>
744+
The center of mass of the path.
745+
746+
Raises
747+
------
748+
ValueError
749+
An empty path has no well-defined center of mass.
750+
751+
In addition, if a specific *dimension* is requested and that
752+
dimension is not well-defined, an error is raised. This can happen
753+
if::
754+
755+
1) 2D expected value was requested but the path has zero area
756+
2) 1D expected value was requested but the path has only
757+
`Path.MOVETO` directives
758+
3) 0D expected value was requested but the path has NO
759+
subsequent `Path.MOVETO` directives.
760+
761+
This error cannot be raised if the function is allowed to infer
762+
what *dimension* to use.
763+
"""
764+
area = None
765+
cleaned = self.cleaned(**kwargs)
766+
move_codes = cleaned.codes == Path.MOVETO
767+
if len(cleaned.codes) == 0:
768+
raise ValueError("An empty path has no center of mass.")
769+
if dimension is None:
770+
dimension = 2
771+
area = cleaned.signed_area()
772+
if not np.isclose(area, 0):
773+
dimension -= 1
774+
if np.all(move_codes):
775+
dimension = 0
776+
if dimension == 2:
777+
# area computation can be expensive, make sure we don't repeat it
778+
if area is None:
779+
area = cleaned.signed_area()
780+
if np.isclose(area, 0):
781+
raise ValueError("2d expected value over empty area is "
782+
"ill-defined.")
783+
return cleaned._2d_center_of_mass(area)
784+
if dimension == 1:
785+
if np.all(move_codes):
786+
raise ValueError("1d expected value over empty arc-length is "
787+
"ill-defined.")
788+
return cleaned._1d_center_of_mass()
789+
if dimension == 0:
790+
adjacent_moves = (move_codes[1:] + move_codes[:-1]) == 2
791+
if len(move_codes) > 1 and not np.any(adjacent_moves):
792+
raise ValueError("0d expected value with no isolated points "
793+
"is ill-defined.")
794+
return cleaned._0d_center_of_mass()
795+
796+
def _2d_center_of_mass(self, normalization=None):
797+
#TODO: refactor this and signed_area (and maybe others, with
798+
# close= parameter)?
799+
if normalization is None:
800+
normalization = self.signed_area()
801+
r_cm = np.zeros(2)
802+
prev_point = None
803+
prev_code = None
804+
start_point = None
805+
for B, code in self.iter_bezier():
806+
if code == Path.MOVETO:
807+
if prev_code is not None and prev_code is not Path.CLOSEPOLY:
808+
Bclose = BezierSegment(np.array([prev_point, start_point]))
809+
r_cm += Bclose.arc_center_of_mass()
810+
start_point = B.control_points[0]
811+
r_cm += B.arc_center_of_mass()
812+
prev_point = B.control_points[-1]
813+
prev_code = code
814+
# add final implied CLOSEPOLY, if necessary
815+
if start_point is not None \
816+
and not np.all(np.isclose(start_point, prev_point)):
817+
Bclose = BezierSegment(np.array([prev_point, start_point]))
818+
r_cm += Bclose.arc_center_of_mass()
819+
return r_cm / normalization
820+
821+
def _1d_center_of_mass(self):
822+
r_cm = np.zeros(2)
823+
Bs = list(self.iter_bezier())
824+
arc_lengths = np.array([B.arc_length() for B in Bs])
825+
r_cms = np.array([B.center_of_mass() for B in Bs])
826+
total_length = np.sum(arc_lengths)
827+
return np.sum(r_cms*arc_lengths)/total_length
828+
829+
def _0d_center_of_mass(self):
830+
move_verts = self.codes
831+
isolated_verts = move_verts.copy()
832+
if len(move_verts) > 1:
833+
isolated_verts[:-1] = (move_verts[:-1] + move_verts[1:]) == 2
834+
isolated_verts[-1] = move_verts[-1]
835+
num_verts = np.sum(isolated_verts)
836+
return np.sum(self.vertices[isolated_verts], axis=0)/num_verts
837+
701838
def interpolated(self, steps):
702839
"""
703840
Return a new path resampled to length N x steps.

‎lib/matplotlib/tests/test_bezier.py

+38
Original file line numberDiff line numberDiff line change
@@ -20,6 +20,11 @@
2020
]
2121

2222

23+
def _integral_center_of_mass(B):
24+
return np.array([scipy.integrate.quad(lambda t: B(t)[0], 0, 1)[0],
25+
scipy.integrate.quad(lambda t: B(t)[1], 0, 1)[0]])
26+
27+
2328
def _integral_arc_length(B):
2429
dB = B.differentiate(B)
2530
def integrand(t):
@@ -37,6 +42,27 @@ def integrand(t):
3742
return scipy.integrate.quad(integrand, 0, 1)[0]
3843

3944

45+
def _integral_arc_com(B):
46+
dB = B.differentiate(B)
47+
def integrand(t):
48+
dr = dB(t).T
49+
n = np.array([dr[1], -dr[0]])
50+
return B(t).T**2 * n / 2
51+
def integrand_x(t):
52+
return integrand(t)[0]
53+
def integrand_y(t):
54+
return integrand(t)[1]
55+
return np.array([
56+
scipy.integrate.quad(integrand_x, 0, 1)[0],
57+
scipy.integrate.quad(integrand_y, 0, 1)[0]
58+
])
59+
60+
61+
def _integral_com(B):
62+
return np.array([scipy.integrate.quad(lambda t: B(t)[0], 0, 1)[0],
63+
scipy.integrate.quad(lambda t: B(t)[1], 0, 1)[0]])
64+
65+
4066
def test_area_formula():
4167
for B in test_curves:
4268
assert(np.isclose(_integral_arc_area(B), B.arc_area()))
@@ -46,3 +72,15 @@ def test_length_iteration():
4672
for B in test_curves:
4773
assert(np.isclose(_integral_arc_length(B), B.arc_length(
4874
rtol=1e-5, atol=1e-8), rtol=1e-5, atol=1e-8))
75+
76+
77+
def test_center_of_mass_1d():
78+
for B in test_curves:
79+
assert(np.all(np.isclose(B.center_of_mass(),
80+
_integral_center_of_mass(B))))
81+
82+
83+
def test_center_of_mass_2d():
84+
for B in test_curves:
85+
assert(np.all(np.isclose(B.arc_center_of_mass(),
86+
_integral_arc_com(B))))

0 commit comments

Comments
 (0)