20. `geomtools` — Basic geometrical operations.¶

This module defines some basic operations on simple geometrical entities such as points, lines, planes, vectors, segments, triangles, circles.

The operations include intersection, projection, distance computing.

Many of the functions in this module are exported as methods on the Coords and the Geometry derived classes.

Renamed functions:

intersectionPointsLWL    -> intersectLineWithLine
intersectionTimesLWL     -> intersectLineWithLine
intersectionPointsLWP    -> intersectLineWithPlane
intersectionTimesLWP     -> intersectLineWithPlane

intersectionTimesPOP     -> projectPointOnPlane
intersectionPointsPOP    -> projectPointOnPlane
intersectionTimesPOL     -> projectPointOnLine
intersectionPointsPOL    -> projectPointOnLine

distancesPFL             -> distanceFromLine
distancesPFS             -> distanceFromLine

Planned renaming of functions:

areaNormals
degenerate
levelVolumes
smallestDirection
distance                 -> distance
closest
closestPair
projectedArea
polygonNormals
averageNormals
triangleInCircle
triangleCircumCircle
triangleBoundingCircle
triangleObtuse
lineIntersection         -> to be removed (or intersect2DLineWithLine)
displaceLines
segmentOrientation
rotationAngle
anyPerpendicularVector
perpendicularVector
projectionVOV            -> projectVectorOnVector
projectionVOP            -> projectVectorOnPlane


pointsAtLines            -> pointOnLine
pointsAtSegments         -> pointOnSegment
intersectionTimesSWP     -> remove or intersectLineWithLine
intersectionSWP          -> intersectSegmentWithPlane
intersectionPointsSWP    -> intersectSegmentWithPlane
intersectionTimesLWT     -> intersectLineWithTriangle
intersectionPointsLWT    -> intersectLineWithTriangle
intersectionTimesSWT     -> intersectSegmentWithTriangle
intersectionPointsSWT    -> intersectSegmentWithTriangle
intersectionPointsPWP    -> intersectThreePlanes
intersectionLinesPWP     -> intersectTwoPlanes (or intersectPlaneWithPlane)
intersectionSphereSphere -> intersectTwoSpheres (or intersectSphereWithSphere)

faceDistance             ->
edgeDistance             ->
vertexDistance           ->
baryCoords
insideSimplex
insideTriangle

20.1. Classes defined in module geomtools¶

class geomtools.Lines(data)[source]¶

A collection of lines, half-lines(rays) or segments.

This class is intended as a common interface for a collection of lines, half-lines or segments. Many functions in the geomtools module operate on lines (or rays, or segments). Sometimes these are specified by two points, while in other case they need to be specified by one point and a vector. Also the user might store his line data in one of these two ways, and it is annoying and not effective to have to convert between these two multiple times.

Then this class comes as a help. It is a lightweight class that can store lines, rays or segments. The class itself does not discriminate between these: the user is responsible for the interpretation. This is further complicated by the fact that the data structure of a vector is the same as that of a point.

In what follows, a line can mean any of:

a full twosided infinite line,
a ray starting at some point and infinite in one direction,
a finite straight segment between two points.

A Lines instance can be initialized by any of the following data:

a coords_like object with shape (nlines, 2, 3), containing the coordinates of two points on each of the nlines lines. This is the natural way to specify segments, but it can be used for infinite lines and rays as well.
a tuple of two coords_like objects, each with shape (nlines, 2, 3). The first contains a point on each of the lines, the second is a vector along the line. The vectors do not need to be unit length vectors. They typically will not be if this method is used to specify segments.
another Lines instance. In this case a shallow copy of the Lines is created, sharing its data with the input Lines. This is a convenience allowing the user to pass raw data as well as data already converted to a Lines instance in places where line data are expected.

Parameters:

data (coords_like | tuple | Lines) –

A coords_like argument should have a shape (nlines,2,3) and defines nlines lines by the coordinates of two points on each of the lines.

A tuple should contain two (nlines,3) shaped coords_like structures. The first holds a point on the lines, the second a vector along the lines.

If a Lines instance is specified as data, a shallow copy sharing the same data is created.

A Lines instance has the following attributes (onbly the first two are actually stored):

p¶

The first point of the Lines.

Type:	array (nlines,3)

n¶

The vector from the first to the second point.

Type:	array (nlines,3)

coords¶

A Coords containing both points on the lines.

Type:	Coords (nlines, 2, 3)

Examples

Create Lines from two points.

>>> L = Lines([[[2.,0.,0.],[2.,3.,0.]], [[0.,1.,0.],[1.,1.,0.]]])
>>> print(L.p)
[[ 2.  0.  0.]
 [ 0.  1.  0.]]
>>> print(L.n)
[[ 0.  3.  0.]
 [ 1.  0.  0.]]
>>> L.coords
Coords([[[ 2.,  0.,  0.],
         [ 2.,  3.,  0.]],
<BLANKLINE>
        [[ 0.,  1.,  0.],
         [ 1.,  1.,  0.]]])

Create the same Lines from one point and a vector.

>>> M = Lines(([[2.,0.,0.],[0.,1.,0.]], [[0.,3.,0.],[1.,0.,0.]]))
>>> print((M.p == L.p).all() and (M.n == L.n).all())
True
>>> print(id(M.p) == id(L.p) and id(M.n) == id(L.n))
False

Create a shallow copy.

>>> M = Lines(L)
>>> print((M.p == L.p).all() and (M.n == L.n).all())
True
>>> print(id(M.p) == id(L.p) and id(M.n) == id(L.n))
True

toFormex()[source]¶

Convert a Lines to a 2-plex Formex.

Returns:	Formex (nlines, 2, 3) – A plex-2 Formex representing the Lines.

Notes

When drawn, the Lines will always be represented as straight line segments, even if they represent infinitely long lines.

Examples

>>> L = Lines([[[2.,0.,0.],[2.,3.,0.]], [[0.,1.,0.],[1.,1.,0.]]])
>>> print(L.toFormex())
{[2.0,0.0,0.0; 2.0,3.0,0.0], [0.0,1.0,0.0; 1.0,1.0,0.0]}

