240 lines
5.4 KiB
Python
240 lines
5.4 KiB
Python
"""
|
|
Various transforms used for by the 3D code
|
|
"""
|
|
|
|
import numpy as np
|
|
import numpy.linalg as linalg
|
|
|
|
from matplotlib import cbook
|
|
|
|
|
|
@cbook.deprecated("3.1")
|
|
def line2d(p0, p1):
|
|
"""
|
|
Return 2D equation of line in the form ax+by+c = 0
|
|
"""
|
|
# x + x1 = 0
|
|
x0, y0 = p0[:2]
|
|
x1, y1 = p1[:2]
|
|
#
|
|
if x0 == x1:
|
|
a = -1
|
|
b = 0
|
|
c = x1
|
|
elif y0 == y1:
|
|
a = 0
|
|
b = 1
|
|
c = -y1
|
|
else:
|
|
a = y0 - y1
|
|
b = x0 - x1
|
|
c = x0*y1 - x1*y0
|
|
return a, b, c
|
|
|
|
|
|
@cbook.deprecated("3.1")
|
|
def line2d_dist(l, p):
|
|
"""
|
|
Distance from line to point
|
|
line is a tuple of coefficients a, b, c
|
|
"""
|
|
a, b, c = l
|
|
x0, y0 = p
|
|
return abs((a*x0 + b*y0 + c) / np.hypot(a, b))
|
|
|
|
|
|
def _line2d_seg_dist(p1, p2, p0):
|
|
"""distance(s) from line defined by p1 - p2 to point(s) p0
|
|
|
|
p0[0] = x(s)
|
|
p0[1] = y(s)
|
|
|
|
intersection point p = p1 + u*(p2-p1)
|
|
and intersection point lies within segment if u is between 0 and 1
|
|
"""
|
|
|
|
x21 = p2[0] - p1[0]
|
|
y21 = p2[1] - p1[1]
|
|
x01 = np.asarray(p0[0]) - p1[0]
|
|
y01 = np.asarray(p0[1]) - p1[1]
|
|
|
|
u = (x01*x21 + y01*y21) / (x21**2 + y21**2)
|
|
u = np.clip(u, 0, 1)
|
|
d = np.hypot(x01 - u*x21, y01 - u*y21)
|
|
|
|
return d
|
|
|
|
|
|
@cbook.deprecated("3.1")
|
|
def line2d_seg_dist(p1, p2, p0):
|
|
"""distance(s) from line defined by p1 - p2 to point(s) p0
|
|
|
|
p0[0] = x(s)
|
|
p0[1] = y(s)
|
|
|
|
intersection point p = p1 + u*(p2-p1)
|
|
and intersection point lies within segment if u is between 0 and 1
|
|
"""
|
|
return _line2d_seg_dist(p1, p2, p0)
|
|
|
|
|
|
@cbook.deprecated("3.1", alternative="np.linalg.norm")
|
|
def mod(v):
|
|
"""3d vector length"""
|
|
return np.sqrt(v[0]**2+v[1]**2+v[2]**2)
|
|
|
|
|
|
def world_transformation(xmin, xmax,
|
|
ymin, ymax,
|
|
zmin, zmax):
|
|
dx, dy, dz = (xmax-xmin), (ymax-ymin), (zmax-zmin)
|
|
return np.array([[1/dx, 0, 0, -xmin/dx],
|
|
[0, 1/dy, 0, -ymin/dy],
|
|
[0, 0, 1/dz, -zmin/dz],
|
|
[0, 0, 0, 1]])
|
|
|
|
|
|
def view_transformation(E, R, V):
|
|
n = (E - R)
|
|
## new
|
|
# n /= np.linalg.norm(n)
|
|
# u = np.cross(V, n)
|
|
# u /= np.linalg.norm(u)
|
|
# v = np.cross(n, u)
|
|
# Mr = np.diag([1.] * 4)
|
|
# Mt = np.diag([1.] * 4)
|
|
# Mr[:3,:3] = u, v, n
|
|
# Mt[:3,-1] = -E
|
|
## end new
|
|
|
|
## old
|
|
n = n / np.linalg.norm(n)
|
|
u = np.cross(V, n)
|
|
u = u / np.linalg.norm(u)
|
|
v = np.cross(n, u)
|
|
Mr = [[u[0], u[1], u[2], 0],
|
|
[v[0], v[1], v[2], 0],
|
|
[n[0], n[1], n[2], 0],
|
|
[0, 0, 0, 1]]
|
|
#
|
|
Mt = [[1, 0, 0, -E[0]],
|
|
[0, 1, 0, -E[1]],
|
|
[0, 0, 1, -E[2]],
|
|
[0, 0, 0, 1]]
|
|
## end old
|
|
|
|
return np.dot(Mr, Mt)
|
|
|
|
|
|
def persp_transformation(zfront, zback):
|
|
a = (zfront+zback)/(zfront-zback)
|
|
b = -2*(zfront*zback)/(zfront-zback)
|
|
return np.array([[1, 0, 0, 0],
|
|
[0, 1, 0, 0],
|
|
[0, 0, a, b],
|
|
[0, 0, -1, 0]])
|
|
|
|
|
|
def ortho_transformation(zfront, zback):
|
|
# note: w component in the resulting vector will be (zback-zfront), not 1
|
|
a = -(zfront + zback)
|
|
b = -(zfront - zback)
|
|
return np.array([[2, 0, 0, 0],
|
|
[0, 2, 0, 0],
|
|
[0, 0, -2, 0],
|
|
[0, 0, a, b]])
|
|
|
|
|
|
def _proj_transform_vec(vec, M):
|
|
vecw = np.dot(M, vec)
|
|
w = vecw[3]
|
|
# clip here..
|
|
txs, tys, tzs = vecw[0]/w, vecw[1]/w, vecw[2]/w
|
|
return txs, tys, tzs
|
|
|
|
|
|
@cbook.deprecated("3.1")
|
|
def proj_transform_vec(vec, M):
|
|
return _proj_transform_vec(vec, M)
|
|
|
|
|
|
def _proj_transform_vec_clip(vec, M):
|
|
vecw = np.dot(M, vec)
|
|
w = vecw[3]
|
|
# clip here.
|
|
txs, tys, tzs = vecw[0] / w, vecw[1] / w, vecw[2] / w
|
|
tis = (0 <= vecw[0]) & (vecw[0] <= 1) & (0 <= vecw[1]) & (vecw[1] <= 1)
|
|
if np.any(tis):
|
|
tis = vecw[1] < 1
|
|
return txs, tys, tzs, tis
|
|
|
|
|
|
@cbook.deprecated("3.1")
|
|
def proj_transform_vec_clip(vec, M):
|
|
return _proj_transform_vec_clip(vec, M)
|
|
|
|
|
|
def inv_transform(xs, ys, zs, M):
|
|
iM = linalg.inv(M)
|
|
vec = _vec_pad_ones(xs, ys, zs)
|
|
vecr = np.dot(iM, vec)
|
|
try:
|
|
vecr = vecr / vecr[3]
|
|
except OverflowError:
|
|
pass
|
|
return vecr[0], vecr[1], vecr[2]
|
|
|
|
|
|
def _vec_pad_ones(xs, ys, zs):
|
|
return np.array([xs, ys, zs, np.ones_like(xs)])
|
|
|
|
|
|
@cbook.deprecated("3.1")
|
|
def vec_pad_ones(xs, ys, zs):
|
|
return _vec_pad_ones(xs, ys, zs)
|
|
|
|
|
|
def proj_transform(xs, ys, zs, M):
|
|
"""
|
|
Transform the points by the projection matrix
|
|
"""
|
|
vec = _vec_pad_ones(xs, ys, zs)
|
|
return _proj_transform_vec(vec, M)
|
|
|
|
|
|
transform = proj_transform
|
|
|
|
|
|
def proj_transform_clip(xs, ys, zs, M):
|
|
"""
|
|
Transform the points by the projection matrix
|
|
and return the clipping result
|
|
returns txs, tys, tzs, tis
|
|
"""
|
|
vec = _vec_pad_ones(xs, ys, zs)
|
|
return _proj_transform_vec_clip(vec, M)
|
|
|
|
|
|
def proj_points(points, M):
|
|
return np.column_stack(proj_trans_points(points, M))
|
|
|
|
|
|
def proj_trans_points(points, M):
|
|
xs, ys, zs = zip(*points)
|
|
return proj_transform(xs, ys, zs, M)
|
|
|
|
|
|
@cbook.deprecated("3.1")
|
|
def proj_trans_clip_points(points, M):
|
|
xs, ys, zs = zip(*points)
|
|
return proj_transform_clip(xs, ys, zs, M)
|
|
|
|
|
|
def rot_x(V, alpha):
|
|
cosa, sina = np.cos(alpha), np.sin(alpha)
|
|
M1 = np.array([[1, 0, 0, 0],
|
|
[0, cosa, -sina, 0],
|
|
[0, sina, cosa, 0],
|
|
[0, 0, 0, 1]])
|
|
return np.dot(M1, V)
|