Movatterモバイル変換

[0]ホーム
URL:

Skip to content

Navigation Menu

Sign in
Appearance settings

Search code, repositories, users, issues, pull requests...

Provide feedback

We read every piece of feedback, and take your input very seriously.

Saved searches

Use saved searches to filter your results more quickly

 Sign up
Appearance settings

Commit348e7e1

Browse files
committed
add function to get center of mass of path
1 parent8d1e757 commit348e7e1

File tree

3 files changed

+253
-2
lines changed

3 files changed

+253
-2
lines changed

‎lib/matplotlib/bezier.py

Lines changed: 76 additions & 0 deletions
Original file line numberDiff line numberDiff line change
@@ -297,6 +297,73 @@ def arc_length(self, rtol=None, atol=None):
297297
returncurr_est
298298
returncurr_est
299299

300+
defarc_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+
foriinrange(n+1):
361+
forjinrange(n+1):
362+
forkinrange(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+
returnr_cm/6
366+
300367
defarc_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+
defcenter_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+
returnnp.mean(self._cpoints,axis=0)
463+
388464
@classmethod
389465
defdifferentiate(cls,B):
390466
"""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
@@ -694,10 +694,147 @@ def signed_area(self, **kwargs):
694694
# add final implied CLOSEPOLY, if necessary
695695
ifstart_pointisnotNone \
696696
andnotnp.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
returnarea
700700

701+
defcenter_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+
iflen(cleaned.codes)==0:
768+
raiseValueError("An empty path has no center of mass.")
769+
ifdimensionisNone:
770+
dimension=2
771+
area=cleaned.signed_area()
772+
ifnotnp.isclose(area,0):
773+
dimension-=1
774+
ifnp.all(move_codes):
775+
dimension=0
776+
ifdimension==2:
777+
# area computation can be expensive, make sure we don't repeat it
778+
ifareaisNone:
779+
area=cleaned.signed_area()
780+
ifnp.isclose(area,0):
781+
raiseValueError("2d expected value over empty area is "
782+
"ill-defined.")
783+
returncleaned._2d_center_of_mass(area)
784+
ifdimension==1:
785+
ifnp.all(move_codes):
786+
raiseValueError("1d expected value over empty arc-length is "
787+
"ill-defined.")
788+
returncleaned._1d_center_of_mass()
789+
ifdimension==0:
790+
adjacent_moves= (move_codes[1:]+move_codes[:-1])==2
791+
iflen(move_codes)>1andnotnp.any(adjacent_moves):
792+
raiseValueError("0d expected value with no isolated points "
793+
"is ill-defined.")
794+
returncleaned._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+
ifnormalizationisNone:
800+
normalization=self.signed_area()
801+
r_cm=np.zeros(2)
802+
prev_point=None
803+
prev_code=None
804+
start_point=None
805+
forB,codeinself.iter_bezier():
806+
ifcode==Path.MOVETO:
807+
ifprev_codeisnotNoneandprev_codeisnotPath.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+
ifstart_pointisnotNone \
816+
andnotnp.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+
returnr_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()forBinBs])
825+
r_cms=np.array([B.center_of_mass()forBinBs])
826+
total_length=np.sum(arc_lengths)
827+
returnnp.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+
iflen(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+
returnnp.sum(self.vertices[isolated_verts],axis=0)/num_verts
837+
701838
definterpolated(self,steps):
702839
"""
703840
Return a new path resampled to length N x steps.

‎lib/matplotlib/tests/test_bezier.py

Lines changed: 38 additions & 0 deletions
Original file line numberDiff line numberDiff line change
@@ -20,6 +20,11 @@
2020
]
2121

2222

23+
def_integral_center_of_mass(B):
24+
returnnp.array([scipy.integrate.quad(lambdat:B(t)[0],0,1)[0],
25+
scipy.integrate.quad(lambdat:B(t)[1],0,1)[0]])
26+
27+
2328
def_integral_arc_length(B):
2429
dB=B.differentiate(B)
2530
defintegrand(t):
@@ -37,6 +42,27 @@ def integrand(t):
3742
returnscipy.integrate.quad(integrand,0,1)[0]
3843

3944

45+
def_integral_arc_com(B):
46+
dB=B.differentiate(B)
47+
defintegrand(t):
48+
dr=dB(t).T
49+
n=np.array([dr[1],-dr[0]])
50+
returnB(t).T**2*n/2
51+
defintegrand_x(t):
52+
returnintegrand(t)[0]
53+
defintegrand_y(t):
54+
returnintegrand(t)[1]
55+
returnnp.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+
returnnp.array([scipy.integrate.quad(lambdat:B(t)[0],0,1)[0],
63+
scipy.integrate.quad(lambdat:B(t)[1],0,1)[0]])
64+
65+
4066
deftest_area_formula():
4167
forBintest_curves:
4268
assert(np.isclose(_integral_arc_area(B),B.arc_area()))
@@ -46,3 +72,15 @@ def test_length_iteration():
4672
forBintest_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+
deftest_center_of_mass_1d():
78+
forBintest_curves:
79+
assert(np.all(np.isclose(B.center_of_mass(),
80+
_integral_center_of_mass(B))))
81+
82+
83+
deftest_center_of_mass_2d():
84+
forBintest_curves:
85+
assert(np.all(np.isclose(B.arc_center_of_mass(),
86+
_integral_arc_com(B))))

0 commit comments

Comments
 (0)
[8]ページ先頭
©2009-2025 Movatter.jp