20.2. Functions defined in module geomtools¶

geomtools.pointsAtLines(q, m, t)[source]¶

Return the points of lines (q,m) at parameter values t.

Parameters:	q (float array (.., 3)) – Array with (starting) points of a collection of lines. m (float array (.., 3)) – Array with vectors along the lines, broadcast compatible with q. The vectors do not need to have unit length. t (float array) – Array with parameter values, broadcast compatible with q[…] and m[…]. Parametric value 0 is at point q, parametric value 1 is at q + m.
Returns:	Coords – A Coords array with the points at parameter values t.

geomtools.pointsAtSegments(S, t)[source]¶

Return the points of line segments S at parameter values t.

Parameters:	S (float array (.., 2, 3)) – A collection of line segments defined by two points. t (float array) – Array with parameter values, broadcast compatible with S[…]. Parametric value 0 is at point 0, parametric value 1 is at point 1.
Returns:	Coords – A Coords array with the points at parameter values t.

geomtools.intersectLineWithLine(q1, m1, q2, m2, mode='all', times=False)[source]¶

Find the common perpendicular of lines (q1,m1) and lines (q2,m2).

Return the intersection points of lines (q1,m1) and lines (q2,m2) with the perpendiculars between them. For intersecting lines, the corresponding points will coincide.

Parameters:

q1 (float array (nq1, 3)) – Points on the first set of lines.
m1 (float array (nq1, 3)) – Direction vectors of the first set of lines.
q2 (float array (nq2, 3)) – Points on the second set of lines.
m2 (float array (nq2, 3)) – Direction vectors of the second set of lines.
mode ('all' | 'pair') – If ‘all’, the intersection of all lines (q1,m1) with all lines (q2,m2) is computed; nq1 and nq2 can be different. If ‘pair’, nq1 and nq2 should be equal (or 1) and the intersection of pairs of lines is computed (using broadcasting for length 1 data).
times (bool) – If True, return the parameter values of the intersection points instead of the coordinates (see Notes).

Returns:

ar1 (float array) – The intersection points of the common perpendiculars with the first set of lines. See Notes.
ar2 (float array) – The intersection points of the common perpendiculars with the second set of lines. See Notes.

Notes

By default (times=False) the returned arrays are the coordinates of the points, with shape (nq1,nq2,3) for mode ‘all’ and (nq1,3) for mode ‘pair’. With times=True the return values are the parameter values of the intersection points, and the size opf the arrays is (nq1,nq2) for mode ‘all’ and (nq1,) for mode ‘pair’.

The coordinates of a point with parameter value t on a line (q,m) are given by q + t * m.

The intersection of two parallel lines results in NAN-values. These are not removed from the result. The user has to do it himself if needed.

Examples

>>> q,m = [[0,0,0],[0,0,1],[0,0,3]], [[1,0,0],[1,1,0],[0,1,0]]
>>> p,n = [[2.,0.,0.],[0.,0.,0.]], [[0.,1.,0.],[0.,0.,1.]]
>>> x1,x2 = intersectLineWithLine(q,m,p,n)
>>> print(x1)
[[[  2.  0.  0.]
  [  0.  0.  0.]]
<BLANKLINE>
 [[  2.  2.  1.]
  [  0.  0.  1.]]
<BLANKLINE>
 [[ nan nan nan]
  [  0.  0.  3.]]]
>>> print(x2)
[[[  2.  0.  0.]
  [  0.  0.  0.]]
<BLANKLINE>
 [[  2.  2.  0.]
  [  0.  0.  1.]]
<BLANKLINE>
 [[ nan nan nan]
  [  0.  0.  3.]]]

>>> x1,x2 = intersectLineWithLine(q[:2],m[:2],p,n,mode='pair')
>>> print(x1)
[[ 2. 0. 0.]
 [ 0. 0. 1.]]
>>> print(x2)
[[ 2. 0. 0.]
 [ 0. 0. 1.]]

>>> t1,t2 = intersectLineWithLine(q,m,p,n,times=True)
>>> print(t1)
[[  2. -0.]
 [  2. -0.]
 [ nan -0.]]
>>> print(t2)
[[ -0. -0.]
 [  2.  1.]
 [ nan  3.]]

>>> t1,t2 = intersectLineWithLine(q[:2],m[:2],p,n,mode='pair',times=True)
>>> print(t1)
[ 2. -0.]
>>> print(t2)
[-0.  1.]

geomtools.intersectLineWithPlane(q, m, p, n, mode='all', times=False)[source]¶

Find the intersection of lines (q,m) and planes (p,n).

Return the intersection points of lines (q,m) and planes (p,n).

Parameters:

q,`m`: (nq,3) shaped arrays of points and vectors (mode=all) or broadcast compatible arrays (mode=pair), defining a single line or a set of lines.
p,`n`: (np,3) shaped arrays of points and normals (mode=all) or broadcast compatible arrays (mode=pair), defining a single plane or a set of planes.
mode: all to calculate the intersection of each line (q,m) with all planes (p,n) or pair for pairwise intersections.

Returns a (nq,np) shaped (mode=all) array of parameter values t, such that the intersection points are given by q+t*m.

Notice that the result will contain an INF value for lines that are parallel to the plane.

Example:

>>> q,m = [[0,0,0],[0,1,0],[0,0,3]], [[1,0,0],[0,1,0],[0,0,1]]
>>> p,n = [[1.,1.,1.],[1.,1.,1.]], [[1.,1.,0.],[1.,1.,1.]]

>>> t = intersectLineWithPlane(q,m,p,n,times=True)
>>> print(t)
[[  2.  3.]
 [  1.  2.]
 [ inf  0.]]
>>> x = intersectLineWithPlane(q,m,p,n)
>>> print(x)
[[[  2.  0.  0.]
  [  3.  0.  0.]]
<BLANKLINE>
 [[  0.  2.  0.]
  [  0.  3.  0.]]
<BLANKLINE>
 [[ nan nan inf]
  [  0.  0.  3.]]]
