8000 add function to get center of mass of path · matplotlib/matplotlib@4b5fc44 · GitHub
[go: up one dir, main page]
Skip to content

Commit 4b5fc44

Browse files
brunobeltranbrunobeltran
committed
add function to get center of mass of path
1 parent a76069b commit 4b5fc44

File tree

2 files changed

+215
-2
lines changed

2 files changed

+215
-2
lines changed

‎lib/matplotlib/bezier.py

Lines changed: 76 additions & 0 deletions
Original file line numberDiff line numberDiff line change
@@ -306,6 +306,73 @@ def arc_length(self, rtol=None, atol=None):
306306
return curr_est
307307
return curr_est
308308

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

464+
def center_of_mass(self):
465+
"""Return the center of mass of the curve (not the filled curve!)
466+
467+
Notes
468+
-----
469+
Computed as the mean of the control points.
470+
"""
471+
return np.mean(self._cpoints, axis=0)
472+
397473
@classmethod
398474
def differentiate(cls, B):
399475
"""Return the derivative of a BezierSegment, itself a BezierSegment"""

‎lib/matplotlib/path.py

Lines changed: 139 additions & 2 deletions
Original file line numberDiff line numberDiff line change
@@ -706,10 +706,147 @@ def signed_area(self, **kwargs):
706706
# add final implied CLOSEPOLY, if necessary
707707
if start_point is not None \
708708
and not np.all(np.isclose(start_point, prev_point)):
709-
B = BezierSegment(np.array([prev_point, start_point]))
710-
area += B.arc_area()
709+
Bclose = BezierSegment(np.array([prev_point, start_point]))
710+
area += Bclose.arc_area()
711711
return area
712712

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

0 commit comments

Comments
 (0)
0