mathutils -- All the basic math types and functions

Most of the names here are coming from the glm module. But the goal is to harvest at one point all the basic math functions and objects that are used all around madcad.

Tip

All the present functions and types are present in the root madcad module.

Most common glm types

All following objects implement the common operators + - * / <> ==

vec3 = dvec3 module-attribute

mat3 = dmat3 module-attribute

mat4 = dmat4 module-attribute

quat = dquat module-attribute

All glm types exists with several element types and in several precision:

Prefix	Precision
d	f64 aka double precision floating point
f	f32 aka simple precision floating point
i	i32 aka integer
i16	16 bits integer
i8	8 bits integer aka byte
u	unsigned integer (also declines in u16 and u32)
b	bit aka boolean

• Precision specification is put as prefix: dvec3, fvec3, imat4. Notation without prefix refers to the madcad implementation precision: float64 (prefix 'd').
• Object dimension is put as suffix.

In this documentation, when we refer to a 'vector' without explicit type, we obviously mean a vec3 aka. dvec3.

Note

The default glm_ precision type is float32 (prefix 'f'). For convenience, these are overriden in madcad to use float64 for a better precision.

Common vector operations

dot() builtin

dot(x: number, y: number) -> float Returns the dot product of x and y, i.e., result = x * y. dot(x: vecN, y: vecN) -> number Returns the dot product of x and y, i.e., result = x[0] * y[0] + x[1] * y[1] + ... dot(x: quat, y: quat) -> float Returns dot product of x and y, i.e., x[0] * y[0] + x[1] * y[1] + ...

cross() builtin

cross(x: vec3, y: vec3) -> vec3 Returns the cross product of x and y. cross(x: quat, y: quat) -> quat Compute a cross product.

length() builtin

length(x: float) -> float Returns the length of x, i.e., abs(x). length(x: vecN) -> float Returns the length of x, i.e., sqrt(x * x). length(x: quat) -> float Returns the norm of a quaternion.

distance() builtin

distance(p0: float, p1: float) -> float Returns the distance between p0 and p1, i.e., length(p0 - p1). distance(p0: vecN, p1: vecN) -> float Returns the distance between p0 and p1, i.e., length(p0 - p1).

anglebt(x, y)

Angle between two vectors

The result is not sensitive to the lengths of x and y

Source code in madcad/mathutils.py

def anglebt(x,y) -> float:
	''' Angle between two vectors

		![anglebt](../schemes/mathutils-anglebt.svg)

		The result is not sensitive to the lengths of x and y
	'''
	n = length(x)*length(y)
	return acos(min(1,max(-1, dot(x,y)/n)))	if n else 0

arclength(p1, p2, n1, n2)

Length of an arc between p1 and p2, normal to the given vectors in respective points

Source code in madcad/mathutils.py

def arclength(p1, p2, n1, n2):
	''' Length of an arc between p1 and p2, normal to the given vectors in respective points '''
	c = max(0, dot(n1,n2))
	if abs(c-1) < NUMPREC:	return 0
	v = p1-p2
	return sqrt(dot(v,v) / (2-2*c)) * acos(c)

project(vec, dir)

Component of vec along dir, equivalent to :code:dot(vec,dir) / dot(dir,dir) * dir

The result is not sensitive to the length of dir

Source code in madcad/mathutils.py

def project(vec, dir) -> vec3:
	''' Component of `vec` along `dir`, equivalent to :code:`dot(vec,dir) / dot(dir,dir) * dir`

		![project](../schemes/mathutils-project.svg)

		The result is not sensitive to the length of `dir`
	'''
	try:	return dot(vec,dir) / dot(dir,dir) * dir
	except ZeroDivisionError:	
		if dot(vec,vec):		return vec3(nan)
		else:					return vec3(0)

noproject(vec, dir)

Components of vec not along dir, equivalent to :code:vec - project(vec,dir)

The result is not sensitive to the length of dir

Source code in madcad/mathutils.py

def noproject(vec, dir) -> vec3:
	''' Components of `vec` not along `dir`, equivalent to :code:`vec - project(vec,dir)`

		![noproject](../schemes/mathutils-noproject.svg)

		The result is not sensitive to the length of `dir`
	'''
	return vec - project(vec,dir)

unproject(vec, dir)

Return the vector in the given direction as if vec was its projection on it, equivalent to :code:dot(vec,vec) / dot(vec,dir) * dir

The result is not sensitive to the length of dir

Source code in madcad/mathutils.py

def unproject(vec, dir) -> vec3:
	''' Return the vector in the given direction as if `vec` was its projection on it, equivalent to :code:`dot(vec,vec) / dot(vec,dir) * dir`

		![unproject](../schemes/mathutils-unproject.svg)

		The result is not sensitive to the length of `dir`
	'''
	try:	return dot(vec,vec) / dot(vec,dir) * dir
	except ZeroDivisionError:	
		if dot(vec,vec):		return vec3(nan)
		else:					return vec3(0)

reflect() builtin

reflect(I: float, N: float) -> float For the incident vector I and surface orientation N, returns the reflection direction: result = I - 2.0 * dot(N, I) * N. reflect(I: vecN, N: vecN) -> vecN For the incident vector I and surface orientation N, returns the reflection direction: result = I - 2.0 * dot(N, I) * N.

perp(v)

Perpendicular vector to the given vector

Source code in madcad/mathutils.py

def perp(v:vec2) -> vec2:
	''' Perpendicular vector to the given vector

		![perp](../schemes/mathutils-perp.svg)
	'''
	return vec2(-v[1], v[0])

See the glm complete reference

Transformations

transform(*args)

Create an affine transformation matrix.

Supported inputs

:mat4: obviously return it unmodified :float: scale using the given ratio :vec3: translation only :quat, mat3, mat4: rotation only :(vec3,vec3), (vec3,mat3), (vec3,quat): (o,T) translation and rotation :(vec3,vec3,vec3): (x,y,z) base of vectors for rotation :(vec3,vec3,vec3,vec3): (o,x,y,z) translation and base of vectors for rotation