>>> x = intersectLineWithPlane(q[:2],m[:2],p,n,mode='pair')
>>> print(x)
[[ 2. 0. 0.]
 [ 0. 3. 0.]]

geomtools.intersectLineWithTriangle(T, p, p1, method='line', atol=1e-05)[source]¶

Compute the intersection points with a set of lines.

Parameters:

T (coords_like (ntri,3,3)) – The coordinates of the three vertices of ntri triangles.
p (coords_like (nlines,3)) – A first point for each of the lines to intersect.
p1 (coords_like (nlines,3)) – The second point defining the lines to intersect.
method ('line' | 'segment' | ' ray') – Define whether the points p and p1 define an infinitely long line, a finite segment p-p1 or a half infinite ray (p->p1).
atol (float) – Tolerance to be used in deciding whether an intersection point on a border edge is inside the surface or not.

Returns:

X (Coords (nipts, 3)) – A fused set of intersection points
ind (int array (nipts, 3)) – An array identifying the related intersection points, lines and triangle.

Notes

A line laying in the plane of a triangle does not generate intersections.

This method is faster than the similar function intersectionPointsLWT().

Examples

>>> T = Formex('3:.12.34').coords
>>> L = Coords([[[0.,0.,0.], [0.,0.,1.]],
...             [[0.5,0.5,0.5], [0.5,0.5,1.]],
...             [[0.2,0.7,0.5], [0.2,0.8,0.5]],
...             [[0.2,0.7,0.5], [0.2,0.7,0.2]],
...             ])
>>> P, ind = intersectLineWithTriangle(T, L[:,0,:], L[:,1,:])
>>> print(P)
[[ 0.   0.   0. ]
 [ 0.5  0.5  0. ]
 [ 0.2  0.7  0. ]]
>>> print(ind)
[[0 0 0]
 [0 0 1]
 [1 1 0]
 [1 1 1]
 [2 3 1]]
>>> P, ind = intersectLineWithTriangle(T, L[:,0,:], L[:,1,:],
...          method='ray')
>>> print(P)
[[ 0.   0.   0. ]
 [ 0.2  0.7  0. ]]
>>> print(ind)
[[0 0 0]
 [0 0 1]
 [1 3 1]]
>>> P, ind = intersectLineWithTriangle(T, L[:,0,:], L[:,1,:],
...          method='segment')
>>> print(P)
[[ 0.  0.  0.]]
>>> print(ind)
[[0 0 0]
 [0 0 1]]

geomtools.intersectionTimesSWP(S, p, n, mode='all')[source]¶

Return the intersection of line segments S with planes (p,n).

This is like intersectionTimesLWP, but the lines are defined by two points instead of by a point and a vector. Parameters:

S: (nS,2,3) shaped array (mode=all) or broadcast compatible array (mode=pair), defining one or more line segments.
p,`n`: (np,3) shaped arrays of points and normals (mode=all) or broadcast compatible arrays (mode=pair), defining a single plane or a set of planes.
mode: all to calculate the intersection of each line segment S with all planes (p,n) or pair for pairwise intersections.

Returns a (nS,np) shaped (mode=all) array of parameter values t, such that the intersection points are given by (1-t)*S[…,0,:] + t*S[…,1,:].

geomtools.intersectionSWP(S, p, n, mode='all', return_all=False, atol=0.0)[source]¶

Return the intersection points of line segments S with planes (p,n).

Parameters:

S: (nS,2,3) shaped array, defining a single line segment or a set of line segments.
p,`n`: (np,3) shaped arrays of points and normals, defining a single plane or a set of planes.
mode: all to calculate the intersection of each line segment S with all planes (p,n) or pair for pairwise intersections.
return_all: if True, all intersection points of the lines along the segments are returned. Default is to return only the points that lie on the segments.
atol: float tolerance of the points inside the line segments.

Return values if return_all==True:

t: (nS,NP) parametric values of the intersection points along the line segments.
x: the intersection points themselves (nS,nP,3).

Return values if return_all==False:

t: (n,) parametric values of the intersection points along the line segments (n <= nS*nP)
x: the intersection points themselves (n,3).
wl: (n,) line indices corresponding with the returned intersections.
wp: (n,) plane indices corresponding with the returned intersections

geomtools.intersectionPointsSWP(S, p, n, mode='all', return_all=False, atol=0.0)[source]¶

Return the intersection points of line segments S with planes (p,n).

This is like intersectionSWP() but does not return the parameter values. It is equivalent to:

intersectionSWP(S,p,n,mode,return_all)[1:]

geomtools.intersectionTimesLWT(q, m, F, mode='all')[source]¶

Return the intersection of lines (q,m) with triangles F.

Parameters:

q,`m`: (nq,3) shaped arrays of points and vectors (mode=all) or broadcast compatible arrays (mode=pair), defining a single line or a set of lines.
F: (nF,3,3) shaped array (mode=all) or broadcast compatible array (mode=pair), defining a single triangle or a set of triangles.
mode: all to calculate the intersection of each line (q,m) with all triangles F or pair for pairwise intersections.

Returns a (nq,nF) shaped (mode=all) array of parameter values t,: such that the intersection points are given q+tm.

geomtools.intersectionPointsLWT(q, m, F, mode='all', return_all=False)[source]¶

Return the intersection points of lines (q,m) with triangles F.

Parameters:

q,`m`: (nq,3) shaped arrays of points and vectors, defining a single line or a set of lines.
F: (nF,3,3) shaped array, defining a single triangle or a set of triangles.
mode: all to calculate the intersection points of each line (q,m) with all triangles F or pair for pairwise intersections.
return_all: if True, all intersection points are returned. Default is to return only the points that lie inside the triangles.

Returns:	If return_all==True, a (nq,nF,3) shaped (mode=all) array of intersection points, else, a tuple of intersection points with shape (n,3) and line and plane indices with shape (n), where n <= nq*nF.

