Skip to content

Instantly share code, notes, and snippets.

Embed
What would you like to do?
Trying to debug an issue with point projection.
import numpy
import math
def make_sphere(steps=512):
points = list()
colors = list()
for s in range(steps):
y = math.sin((8 * s / float(steps)) * 2 * math.pi)
x = math.cos((32 * s / float(steps)) * 2 * math.pi)
z = math.sin((4 * s / float(steps)) * 2 * math.pi)
r = 128 + (128 * math.sin(s))
g = 128 + (128 * math.sin(3 * s))
b = 128 + (128 * math.sin(5 * s))
points.append([x, y, z])
colors.append([r, g, b])
return points, colors
def make_image(filename, num_pts=1000, scale=1):
from PIL import Image
sphere = numpy.zeros((num_pts, 3))
for i in range(num_pts):
sphere[i, :] = random_on_unit_circle()
proj = project_points(sphere*scale, numpy.asarray([-10, 1, 0]), numpy.asarray([0, 1.0, 0]), 0, (0, 0), 0.1, 0)
img = Image.new('RGB', (640,480))
for p in range(proj.shape[0]):
if proj[p,0] < -320 or proj[p,0] >= 320 or proj[p,1] < -240 or proj[p,1] >= 240:
continue
color = tuple(int(255*x) for x in list(sphere[p,:]))
img.putpixel((320+int(proj[p,0]), 240+int(proj[p,1])), color)
img.save(filename)
def test_render(filename, frames):
import numpngw
numpngw.write_apng(filename, frames, delay=250, use_palette=True)
def random_on_unit_circle():
pt = numpy.random.uniform(-1, 1, (3,))
return pt/numpy.linalg.norm(pt)
def points_in_front(points, camera_pos, camera_target):
indices = 1*(numpy.inner((points-camera_pos), (camera_target-camera_pos)) > 0)
return points[indices]
# Projective geometry style:
def build_intrinsic_translation(x, y):
# Principle point offset
return numpy.asarray([
[1, 0, x],
[0, 1, y],
[0, 0, 1]
])
def build_intrinsic_scaling_matrix(fx, fy):
# Focal length
return numpy.asarray([
[fx, 0, 0],
[0, fy, 0],
[0, 0, 1]
])
def build_intrinsic_shear(shear, fx):
# Axis skew
return numpy.asarray([
[1, shear/fx, 0],
[0, 1, 0],
[0, 0, 1]
])
def build_intrinsic_matrix(principle_point_offset, focal_length, shear):
# 2d translation * 2d scaling * 2d shear
return numpy.asarray([
[focal_length, shear, principle_point_offset[0]],
[0, focal_length, principle_point_offset[1]],
[0, 0, 1]
], dtype=numpy.float)
def build_extrinsic_transform(dx, dy, dz):
return numpy.asarray([
[1, 0, 0, dx],
[0, 1, 0, dy],
[0, 0, 1, dz]
])
def build_extrinsic_rotation_axis_angle(axis, angle):
x = axis[0]
y = axis[1]
z = axis[2]
c = math.cos(angle)
s = math.sin(angle)
n = 1-c
return numpy.asarray([
[x*x*n + c, x*y*n - z*s, x*z*n + y*s],
[y*c*n + z*s, y*y*n + c, y*z*n - x*s],
[z*x*n - y*s, z*y*n + x*s, z*z*n + c]
])
def build_extrinsic_rotation_look_at(position, target, up):
look = target - position
look /= numpy.linalg.norm(look)
s = numpy.cross(look, up)
s /= numpy.linalg.norm(s)
u = numpy.cross(s, look)
return numpy.vstack([
s,
u,
-look
])
def build_extrinsic(camera_pos, camera_target, up):
# Can do ex = [R|t] for a 3x4 matrix or can combine into rot and trans in one op with
# [R t] via [I|t] x [R|0]
# [0 1] [0|1] [ |1]
extrinsic = numpy.zeros((3,4))
#transform = build_extrinsic_transform(camera_pos[0], camera_pos[1], camera_pos[2])
rotation = build_extrinsic_rotation_look_at(camera_pos, camera_target, up)
#extrinsic[0:3,3] = transform
extrinsic[0,3] = -camera_pos[0]
extrinsic[1,3] = -camera_pos[1]
extrinsic[2,3] = -camera_pos[2]
extrinsic[0:3,0:3] = rotation
return extrinsic
def project_points(points, camera_pos, camera_target, rotation_around_target, principle_point, focal_length, skew):
intrinsic = build_intrinsic_matrix(principle_point, focal_length, skew)
up = numpy.dot(numpy.asarray([0,1.0,0]), build_extrinsic_rotation_axis_angle(camera_target - camera_pos, rotation_around_target))
extrinsic = build_extrinsic(camera_pos, camera_target, up)
final = numpy.dot(intrinsic, extrinsic) # Should be 3x4
# Augment points -> homogeneous coordinates, nx4
print(points[0:5,:])
pts = numpy.hstack((points, numpy.ones((points.shape[0], 1)))).T # 4xn
print(pts[:,0:5])
pts = numpy.dot(final, pts) # 3xn
print(pts[:,0:5])
return (pts / pts[2,:]).T[:,0:2]
# OpenGL Style:
# Ignore all these.
# View Matrix transforms points from world space to view/camera space.
# v' = Projection * View * Model * p
def view_matrix(camera_position, camera_target, camera_up):
"""Let camera_target be in world coords and camera_up be in unit coords."""
zaxis = camera_position - camera_target # Forward
zaxis /= numpy.linalg.norm(zaxis)
xaxis = numpy.cross(camera_up, zaxis) # Right
xaxis /= numpy.linalg.norm(xaxis)
yaxis = numpy.cross(zaxis, xaxis) # Up
mat = numpy.zeros((4,4))
mat[0:3,0] = xaxis
mat[0:3,1] = yaxis
mat[0:3,2] = zaxis
mat[3,0] = -numpy.inner(xaxis, camera_position)
mat[3,1] = -numpy.inner(yaxis, camera_position)
mat[3,2] = -numpy.inner(zaxis, camera_position)
mat[3,3] = 1.0
return mat
def projection_matrix(near, far, right, left, top, bottom):
mat = numpy.zeros((4,4))
mat[0,0] = (2.0*near)/(right-left)
mat[0,2] = (right+left)/(right-left)
mat[1,1] = (2.0*near)/(top-bottom)
mat[1,2] = (top+bottom)/(top-bottom)
mat[2,2] = -(far+near)/(far-near)
mat[2,3] = (-2.0*far*near)/(far-near)
mat[3,2] = -1
return mat
def project_points_to_screen(points, camera_position, camera_target, camera_up, near, far, resolution):
"""Given an array of points in 3D and a camera, calculate the 2D screen points for the items in front."""
vm = view_matrix(camera_position, camera_target, camera_up)
pm = projection_matrix(near, far, (resolution[0]/2), -(resolution[0]/2), resolution[1]/2, -(resolution[1]/2))
final_matrix = numpy.dot(vm,pm)
pts = numpy.dot(numpy.hstack([points, numpy.ones([512,1])]), final_matrix)
# Renorm
pts = numpy.divide(pts.T, pts[:,3]).T
pts = numpy.divide(pts.T, pts[:,2]).T
return pts
Sign up for free to join this conversation on GitHub. Already have an account? Sign in to comment
You can’t perform that action at this time.