Source code in madcad/mathutils.py

def transform(*args) -> mat4:
	''' Create an affine transformation matrix.

		Supported inputs:
			:mat4:                                    obviously return it unmodified
			:float:                                   scale using the given ratio 
			:vec3:                                    translation only
			:quat, mat3, mat4:                        rotation only
			:(vec3,vec3), (vec3,mat3), (vec3,quat):   `(o,T)` translation and rotation
			:(vec3,vec3,vec3):                        `(x,y,z)` base of vectors for rotation
			:(vec3,vec3,vec3,vec3):                   `(o,x,y,z)` translation and base of vectors for rotation
	'''
	if len(args) == 1 and isinstance(args[0], (tuple,list)):
		args = args[0]
	if len(args) == 1:
		if isinstance(args[0], mat4):	return args[0]
		elif isinstance(args[0], (mat3, quat, int, float)):	return mat4(args[0])
		elif isinstance(args[0], vec3):	return translate(args[0])
	elif len(args) == 2:
		if isinstance(args[0], vec3):
			if   isinstance(args[1], (mat3, quat, int, float)):		m = mat4(args[1])
			elif isinstance(args[1], vec3):		m = mat4(quat(args[1]))
			m[3] = vec4(args[0], 1)
			return m
	elif isinstance(args[0], vec3) and len(args) == 3:			
		return mat4(mat3(*args))
	elif isinstance(args[0], vec3) and len(args) == 4:
		m = mat4(mat3(*args[1:]))
		m[3] = vec4(args[0], 1)
		return m

	raise TypeError('a transformation must be a  mat3, mat4, quat, (O,mat3), (O,quat), (0,x,y,z), not {}'.format(args))

translate() builtin

translate(v: vec3) -> mat4x4 Builds a translation 4 x 4 matrix created from a vector of 3 components. translate(v: vec2) -> mat3x3 Builds a translation 3 x 3 matrix created from a vector of 2 components. translate(m: mat4x4, v: vec3) -> mat4x4 Builds a translation 4 x 4 matrix created from a vector of 3 components. m is the input matrix multiplied by this translation matrix translate(m: mat3x3, v: vec2) -> mat3x3 Builds a translation 3 x 3 matrix created from a vector of 2 components. m is the input matrix multiplied by this translation matrix

rotate() builtin

rotate(angle: number, axis: vec3) -> mat4x4 Builds a rotation 4 x 4 matrix created from an axis vector and an angle. rotate(angle: number) -> mat3x3 Builds a rotation 3 x 3 matrix created from an angle. rotate(m: mat4x4, angle: number, axis: vec3) -> mat4x4 Builds a rotation 4 x 4 matrix created from an axis vector and an angle. m is the input matrix multiplied by this translation matrix rotate(m: mat3x3, angle: number) -> mat3x3 Builds a rotation 3 x 3 matrix created from an angle. m is the input matrix multiplied by this translation matrix rotate(v: vec2, angle: float) -> vec2 Rotate a two dimensional vector. rotate(v: vec3, angle: float, normal: vec3) -> vec3 Rotate a three dimensional vector around an axis. rotate(v: vec4, angle: float, normal: vec3) -> vec4 Rotate a four dimensional vector around an axis. rotate(q: quat, angle: float, axis: vec3) -> quat Rotates a quaternion from a vector of 3 components axis and an angle.

scaledir(dir, factor=None)

Return a mat3 scaling in the given direction, with the given factor (1 means original scale)

If factor is None, the length of dir is used, but it can leads to precision loss on direction when too small.

Source code in madcad/mathutils.py

def scaledir(dir, factor=None) -> mat3:
	''' Return a mat3 scaling in the given direction, with the given factor (1 means original scale)

		![scaledir](../schemes/mathutils-scaledir.svg)

		If factor is None, the length of dir is used, but it can leads to precision loss on direction when too small.
	'''
	if factor is None:
		factor = length(dir)
		dir = dir / factor
	return mat3(1) + (factor-1)*mat3(dir[0]*dir, dir[1]*dir, dir[2]*dir)

scale() builtin

scale(v: vec3) -> mat4x4 Builds a scale 4 x 4 matrix created from 3 scalars. scale(v: vec2) -> mat3x3 Builds a scale 3 x 3 matrix created from a vector of 2 components. scale(m: mat4x4, v: vec3) -> mat4x4 Builds a scale 4 x 4 matrix created from 3 scalars. m is the input matrix multiplied by this translation matrix scale(m: mat3x3, v: vec2) -> mat3x3 Builds a scale 3 x 3 matrix created from a vector of 2 components. m is the input matrix multiplied by this translation matrix

rotatearound(angle, *args)

Return a transformation matrix for a rotation around an axis

rotatearound(angle, axis) rotatearound(angle, origin, dir)

Source code in madcad/mathutils.py

def rotatearound(angle, *args) -> mat4:
	''' Return a transformation matrix for a rotation around an axis

		rotatearound(angle, axis)
		rotatearound(angle, origin, dir)
	'''
	if len(args) == 1:		origin, dir = args[0]
	elif len(args) == 2:	origin, dir = args
	else:
		raise TypeError('invalid use of rotatearound')

	r = mat3_cast(angleAxis(angle, dir))
	m = mat4(r)
	m[3] = vec4(origin - r*origin, 1)
	return m

dirbase(dir, align=vec3(1, 0, 0))

Return a base using the given direction as z axis (and the nearer vector to align as x)

Source code in madcad/mathutils.py