geomtools.intersectionTimesSWT(S, F, mode='all')[source]¶

Return the intersection of lines segments S with triangles F.

Parameters:

S: (nS,2,3) shaped array (mode=all) or broadcast compatible array (mode=pair), defining a single line segment or a set of line segments.
F: (nF,3,3) shaped array (mode=all) or broadcast compatible array (mode=pair), defining a single triangle or a set of triangles.
mode: all to calculate the intersection of each line segment S with all triangles F or pair for pairwise intersections.

Returns a (nS,nF) shaped (mode=all) array of parameter values t, such that the intersection points are given by (1-t)*S[…,0,:] + t*S[…,1,:].

geomtools.intersectionPointsSWT(S, F, mode='all', return_all=False)[source]¶

Return the intersection points of lines segments S with triangles F.

Parameters:

S: (nS,2,3) shaped array, defining a single line segment or a set of line segments.
F: (nF,3,3) shaped array, defining a single triangle or a set of triangles.
mode: all to calculate the intersection points of each line segment S with all triangles F or pair for pairwise intersections.
return_all: if True, all intersection points are returned. Default is to return only the points that lie on the segments and inside the triangles.

Returns:	If return_all==True, a (nS,nF,3) shaped (mode=all) array of intersection points, else, a tuple of intersection points with shape (n,3) and line and plane indices with shape (n), where n <= nS*nF.

geomtools.intersectionPointsPWP(p1, n1, p2, n2, p3, n3, mode='all')[source]¶

Return the intersection points of planes (p1,n1), (p2,n2) and (p3,n3).

Parameters:

pi,`ni` (i=1…3): (npi,3) shaped arrays of points and normals (mode=all) or broadcast compatible arrays (mode=pair), defining a single plane or a set of planes.
mode: all to calculate the intersection of each plane (p1,n1) with all planes (p2,n2) and (p3,n3) or pair for pairwise intersections.

Returns a (np1,np2,np3,3) shaped (mode=all) array of intersection points.

geomtools.intersectionLinesPWP(p1, n1, p2, n2, mode='all')[source]¶

Return the intersection lines of planes (p1,n1) and (p2,n2).

Parameters:

pi,`ni` (i=1…2): (npi,3) shaped arrays of points and normals (mode=all) or broadcast compatible arrays (mode=pair), defining a single plane or a set of planes.
mode: all to calculate the intersection of each plane (p1,n1) with all planes (p2,n2) or pair for pairwise intersections.

Returns a tuple of (np1,np2,3) shaped (mode=all) arrays of intersection points q and vectors m, such that the intersection lines are given by q+t*m.

geomtools.intersectionSphereSphere(R, r, d)[source]¶

Intersection of two spheres (or two circles in the x,y plane).

Computes the intersection of two spheres with radii R, resp. r, having their centres at distance d <= R+r. The intersection is a circle with its center on the segment connecting the two sphere centers at a distance x from the first sphere, and having a radius y. The return value is a tuple x,y.

geomtools.projectPointOnPlane(X, p, n, mode='all')[source]¶

Return the projection of points X on planes (p,n).

Parameters:

X: a (nx,3) shaped array of points.
p, n: (np,3) shaped arrays of points and normals defining np planes.
mode: ‘all’ or ‘pair:
- if ‘all’, the projection of all points on all planes is computed; nx and np can be different.
- if ‘pair’: nx and np should be equal (or 1) and the projection of pairs of point and plane are computed (using broadcasting for length 1 data).

Returns a float array of size (nx,np,3) for mode ‘all’, or size (nx,3) for mode ‘pair’.

Example:

>>> X = Coords([[0.,1.,0.],[3.,0.,0.],[4.,3.,0.]])
>>> p,n = [[2.,0.,0.],[0.,1.,0.]], [[1.,0.,0.],[0.,1.,0.]]
>>> print(projectPointOnPlane(X,p,n))
[[[ 2. 1. 0.]
  [ 0. 1. 0.]]
<BLANKLINE>
 [[ 2. 0. 0.]
  [ 3. 1. 0.]]
<BLANKLINE>
 [[ 2. 3. 0.]
  [ 4. 1. 0.]]]
>>> print(projectPointOnPlane(X[:2],p,n,mode='pair'))
[[ 2. 1. 0.]
 [ 3. 1. 0.]]

geomtools.projectPointOnPlaneTimes(X, p, n, mode='all')[source]¶

Return the projection of points X on planes (p,n).

This is like projectPointOnPlane() but instead of returning the projected points, returns the parametric values t along the lines (X,n), such that the projection points can be computed from X+t*n.

Parameters: see projectPointOnPlane().

Returns a float array of size (nx,np) for mode ‘all’, or size (nx,) for mode ‘pair’.

Example:

>>> X = Coords([[0.,1.,0.],[3.,0.,0.],[4.,3.,0.]])
>>> p,n = [[2.,0.,0.],[0.,1.,0.]], [[1.,0.,0.],[0.,1.,0.]]
>>> print(projectPointOnPlaneTimes(X,p,n))
[[ 2.  0.]
 [-1.  1.]
 [-2. -2.]]
>>> print(projectPointOnPlaneTimes(X[:2],p,n,mode='pair'))
[ 2. 1.]

geomtools.projectPointOnLine(X, p, n, mode='all')[source]¶

Return the projection of points X on lines (p,n).

Parameters:

X: a (nx,3) shaped array of points.
p, n: (np,3) shaped arrays of points and vectors defining np lines.
mode: ‘all’ or ‘pair:
- if ‘all’, the projection of all points on all lines is computed; nx and np can be different.
- if ‘pair’: nx and np should be equal (or 1) and the projection of pairs of point and line are computed (using broadcasting for length 1 data).

Returns a float array of size (nx,np,3) for mode ‘all’, or size (nx,3) for mode ‘pair’.

Example:

