Mapping the Surface of the Moon
In the past, one needed to have advanced equipment and complex technology to be able to map the surface of the Moon. However, in the modern day, anyone can create a map of the Moon. The purpose of this project is to show that anyone can make a map of the surface (or at least the visible side)—no matter who you are or where you come from—as long as you have the right knowledge!1. Data Collection
Initially, we planned to capture our own images of the Moon but ultimately opted for high-quality images from NASA's SVS Dial-A-Moon database. To correct for libration and align the features consistently, we manually adjusted the images and used MATLAB for refinement.
Here's an example of some photos we took though:
2. Photometric Stereo
To create a detailed normal map, we used the photometric stereo algorithm, introduced by Robert J. Woodham in 1980. This method estimates surface normals using multiple images taken under varying lighting conditions.
Mathematical Formulation
The photometric stereo technique relies on the equation:
$$ I = \frac{\rho}{\pi} \hat{s} \cdot \hat{n} $$
where:
- \( I \) is the pixel intensity
- \( \rho \) is the albedo (reflectance)
- \( \hat{s} \) is the light source direction vector
- \( \hat{n} \) is the unit surface normal vector
Solving for Surface Normals
When using three images, we form a \( 3\times 3 \) matrix that can be inverted:
This formula applies to each individual pixel across the 3 images (i.e., \( I_1 \), \( I_2 \), and \( I_3 \) are the intensities of the pixel in each of the three images).
For more than three images, we use the Least Squares method to estimate the normal vector:
$$ \begin{bmatrix} I_1 \\ I_2 \\ I_3 \\ \vdots \\ I_{24} \end{bmatrix} = \frac{\rho}{\pi} \begin{bmatrix} s_{x,1} & s_{y,1} & s_{z,1} \\ s_{x,2} & s_{y,2} & s_{z,2} \\ s_{x,3} & s_{y,3} & s_{z,3} \\ \vdots & \vdots & \vdots \\ s_{x,24} & s_{y,24} & s_{z,24} \end{bmatrix} \mathbf{n} $$ $$\mathbf{I} = \dfrac{\rho}{\pi}\mathbf{S}\mathbf{n}$$ $$\dfrac{\rho}{\pi}\mathbf{S}^T \mathbf{S} \mathbf{n} = \mathbf{S}^T \mathbf{I}$$ $$\mathbf{n} = \dfrac{\pi}{\rho}\left(\mathbf{S}^T \mathbf{S}\right)^{-1} \mathbf{S}^T \mathbf{I}$$where \( \mathbf{S} \) is the light source matrix and \( \mathbf{I} \) is the intensity vector. The reason for this is that using more than 3 images results in a \( n\times 3 \) matrix with \( n \) rows for \( n \) images and \( 3 \) columns (one for each of the \( x \), \( y \), and \( z \) directions of the light source vector in space).
The reason why doing this works is because multiplying a matrix \( \mathbf{S} \) with its transpose \( \mathbf{S}^T \) gives a square, invertible matrix.
3. Heightmap Rendering
To reconstruct a 3D model from the normal maps, we integrated the normal components similarly to a Riemann sum.
The height variation in the x-direction is computed as:
\( h_x = \dfrac{n_x}{n_z} d \)
where \( d \) is the pixel size. A similar approach is used for the y-direction.
So, for a pixel at position \( \left(x_0, y_0\right) \), with the origin defined to be at the top left corner, the height at that pixel would be:
\( \displaystyle h\left(x_0, y_0\right) = \dfrac{d}{n_z}\left(\sum_{k=0}^{x_0} n_x + \sum_{k=0}^{y_0} n_y\right) \)