def dirbase(dir, align=vec3(1,0,0)):
	''' Return a base using the given direction as z axis (and the nearer vector to align as x) '''
	x = noproject(align, dir)
	if not length2(x) > NUMPREC**2:
		align = vec3(align[2],-align[0],align[1])
		x = noproject(align, dir)
	if not length2(x) > NUMPREC**2:
		align = vec3(align[1],-align[2],align[0])
		x = noproject(align, dir)
	x = normalize(x)
	y = cross(dir, x)
	return x,y,dir

transpose() builtin

transpose(x: matNxM) -> matMxN Returns the transposed matrix of x.

inverse() builtin

inverse(m: matSxS) -> matSxS Return the inverse of a squared matrix. inverse(q: quat) -> quat Return the inverse of a quaternion.

affineInverse() builtin

affineInverse(m: matSxS) -> matSxS Fast matrix inverse for affine matrix.

angleAxis() builtin

angleAxis(angle: float, axis: vec3) -> quat Build a quaternion from an angle and a normalized axis.

Interpolation functions

mix() builtin

mix(x: number, y: number, a: float) -> number Returns x * (1.0 - a) + y * a, i.e., the linear blend of x and y using the floating-point value a. The value for a is not restricted to the range [0, 1]. mix(x: number, y: number, a: bool) -> number Returns y if a is True and x otherwise. mix(x: vecN, y: vecN, a: fvecN) -> vecN Returns x * (1.0 - a) + y * a, i.e., the linear blend of x and y using the floating-point value a. The value for a is not restricted to the range [0, 1]. mix(x: vecN, y: vecN, a: bvecN) -> vecN For each component index i: Returns y[i] if a[i] is True and x[i] otherwise. mix(x: matNxM, y: matNxM, a: fmatNxM) -> matNxM Returns x * (1.0 - a) + y * a, i.e., the linear blend of x and y using the floating-point value a for each component. The value for a is not restricted to the range [0, 1]. mix(x: matNxM, y: matNxM, a: float) -> matNxM Returns x * (1.0 - a) + y * a, i.e., the linear blend of x and y using the floating-point value a for each component. The value for a is not restricted to the range [0, 1]. mix(x: quat, y: quat, a: float) -> quat Spherical linear interpolation of two quaternions. The interpolation is oriented and the rotation is performed at constant speed. For short path spherical linear interpolation, use the slerp function.

hermite = interpol2 module-attribute

step() builtin

step(edge: number, x: number) -> float Returns 0.0 if x < edge, otherwise it returns 1.0. step(edge: number, x: vecN) -> vecN For every component c of x: Returns 0.0 if c < edge, otherwise it returns 1.0. step(edge: vecN, x: vecN) -> vecN For every index i: Returns 0.0 if x[i] < edge[i], otherwise it returns 1.0.

linstep(start, stop, x)

like smoothstep but with a linear ramp between start and stop

Source code in madcad/mathutils.py

def linstep(start, stop, x):
	''' like smoothstep but with a linear ramp between `start` and `stop`

		![linstep](../schemes/mathutils-step.svg)
	'''
	if x <= start:	return 0
	if x >= stop:	return 1
	return (x-start)/(stop-start)

smoothstep() builtin

smoothstep(edge0: number, edge1: number, x: number) -> float Returns 0.0 if x <= edge0 and 1.0 if x >= edge1 and performs smooth Hermite interpolation between 0 and 1 when edge0 < x < edge1. This is useful in cases where you would want a threshold function with a smooth transition. This is equivalent to : t = clamp((x - edge0) / (edge1 - edge0), 0, 1) return t * t * (3 - 2 * t) Results are undefined if edge0 >= edge1. smoothstep(edge0: number, edge1: number, x: vecN) -> vecN Returns 0.0 if x <= edge0 and 1.0 if x >= edge1 and performs smooth Hermite interpolation between 0 and 1 when edge0 < x < edge1. This is useful in cases where you would want a threshold function with a smooth transition. This is equivalent to : t = clamp((x - edge0) / (edge1 - edge0), 0, 1) return t * t * (3 - 2 * t) Results are undefined if edge0 >= edge1. smoothstep(edge0: vecN, edge1: vecN, x: vecN) -> vecN Returns 0.0 if x <= edge0 and 1.0 if x >= edge1 and performs smooth Hermite interpolation between 0 and 1 when edge0 < x < edge1. This is useful in cases where you would want a threshold function with a smooth transition. This is equivalent to : t = clamp((x - edge0) / (edge1 - edge0), 0, 1) return t * t * (3 - 2 * t) Results are undefined if edge0 >= edge1.

clamp() builtin

clamp(x: number, minVal: number, maxVal: number) -> number Returns min(max(x, minVal), maxVal). clamp(x: vecN, minVal: number, maxVal: number) -> vecN Returns min(max(x, minVal), maxVal) for each component in x using the floating-point values minVal and maxVal. clamp(x: vecN, minVal: vecN, maxVal: vecN) -> vecN Returns min(max(x, minVal), maxVal) for each component in x using the floating-point values minVal and maxVal.

lerp() builtin

lerp(x: float, y: float, a: float) -> float Returns x * (1.0 - a) + y * a, i.e., the linear blend of x and y using the floating-point value a. The value for a is not restricted to the range [0, 1]. lerp(x: vecN, y: vecN, a: float) -> vecN Returns x * (1.0 - a) + y * a, i.e., the linear blend of x and y using the floating-point value a. The value for a is not restricted to the range [0, 1]. lerp(x: vecN, y: vecN, a: vecN) -> vecN Returns x * (1.0 - a) + y * a, i.e., the linear blend of x and y using the vector a. The value for a is not restricted to the range [0, 1]. lerp(x: quat, y: quat, a: float) -> quat Linear interpolation of two quaternions. The interpolation is oriented.

slerp() builtin

slerp(x: quat, y: quat, a: float) -> quat Spherical linear interpolation of two quaternions. The interpolation always take the short path and the rotation is performed at constant speed. slerp(x: vec3, y: vec3, a: float) -> vec3 Returns Spherical interpolation between two vectors.