>>> X = Coords([[0.,1.,0.],[3.,0.,0.],[4.,3.,0.]])
>>> p,n = [[2.,0.,0.],[0.,1.,0.]], [[0.,2.,0.],[1.,0.,0.]]
>>> print(projectPointOnLine(X,p,n))
[[[ 2. 1. 0.]
  [ 0. 1. 0.]]
<BLANKLINE>
 [[ 2. 0. 0.]
  [ 3. 1. 0.]]
<BLANKLINE>
 [[ 2. 3. 0.]
  [ 4. 1. 0.]]]
>>> print(projectPointOnLine(X[:2],p,n,mode='pair'))
[[ 2. 1. 0.]
 [ 3. 1. 0.]]

geomtools.projectPointOnLineTimes(X, p, n, mode='all')[source]¶

Return the projection of points X on lines (p,n).

This is like projectPointOnLine() but instead of returning the projected points, returns the parametric values t along the lines (X,n), such that the projection points can be computed from p+t*n.

Parameters: see projectPointOnLine().

Returns a float array of size (nx,np) for mode ‘all’, or size (nx,) for mode ‘pair’.

Example:

>>> X = Coords([[0.,1.,0.],[3.,0.,0.],[4.,3.,0.]])
>>> p,n = [[2.,0.,0.],[0.,1.,0.]], [[0.,1.,0.],[1.,0.,0.]]
>>> print(projectPointOnLineTimes(X,p,n))
[[ 1. 0.]
 [ 0. 3.]
 [ 3. 4.]]
>>> print(projectPointOnLineTimes(X[:2],p,n,mode='pair'))
[ 1. 3.]

geomtools.distanceFromLine(X, lines, mode='all')[source]¶

Return the distance of points X from lines (p,n).

Parameters:

