Pythonで、SpaceCarvingのためのボクセルデータ投影をやってみました。
Pythonだと本当にすっきりと書けていいですね。
#coding:utf-8
import numpy as np
import math
from PIL import Image

# ユーザー定義変数
voxelWidth = 50 # ボクセルの各次元の個数
spaceSize = 1000.0 # ボクセル群が埋める空間の大きさ[mm]
radius = 250.0 # 球の半径[mm]
imageShape = np.array((640, 640)) # 画像解像度
viewAngleDeg = 90.0 # 画角[deg]
cameraPos = np.array([0.0, 0.0, -500.0]) # カメラ位置[mm]

# 自動算出変数
voxelSize = spaceSize / voxelWidth # ボクセル1つの大きさ[mm]
principlePoint = imageShape / 2.0 # 主点
viewAngleRad = viewAngleDeg / 180.0 * math.pi # 画角[rad]
focalLengthPix = imageShape / 2 / math.tan(viewAngleRad / 2.0) # 焦点距離[pix]

# 中心からの距離を要素として持つボクセルデータを生成
print 'calc distance'
voxelNo = np.arange(voxelWidth**3) # ボクセル番号。座標計算に用いる。
coords3d = np.zeros((3, voxelWidth**3)) # 三次元座標。各列がボクセル、各行がX, Y, Z。
coords3d[0,:] = (voxelNo % voxelWidth - voxelWidth / 2 + 0.5) * voxelSize # X座標
coords3d[1,:] = (voxelNo / voxelWidth % voxelWidth - voxelWidth / 2 + 0.5) * voxelSize # Y座標
coords3d[2,:] = (voxelNo / voxelWidth / voxelWidth - voxelWidth / 2 + 0.5) * voxelSize # Z座標
distance = np.sqrt(np.sum((coords3d * coords3d), axis=0)) # 各ボクセルの中心からの距離

# 球のボクセルデータを生成。内部は1。外部は0
print 'make voxel data'
voxels = np.zeros(voxelWidth**3)
voxels[distance < radius] = 1
del distance

# 射影行列の生成
print 'make projection matrix'
intrinsic = np.array([[focalLengthPix[0], 0.0, principlePoint[0]],\
[0.0, focalLengthPix[1], principlePoint[1]],\
[0.0, 0.0, 1.0]])
extrinsic = np.hstack((np.eye(3), -cameraPos.reshape(1, len(cameraPos)).transpose()))
proj = np.dot(intrinsic, extrinsic)

# 投影
print 'projection'
coords3d = np.vstack((coords3d, np.ones(voxelWidth**3)))
coords2d = np.dot(proj, coords3d)
coords2d[:2,:] /= coords2d[2,:]
del coords3d

# 球内部の投影座標を抽出
print 'extract coord'
intersectCoord = np.array(coords2d[:2, (voxels == 1) *\
(coords2d[0,:] >= 0) *\
(coords2d[0,:] < imageShape[0]) *\
(coords2d[1,:] >= 0) *\
(coords2d[1,:] < imageShape[1])], dtype=np.int64)
del coords2d, voxels

# 画像生成
print 'make image'
image = np.zeros(imageShape)
for i in range(intersectCoord.shape[1]):
image[intersectCoord[0,i], intersectCoord[1,i]] = 255

# 画像表示
Image.fromarray(image).show()
以下、投影結果の画像です。
voxelWidth=400で、メモリ6GBくらい使います。
全ボクセルの三次元座標を行列として確保してるのが痛い。

・voxelWidth=100
100

・voxelWidth=200
200


・voxelWidth=300
300

・voxelWidth=400
400