Trying to debug an issue with point projection.

Projection Issue

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