X (coords_like (nx,3)) – A collection of points.
lines (line_like) –
One of the following definitions of the line(s):
- a tuple (p,n), where both p and n are (np,3) shaped arrays of respectively points and vectors defining np lines;
- an (np,2,3) shaped array containing two points of each line.
mode ('all' or 'pair:) – If ‘all’, the distance of all points to all lines is computed; nx and np can be different. If ‘pair’: nx and np should be equal (or 1) and the distance of pairs of point and line are computed (using broadcasting for length 1 data).

Returns:

float array – A float array of size (nx,np) for mode ‘all’, or size (nx) for mode ‘pair’, with the distances between the points and the lines.

Examples

>>> X = Coords([[0.,1.,0.],[3.,0.,0.],[4.,3.,0.]])
>>> L = Lines(([[2.,0.,0.],[0.,1.,0.]], [[0.,3.,0.],[1.,0.,0.]]))
>>> print(distanceFromLine(X,L))
[[ 2. 0.]
 [ 1. 1.]
 [ 2. 2.]]
>>> print(distanceFromLine(X[:2],L,mode='pair'))
[ 2. 1.]
>>> L = Lines(([[[2.,0.,0.],[2.,2.,0.]], [[0.,1.,0.],[1.,1.,0.]]]))
>>> print(distanceFromLine(X,L))
[[ 2. 0.]
 [ 1. 1.]
 [ 2. 2.]]
>>> print(distanceFromLine(X[:2],L,mode='pair'))
[ 2. 1.]

geomtools.pointNearLine(X, p, n, atol, nproc=1)[source]¶

Find the points from X that are near to lines (p,n).

Finds the points from X that are closer than atol to any of the lines (p,n).

Parameters:	X (coords_like (npts,3)) – An array of points. p (coords_like (nlines,3)) – An array of points which together with n define the lines. n (coords_like (nlines,3)) – An array of direction vectors which together with p define the lines. atol (float or float array (npts,)) – A global or pointwise tolerance to be used in determining close points.
Returns:	list of int arrays – An index array for each of the lines, holding the indices of the points that are close to that line.

Examples

>>> X = Coords([[0.,1.,0.],[3.,0.,0.],[4.,3.,0.]])
>>> p,n = [[2.,0.,0.],[0.,1.,0.]], [[0.,3.,0.],[1.,0.,0.]]
>>> print(pointNearLine(X,p,n,1.5))
[array([1]), array([0, 1])]
>>> print(pointNearLine(X,p,n,1.5,2))
[array([1]), array([0, 1])]

geomtools.faceDistance(X, Fp, return_points=False)[source]¶

Compute the closest perpendicular distance to a set of triangles.

X is a (nX,3) shaped array of points. Fp is a (nF,3,3) shaped array of triangles.

Note that some points may not have a normal with footpoint inside any of the facets.

The return value is a tuple OKpid,OKdist,OKpoints where:

OKpid is an array with the point numbers having a normal distance;
OKdist is an array with the shortest distances for these points;
OKpoints is an array with the closest footpoints for these points and is only returned if return_points = True.

geomtools.edgeDistance(X, Ep, return_points=False)[source]¶

Compute the closest perpendicular distance of points X to a set of edges.

X is a (nX,3) shaped array of points. Ep is a (nE,2,3) shaped array of edge vertices.

Note that some points may not have a normal with footpoint inside any of the edges.

The return value is a tuple OKpid,OKdist,OKpoints where:

OKpid is an array with the point numbers having a normal distance;
OKdist is an array with the shortest distances for these points;
OKpoints is an array with the closest footpoints for these points and is only returned if return_points = True.

geomtools.vertexDistance(X, Vp, return_points=False)[source]¶

Compute the closest distance of points X to a set of vertices.

X is a (nX,3) shaped array of points. Vp is a (nV,3) shaped array of vertices.

The return value is a tuple OKdist,OKpoints where:

OKdist is an array with the shortest distances for the points;
OKpoints is an array with the closest vertices for the points and is only returned if return_points = True.

geomtools.areaNormals(x)[source]¶

Compute the area and normal vectors of a collection of triangles.

x is an (ntri,3,3) array with the coordinates of the vertices of ntri triangles.

Returns a tuple (areas,normals) with the areas and the normals of the triangles. The area is always positive. The normal vectors are normalized.

geomtools.degenerate(area, normals)[source]¶

Return a list of the degenerate faces according to area and normals.

area,normals are equal sized arrays with the areas and normals of a list of faces, such as the output of the areaNormals() function.

A face is degenerate if its area is less or equal to zero or the normal has a nan (not-a-number) value.

Returns a list of the degenerate element numbers as a sorted array.

geomtools.hexVolume(x)[source]¶

Compute the volume of hexahedrons.

Parameters:

x: float array (nelems,8,3)

Returns a float array (nelems) withe the approximate volume of the hexahedrons formed by each 8-tuple of vertices. The volume is obained by dividing the hexahedron in 24 tetrahedrons and using the formulas from http://www.osti.gov/scitech/servlets/purl/632793

Example:

>>> from pyformex.elements import Hex8
>>> X = Coords(Hex8.vertices).reshape(-1,8,3)
>>> print(hexVolume(X))
[ 1.]

geomtools.levelVolumes(x)[source]¶

Compute the level volumes of a collection of elements.

x is an (nelems,nplex,3) array with the coordinates of the nplex vertices of nelems elements, with nplex equal to 2, 3 or 4.

If nplex == 2, returns the lengths of the straight line segments. If nplex == 3, returns the areas of the triangle elements. If nplex == 4, returns the signed volumes of the tetrahedron elements. Positive values result if vertex 3 is at the positive side of the plane defined by the vertices (0,1,2). Negative volumes are reported for tetrahedra having reversed vertex order.

For any other value of nplex, raises an error. If successful, returns an (nelems,) shaped float array.

geomtools.inertialDirections(x)[source]¶

Return the directions and dimension of a Coords based of inertia.

x: a Coords-like array

Returns a tuple of the principal direction vectors and the sizes along these directions, ordered from the smallest to the largest direction.

geomtools.smallestDirection(x, method='inertia', return_size=False)[source]¶

Return the direction of the smallest dimension of a Coords.

x: a Coords-like array
method: one of ‘inertia’ or ‘random’
return_size: if True and method is ‘inertia’, a tuple of a direction vector and the size along that direction and the cross directions; else, only return the direction vector.

geomtools.largestDirection(x, return_size=False)[source]¶

Return the direction of the largest dimension of a Coords.

x: a Coords-like array
return_size: if True and method is ‘inertia’, a tuple of a direction vector and the size along that direction and the cross directions; else, only return the direction vector.

geomtools.distance(X, Y)[source]¶

Return the distance of all points of X to those of Y.

Parameters:

X: (nX,3) shaped array of points.
Y: (nY,3) shaped array of points.

Returns an (nX,nT) shaped array with the distances between all points of X and Y.

geomtools.closest(X, Y=None, return_dist=False)[source]¶

Find the point of Y closest to each of the points of X.

Parameters:

X: (nX,3) shaped array of points
Y: (nY,3) shaped array of points. If None, Y is taken equal to X, allowing to search for the closest point in a single set. In the latter case, the point itself is excluded from the search (as otherwise that would obviously be the closest one).
return_dist: bool. If True, also returns the distances of the closest points.

Returns:

ind: (nX,) int array with the index of the closest point in Y to the points of X
dist: (nX,) float array with the distance of the closest point. This is equal to length(X-Y[ind]). It is only returned if return_dist is True.

geomtools.closestPair(X, Y)[source]¶

Find the closest pair of points from X and Y.

Parameters:

X: (nX,3) shaped array of points
Y: (nY,3) shaped array of points

Returns a tuple (i,j,d) where i,j are the indices in X,Y identifying the closest points, and d is the distance between them.

geomtools.projectedArea(x, dir)[source]¶

Compute projected area inside a polygon.

Parameters:

x: (npoints,3) Coords with the ordered vertices of a (possibly nonplanar) polygonal contour.
dir: either a global axis number (0, 1 or 2) or a direction vector consisting of 3 floats, specifying the projection direction.

Returns a single float value with the area inside the polygon projected in the specified direction.

Note that if the polygon is planar and the specified direction is that of the normal on its plane, the returned area is that of the planar figure inside the polygon. If the polygon is nonplanar however, the area inside the polygon is not defined. The projected area in a specified direction is, since the projected polygon is a planar one.

geomtools.polygonNormals(x)[source]¶

Compute normals in all points of polygons in x.

x is an (nel,nplex,3) coordinate array representing nel (possibly nonplanar) polygons.

The return value is an (nel,nplex,3) array with the unit normals on the two edges ending in each point.

geomtools.averageNormals(coords, elems, atNodes=False, treshold=None)[source]¶

Compute average normals at all points of elems.

coords is a (ncoords,3) array of nodal coordinates. elems is an (nel,nplex) array of element connectivity.

The default return value is an (nel,nplex,3) array with the averaged unit normals in all points of all elements. If atNodes == True, a more compact array with the unique averages at the nodes is returned.

geomtools.triangleInCircle(x)[source]¶

Compute the incircles of the triangles x.

The incircle of a triangle is the largest circle that can be inscribed in the triangle.

x is a Coords array with shape (ntri,3,3) representing ntri triangles.

Returns a tuple r,C,n with the radii, Center and unit normals of the incircles.

Example:

>>> X = Formex(Coords([1.,0.,0.])).rosette(3,120.)
>>> print(X)
{[1.0,0.0,0.0], [-0.5,0.866025,0.0], [-0.5,-0.866025,0.0]}
>>> radius, center, normal = triangleInCircle(X.coords.reshape(-1,3,3))
>>> print(radius)
[ 0.5]
>>> print(center)
[[ 0. 0. 0.]]

geomtools.triangleCircumCircle(x, bounding=False)[source]¶

Compute the circumcircles of the triangles x.

x is a Coords array with shape (ntri,3,3) representing ntri triangles.

Returns a tuple r,C,n with the radii, Center and unit normals of the circles going through the vertices of each triangle.

If bounding=True, this returns the triangle bounding circle.

geomtools.triangleBoundingCircle(x)[source]¶

Compute the bounding circles of the triangles x.

The bounding circle is the smallest circle in the plane of the triangle such that all vertices of the triangle are on or inside the circle. If the triangle is acute, this is equivalent to the triangle’s circumcircle. It the triangle is obtuse, the longest edge is the diameter of the bounding circle.

x is a Coords array with shape (ntri,3,3) representing ntri triangles.

Returns a tuple r,C,n with the radii, Center and unit normals of the bounding circles.

geomtools.triangleObtuse(x)[source]¶

Check for obtuse triangles.

x is a Coords array with shape (ntri,3,3) representing ntri triangles.

Returns an (ntri) array of True/False values indicating whether the triangles are obtuse.

geomtools.lineIntersection(P0, N0, P1, N1)[source]¶

Find the intersection of 2 (sets of) lines.

This relies on the lines being pairwise coplanar.

geomtools.displaceLines(A, N, C, d)[source]¶

Move all lines (A,N) over a distance a in the direction of point C.

A,N are arrays with points and directions defining the lines. C is a point. d is a scalar or a list of scalars. All line elements of F are translated in the plane (line,C) over a distance d in the direction of the point C. Returns a new set of lines (A,N).

geomtools.segmentOrientation(vertices, vertices2=None, point=None)[source]¶

Determine the orientation of a set of line segments.

vertices and vertices2 are matching sets of points. point is a single point. All arguments are Coords objects.

Line segments run between corresponding points of vertices and vertices2. If vertices2 is None, it is obtained by rolling the vertices one position foreward, thus corresponding to a closed polygon through the vertices). If point is None, it is taken as the center of vertices.

