A script that demonstrate how to compute the shadow of a 3D object

compute_shadow.py

#!/usr/bin/env python
"""
        DO WHAT THE FUCK YOU WANT TO PUBLIC LICENSE
                    Version 2, December 2004

 Copyright (C) 2017 Ruben Di Battista <[email protected]>

 Everyone is permitted to copy and distribute verbatim or modified
 copies of this license document, and changing it is allowed as long
 as the name is changed.

            DO WHAT THE FUCK YOU WANT TO PUBLIC LICENSE
   TERMS AND CONDITIONS FOR COPYING, DISTRIBUTION AND MODIFICATION

  0. You just DO WHAT THE FUCK YOU WANT TO.

NOTES:
    - Theory discussion: http://math.stackexchange.com/q/2161369/41944
    - Works only with Python2 due to np.genfromtxt handling of unicode strings
        in Python3.
    - Test mesh 'Male.OBJ' here: http://tf3dm.com/3d-model/male-base-88907.html
"""

import numpy as np
import re
import matplotlib.pyplot as plt

from mpl_toolkits.mplot3d import Axes3D


def rotation_matrix(axis, theta):
    """
    Return the rotation matrix associated with counterclockwise rotation about
    the given axis by theta radians.
    From: http://stackoverflow.com/a/6802723/1334711
    """
    axis = np.asarray(axis)
    axis = axis/np.linalg.norm(axis)
    a = np.cos(theta/2.0)
    b, c, d = -axis*np.sin(theta/2.0)
    aa, bb, cc, dd = a*a, b*b, c*c, d*d
    bc, ad, ac, ab, bd, cd = b*c, a*d, a*c, a*b, b*d, c*d
    return np.matrix(np.array([[aa+bb-cc-dd, 2*(bc+ad), 2*(bd-ac)],
                     [2*(bc-ad), aa+cc-bb-dd, 2*(cd+ab)],
                     [2*(bd+ac), 2*(cd-ab), aa+dd-bb-cc]]))


# ::: PLANE OF PROJECTION :::
# Plane to project
n = np.matrix([0, 0, 1]).T

# A point on the plane
r_0 = np.matrix(np.zeros(3)).T

# ::: SUN DIRECTION :::
azimuth = 0*np.pi/180
elevation = 20*np.pi/180
I = np.matrix(np.eye(3))

# Rotate first azimuth (z-axis)
I_rot = rotation_matrix(I[:, 2], azimuth)*I

# Rotate then elevation (y-axis)
I_rot_rot = rotation_matrix(I[:, 1], elevation)*I_rot

# The direction vector of the sun is the first column of the last matrix
s = I_rot_rot[:, 0]

# Only vertices (line start by 'v')
lines = (line for line in open('Male.OBJ', 'r')
         if re.match(r'^v', line))

verts = np.genfromtxt(lines, usecols=[1, 2, 3]).T

# Rotate guy 90 deg around x as standing along z
verts = rotation_matrix(np.array([1, 0, 0]), np.pi/2)*verts

# Rotate him now (around z) to face the sun on axis x
verts = rotation_matrix(np.array([0, 0, 1]), -np.pi/2)*verts

# (Standing on the ground) Feet are put on the plane z=0
# Get lowest vertex
feet_coors = np.abs(np.min(verts[2, :]))
verts[2, :] = verts[2, :] + feet_coors*np.ones(verts.shape[1])

# Relative vector
rr = verts-r_0

# Projection Matrix
k = I[:, 2]
k_proj = k - np.multiply(1/(k.T*s), s)
P = np.hstack((I[:, 0:2], k_proj))

# Compute projected points
proj = P*verts

fig = plt.figure()
ax = Axes3D(fig)

# 2D plot of shadow

# Plots only a subset of points
# Back to array since Matplotlib plots strange things if 
# using np.matrix types (http://stackoverflow.com/q/42475393/1334711)
verts = np.array(verts[:, ::100])
proj = np.array(proj[:, ::100])


ax.plot(proj[0, :], proj[1, :], proj[2, :], '+')
ax.plot(verts[0, :], verts[1, :], verts[2, :], 'o')

plt.show()