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

Commit a933565

Browse files
brunobeltranbrunobeltran
committed
add function to get center of mass of path
1 parent aca0710 commit a933565

File tree

2 files changed

+213
-4
lines changed

2 files changed

+213
-4
lines changed

‎lib/matplotlib/bezier.py

Lines changed: 72 additions & 0 deletions
Original file line numberDiff line numberDiff line change
@@ -218,6 +218,73 @@ def point_at_t(self, t):
218218
def arc_length(self, tol=1e-6):
219219
pass
220220

221+
def arc_center_of_mass(self):
222+
r"""
223+
Center of mass of the (even-odd-rendered) area swept out by the ray
224+
from the origin to the path.
225+
226+
Summing this vector for each segment along a closed path will produce
227+
that area's center of mass.
228+
229+
Returns
230+
-------
231+
r_cm : (2,) np.array<float>
232+
the "arc's center of mass"
233+
234+
Notes
235+
-----
236+
A simple analytical form can be derived for general Bezier curves.
237+
Suppose the curve was closed, so :math:`B(0) = B(1)`. Call the area
238+
enclosed by :math:`B(t)` :math:`B_\text{int}`. The center of mass of
239+
:math:`B_\text{int}` is defined by the expected value of the position
240+
vector `\vec{r}`
241+
242+
.. math::
243+
244+
\vec{R}_\text{cm} = \int_{B_\text{int}} \vec{r} \left( \frac{1}{
245+
\int_{B_\text{int}}} d\vec{r} \right) d\vec{r}
246+
247+
where :math:`(1/\text{Area}(B_\text{int})` can be interpreted as a
248+
probability density.
249+
250+
In order to compute this integral, we choose two functions
251+
:math:`F_0(x,y) = [x^2/2, 0]` and :math:`F_1(x,y) = [0, y^2/2]` such
252+
that :math:`[\div \cdot F_0, \div \cdot F_1] = \vec{r}`. Then, applying
253+
the divergence integral (componentwise), we get that
254+
255+
.. math::
256+
\vec{R}_\text{cm} &= \oint_{B(t)} F \cdot \vec{n} dt \\
257+
&= \int_0^1 \left[ \begin{array}{1}
258+
B^{(0)}(t) \frac{dB^{(1)}(t)}{dt} \\
259+
- B^{(1)}(t) \frac{dB^{(0)}(t)}{dt} \end{array} \right] dt
260+
261+
After expanding in Berstein polynomials and moving the integral inside
262+
all the sums, we get that
263+
264+
.. math::
265+
\vec{R}_\text{cm} = \frac{1}{6} \sum_{i,j=0}^n\sum_{k=0}^{n-1}
266+
\frac{{n \choose i}{n \choose j}{{n-1} \choose k}}
267+
{{3n - 1} \choose {i + j + k}}
268+
\left(\begin{array}{1}
269+
P^{(0)}_i P^{(0)}_j (P^{(1)}_{k+1} - P^{(1)}_k)
270+
- P^{(1)}_i P^{(1)}_j (P^{(0)}_{k+1} - P^{(0)}_k)
271+
\right) \end{array}
272+
273+
where :math:`P_i = [P^{(0)}_i, P^{(1)}_i]` is the :math:`i`'th control
274+
point of the curve and :math:`n` is the degree of the curve.
275+
"""
276+
n = self.degree
277+
r_cm = np.zeros(2)
278+
P = self.control_points
279+
dP = np.diff(P, axis=0)
280+
Pn = np.array([[1, -1]])*dP[:, ::-1] # n = [y, -x]
281+
for i in range(n + 1):
282+
for j in range(n + 1):
283+
for k in range(n):
284+
r_cm += _comb(n, i) * _comb(n, j) * _comb(n - 1, k) \
285+
* P[i]*P[j]*Pn[k] / _comb(3*n - 1, i + j + k)
286+
return r_cm/6
287+
221288
def arc_area(self):
222289
r"""
223290
(Signed) area swept out by ray from origin to curve.
@@ -306,6 +373,11 @@ def arc_area(self):
306373
* (P[j, 0]*dP[k, 1] - P[j, 1]*dP[k, 0])
307374
return (1/4)*area
308375

376+
def center_of_mass(self):
377+
# return np.array([scipy.integrate.quad(lambda t: B(t)[0], 0, 1)[0],
378+
# scipy.integrate.quad(lambda t: B(t)[1], 0, 1)[0]])
379+
pass
380+
309381
@classmethod
310382
def differentiate(cls, B):
311383
"""Return the derivative of a BezierSegment, itself a BezierSegment"""

‎lib/matplotlib/path.py

Lines changed: 141 additions & 4 deletions
Original file line numberDiff line numberDiff line change
@@ -467,10 +467,10 @@ def iter_bezier(self, **kwargs):
467467
yield BezierSegment(np.array([prev_vert, vertices])), code
468468
elif code == Path.CURVE3:
469469
yield BezierSegment(np.array([prev_vert, vertices[:2],
470-
vertices[2:]])), code
470+
vertices[2:]])), code
471471
elif code == Path.CURVE4:
472472
yield BezierSegment(np.array([prev_vert, vertices[:2],
473-
vertices[2:4], vertices[4:]])), code
473+
vertices[2:4], vertices[4:]])), code
474474
elif code == Path.CLOSEPOLY:
475475
yield BezierSegment(np.array([prev_vert, first_vert])), code
476476
elif code == Path.STOP:
@@ -690,10 +690,147 @@ def signed_area(self, **kwargs):
690690
# add final implied CLOSEPOLY, if necessary
691691
if start_point is not None \
692692
and not np.all(np.isclose(start_point, prev_point)):
693-
B = BezierSegment(np.array([prev_point, start_point]))
694-
area += B.arc_area()
693+
Bclose = BezierSegment(np.array([prev_point, start_point]))
694+
area += Bclose.arc_area()
695695
return area
696696

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

0 commit comments

Comments
 (0)
0