The orientation algorithm checks whether the line segments turn positively around the point.

Returns an array with +1/-1 for positive/negative oriented segments.

geomtools.rotationAngle(A, B, m=None, angle_spec=0.017453292519943295)[source]¶

Return rotation angles and vectors for rotations of A to B.

A and B are (n,3) shaped arrays where each line represents a vector. This function computes the rotation from each vector of A to the corresponding vector of B. If m is None, the return value is a tuple of an (n,) shaped array with rotation angles (by default in degrees) and an (n,3) shaped array with unit vectors along the rotation axis. If m is a (n,3) shaped array with vectors along the rotation axis, the return value is a (n,) shaped array with rotation angles. The returned angle is then the angle between the planes formed by the axis and the vectors. Specify angle_spec=RAD to get the angles in radians.

geomtools.anyPerpendicularVector(A)[source]¶

Return arbitrary vectors perpendicular to vectors of A.

A is a (n,3) shaped array of vectors. The return value is a (n,3) shaped array of perpendicular vectors.

The returned vector is always a vector in the x,y plane. If the original is the z-axis, the result is the x-axis.

geomtools.perpendicularVector(A, B)[source]¶: Return vectors perpendicular on both A and B.

geomtools.projectionVOV(A, B)[source]¶: Return the projection of vector of A on vector of B.

geomtools.projectionVOP(A, n)[source]¶: Return the projection of vector of A on a plane with normal n.

geomtools.baryCoords(S, P)[source]¶

Compute the barycentric coordinates of points wrt. simplexes.

An n-simplex is a geometrical structure defined by n+1 vertices and bordered by the convex hull of those vertices. In practice it is either:

1-simplex: line segment (nplex=2)
2-simplex: triangle (nplex=3)
3-simplex: tetrahedron (nplex=4)

The barycentric coordinates of a 3d point with respect to a simplex of a lower order are the barycentric coordinates of the projection of that point on the simplex.

Parameters:	S (coords_like (nel, nplex, 3)) – A set of nel n-simplexes (n=nplex-1). P (coords_like (nel, npts, 3) or (1, npts, 3)) – A set of npts points (for each of the simplexes in S) for which the barycentric coordinates are to be computed. If the shape is (1,npts,3), the same points are used with each of the simplexes.
Returns:	float array (nel, npts, nplex) – The barycentric coordinates of the points with respect to the simplexes.

See also

insideSimplex(): test if a (projection) point falls within a simplex

Examples

>>> S = Coords('.1.6.4').reshape(3,2,3)
>>> P = Coords([[[0,0,0],[0.2,0.2,0],[0.5,0.5,0],[0.5,0.7,0]]])
>>> baryCoords(S,P)
array([[[ 1. ,  0. ],
        [ 0.8,  0.2],
        [ 0.5,  0.5],
        [ 0.5,  0.5]],
<BLANKLINE>
       [[ 0.5,  0.5],
        [ 0.5,  0.5],
        [ 0.5,  0.5],
        [ 0.4,  0.6]],
<BLANKLINE>
       [[ 0. ,  1. ],
        [ 0.2,  0.8],
        [ 0.5,  0.5],
        [ 0.7,  0.3]]])
>>> S1 = Coords('016').reshape(1,3,3)
>>> baryCoords(S1,P)
array([[[ 1. ,  0. ,  0. ],
        [ 0.6,  0.2,  0.2],
        [ 0. ,  0.5,  0.5],
        [-0.2,  0.5,  0.7]]])

geomtools.insideSimplex(S, P, atol=0.0)[source]¶

Check which points P are inside the simplexes S.

An n-simplex is a geometrical structure defined by n+1 vertices and bordered by the convex hull of those vertices. In practice it is either:

1-simplex: line segment (nplex=2)
2-simplex: triangle (nplex=3)
3-simplex: tetrahedron (nplex=4)

A 3d point is ‘inside’ a simplex of a lower order if its projection on that simplex falls within the simplex. This is satisfied iff all the barycentric coordinates of the point with respect to the simplex are in the range [0.0, 1.0].

Parameters:

S (coords_like (nel, nplex, 3)) – A set of nel n-simplexes (n=nplex-1).
P (coords_like (nel, npts, 3) or (1, npts, 3)) – A set of npts points that have to be tested against the nel simplexes. If the shape is (1,npts,3), the same points are tested against each of the simplexes. With a shape (nel, npts, 3) each simplex has its own set of points to be tested.
atol (float, optional) – If provided, points which are falling outside the simplex within this tolerance (in parametric space), are also reported as inside.

