Calculate the number of steradians contained in a spherical wedge with radius $r = 1$, defined by $\theta=\frac{\pi}{6}$, $\phi=\frac{\pi}{6}$ centered around vector $(\frac{\sqrt{2}}{4}, \frac{\sqrt{2}}{4}, \frac{\sqrt{3}}{2})$.
Consider a camera looking at Lambertian wall with constant albedo, illuminated by a light source at infinity such that the radiance emitted by the wall is $L$ in all directions. The angle between the optical axis and the wall’s surface normal is 60 degrees. The focal length of the camera is 50mm, and the pixels are 1mm by 1mm.
We will solve this for a pixel centered at the optical axis and also consider three distances along the line of sight between the pinhole and wall (1000mm, 2000mm, 4000mm). We will also walk through four steps and take advantage of small area approximations.
We see a diffuse torus centered at the origin in an orthographic camera, looking down the z-axis. The parameters of the torus are shown in the Figure "Problem3 torus" and the albedo is $\rho$.
This torus is illuminated by a distant point light source whose direction is $(0,0,1)$. There is no other illumination.
What is the brightness at a point $(x, y)$ on the surface?
We have a square area source and a square occluder, both parallel to a plane.
The edge lengths of the source and occluder are 2 and 4, respectively, and they are vertically above one another with their centers aligned. The distances from the occluder to the source and plane are both 3.
What is the area of the umbra on the plane?
What is the area of the penumbra on the plane?
The goal of this problem is to implement a couple of different algorithms that reconstruct a surface using the concept of photometric stereo.
Additionally, you will implement the specular removal technique of Mallick et al., which enables photometric stereo reconstruction of certain non-Lambertian materials.
You can assume a Lambertian reflectance function once specularties are removed, but the albedo is unknown and non-constant in the images.
Your program will take in multiple images as input along with the light source direction (and color when necessary) for each image.
Synthetic Images, Specular Sphere Images, Pear Images for Part 1, 2, 3: Available in *.pickle
files (graciously provided by Satya Mallick) which contain
im1
, im2
, im3
, im4
... images.
l1
, l2
, l3
, l4
... light source directions.
c
(when required) color of light source.
Implement the photometric stereo technique described in Forsyth and Ponce 2.2.4 (Photometric Stereo: Shape from Multiple Shaded Images) and the lecture notes.
Your program should have two parts:
Read in the images and corresponding light source directions, and estimate the surface normals and albedo map.
Reconstruct the depth map from the surface normals. You can first try the naive scanline-based shape by integration method described in the book. If this does not work well on real images, you can use the implementation of the Horn integration technique given below in horn_integrate
function.
Try using only im1
, im2
and im4
first. Display your outputs as mentioned below.
Then use all four images. (Most accurate).
For each of the above cases you must output:
The estimated albedo map.
The estimated surface normals by showing both
A wireframe of depth map.
An example of outputs is shown in the Figure "Problem5 example".
Note: You will find all the data for this part in synthetic_data.pickle
.
## Example: How to read and access data from a pickle
import pickle
import matplotlib.pyplot as plt
%matplotlib inline
pickle_in = open("synthetic_data.pickle", "rb")
#data = pickle.load(pickle_in)
data = pickle.load(pickle_in, encoding="latin1")
# data is a dict which stores each element as a key-value pair.
print("Keys: " + str(data.keys()))
# To access the value of an entity, refer it by its key.
print("Image:")
plt.imshow(data["im1"], cmap = "gray")
plt.show()
print("Light source direction: " + str(data["l1"]))
import numpy as np
from scipy.signal import convolve
from numpy import linalg
def horn_integrate(gx, gy, mask, niter):
'''
horn_integrate recovers the function g from its partial
derivatives gx and gy.
mask is a binary image which tells which pixels are
involved in integration.
niter is the number of iterations.
typically 100,000 or 200,000,
although the trend can be seen even after 1000 iterations.
'''
g = np.ones(np.shape(gx))
gx = np.multiply(gx, mask)
gy = np.multiply(gy, mask)
A = np.array([[0,1,0],[0,0,0],[0,0,0]]) #y-1
B = np.array([[0,0,0],[1,0,0],[0,0,0]]) #x-1
C = np.array([[0,0,0],[0,0,1],[0,0,0]]) #x+1
D = np.array([[0,0,0],[0,0,0],[0,1,0]]) #y+1
d_mask = A + B + C + D
den = np.multiply(convolve(mask,d_mask,mode="same"),mask)
den[den == 0] = 1
rden = 1.0 / den
mask2 = np.multiply(rden, mask)
m_a = convolve(mask, A, mode="same")
m_b = convolve(mask, B, mode="same")
m_c = convolve(mask, C, mode="same")
m_d = convolve(mask, D, mode="same")
term_right = np.multiply(m_c, gx) + np.multiply(m_d, gy)
t_a = -1.0 * convolve(gx, B, mode="same")
t_b = -1.0 * convolve(gy, A, mode="same")
term_right = term_right + t_a + t_b
term_right = np.multiply(mask2, term_right)
for k in range(niter):
g = np.multiply(mask2, convolve(g, d_mask, mode="same")) + term_right
return g
def photometric_stereo(images, lights, mask):
'''
your implementaion
'''
b = np.zeros((images.shape[1],images.shape[2],3))
albedo = np.ones(images[0].shape)
normals = np.dstack((np.zeros(images[0].shape),
np.zeros(images[0].shape),
np.ones(images[0].shape)))
p = np.zeros(images[0].shape)
q = np.zeros(images[0].shape)
for i in range(images.shape[1]):
for j in range(images.shape[2]):
b[i,j,:] = np.dot(np.dot(np.linalg.inv(np.dot(lights.T,lights)), lights.T),images[:,i,j])
albedo[i,j] = np.linalg.norm(b[i,j,:])
normals[i,j,:] = b[i,j,:]/albedo[i,j]
p[i,j] = normals[i,j,0]/normals[i,j,2]
q[i,j] = normals[i,j,1]/normals[i,j,2]
H = np.ones(images[0].shape)
H_horn = horn_integrate(p,q,mask,10000)
p = np.multiply(mask,p)
q = np.multiply(mask,q)
for i in range(1,images[0].shape[0]):
H[i,0] = H[i-1,0] + q[i,0]
for i in range(images[0].shape[0]):
for j in range(1,images[0].shape[1]):
H[i,j] = H[i,j-1] + p[i,j]
return albedo, normals, H, H_horn
from mpl_toolkits.mplot3d import Axes3D
pickle_in = open("synthetic_data.pickle", "rb")
#data = pickle.load(pickle_in)
data = pickle.load(pickle_in, encoding="latin1")
lights = np.vstack((data["l1"], data["l2"], data["l4"]))
# lights = np.vstack((data["l1"], data["l2"], data["l3"], data["l4"]))
images = []
images.append(data["im1"])
images.append(data["im2"])
# images.append(data["im3"])
images.append(data["im4"])
images = np.array(images)
mask = np.ones(data["im1"].shape)
albedo, normals, depth, horn = photometric_stereo(images, lights, mask)
# --------------------------------------------------------------------------
# Following code is just a working example so you don't get stuck with any
# of the graphs required. You may want to write your own code to align the
# results in a better layout.
# --------------------------------------------------------------------------
# Stride in the plot, you may want to adjust it to different images
stride = 5
# showing albedo map
fig = plt.figure()
albedo_max = albedo.max()
albedo = albedo / albedo_max
plt.imshow(albedo, cmap="gray")
plt.show()
# showing normals as three separate channels
figure = plt.figure()
ax1 = figure.add_subplot(131)
ax1.imshow(normals[..., 0])
ax2 = figure.add_subplot(132)
ax2.imshow(normals[..., 1])
ax3 = figure.add_subplot(133)
ax3.imshow(normals[..., 2])
plt.show()
# showing normals as quiver
X, Y, _ = np.meshgrid(np.arange(0,np.shape(normals)[0], stride),
np.arange(0,np.shape(normals)[1], stride),
np.arange(1))
X = X[..., 0]
Y = Y[..., 0]
Z = depth[::stride,::stride].T
NX = normals[..., 0][::stride,::-stride].T
NY = normals[..., 1][::-stride,::stride].T
NZ = normals[..., 2][::stride,::stride].T
fig = plt.figure(figsize=(5, 5))
ax = fig.gca(projection='3d')
plt.quiver(X,Y,Z,NX,NY,NZ, length=5)
plt.show()
# plotting wireframe depth map
H = depth[::stride,::stride]
fig = plt.figure()
ax = fig.gca(projection='3d')
ax.plot_surface(X,Y, H.T)
plt.show()
H = horn[::stride,::stride]
fig = plt.figure()
ax = fig.gca(projection='3d')
ax.plot_surface(X,Y, H.T)
plt.show()
from mpl_toolkits.mplot3d import Axes3D
pickle_in = open("synthetic_data.pickle", "rb")
#data = pickle.load(pickle_in)
data = pickle.load(pickle_in, encoding="latin1")
#lights = np.vstack((data["l1"], data["l2"], data["l4"]))
lights = np.vstack((data["l1"], data["l2"], data["l3"], data["l4"]))
images = []
images.append(data["im1"])
images.append(data["im2"])
images.append(data["im3"])
images.append(data["im4"])
images = np.array(images)
mask = np.ones(data["im1"].shape)
albedo, normals, depth, horn = photometric_stereo(images, lights, mask)
# --------------------------------------------------------------------------
# Following code is just a working example so you don't get stuck with any
# of the graphs required. You may want to write your own code to align the
# results in a better layout.
# --------------------------------------------------------------------------
# Stride in the plot, you may want to adjust it to different images
stride = 5
# showing albedo map
fig = plt.figure()
albedo_max = albedo.max()
albedo = albedo / albedo_max
plt.imshow(albedo, cmap="gray")
plt.show()
# showing normals as three separate channels
figure = plt.figure()
ax1 = figure.add_subplot(131)
ax1.imshow(normals[..., 0])
ax2 = figure.add_subplot(132)
ax2.imshow(normals[..., 1])
ax3 = figure.add_subplot(133)
ax3.imshow(normals[..., 2])
plt.show()
# showing normals as quiver
X, Y, _ = np.meshgrid(np.arange(0,np.shape(normals)[0], stride),
np.arange(0,np.shape(normals)[1], stride),
np.arange(1))
X = X[..., 0]
Y = Y[..., 0]
Z = depth[::stride,::stride].T
NX = normals[..., 0][::stride,::-stride].T
NY = normals[..., 1][::-stride,::stride].T
NZ = normals[..., 2][::stride,::stride].T
fig = plt.figure(figsize=(5, 5))
ax = fig.gca(projection='3d')
plt.quiver(X,Y,Z,NX,NY,NZ, length=5)
plt.show()
# plotting wireframe depth map
H = depth[::stride,::stride]
fig = plt.figure()
ax = fig.gca(projection='3d')
ax.plot_surface(X,Y, H.T)
plt.show()
H = horn[::stride,::stride]
fig = plt.figure()
ax = fig.gca(projection='3d')
ax.plot_surface(X,Y, H.T)
plt.show()
Implement the specularity removal technique described in Beyond Lambert: Reconstructing Specular Surfaces Using Color (by Mallick, Zickler, Kriegman, and Belhumeur; CVPR 2005).
Your program should input an RGB image and light source color and output the corresponding SUV image.
Try this out first with the specular sphere images and then with the pear images.
For each specular sphere and pear images, include
The original image (in RGB colorspace).
The recovered $S$ channel of the image.
The recovered diffuse part of the image - Use $G = \sqrt{U^2+V^2}$ to represent the diffuse part.
Note: You will find all the data for this part in specular_sphere.pickle
and specular_pear.pickle
.
def get_rot_mat(rot_v, unit=None):
'''
Takes a vector and returns the rotation matrix required to align the
unit vector(2nd arg) to it.
'''
if unit is None:
unit = [1.0, 0.0, 0.0]
rot_v = rot_v/np.linalg.norm(rot_v)
uvw = np.cross(rot_v, unit) #axis of rotation
rcos = np.dot(rot_v, unit) #cos by dot product
rsin = np.linalg.norm(uvw) #sin by magnitude of cross product
#normalize and unpack axis
if not np.isclose(rsin, 0):
uvw = uvw/rsin
u, v, w = uvw
# Compute rotation matrix
R = (
rcos * np.eye(3) +
rsin * np.array([
[ 0, -w, v],
[ w, 0, -u],
[-v, u, 0]
]) +
(1.0 - rcos) * uvw[:,None] * uvw[None,:]
)
return R
def RGBToSUV(I_rgb, rot_vec):
'''
your implementation which takes an RGB image and a vector encoding
the orientation of S channel wrt to RGB
'''
R = get_rot_mat(rot_vec)
S = np.ones(I_rgb.shape[:2])
G = np.ones(I_rgb.shape[:2])
SUV = np.zeros(I_rgb.shape)
for i in range(I_rgb.shape[0]):
for j in range(I_rgb.shape[1]):
SUV[i][j] = np.dot(R,I_rgb[i][j])
S[i][j] = SUV[i][j][0]
G[i][j] = np.linalg.norm((SUV[i][j][1],SUV[i][j][2]))
return S, G
pickle_in = open("specular_sphere.pickle", "rb")
#data = pickle.load(pickle_in)
data = pickle.load(pickle_in, encoding="latin1")
# For image 1
im = 'im1'
S, G = RGBToSUV(data[im], np.hstack((data["c"][0][0],
data["c"][1][0],
data["c"][2][0])))
plt.subplot(131)
plt.imshow((data[im]-np.min(data[im]))/(np.max(data[im])-np.min(data[im])))
plt.title('original image')
plt.subplot(132)
plt.imshow(S,cmap='gray')
plt.title('S channel')
plt.subplot(133)
plt.imshow(G,cmap='gray')
plt.title('diffuse part of the image')
plt.show()
# For image 2
im = 'im2'
S, G = RGBToSUV(data[im], np.hstack((data["c"][0][0],
data["c"][1][0],
data["c"][2][0])))
plt.subplot(131)
plt.imshow((data[im]-np.min(data[im]))/(np.max(data[im])-np.min(data[im])))
plt.title('original image')
plt.subplot(132)
plt.imshow(S,cmap='gray')
plt.title('S channel')
plt.subplot(133)
plt.imshow(G,cmap='gray')
plt.title('diffuse part of the image')
plt.show()
# For image 3
im = 'im3'
S, G = RGBToSUV(data[im], np.hstack((data["c"][0][0],
data["c"][1][0],
data["c"][2][0])))
plt.subplot(131)
plt.imshow((data[im]-np.min(data[im]))/(np.max(data[im])-np.min(data[im])))
plt.title('original image')
plt.subplot(132)
plt.imshow(S,cmap='gray')
plt.title('S channel')
plt.subplot(133)
plt.imshow(G,cmap='gray')
plt.title('diffuse part of the image')
plt.show()
# For image 4
im = 'im4'
S, G = RGBToSUV(data[im], np.hstack((data["c"][0][0],
data["c"][1][0],
data["c"][2][0])))
plt.subplot(131)
plt.imshow((data[im]-np.min(data[im]))/(np.max(data[im])-np.min(data[im])))
plt.title('original image')
plt.subplot(132)
plt.imshow(S,cmap='gray')
plt.title('S channel')
plt.subplot(133)
plt.imshow(G,cmap='gray')
plt.title('diffuse part of the image')
plt.show()
pickle_in = open("specular_pear.pickle", "rb")
#data = pickle.load(pickle_in)
data = pickle.load(pickle_in, encoding="latin1")
# For image 1
im = 'im1'
S, G = RGBToSUV(data[im], np.hstack((data["c"][0][0],
data["c"][1][0],
data["c"][2][0])))
plt.subplot(131)
plt.imshow((data[im]-np.min(data[im]))/(np.max(data[im])-np.min(data[im])))
plt.title('original image')
plt.subplot(132)
plt.imshow(S,cmap='gray')
plt.title('S channel')
plt.subplot(133)
plt.imshow(G,cmap='gray')
plt.title('diffuse part of the image')
plt.show()
# For image 2
im = 'im2'
S, G = RGBToSUV(data[im], np.hstack((data["c"][0][0],
data["c"][1][0],
data["c"][2][0])))
plt.subplot(131)
plt.imshow((data[im]-np.min(data[im]))/(np.max(data[im])-np.min(data[im])))
plt.title('original image')
plt.subplot(132)
plt.imshow(S,cmap='gray')
plt.title('S channel')
plt.subplot(133)
plt.imshow(G,cmap='gray')
plt.title('diffuse part of the image')
plt.show()
# For image 3
im = 'im3'
S, G = RGBToSUV(data[im], np.hstack((data["c"][0][0],
data["c"][1][0],
data["c"][2][0])))
plt.subplot(131)
plt.imshow((data[im]-np.min(data[im]))/(np.max(data[im])-np.min(data[im])))
plt.title('original image')
plt.subplot(132)
plt.imshow(S,cmap='gray')
plt.title('S channel')
plt.subplot(133)
plt.imshow(G,cmap='gray')
plt.title('diffuse part of the image')
plt.show()
# For image 4
im = 'im4'
S, G = RGBToSUV(data[im], np.hstack((data["c"][0][0],
data["c"][1][0],
data["c"][2][0])))
plt.subplot(131)
plt.imshow((data[im]-np.min(data[im]))/(np.max(data[im])-np.min(data[im])))
plt.title('original image')
plt.subplot(132)
plt.imshow(S,cmap='gray')
plt.title('S channel')
plt.subplot(133)
plt.imshow(G,cmap='gray')
plt.title('diffuse part of the image')
plt.show()
Combine parts 1 and 2 by running your photometric stereo code on the diffuse components of the specular sphere and pear images.
For comparison, run your photometric stereo code on the original images (converted to grayscale) as well. You should notice erroneous "bumps" in the resulting reconstructions, as a result of violating the Lambertian assumption.
For each specular sphere and pear image sets, using all the four images, include:
The estimated albedo map (original and diffuse)
The estimated surface normals (original and diffuse) by showing both
A wireframe of depth map (original and diffuse)
def rgb2gray(rgb):
return np.dot(rgb[...,:3], [0.299, 0.587, 0.114])
# ---------------------------------------------------------------------------
# You may reuse the code for photometric_stereo here.
# Write your code below to process the data and send it to photometric_stereo
# and display the albedo, normals and depth maps.
# ---------------------------------------------------------------------------
from mpl_toolkits.mplot3d import Axes3D
pickle_in = open("specular_sphere.pickle", "rb")
#data = pickle.load(pickle_in)
data = pickle.load(pickle_in, encoding="latin1")
#lights = np.vstack((data["l1"], data["l2"], data["l4"]))
lights = np.vstack((data["l1"], data["l2"], data["l3"], data["l4"]))
images = []
images.append(rgb2gray(data["im1"]))
images.append(rgb2gray(data["im2"]))
images.append(rgb2gray(data["im3"]))
images.append(rgb2gray(data["im4"]))
images = np.array(images)
mask = np.ones(rgb2gray(data["im1"]).shape)
a = rgb2gray(data["im4"])+rgb2gray(data["im3"])+rgb2gray(data["im2"])+rgb2gray(data["im1"])
thresh = np.min(a)+0.1*(np.max(a)-np.min(a))
mask[a < thresh] = 0
albedo, normals, depth, horn = photometric_stereo(images, lights, mask)
# --------------------------------------------------------------------------
# Following code is just a working example so you don't get stuck with any
# of the graphs required. You may want to write your own code to align the
# results in a better layout.
# --------------------------------------------------------------------------
# Stride in the plot, you may want to adjust it to different images
stride = 5
# showing albedo map
fig = plt.figure()
albedo_max = albedo.max()
albedo = albedo / albedo_max
plt.imshow(albedo, cmap="gray")
plt.show()
# showing normals as three separate channels
figure = plt.figure()
ax1 = figure.add_subplot(131)
ax1.imshow(normals[..., 0])
ax2 = figure.add_subplot(132)
ax2.imshow(normals[..., 1])
ax3 = figure.add_subplot(133)
ax3.imshow(normals[..., 2])
plt.show()
# showing normals as quiver
X, Y, _ = np.meshgrid(np.arange(0,np.shape(normals)[0], stride),
np.arange(0,np.shape(normals)[1], stride),
np.arange(1))
X = X[..., 0]
Y = Y[..., 0]
Z = depth[::stride,::stride].T
NX = normals[..., 0][::stride,::-stride].T
NY = normals[..., 1][::-stride,::stride].T
NZ = normals[..., 2][::stride,::stride].T
fig = plt.figure(figsize=(5, 5))
ax = fig.gca(projection='3d')
plt.quiver(X,Y,Z,NX,NY,NZ, length=5)
plt.show()
# plotting wireframe depth map
H = depth[::stride,::stride]
fig = plt.figure()
ax = fig.gca(projection='3d')
ax.plot_surface(X,Y, H.T)
plt.show()
H = horn[::stride,::stride]
fig = plt.figure()
ax = fig.gca(projection='3d')
ax.plot_surface(X,Y, H.T)
plt.show()
S1, G1 = RGBToSUV(data['im1'], np.hstack((data["c"][0][0],
data["c"][1][0],
data["c"][2][0])))
S2, G2 = RGBToSUV(data['im2'], np.hstack((data["c"][0][0],
data["c"][1][0],
data["c"][2][0])))
S3, G3 = RGBToSUV(data['im3'], np.hstack((data["c"][0][0],
data["c"][1][0],
data["c"][2][0])))
S4, G4 = RGBToSUV(data['im4'], np.hstack((data["c"][0][0],
data["c"][1][0],
data["c"][2][0])))
images = []
images.append(G1)
images.append(G2)
images.append(G3)
images.append(G4)
images = np.array(images)
albedo, normals, depth, horn = photometric_stereo(images, lights, mask)
# --------------------------------------------------------------------------
# Following code is just a working example so you don't get stuck with any
# of the graphs required. You may want to write your own code to align the
# results in a better layout.
# --------------------------------------------------------------------------
# Stride in the plot, you may want to adjust it to different images
stride = 5
# showing albedo map
fig = plt.figure()
albedo_max = albedo.max()
albedo = albedo / albedo_max
plt.imshow(albedo, cmap="gray")
plt.show()
# showing normals as three separate channels
figure = plt.figure()
ax1 = figure.add_subplot(131)
ax1.imshow(normals[..., 0])
ax2 = figure.add_subplot(132)
ax2.imshow(normals[..., 1])
ax3 = figure.add_subplot(133)
ax3.imshow(normals[..., 2])
plt.show()
# showing normals as quiver
X, Y, _ = np.meshgrid(np.arange(0,np.shape(normals)[0], stride),
np.arange(0,np.shape(normals)[1], stride),
np.arange(1))
X = X[..., 0]
Y = Y[..., 0]
Z = depth[::stride,::stride].T
NX = normals[..., 0][::stride,::-stride].T
NY = normals[..., 1][::-stride,::stride].T
NZ = normals[..., 2][::stride,::stride].T
fig = plt.figure(figsize=(5, 5))
ax = fig.gca(projection='3d')
plt.quiver(X,Y,Z,NX,NY,NZ, length=5)
plt.show()
# plotting wireframe depth map
H = depth[::stride,::stride]
fig = plt.figure()
ax = fig.gca(projection='3d')
ax.plot_surface(X,Y, H.T)
plt.show()
H = horn[::stride,::stride]
fig = plt.figure()
ax = fig.gca(projection='3d')
ax.plot_surface(X,Y, H.T)
plt.show()
# ---------------------------------------------------------------------------
# You may reuse the code for photometric_stereo here.
# Write your code below to process the data and send it to photometric_stereo
# and display the albedo, normals and depth maps.
# ---------------------------------------------------------------------------
from mpl_toolkits.mplot3d import Axes3D
pickle_in = open("specular_pear.pickle", "rb")
#data = pickle.load(pickle_in)
data = pickle.load(pickle_in, encoding="latin1")
#lights = np.vstack((data["l1"], data["l2"], data["l4"]))
lights = np.vstack((data["l1"], data["l2"], data["l3"], data["l4"]))
images = []
images.append(rgb2gray(data["im1"]))
images.append(rgb2gray(data["im2"]))
images.append(rgb2gray(data["im3"]))
images.append(rgb2gray(data["im4"]))
images = np.array(images)
mask = np.ones(rgb2gray(data["im1"]).shape)
a = rgb2gray(data["im4"])+rgb2gray(data["im3"])+rgb2gray(data["im2"])+rgb2gray(data["im1"])
thresh = np.min(a)+0.03*(np.max(a)-np.min(a))
mask[a < thresh] = 0
albedo, normals, depth, horn = photometric_stereo(images, lights, mask)
# --------------------------------------------------------------------------
# Following code is just a working example so you don't get stuck with any
# of the graphs required. You may want to write your own code to align the
# results in a better layout.
# --------------------------------------------------------------------------
# Stride in the plot, you may want to adjust it to different images
stride = 5
# showing albedo map
fig = plt.figure()
albedo_max = albedo.max()
albedo = albedo / albedo_max
plt.imshow(albedo, cmap="gray")
plt.show()
# showing normals as three separate channels
figure = plt.figure()
ax1 = figure.add_subplot(131)
ax1.imshow(normals[..., 0])
ax2 = figure.add_subplot(132)
ax2.imshow(normals[..., 1])
ax3 = figure.add_subplot(133)
ax3.imshow(normals[..., 2])
plt.show()
# showing normals as quiver
X, Y, _ = np.meshgrid(np.arange(0,np.shape(normals)[0], stride),
np.arange(0,np.shape(normals)[1], stride),
np.arange(1))
X = X[..., 0]
Y = Y[..., 0]
Z = depth[::stride,::stride].T
NX = normals[..., 0][::stride,::-stride].T
NY = normals[..., 1][::-stride,::stride].T
NZ = normals[..., 2][::stride,::stride].T
fig = plt.figure(figsize=(5, 5))
ax = fig.gca(projection='3d')
plt.quiver(X,Y,Z,NX,NY,NZ, length=5)
plt.show()
# plotting wireframe depth map
H = depth[::stride,::stride]
fig = plt.figure()
ax = fig.gca(projection='3d')
ax.plot_surface(X,Y, H.T)
plt.show()
H = horn[::stride,::stride]
fig = plt.figure()
ax = fig.gca(projection='3d')
ax.plot_surface(X,Y, H.T)
plt.show()
S1, G1 = RGBToSUV(data['im1'], np.hstack((data["c"][0][0],
data["c"][1][0],
data["c"][2][0])))
S2, G2 = RGBToSUV(data['im2'], np.hstack((data["c"][0][0],
data["c"][1][0],
data["c"][2][0])))
S3, G3 = RGBToSUV(data['im3'], np.hstack((data["c"][0][0],
data["c"][1][0],
data["c"][2][0])))
S4, G4 = RGBToSUV(data['im4'], np.hstack((data["c"][0][0],
data["c"][1][0],
data["c"][2][0])))
images = []
images.append(G1)
images.append(G2)
images.append(G3)
images.append(G4)
images = np.array(images)
albedo, normals, depth, horn = photometric_stereo(images, lights, mask)
# --------------------------------------------------------------------------
# Following code is just a working example so you don't get stuck with any
# of the graphs required. You may want to write your own code to align the
# results in a better layout.
# --------------------------------------------------------------------------
# Stride in the plot, you may want to adjust it to different images
stride = 5
# showing albedo map
fig = plt.figure()
albedo_max = albedo.max()
albedo = albedo / albedo_max
plt.imshow(albedo, cmap="gray")
plt.show()
# showing normals as three separate channels
figure = plt.figure()
ax1 = figure.add_subplot(131)
ax1.imshow(normals[..., 0])
ax2 = figure.add_subplot(132)
ax2.imshow(normals[..., 1])
ax3 = figure.add_subplot(133)
ax3.imshow(normals[..., 2])
plt.show()
# showing normals as quiver
X, Y, _ = np.meshgrid(np.arange(0,np.shape(normals)[0], stride),
np.arange(0,np.shape(normals)[1], stride),
np.arange(1))
X = X[..., 0]
Y = Y[..., 0]
Z = depth[::stride,::stride].T
NX = normals[..., 0][::stride,::-stride].T
NY = normals[..., 1][::-stride,::stride].T
NZ = normals[..., 2][::stride,::stride].T
fig = plt.figure(figsize=(5, 5))
ax = fig.gca(projection='3d')
plt.quiver(X,Y,Z,NX,NY,NZ, length=5)
plt.show()
# plotting wireframe depth map
H = depth[::stride,::stride]
fig = plt.figure()
ax = fig.gca(projection='3d')
ax.plot_surface(X,Y, H.T)
plt.show()
H = horn[::stride,::stride]
fig = plt.figure()
ax = fig.gca(projection='3d')
ax.plot_surface(X,Y, H.T)
plt.show()
In this portion of the assignment we will be exploring different methods of approximating local illumination of objects in a scene. As discovered in the photometeric stereo portion of this homework, we know that different light models work better with different view, illumination sources and materials. This last section of the homework will be an exercise in rendering surfaces. Here, you need use the surface normals from Part 3 of Problem 5 to calculate the image intensity of the specular sphere and pear, with various light sources, different materials, and using a number of illumination models. For the sake of simplicity, multiple reflections of light rays, and occlusion of light rays due to object/scene can be ignored.
The surface normals of the specular sphere and the pear from Part 3 of Problem 5. For comparison, You should display the rendering results for both normals calculated from the original image and the diffuse components.
Assume that the albedo map is uniform.
One of the simplest models available to render 3D objections with illumination is the Lambertian model. This model finds the apparent brightness to an observer using the direction of the light source $\mathbf{L}$ and the normal vector on the surface of the object $\mathbf{N}$. The brightness intensity at a given point on an object’s surface, $\mathbf{I_d}$, with a single light source is found using the following relationship:
$$\mathbf{I_d} = \mathbf{L} \cdot \mathbf{N} (I_l\mathbf{C})$$where, $\mathbf{C}$ and $I_l$ are the the color and intensity of the light source respectively.
One major drawback of Lambertian illumination is that it only considers the diffuse light in its calculation of brightness intensity. One other major component to illumination rendering is the specular component. The specular reflectance is the component of light that is reflected in a single direction, as opposed to all directions, which is the case in diffuse reflectance. One of the most used models to compute surface brightness with specular components is the Phong illumination model. This model combines ambient lighting, diffused reflectance as well as specular reflectance to find the brightness on a surface. Phong shading also considers the material in the scene which is characterized by four values: the ambient reflection constant ($k_a$), the diffuse reflection constant ($k_d$), the specular reflection constant ($k_s$) and $\alpha$ the Phong constant, which is the ‘shininess’ of an object. Furthermore, since the specular component produces ‘rays’, only some of which would be observed by a single observer, the observer’s viewing direction ($\mathbf{V}$) must also be known. For some scene with known material parameters with $M$ light sources the light intensity $\mathbf{I}_{phong}$ on a surface with normal vector $\mathbf{N}$ seen from viewing direction $\mathbf{V}$ can be computed by:
$$\mathbf{I}_{phong} = k_{a}\mathbf{I}_{a} + \sum_{m\in M}\left\{k_d(\mathbf{L}_{m}\cdot\mathbf{N})\mathbf{I}_{m,d} + k_{s}(\mathbf{R}_{m}\cdot\mathbf{V})^{\alpha}\mathbf{I}_{m,s}\right\}\text{,}$$$$\mathbf{R}_{m} = 2\mathbf{N}(\mathbf{L}_{m}\cdot\mathbf{N}) - \mathbf{L}_{m}\text{,}$$where $\mathbf{I}_{a}$, is the color and intensity of the ambient lighting, $\mathbf{I}_{m,d}$ and $\mathbf{I}_{m,s}$ are the color values for the diffuse and specular light of the $m$th light source.
Please complete the following:
Write the function lambertian()
that calculates the Lambertian light intensity given the light direction $\mathbf{L}$ with color and intensity $\mathbf{C}$ and $I_l = 1$, and normal vector $\mathbf{N}$. Then use this function in a program that calculates and displays the specular sphere and the pear using each of the two lighting sources found in Table 1. Note: You do not need to worry about material coefficients in this model.
Write the function phong()
that calculates the Phong light intensity given the material constants $(k_a, k_d, k_s, \alpha)$, $\mathbf{V} = (0, 0, 1)^\top$, $\mathbf{N}$ and some number of $M$ light sources. Then use this function in a program that calculates and displays the specular sphere and the pear using each of the sets of coefficients found in Table 2 with each light source individually, and both light sources combined.
Hint: To avoid artifacts due to shadows, ensure that any negative intensities found are set to zero.
Table 1: Light Sources
$m$ | Location | Color (RGB) |
---|---|---|
1 | $(-\tfrac{1}{3},\tfrac{1}{3},\tfrac{1}{3})^{\top}$ | $(1,1,1)$ |
2 | $(1,0,0)^{\top}$ | $(1,.5,.5)$ |
Table 2: Material Coefficients
Mat. | $k_a$ | $k_d$ | $k_s$ | $\alpha$ |
---|---|---|---|---|
1 | $0$ | $0.1$ | $0.75$ | $5$ |
2 | $0$ | $0.5$ | $0.1$ | $5$ |
3 | $0$ | $0.5$ | $0.5$ | $10$ |
def lambertian(normals, lights, color, intensity, mask):
'''Your implementation'''
image = np.ones((normals.shape[0], normals.shape[1], 3))
for i in range(normals.shape[0]):
for j in range(normals.shape[1]):
image[i][j] = np.dot(np.dot(lights,normals[i][j]),np.dot(intensity,color))
image[ image < 0 ] = 0
image = np.multiply(image,np.stack((mask,mask,mask),axis = 2))
return image
# Output the rendering results
from mpl_toolkits.mplot3d import Axes3D
pickle_in = open("specular_sphere.pickle", "rb")
#data = pickle.load(pickle_in)
data = pickle.load(pickle_in, encoding="latin1")
#lights = np.vstack((data["l1"], data["l2"], data["l4"]))
lights = np.vstack((data["l1"], data["l2"], data["l3"], data["l4"]))
images = []
images.append(rgb2gray(data["im1"]))
images.append(rgb2gray(data["im2"]))
images.append(rgb2gray(data["im3"]))
images.append(rgb2gray(data["im4"]))
images = np.array(images)
mask = np.ones(rgb2gray(data["im1"]).shape)
a = rgb2gray(data["im4"])+rgb2gray(data["im3"])+rgb2gray(data["im2"])+rgb2gray(data["im1"])
thresh = np.min(a)+0.1*(np.max(a)-np.min(a))
mask[a < thresh] = 0
albedo, normals, depth, horn = photometric_stereo(images, lights, mask)
print('original sphere images')
lights1 = (-1/3,1/3,1/3)
color1 = (1,1,1)
intensity = 1
image = lambertian(normals, lights1, color1, intensity, mask)
plt.subplot(121)
#image = (image-np.min(image))/(np.max(image)-np.min(image))
plt.imshow(image)
plt.title('result for light source 1')
lights2 = (1,0,0)
color2 = (1,0.5,0.5)
intensity = 1
image = lambertian(normals, lights2, color2, intensity, mask)
plt.subplot(122)
#image = (image-np.min(image))/(np.max(image)-np.min(image))
plt.imshow(image)
plt.title('result for light source 2')
plt.show()
S1, G1 = RGBToSUV(data['im1'], np.hstack((data["c"][0][0],
data["c"][1][0],
data["c"][2][0])))
S2, G2 = RGBToSUV(data['im2'], np.hstack((data["c"][0][0],
data["c"][1][0],
data["c"][2][0])))
S3, G3 = RGBToSUV(data['im3'], np.hstack((data["c"][0][0],
data["c"][1][0],
data["c"][2][0])))
S4, G4 = RGBToSUV(data['im4'], np.hstack((data["c"][0][0],
data["c"][1][0],
data["c"][2][0])))
images = []
images.append(G1)
images.append(G2)
images.append(G3)
images.append(G4)
images = np.array(images)
albedo, normals, depth, horn = photometric_stereo(images, lights, mask)
print('diffuse part')
image = lambertian(normals, lights1, color1, intensity, mask)
plt.subplot(121)
#image = (image-np.min(image))/(np.max(image)-np.min(image))
plt.imshow(image)
plt.title('result for light source 1')
image = lambertian(normals, lights2, color2, intensity, mask)
plt.subplot(122)
#image = (image-np.min(image))/(np.max(image)-np.min(image))
plt.imshow(image)
plt.title('result for light source 2')
plt.show()
# Output the rendering results
from mpl_toolkits.mplot3d import Axes3D
pickle_in = open("specular_pear.pickle", "rb")
#data = pickle.load(pickle_in)
data = pickle.load(pickle_in, encoding="latin1")
#lights = np.vstack((data["l1"], data["l2"], data["l4"]))
lights = np.vstack((data["l1"], data["l2"], data["l3"], data["l4"]))
images = []
images.append(rgb2gray(data["im1"]))
images.append(rgb2gray(data["im2"]))
images.append(rgb2gray(data["im3"]))
images.append(rgb2gray(data["im4"]))
images = np.array(images)
mask = np.ones(rgb2gray(data["im1"]).shape)
a = rgb2gray(data["im4"])+rgb2gray(data["im3"])+rgb2gray(data["im2"])+rgb2gray(data["im1"])
thresh = np.min(a)+0.03*(np.max(a)-np.min(a))
mask[a < thresh] = 0
albedo, normals, depth, horn = photometric_stereo(images, lights, mask)
print('original sphere images')
lights1 = (-1/3,1/3,1/3)
color1 = (1,1,1)
intensity = 1
image = lambertian(normals, lights1, color1, intensity, mask)
plt.subplot(121)
#image = (image-np.min(image))/(np.max(image)-np.min(image))
plt.imshow(image)
plt.title('result for light source 1')
lights2 = (1,0,0)
color2 = (1,0.5,0.5)
intensity = 1
image = lambertian(normals, lights2, color2, intensity, mask)
plt.subplot(122)
#image = (image-np.min(image))/(np.max(image)-np.min(image))
plt.imshow(image)
plt.title('result for light source 2')
plt.show()
S1, G1 = RGBToSUV(data['im1'], np.hstack((data["c"][0][0],
data["c"][1][0],
data["c"][2][0])))
S2, G2 = RGBToSUV(data['im2'], np.hstack((data["c"][0][0],
data["c"][1][0],
data["c"][2][0])))
S3, G3 = RGBToSUV(data['im3'], np.hstack((data["c"][0][0],
data["c"][1][0],
data["c"][2][0])))
S4, G4 = RGBToSUV(data['im4'], np.hstack((data["c"][0][0],
data["c"][1][0],
data["c"][2][0])))
images = []
images.append(G1)
images.append(G2)
images.append(G3)
images.append(G4)
images = np.array(images)
albedo, normals, depth, horn = photometric_stereo(images, lights, mask)
print('diffuse part')
image = lambertian(normals, lights1, color1, intensity, mask)
plt.subplot(121)
#image = (image-np.min(image))/(np.max(image)-np.min(image))
plt.imshow(image)
plt.title('result for light source 1')
image = lambertian(normals, lights2, color2, intensity, mask)
plt.subplot(122)
#image = (image-np.min(image))/(np.max(image)-np.min(image))
plt.imshow(image)
plt.title('result for light source 2')
plt.show()
def phong(normals, lights, color, material, view, mask):
'''Your implementation'''
image = np.ones((normals.shape[0], normals.shape[1], 3))
for i in range(normals.shape[0]):
for j in range(normals.shape[1]):
Rm = 2*np.dot(normals[i][j],np.dot(lights,normals[i][j]))-np.transpose(lights)
image[i][j] = np.dot((material[1]*np.dot(lights,normals[i][j])),color)+np.dot((material[2]*np.power(np.dot(Rm,np.transpose(view)),material[3])),color)
image[image < 0] = 0
image = np.multiply(image,np.stack((mask,mask,mask),axis = 2))
return image
def phong2(normals, lights1, color1, lights2, color2, material, view, mask):
'''Your implementation'''
image = np.ones((normals.shape[0], normals.shape[1], 3))
for i in range(normals.shape[0]):
for j in range(normals.shape[1]):
Rm1 = 2*np.dot(normals[i][j],np.dot(lights1,normals[i][j]))-np.transpose(lights1)
I1 = np.dot((material[1]*np.dot(lights1,normals[i][j])),color1)+np.dot((material[2]*np.power(np.dot(Rm1,np.transpose(view)),material[3])),color1)
Rm2 = 2*np.dot(normals[i][j],np.dot(lights2,normals[i][j]))-np.transpose(lights2)
I2 = np.dot((material[1]*np.dot(lights2,normals[i][j])),color2)+np.dot((material[2]*np.power(np.dot(Rm2,np.transpose(view)),material[3])),color2)
image[i][j] = I1 + I2
#image = (image-np.min(image))/(np.max(image)-np.min(image))
image[image < 0] = 0
image = np.multiply(image,np.stack((mask,mask,mask),axis = 2))
return image
# Output the rendering results
from mpl_toolkits.mplot3d import Axes3D
pickle_in = open("specular_sphere.pickle", "rb")
#data = pickle.load(pickle_in)
data = pickle.load(pickle_in, encoding="latin1")
#lights = np.vstack((data["l1"], data["l2"], data["l4"]))
lights = np.vstack((data["l1"], data["l2"], data["l3"], data["l4"]))
images = []
images.append(rgb2gray(data["im1"]))
images.append(rgb2gray(data["im2"]))
images.append(rgb2gray(data["im3"]))
images.append(rgb2gray(data["im4"]))
images = np.array(images)
mask = np.ones(rgb2gray(data["im1"]).shape)
a = rgb2gray(data["im4"])+rgb2gray(data["im3"])+rgb2gray(data["im2"])+rgb2gray(data["im1"])
thresh = np.min(a)+0.1*(np.max(a)-np.min(a))
mask[a < thresh] = 0
albedo, normals, depth, horn = photometric_stereo(images, lights, mask)
# Output the rendering results
material = [(0,0.1,0.75,5),(0,0.5,0.1,5),(0,0.5,0.5,10)]
lights1 = (-1/3,1/3,1/3)
color1 = (1,1,1)
lights2 = (1,0,0)
color2 = (1,0.5,0.5)
view = (0,0,1)
def plot_phong(normals,lights1,color1,lights2, color2, material, view, mask):
fig=plt.figure(1,figsize=(10,10))
fig.add_subplot(331)
image = phong(normals, lights1, color1, material[0], view, mask)
#image = (image-np.min(image))/(np.max(image)-np.min(image))
plt.imshow(image)
plt.title('material 1 light source 1')
image = phong(normals, lights1, color1, material[1], view, mask)
fig.add_subplot(332)
#image = (image-np.min(image))/(np.max(image)-np.min(image))
plt.imshow(image)
plt.title('material 2 light source 1')
image = phong(normals, lights1, color1, material[2], view, mask)
fig.add_subplot(333)
#image = (image-np.min(image))/(np.max(image)-np.min(image))
plt.imshow(image)
plt.title('material 3 light source 1')
image = phong(normals, lights2, color2, material[0], view, mask)
fig.add_subplot(334)
#image = (image-np.min(image))/(np.max(image)-np.min(image))
plt.imshow(image)
plt.title('material 1 light source 2')
image = phong(normals, lights2, color2, material[1], view, mask)
fig.add_subplot(335)
#image = (image-np.min(image))/(np.max(image)-np.min(image))
plt.imshow(image)
plt.title('material 2 light source 2')
image = phong(normals, lights2, color2, material[2], view, mask)
fig.add_subplot(336)
#image = (image-np.min(image))/(np.max(image)-np.min(image))
plt.imshow(image)
plt.title('material 3 light source 2')
image = phong2(normals,lights1,color1,lights2, color2, material[0], view, mask)
fig.add_subplot(337)
#image = (image-np.min(image))/(np.max(image)-np.min(image))
plt.imshow(image)
plt.title('material 1 both lights')
image = phong2(normals,lights1,color1,lights2, color2, material[1], view, mask)
fig.add_subplot(338)
#image = (image-np.min(image))/(np.max(image)-np.min(image))
plt.imshow(image)
plt.title('material 2 both lights')
image = phong2(normals,lights1,color1,lights2, color2, material[2], view, mask)
fig.add_subplot(339)
#image = (image-np.min(image))/(np.max(image)-np.min(image))
plt.imshow(image)
plt.title('material 3 both lights')
plt.show()
plot_phong(normals,lights1,color1,lights2, color2, material, view, mask)
S1, G1 = RGBToSUV(data['im1'], np.hstack((data["c"][0][0],
data["c"][1][0],
data["c"][2][0])))
S2, G2 = RGBToSUV(data['im2'], np.hstack((data["c"][0][0],
data["c"][1][0],
data["c"][2][0])))
S3, G3 = RGBToSUV(data['im3'], np.hstack((data["c"][0][0],
data["c"][1][0],
data["c"][2][0])))
S4, G4 = RGBToSUV(data['im4'], np.hstack((data["c"][0][0],
data["c"][1][0],
data["c"][2][0])))
images = []
images.append(G1)
images.append(G2)
images.append(G3)
images.append(G4)
images = np.array(images)
albedo, normals, depth, horn = photometric_stereo(images, lights, mask)
plot_phong(normals,lights1,color1,lights2, color2, material, view, mask)
from mpl_toolkits.mplot3d import Axes3D
pickle_in = open("specular_pear.pickle", "rb")
#data = pickle.load(pickle_in)
data = pickle.load(pickle_in, encoding="latin1")
#lights = np.vstack((data["l1"], data["l2"], data["l4"]))
lights = np.vstack((data["l1"], data["l2"], data["l3"], data["l4"]))
images = []
images.append(rgb2gray(data["im1"]))
images.append(rgb2gray(data["im2"]))
images.append(rgb2gray(data["im3"]))
images.append(rgb2gray(data["im4"]))
images = np.array(images)
mask = np.ones(rgb2gray(data["im1"]).shape)
a = rgb2gray(data["im4"])+rgb2gray(data["im3"])+rgb2gray(data["im2"])+rgb2gray(data["im1"])
thresh = np.min(a)+0.03*(np.max(a)-np.min(a))
mask[a < thresh] = 0
albedo, normals, depth, horn = photometric_stereo(images, lights, mask)
plot_phong(normals,lights1,color1,lights2, color2, material, view, mask)
S1, G1 = RGBToSUV(data['im1'], np.hstack((data["c"][0][0],
data["c"][1][0],
data["c"][2][0])))
S2, G2 = RGBToSUV(data['im2'], np.hstack((data["c"][0][0],
data["c"][1][0],
data["c"][2][0])))
S3, G3 = RGBToSUV(data['im3'], np.hstack((data["c"][0][0],
data["c"][1][0],
data["c"][2][0])))
S4, G4 = RGBToSUV(data['im4'], np.hstack((data["c"][0][0],
data["c"][1][0],
data["c"][2][0])))
images = []
images.append(G1)
images.append(G2)
images.append(G3)
images.append(G4)
images = np.array(images)
albedo, normals, depth, horn = photometric_stereo(images, lights, mask)
plot_phong(normals,lights1,color1,lights2, color2, material, view, mask)