linrange(start, stop=None, step=None, div=0, end=True)

Yield successive intermediate values between start and stop

    stepping:

    - if `step` is given, it will be the amount between raised value until it gets over `stop`
    - if `div` is given, it will be the number of intermediate steps between `start` and `stop` (with linear spacing)

    ending:

    - if `end` is True, it will stop iterate with value `stop` (or just before)
    - if `end` is False, it will stop iterating just before `stop` and never with `stop`

Examples:

>>> list(linrange(5, -5, div=1))
[5, 0, -5]

>>> list(linrange(5, -5, div=10)

NOTE

If step is given and is not a multiple of stop-start then end has no influence

Source code in madcad/mathutils.py

def linrange(start, stop=None, step=None, div=0, end=True):
	''' Yield successive intermediate values between start and stop 

		stepping:

		- if `step` is given, it will be the amount between raised value until it gets over `stop`
		- if `div` is given, it will be the number of intermediate steps between `start` and `stop` (with linear spacing)

		ending:

		- if `end` is True, it will stop iterate with value `stop` (or just before)
		- if `end` is False, it will stop iterating just before `stop` and never with `stop`

	Examples:
		>>> list(linrange(5, -5, div=1))
		[5, 0, -5]

		>>> list(linrange(5, -5, div=10)


	NOTE:  
		If step is given and is not a multiple of `stop-start` then `end` has no influence
	'''
	if stop is None:	start, stop = 0, start
	if step is None:	step = (stop-start)/(div+1)
	elif step * (stop-start) < 0:	step = -step
	if not end:			stop -= step
	stop += NUMPREC*stop

	t = start
	while (stop-t)*step >= 0:
		yield t
		t += step

intri_smooth(pts, ptangents, a, b)

Cubic interpolation over a triangle, edges are guaranteed to fit an interpol2 curve using the edge tangents

Note

If the tangents lengths are set to the edge lengths, that version gives a result that only blends between the curved edges, a less bulky result than intri_sphere

Source code in madcad/mathutils.py

def intri_smooth(pts, ptangents, a,b):
	''' Cubic interpolation over a triangle, edges are guaranteed to fit an interpol2 curve using the edge tangents

	Note:
		If the tangents lengths are set to the edge lengths, that version gives a result that only blends between the curved edges, a less bulky result than `intri_sphere`
	'''
	A,B,C = pts
	ta,tb,tc = ptangents
	c = 1-a-b
	return (	0
			+	a**2 * (A + b*ta[0] + c*ta[1] + b*c*(ta[0]+ta[1]))
			+	b**2 * (B + c*tb[0] + a*tb[1] + c*a*(tb[0]+tb[1]))
			+	c**2 * (C + a*tc[0] + b*tc[1] + a*b*(tc[0]+tc[1]))
			+	2*(b*c + c*a + a*b) * (a*A + b*B + c*C)
			)

intri_sphere(pts, ptangents, a, b, etangents=None)

Cubic interpolation over a triangle (2 dimension space), edges are guaranteed to fit an interpol2 curve using the edge tangents

Note

If the tangents lengths are set to the edge lengths, that version gives a result close to a sphere surface

Source code in madcad/mathutils.py

def intri_sphere(pts, ptangents, a,b, etangents=None):
	''' Cubic interpolation over a triangle (2 dimension space), edges are guaranteed to fit an interpol2 curve using the edge tangents

	Note:
		If the tangents lengths are set to the edge lengths, that version gives a result close to a sphere surface
	'''
	A,B,C = pts
	ta,tb,tc = ptangents
	c = 1-a-b
	P = a*A + b*B + c*C +  a*b*c * (ta[0] + ta[1] + tb[0] + tb[1] + tc[0] + tc[1])
	return (	0
			+	a**2 * (A + b*ta[0] + c*ta[1])
			+	b**2 * (B + c*tb[0] + a*tb[1])
			+	c**2 * (C + a*tc[0] + b*tc[1])
			+	2*(b*c + c*a + a*b) * P
			)

intri_parabolic(pts, ptangents, a, b, etangents=None)

Quadratic interpolation over a triangle, edges are NOT fitting an interpol2 curve

Source code in madcad/mathutils.py

def intri_parabolic(pts, ptangents, a,b, etangents=None):
	''' Quadratic interpolation over a triangle, edges are NOT fitting an interpol2 curve '''
	A,B,C = pts
	ta,tb,tc = ptangents
	c = 1-a-b
	return (	0
			+	a * (A + b*ta[0] + c*ta[1])
			+	b * (B + c*tb[0] + a*tb[1])
			+	c * (C + a*tc[0] + b*tc[1])
			)

Distances

distance_pa(pt, axis)

Point - axis distance

Source code in madcad/mathutils.py

def distance_pa(pt, axis):
	''' Point - axis distance '''
	return length(noproject(pt-axis[0], axis[1]))

distance_pe(pt, edge)

Point - edge distance

Source code in madcad/mathutils.py

def distance_pe(pt, edge):
	''' Point - edge distance '''
	dir = edge[1]-edge[0]
	l = length2(dir)
	if not l:	return 0
	x = dot(pt-edge[0], dir)/l
	if   x < 0:	return distance(pt,edge[0])
	elif x > 1:	return distance(pt,edge[1])
	else:
		return length(noproject(pt-edge[0], dir))

distance_aa(a1, a2)

Axis - axis distance

Source code in madcad/mathutils.py

def distance_aa(a1, a2):
	''' Axis - axis distance '''
	return length(project(a1[0]-a2[0], cross(a1[1], a2[1])))

distance_ae(axis, edge)

Axis - edge distance

Source code in madcad/mathutils.py

def distance_ae(axis, edge):
	''' Axis - edge distance '''
	x = axis[1]
	z = normalize(cross(x, edge[1]-edge[0]))
	y = cross(z,x)
	s1 = dot(edge[0]-axis[0], y)
	s2 = dot(edge[1]-axis[0], y)
	if s1*s2 < 0:
		return dot(edge[0]-axis[1], z)
	elif abs(s1) < abs(s2):
		return distance_pa(edge[0], axis)
	else:
		return distance_pa(edge[1], axis)

Locally defined data types

Axis(origin, direction=None, interval=None)

Bases: object

A 3D (zeroed) axis with an origin and a direction

    ![axis](../screenshots/primitives-axis.png)

    Mathematically speaking, a 3D axis doesn't necessarily have an origin, since any point on it can be its start, but for implementation and convenience reasons this axis has

Note:

    in previous madcad versions, axis were often tuples and not instances of this class. This is why this class has a `__getitem__` allowing to be used like a tuple. But this class should be used instead now.

Source code in madcad/mathutils.py

def __init__(self, origin, direction=None, interval=None):
	if direction is None:
		origin, direction = vec3(0), origin
	self.origin, self.direction = origin, direction
	self.interval = interval

__slots__ = ('origin', 'direction', 'interval') class-attribute instance-attribute

interval = interval instance-attribute

slvvars = ('origin', 'direction') class-attribute instance-attribute

__getitem__(i)

behave like the axis was a tuple (origin, direction)

Source code in madcad/mathutils.py

def __getitem__(self, i):
	''' behave like the axis was a tuple (origin, direction) '''
	if i==0:	return self.origin
	elif i==1:	return self.direction
	else:		raise IndexError('an axis has only 2 components')

flip()

switch the axis direction

Source code in madcad/mathutils.py

def flip(self) -> 'Axis':
	''' switch the axis direction '''
	return Axis(self.origin, -self.direction, self.interval)

offset(increment)

move the axis origin along its direction

Source code in madcad/mathutils.py

def offset(self, increment) -> 'Axis':
	''' move the axis origin along its direction '''
	return Axis(self.origin + self.direction*increment, self.direction, self.interval)

transform(transform)

move the axis by the given transformation

Source code in madcad/mathutils.py

def transform(self, transform) -> 'Axis':
	''' move the axis by the given transformation '''
	if isinstance(transform, (float,int)):		return Axis(transform*self.origin, self.direction, self.interval)
	elif isinstance(transform, vec3):			return Axis(transform+self.origin, self.direction, self.interval)
	elif isinstance(transform, (mat3, quat)):	return Axis(transform*self.origin, normalize(transform*self.direction), self.interval)
	elif isinstance(transform, (mat4)):			return Axis(transform*self.origin, normalize(mat3(transform)*self.direction), self.interval)
	raise TypeError('transform must be one of float, vec3, mat3, quat, mat4')

slv_tangent(pt)

Source code in madcad/mathutils.py

def slv_tangent(self, pt):
	return self.direction

__eq__(other)

Source code in madcad/mathutils.py

def __eq__(self, other):
	return self is other or isinstance(other, Axis) and (
		self.origin == other.origin and self.direction == other.direction
		)

__repr__()

Source code in madcad/mathutils.py

def __repr__(self):
	return 'Axis({}, {})'.format(self.origin, self.direction)

display(scene)

Source code in madcad/mathutils.py

def display(self, scene):
	from .rendering.d3.marker import AxisDisplay
	return AxisDisplay(scene, (self.origin, self.direction), self.interval)

isaxis(obj)

Return True if the given object is considered to be an axis. An axis can be an instance of Axis or a tuple (vec3, vec3)

Source code in madcad/mathutils.py

def isaxis(obj):
	''' Return True if the given object is considered to be an axis.
		An axis can be an instance of `Axis` or a tuple `(vec3, vec3)`
	'''
	return isinstance(obj, Axis) or isinstance(obj, tuple) and len(obj)==2 and isinstance(obj[0],vec3) and isinstance(obj[1],vec3)

Box(min=None, max=None, center=vec3(0), size=vec3(-inf), width=None, alignment=0.5)

This class describes a box always orthogonal to the base axis, used as convex for area delimitations

box can be created either by poviding:

- min and max corners
- center position and size vector
- center only (box will be empty)
- size only (box will be centered around corner)

This class is independent from the dimension or number precision of the used vectors. You can for instance have a Box of vec2 as well as a box of vec3. However boxes with different vector types cannot interperate.

Attributes:

min:        vector of minimum coordinates of the box (usually bottom left corner)
max:        vector of maximum coordinates of the box (usually top right corner)

Source code in madcad/box.py

def __init__(self, min=None, max=None, center=vec3(0), size=vec3(-inf), width=None, alignment:vec3|float=0.5):
    # DEPRECATED
    if width is not None:
        size = width
    if min and max:			self.min, self.max = min, max
    else:					self.min, self.max = center-alignment*size, center+(1-alignment)*size

__slots__ = ('min', 'max') class-attribute instance-attribute

center property

Mid coordinates of the box

width property

Diagonal vector of the box deprecated

size property

Diagonal vector of the box

__ior__ = merge_update class-attribute instance-attribute

__iand__ = intersection_update class-attribute instance-attribute

__or__ = merge class-attribute instance-attribute

__and__ = intersection class-attribute instance-attribute

from_iter(it) staticmethod

create a box bounding the positions given by the input iterable

Source code in madcad/box.py

@staticmethod
def from_iter(it) -> Box:
    ''' create a box bounding the positions given by the input iterable '''
    it = iter(it)
    first = next(it)
    new = Box(first, first)
    for p in it:
        new |= p
    return new

from_torch(array) staticmethod

create from a tuple (xmin, ymin, xmax, ymax) which is the pytorch convention

Source code in madcad/box.py

@staticmethod
def from_torch(array: tuple) -> Box[vec2]:
    ''' create from a tuple `(xmin, ymin, xmax, ymax)` which is the pytorch convention '''
    return Box(vec2(array[0:2]), vec2(array[2:4]))

from_cv(array) staticmethod

create from a tuple (xmin, ymax, width, height) which is the opencv convention

Source code in madcad/box.py

@staticmethod
def from_cv(array: tuple) -> Box[vec2]:
    ''' create from a tuple `(xmin, ymax, width, height)` which is the opencv convention '''
    min = vec2(array[0:2])
    return Box(min, min + vec2(array[2:4]))

from_matrix(matrix, centered=False) staticmethod

reciprocal of self.to_matrix()

Source code in madcad/box.py

@staticmethod
def from_matrix(matrix:mat4, centered=False) -> Box:
    ''' reciprocal of `self.to_matrix()` '''
    m = matrix
    s = _from_affine_mat[type(m)]
    o = _from_affine_mat[type(m)]
    for i in range(len(v)):
        o[i] = m[-1][i]
        s[i] = m[i][i]
    if centered:
        return Box(o-s, o+s)
    else:
        return Box(o, o+s)

to_torch()

convert to a tuple (xmin, ymin, xmax, ymax)

Source code in madcad/box.py

def to_torch(self) -> tuple:
    ''' convert to a tuple `(xmin, ymin, xmax, ymax)` '''
    return (*self.min, *self.max)

to_cv()

convert to a tuple (xmin, ymax, width, height)

Source code in madcad/box.py

def to_cv(self) -> tuple:
    ''' convert to a tuple `(xmin, ymax, width, height)` '''
    return (*self.min, *self.size())

to_matrix(centered=False)

return a matrix interpolating between corners - not centered: (0,0,0) and (1,1,1) are opposite corners - centered: (-1,-1,-1) and (1,1,1) are opposite corners

Source code in madcad/box.py

def to_matrix(self, centered=False) -> mat4:
    ''' return a matrix interpolating between corners
        - not centered: (0,0,0) and (1,1,1) are opposite corners
        - centered:  (-1,-1,-1) and (1,1,1) are opposite corners
    '''
    if centered:
        s = self.size/2
        o = self.min + self.size/2
    else:
        s = self.size
        o = self.min
    m = _to_affine_mat[type(s)]()
    for i in range(len(s)):
        m[i][i] = s[i]
        m[len(s)][i] = o[i]
    return m

volume()

Volume inside

Source code in madcad/box.py

def volume(self) -> float:
    ''' Volume inside '''
    v = 1
    for edge in self.size:
        v *= edge
    return v

map(alignment)

linearly interpolate between box corners

Source code in madcad/box.py

def map(self, alignment:vec3) -> vec3:
    ''' linearly interpolate between box corners '''
    return glm.mix(self.min, self.max, alignment)

corners()

return an iterable of the corner points of the box

Source code in madcad/box.py

def corners(self) -> Iterator[vec3]:
    ''' return an iterable of the corner points of the box '''
    vec = copy(self.min)
    d = len(self.min)
    for i in range(2**d):
        for j in range(d):
            vec[j] = self.max[j] if (i>>j)&1 else self.min[j]
        yield copy(vec)

slice()

create a tuple of slice, a slice each dimension of the box

Source code in madcad/box.py

def slice(self):
    ''' create a tuple of slice, a slice each dimension of the box '''
    return tuple(slice(start,stop)  for start, stop in zip(self.min, self.max))

transform(transformation)

box bounding the current one in a transformed space

Source code in madcad/box.py

def transform(self, transformation) -> Box:
    ''' box bounding the current one in a transformed space '''
    if not self.isvalid():	return self
    return Box.from_iter(transformation * p   for p in self.corners())

touch_borders(other)

True if the box is touchng one or more border of the other box

Source code in madcad/box.py

def touch_borders(self, other:Box) -> bool:
    ''' True if the box is touchng one or more border of the other box '''
    return any((self.min <= other.min) + (self.max >= other.max))

touch_corners(other)

True if the box is touching one or more corner of the other box

Source code in madcad/box.py

def touch_corners(self, other:Box) -> bool:
    ''' True if the box is touching one or more corner of the other box '''
    return all((self.min <= other.min) + (self.max >= other.max))

contains(other)

return True if the given point is inside or on the surface of the box

Source code in madcad/box.py

def contains(self, other:Box|vec3) -> bool:
    ''' return True if the given point is inside or on the surface of the box '''
    if isinstance(other, Box):
        return all((self.min <= other.min) * (other.max <= self.max))
    else:
        return all((self.min <= other) * (other <= self.max))

inside(other)

return True if the given point is strictly inside the box

Source code in madcad/box.py

def inside(self, other:Box|vec3) -> bool:
    ''' return True if the given point is strictly inside the box '''
    if isinstance(other, Box):
        return all((self.min < other.min) * (other.max < self.max))
    else:
        return all((self.min < other) * (other < self.max))

isvalid()

Return True if the box defines a valid space (min coordinates <= max coordinates)

Source code in madcad/box.py

def isvalid(self) -> bool:
    ''' Return True if the box defines a valid space (min coordinates <= max coordinates) '''
    return any(self.min <= self.max)

isempty()

Return True if the box contains a non null volume

Source code in madcad/box.py

def isempty(self) -> bool:
    ''' Return True if the box contains a non null volume '''
    return not any(self.min < self.max)

__add__(other)

component-wise addition

Source code in madcad/box.py

def __add__(self, other:Box|vec3):
    ''' component-wise addition '''
    if isinstance(other, Box):	
        return Box(self.min + other.min, self.max + other.max)
    else:		
        return Box(self.min + other, self.max + other)

__sub__(other)

component-wise substraction

Source code in madcad/box.py

def __sub__(self, other:Box|vec3):
    ''' component-wise substraction '''
    if isinstance(other, Box):	
        return Box(self.min - other.min, self.max - other.max)
    else:		
        return Box(self.min - other, self.max - other)

__iadd__(other)

Source code in madcad/box.py

def __iadd__(self, other:Box|vec3):
    if isinstance(other, Box):	
        self.min += other.min
        self.max += other.max
    else:
        self.min += other
        self.max += other
    return self

__isub__(other)

Source code in madcad/box.py

def __isub__(self, other:Box|vec3):
    if isinstance(other, Box):	
        self.min -= other.min
        self.max -= other.max
    else:
        self.min -= other
        self.max -= other
    return self

merge_update(other)

Extend the volume of the current box to bound the given point or box

Source code in madcad/box.py

def merge_update(self, other:Box|vec3) -> Box:
    ''' Extend the volume of the current box to bound the given point or box '''
    if isinstance(other, Box):
        self.min = glm.min(self.min, other.min)
        self.max = glm.max(self.max, other.max)
    else:
        self.min = glm.min(self.min, other)
        self.max = glm.max(self.max, other)
    return self

intersection_update(other)

Reduce the volume of the current box to the intersection between the 2 boxes

Source code in madcad/box.py

def intersection_update(self, other:Box|vec3) -> Box:
    ''' Reduce the volume of the current box to the intersection between the 2 boxes '''
    if isinstance(other, Box):
        self.min = glm.max(self.min, other.min)
        self.max = glm.min(self.max, other.max)
    else:
        raise TypeError('expected a Box'.format(type(other)))
    return self

merge(other)

Return a box containing the current and the given box (or point)

Example:

>>> Box(vec2(1,2), vec2(2,3)) .merge(vec2(1,4))
Box(vec2(1,2), vec2(2,4))

>>> Box(vec2(1,2), vec2(2,3)) .merge(Box(vec2(1,-4), vec2(2,8)))
Box(vec2(1,-4), vec2(2,8))

Source code in madcad/box.py

def merge(self, other:Box|vec3) -> Box:
    ''' Return a box containing the current and the given box (or point) 

        Example:

            >>> Box(vec2(1,2), vec2(2,3)) .merge(vec2(1,4))
            Box(vec2(1,2), vec2(2,4))

            >>> Box(vec2(1,2), vec2(2,3)) .merge(Box(vec2(1,-4), vec2(2,8)))
            Box(vec2(1,-4), vec2(2,8))

    '''
    if isinstance(other, Box):
        return Box(	glm.min(self.min, other.min),
                    glm.max(self.max, other.max))
    else:
        return Box(	glm.min(self.min, other),
                    glm.max(self.max, other))

intersection(other)

Return a box for the volume common to the current and the given box

Example:

>>> Box(vec2(-1,2), vec2(2,3)) .intersection(Box(vec2(1,-4), vec2(2,8)))
Box(vec2(1,2), vec2(2,3))

Source code in madcad/box.py

def intersection(self, other) -> Box:
    ''' Return a box for the volume common to the current and the given box 

        Example:

            >>> Box(vec2(-1,2), vec2(2,3)) .intersection(Box(vec2(1,-4), vec2(2,8)))
            Box(vec2(1,2), vec2(2,3))

    '''
    return Box(	glm.max(self.min, other.min),
                glm.min(self.max, other.max))

cast(vec)

convert to a box with an other vector type

Source code in madcad/box.py

def cast(self, vec:type) -> Box:
    ''' convert to a box with an other vector type '''
    return Box(vec(self.min), vec(self.max))

__bool__()

Source code in madcad/box.py

def __bool__(self):
    return self.isvalid()

__array__()

Source code in madcad/box.py

def __array__(self):
    return np.array(self.to_torch())

__repr__()

Source code in madcad/box.py

def __repr__(self):
    return '{}({}, {})'.format(self.__class__.__name__, repr(self.min), repr(self.max))

Screw(location, resultant, moment)

A 3D torsor aka Screw - is a mathematical object defined as follow

• a resultant vector R
• a moment vector field M

The moment M is a vectorial field function, satisfying the relationship: M(A) = M(B) + cross(R, A-B)

Therefore it is possible to represent a localized screw such as

• R = resultant
• M(P) = moment vector at position P

Torsor are useful for generalized solid mechanics to handle multiple variables of the same nature

• wrench (force screw): Screw(position, force, torque)
• twist (velocity screw): Screw(position, angular_velocity, linear_velocity)
• momentum (momentum screw): Screw(position, mass * linear_velocity_com, inertia * angular_velocity_com)

All these screws makes it possible to represent all these values independently from expression location

Source code in madcad/mathutils.py

def __init__(self, location, resultant, moment):
	self.location = location
	self.resultant = resultant
	self.moment = moment

location = location instance-attribute

resultant = resultant instance-attribute

moment = moment instance-attribute

transform(transform)

change the screw from a frame to an other

Source code in madcad/mathutils.py

def transform(self, transform:mat4) -> Screw:
	''' change the screw from a frame to an other '''
	orient = mat3(transform)
	return Screw(
		transform * self.location,
		orient * self.resultant,
		orient * self.moment,
		)

locate(location)

relocate the screw at the given location

    this changes the resultant and moment vectors, but doesn't change the screw itself

Warning

this is different from translating the screw

Source code in madcad/mathutils.py

def locate(self, location:vec3) -> Screw:
	''' relocate the screw at the given location

		this changes the resultant and moment vectors, but doesn't change the screw itself

	Warning:
		this is different from translating the screw 
	'''
	return Screw(
		location,
		self.resultant,
		self.moment + cross(self.resultant, location - self.location),
		)

comoment(other)

screws comoment, performs relocation if necessary

Source code in madcad/mathutils.py

def comoment(self, other:Screw) -> float:
	''' screws comoment, performs relocation if necessary '''
	other = other.locate(self.location)
	return dot(self.moment, other.resultant) + dot(self.resultat * other.moment)

__add__(other)

screws addition, performs relocation if necessary

Source code in madcad/mathutils.py

def __add__(self, other:Screw) -> Screw:
	''' screws addition, performs relocation if necessary '''
	other = other.locate(self.location)
	return Screw(
		self.location,
		self.resultant + other.resultant,
		self.moment + other.moment,
		)

__sub__(other)

screws substraction, performs relocation if necessary

Source code in madcad/mathutils.py

def __sub__(self, other:Screw) -> Screw:
	''' screws substraction, performs relocation if necessary '''
	other = other.locate(self.location)
	return Screw(
		self.location,
		self.resultant - other.resultant,
		self.moment - other.moment,
		)

__neg__()

Source code in madcad/mathutils.py

def __neg__(self) -> Screw:
	return Screw(
		self.location,
		-self.resultant,
		-self.moment,
		)

__mul__(other)

Source code in madcad/mathutils.py

def __mul__(self, other:float) -> Screw:
	return Screw(
		self.location,
		self.resultant * other,
		self.moment * other,
		)	

__div__(other)

Source code in madcad/mathutils.py

def __div__(self, other:float) -> Screw:
	return Screw(
		self.location,
		self.resultant / other,
		self.moment / other,
		)

__eq__(other)

screws equality, performs relocation if necessary

Source code in madcad/mathutils.py

def __eq__(self, other:Screw) -> bool:
	''' screws equality, performs relocation if necessary '''
	other = other.locate(self.location)
	return glm.all(self.moment == other.moment) and glm.all(self.resultant == other.resultant)

to_matrix(location=None)

convert the screw to its matrix form

represents the relocate operation for any screw and the derivative of the rotation matrix for velocity screws

Source code in madcad/mathutils.py

def to_matrix(self, location=None) -> mat4:
	''' convert the screw to its matrix form

		represents the relocate operation for any screw and the derivative of the rotation matrix for velocity screws
	'''
	base = self.locate(location or vec3(0))
	return skew(base.resultant, base.moment)

from_matrix(mat, location=None) staticmethod

convert a velocity matrix into a screw at world location 0

Source code in madcad/mathutils.py

@staticmethod
def from_matrix(mat: mat3|mat4, location=None) -> Screw:
	''' convert a velocity matrix into a screw at world location 0 '''
	return Screw(location or vec3(0), *unskew(mat))

from_rate(f, t, dt=1e-06) staticmethod

compute a Screw of a frame by rating the given function f at the given instant t

f is supposed to - take a time instant and return a frame matrix - be continuous and f(t-dt) and f(t+dt) to exist

Source code in madcad/mathutils.py

@staticmethod
def from_rate(f:Callable[[float],mat3|mat4], t:float, dt=1e-6) -> Screw:
	''' compute a Screw of a frame by rating the given function `f` at the given instant `t`

		`f` is supposed to
		- take a time instant and return a frame matrix
		- be continuous and `f(t-dt)` and `f(t+dt)` to exist
	'''
	return Screw.from_matrix((f(t+dt) - f(t-dt)) / (2*dt) @ affineInverse(f(t)))

__repr__()

Source code in madcad/mathutils.py

def __repr__(self):
	return '{}(\n\t{}, \n\t{}, \n\t{})'.format(self.__class__.__name__, 
		repr(self.location), repr(self.resultant), repr(self.moment))

display(scene)

Source code in madcad/mathutils.py

def display(self, scene):
	# TODO draw resultant and moment vectors at the current location
	raise NotImplementedError("In developement")

comoment(t1, t2)

Comomentum of screws: dot(M1, R2) + dot(M2, R1)

The result is independent of torsors location

Source code in madcad/mathutils.py

def comoment(t1:Screw, t2:Screw) -> float:
	''' Comomentum of screws:   `dot(M1, R2)  +  dot(M2, R1)`

		The result is independent of torsors location
	'''
	t2 = t2.locate(t1.location)
	return dot(t1.moment, t2.resultant) + dot(t2.moment, t1.resultant)

skew(r, t=None)

skew matrix of 3D vector v (pre cross product matrix). it t is given, it provides a fourth column (usefull for affine transforms)

r = vec3(...) a = vec3(...) cross(r, a) == skew(r) * a

Source code in madcad/mathutils.py

def skew(r:vec3, t:vec3=None) -> mat3|mat4:
	''' skew matrix of 3D vector `v` (pre cross product matrix). it `t` is given, it provides a fourth column (usefull for affine transforms) 

		>>> r = vec3(...)
		>>> a = vec3(...)
		>>> cross(r, a) == skew(r) * a
	'''
	if t is None:
		return mat3(
			 0,   +r.z, -r.y,
			-r.z,  0,   +r.x,
			+r.y, -r.x,  0,
			)
	else:
		return mat4(
			 0,   +r.z, -r.y, 0,
			-r.z,  0,   +r.x, 0,
			+r.y, -r.x,  0,   0,
			 t.x,  t.y,  t.z, 0,
			)

unskew(m)

vector(s) from which the given matrix is a skew matrix (approximated it it is not actually a skew matrix)

r = mat3(...) unskew(skew(r)) == r

Source code in madcad/mathutils.py

def unskew(m: mat3|mat4) -> vec3|tuple[vec3, vec3]:
	''' vector(s) from which the given matrix is a skew matrix (approximated it it is not actually a skew matrix) 

		>>> r = mat3(...)
		>>> unskew(skew(r)) == r
	'''
	if isinstance(m, mat3):
		return vec3(
			m[1][2] - m[2][1], 
			m[2][0] - m[0][2],
			m[0][1] - m[1][0],
			) * 0.5
	elif isinstance(m, mat4):
		return (
			vec3(
				m[1][2] - m[2][1], 
				m[2][0] - m[0][2],
				m[0][1] - m[1][0],
				) * 0.5,
			vec3(m[3]),
			)
	else:
		raise TypeError('unskew only supports mat3 and mat4')