Returns:

bool array (nel, npts) – The array holds True where the (projection of the) points fall inside the simplex, to the accuracy atol.

See also

baryCoords(): compute the barycentric coordiantes of the points

Examples

>>> S = Coords('.1.6.4').reshape(3,2,3)
>>> P = Coords([[[0.2,0.2,0],[0.5,0.5,0],[0.5,0.7,0],[0.5,1.2,0]]])
>>> insideSimplex(S,P)
array([[ True,  True,  True,  True],
       [ True,  True,  True,  True],
       [ True,  True,  True, False]])
>>> T = Coords('132.14').reshape(2,3,3)
>>> insideSimplex(T,P)
array([[ True,  True, False, False],
       [False,  True,  True, False]])

geomtools.insideTriangle(S, P, atol=0.0)¶

Check which points P are inside the simplexes S.

An n-simplex is a geometrical structure defined by n+1 vertices and bordered by the convex hull of those vertices. In practice it is either:

1-simplex: line segment (nplex=2)
2-simplex: triangle (nplex=3)
3-simplex: tetrahedron (nplex=4)

A 3d point is ‘inside’ a simplex of a lower order if its projection on that simplex falls within the simplex. This is satisfied iff all the barycentric coordinates of the point with respect to the simplex are in the range [0.0, 1.0].

Parameters:

S (coords_like (nel, nplex, 3)) – A set of nel n-simplexes (n=nplex-1).
P (coords_like (nel, npts, 3) or (1, npts, 3)) – A set of npts points that have to be tested against the nel simplexes. If the shape is (1,npts,3), the same points are tested against each of the simplexes. With a shape (nel, npts, 3) each simplex has its own set of points to be tested.
atol (float, optional) – If provided, points which are falling outside the simplex within this tolerance (in parametric space), are also reported as inside.

Returns:

bool array (nel, npts) – The array holds True where the (projection of the) points fall inside the simplex, to the accuracy atol.

See also

baryCoords(): compute the barycentric coordiantes of the points

Examples

>>> S = Coords('.1.6.4').reshape(3,2,3)
>>> P = Coords([[[0.2,0.2,0],[0.5,0.5,0],[0.5,0.7,0],[0.5,1.2,0]]])
>>> insideSimplex(S,P)
array([[ True,  True,  True,  True],
       [ True,  True,  True,  True],
       [ True,  True,  True, False]])
>>> T = Coords('132.14').reshape(2,3,3)
>>> insideSimplex(T,P)
array([[ True,  True, False, False],
       [False,  True,  True, False]])

geomtools.insideSegment(Q, Q1, P, atol=0.0)[source]¶

Check if projections of points are inside line segments.

Checks which projections of points P on the line segments Q-Q1 are within the line segments.

Parameters:

Q (array_like (nseg, 3)) – First point of nseg line segments.
Q1 (array_like (nseg, 3)) – Second point of nseg line segments.
P (array_like (nseg, 3)) – The points to test against the segments. The testing is done by pair: one point for each segment. The test pertains to the projection of the point on the line containing the segment.
atol (float, optional) – If provided, points which are falling outside the segment but within this tolerance are also reported as inside.

Returns:

bool array (nseg) – The array holds True where the (projection of the) point falls inside the corresponding segemnt, within the accuracy atol.

Notes

On large data sets this is (slightly) more efficient than the equivalent:

S = np.stack([Q, Q1], axis=1)
P = P[:,np.newaxis]
t = geomtools.insideSimplex(S,P).reshape(-1)

See also

insideSimplex(): a general test of (projection of) points inside simplexes
insideRay(): a similar test for points on a half-line (ray)

Examples

>>> nseg = 13
>>> Q = Coords(np.random.rand(nseg,3))
>>> Q1 = Coords(np.random.rand(nseg,3))
>>> u = at.uniformParamValues(nseg-1, -0.1, 1.1).reshape(-1,1,1)
>>> P = Q[:,np.newaxis] * (1.-u) + Q1[:, np.newaxis] * u
>>> P = P.reshape(nseg,3)
>>> insideSegment(Q, Q1, P, atol=1.e-5)
array([False,  True,  True,  True,  True,  True,  True,  True,  True,
        True,  True,  True, False])
>>> insideRay(Q, Q1-Q, P, atol=1.e-5)
array([False,  True,  True,  True,  True,  True,  True,  True,  True,
        True,  True,  True,  True])

geomtools.insideRay(Q, v, P, atol=0.0)[source]¶

Check if projections of points are inside rays (half-lines).

Checks which projections of points P on the line (Q,v) are within the positive half-lines (rays).

Parameters:

Q (array_like (nlines, 3)) – Starting point of nlines rays (half-lines).
v (array_like (nlines, 3)) – Direction vectors of nlines rays (half-lines).
P (array_like (nlines, 3)) – The points to test against the rays. The testing is done by pair: one point for each segment. The test pertains to the projection of the point on the line containing the ray.
atol (float, optional) – If provided, points which are falling outside the ray but within this tolerance are also reported as inside.

Returns:

bool array (nseg) – The array holds True where the (projection of the) point falls inside the corresponding ray, within the accuracy atol.

Notes

On large data sets this is (slightly) more efficient than the equivalent:

S = np.stack([Q, Q+v], axis=1)
P = P[:,np.newaxis]
B = geomtools.baryCoords(S,P).reshape(-1,2)
t = where(B[:,1] >= atol)

See also

insideSimplex(): a general test of (projection of) points inside simplexes
insideRay(): a similar test for points on a half-line (ray)

Table of Contents

Previous topic

Next topic

20. `geomtools` — Basic geometrical operations.¶

20.1. Classes defined in module geomtools¶

20.2. Functions defined in module geomtools¶

20. geomtools — Basic geometrical operations.¶

20.1. Classes defined in module geomtools¶

20.2. Functions defined in module geomtools¶

20. `geomtools` — Basic geometrical operations.¶