In [4]:
maple_leaf = Tracer2D("data/maple2d.json")
The word epicycle derives from the Ancient Greek word ἐπίκυκλος, literally meaning "circle moving on another circle". While it was originally proposed to discuss retrograde motion of planetary bodies, it is now Youtube folklore that any curve can be traced as an epicycle, provided enough circles. This was made popular by 3Blue1Brown (you should watch it if you haven't before). These animations are incredible in both their visual beauty and their mathematical elegance, as they arise from the Fourier transform, in an a manner that seems out of thin air.
Say we have $N$ spaced out samples from the curve. Assuming we're drawing in 2D, we can encode these as complex numbers $z_0, z_1, \ldots, z_{N-1} \in \mathbb{C}$. Applying the discrete Fourier transform gives us coefficients $c_0, c_1, \ldots, c_{N-1} \in \mathbb{C}$. Then the inverse transform (under the typical convention) tells us $$z_t = \frac{1}{N} \sum_{k=0}^{N-1} c_k e^{\frac{2\pi i}{N} kt}.$$ In particular, each of the terms $e^{\frac{2\pi i}{N} kt}$ (representing the orthodox orthogonal basis of periodic functions) exactly corresponds to a rotation with some frequency. This is the key fact that causes the epicycle animation to pop out. Although it may seem like a coincidence initially, it is actually due to the deep connection between the Fourier transform and $SO(2)$.
The typical Fourier series is only able to represent curves in space, as the frequency space is only parameterized by a single parameter. In this project, I will demonstrate a generalization of epicycles in order to model 2-manifolds embedded in $\mathbb{R}^3$.
In the code below, we generate the epicycles to trace out a maple leaf. The maple leaf is represented as a curve in 2D, which is in turn sampled 100 times to give a discrete series to perform DFT on.
import numpy as np
import json
from scipy.fft import fft, ifft
import matplotlib as mpl
import matplotlib.pyplot as plt
from tqdm.notebook import tqdm
import meshplot as mp
import open3d as o3d
import ipywidgets as widgets
from IPython.display import display
%matplotlib notebook
VERBOSE = True
Jupyter environment detected. Enabling Open3D WebVisualizer. [Open3D INFO] WebRTC GUI backend enabled. [Open3D INFO] WebRTCWindowSystem: HTTP handshake server disabled.
class Tracer2D:
def __init__(self, fp):
"""
Initialize structure for tracing known 2d points
"""
self.fp = fp
with open(fp, "r") as f:
self.raw_data = np.array(list(json.load(f)))
self.data = self.raw_data[:,0] - 1j * self.raw_data[:,1]
self.data -= np.mean(self.data)
self.n = self.data.shape[0]
self.spectral_data = fft(self.data)
self.order = np.argsort(-np.abs(self.spectral_data))
self.artists = []
self.init_plot()
def __length__(self):
return self.n
def init_plot(self):
"""
Initialize plot and save artists
"""
self.ax = plt.gca()
self.ax.cla()
self.ax.set_aspect("equal")
self.ax.set_ylim(-250, 250)
self.ax.set_xlim(-250, 250)
cur_pos = 0
for i in range(self.n):
point = plt.scatter([cur_pos.real], [cur_pos.imag], color="green", s=2)
self.ax.add_artist(point)
self.artists.append(point)
rot = 1/self.n * self.spectral_data[self.order[i]]
circle = plt.Circle((cur_pos.real, cur_pos.imag), np.abs(rot), color="black", fill=False)
self.ax.add_artist(circle)
self.artists.append(circle)
cur_pos += rot
point = plt.scatter([cur_pos.real], [cur_pos.imag], color="red", s=4)
self.ax.add_artist(point)
self.artists.append(point)
prefix = np.arange(self.n)/self.n < 1
trace, = plt.plot(self.data[prefix].real, self.data[prefix].imag, color="red")
self.ax.add_artist(trace)
self.artists.append(trace)
def get_plot(self, t):
"""
Update artists
"""
cur_pos = 0
for i in range(self.n):
rot = 1/self.n * self.spectral_data[self.order[i]] * np.exp(self.order[i] * t * 2 * np.pi * 1j)
self.artists[2*i].set_offsets([[cur_pos.real,cur_pos.imag]])
self.artists[2*i+1].center = (cur_pos.real, cur_pos.imag)
cur_pos += rot
self.artists[-2].set_offsets([[cur_pos.real,cur_pos.imag]])
prefix = np.arange(self.n)/self.n < t
self.artists[-1].set_data(self.data[prefix].real, self.data[prefix].imag)
freq_slider = widgets.FloatSlider(
value=0.0,
min=0.0,
max=1.0,
step=0.01,
description="t",
readout_format=".01f",
)
def plot(tracer):
t = freq_slider.value
tracer.get_plot(t)
We can obtain a set of equally-spaced samples from the outline of a maple leaf using coordinator. These samples are stored in maple2d.json.
maple_leaf = Tracer2D("data/maple2d.json")
def update_plot(event):
plot(maple_leaf)
freq_slider.observe(update_plot, names="value")
freq_slider
FloatSlider(value=0.0, description='t', max=1.0, readout_format='.01f', step=0.01)
The most direct extension of FFT to spheres is the spherical harmonics, which can be used as a basis to express any complex function over the sphere. In particular, $$F(\phi, \theta) = \sum_{l=0}^L \sum_{m=-l}^l f_l^m Y_l^m(\theta, \phi).$$ However, it is more important for our purposes that the image of the basis functions are on the sphere rather than the domain.
Another proposed decomposition is the Quaternion Fourier Transform (1). In the discrete sense, this transforms a 2D grid of quaternions $f(n, m)$ for $0 \le n < N, 0 \le m < M$ and produces $$F(u, v) = \frac{1}{NM} \sum_{m=0}^{M-1} \sum_{n=0}^{N-1} e^{-2\pi \mu \left(\frac{mv}{M} + \frac{nu}{N} \right)} f(n,m)$$ where $\mu$ is any unit pure quaternion. The inverse is similar. Just as with the Fourier transform, we note that $e^{-2\pi \mu \left(\cdots \right)}$ is in fact a unit quaternion, and so corresponds to an element of $SO(3)$. The Quaternion Fourier Transform was initially proposed for image processing purposes (2), hence the 2D grid input. However, the quaternions do not perform conjugation in the formula which makes it difficult to extract spheres from this.
As there has not been much work for this particular problem and current options are lacking, we must create our own method.
There is no lack of representations for 3D shapes, whether by their surfaces or by their volumes. The most common surface representation is the triangle mesh. In our setting, we assume that we have a triangle mesh already parameterized by two variables. Then we can sample a point cloud representation. From there, we will fit a sum of spheres to obtain our novel representation. Implicit in our goal of parameterization with spheres is that we are working with a closed surface of genus zero, meaning it may be deformed to a sphere.
def display_point_cloud(points, size=0.05):
"""
Plot point cloud with random colours
"""
colors = np.array(mpl.cm.get_cmap("tab20")(np.random.uniform(0, 1, (points.shape[0],))))[:,:3]
mp.plot(points, c=colors, shading={"point_size": size})
def display_mesh(mesh):
"""
Plot mesh with random colours
"""
vertices = np.array(mesh.vertices)
faces = np.array(mesh.triangles)
colors = np.array(mpl.cm.get_cmap("tab20")(np.random.uniform(0, 1, (faces.shape[0],))))[:,:3]
mp.plot(vertices, faces, c=colors)
Our problem will be to find a function that maps a sphere to a given shape as a sum of spheres. We represent the sphere with spherical coordinates $\theta \in [0, \pi], \phi \in [0, 2\pi)$. So our algorithm will take as input samples $z(\theta_i, \phi_i)$ on the shape's surface for some chosen $\theta_i, \phi_i$ and expand it to a function $Z(\theta, \phi)$ for all $\theta, \phi$ which coincides to the known samples for the specified $\theta_i, \phi_i$. The use of spherical coordinates allows for a natural mapping from a point on a sphere to a point on the reconstructed shape.
def spherical_to_cartesian(theta, phi):
"""
Convert spherical coordinates to Cartesian coordinates
"""
return np.array([
np.sin(theta) * np.cos(phi),
np.sin(theta) * np.sin(phi),
np.cos(theta)
])
def evenly_spaced_angles(J, K):
js, ks = np.meshgrid(np.arange(J), np.arange(K), indexing="ij")
return 2 * np.pi * js / J, 2 * np.pi * ks / K
if VERBOSE:
thetas, phis = evenly_spaced_angles(40, 40)
points = spherical_to_cartesian(thetas, phis).transpose(1, 2, 0).reshape(-1, 3)
print("Displaying evenly sampled spherical coordinates")
display_point_cloud(points)
Displaying evenly sampled spherical coordinates
Renderer(camera=PerspectiveCamera(children=(DirectionalLight(color='white', intensity=0.6, position=(0.0, 0.0,…
As a test case, we will reconstruct the Stanford bunny. The algorithm described in (4) implicitly maps genus-zero closed surfaces to spheres by fairing them. We use the meshes below which were obtained as the output of this algorithm. They share mesh topologies, but have different vertex positions, which gives a pointwise mapping from the sphere to the Stanford bunny. So there is a correspondence between the vertices of the sphere $v^s_i$ and the vertices of the target shape $v^t_i$.
target_file = "data/flow00000000.obj"
sphere_file = "data/flow00000217.obj"
target = o3d.io.read_triangle_mesh(target_file)
sphere = o3d.io.read_triangle_mesh(sphere_file)
if VERBOSE:
print("Displaying target mesh")
display_mesh(target)
print("Displaying sphere mesh")
display_mesh(sphere)
Displaying target mesh
Renderer(camera=PerspectiveCamera(children=(DirectionalLight(color='white', intensity=0.6, position=(0.0833379…
Displaying sphere mesh
Renderer(camera=PerspectiveCamera(children=(DirectionalLight(color='white', intensity=0.6, position=(-1.999735…
We can interpolate these pointwise mappings to sample from arbitrary $\theta, \phi$. In particular, for a given spherical coordinate $\theta, \phi$, we consider the ray $(r, \theta, \phi)$, $r > 0$. This intersects the sphere mesh at a face with vertices $v^s_0, v^s_1, v^s_2$. In particular, we can identify the point of intersection $$\texttt{SphericalToCartesian}(t_{\text{hit}}, \theta, \phi) = a v^s_0 + b v^s_1 + c v^s_2$$ for barycentric coordinates $a, b, c \ge 0,\, a+b+c=1$. Using this, we identify our sample $z(\theta, \phi)$ as $a v^t_0 + b v^t_1 + c v^t_2$.
def generate_samples(thetas, phis, target, sphere):
"""
Generate samples at given thetas and phis
target is the target mesh (e.g. Stanford bunny) and sphere is the corresponding sphere mesh
"""
target_vertices = np.array(target.vertices)
sphere_vertices = np.array(sphere.vertices)
sphere_triangles = np.array(sphere.triangles)
scene = o3d.t.geometry.RaycastingScene()
scene.add_triangles(o3d.t.geometry.TriangleMesh.from_legacy(sphere))
size = thetas.shape[0]
vals = np.zeros((size, 3))
rays = []
for i in range(size):
ray = spherical_to_cartesian(thetas[i], phis[i])
rays.append([0, 0, 0, ray[0], ray[1], ray[2]])
orays = o3d.core.Tensor(rays, dtype=o3d.core.Dtype.Float32)
result = scene.cast_rays(orays)
ts = result['t_hit'].numpy()
ids = result['primitive_ids'].numpy()
count = 0
for i in range(size):
ray = spherical_to_cartesian(thetas[i], phis[i])
ai = sphere_triangles[ids[count],0]
bi = sphere_triangles[ids[count],1]
ci = sphere_triangles[ids[count],2]
p = ts[count] * ray
area = np.linalg.norm(np.cross(sphere_vertices[ai,:]-sphere_vertices[bi,:], sphere_vertices[ai,:]-sphere_vertices[ci,:]))
bary_a = np.linalg.norm(np.cross(p-sphere_vertices[bi,:], p-sphere_vertices[ci,:])) / area
bary_b = np.linalg.norm(np.cross(p-sphere_vertices[ci,:], p-sphere_vertices[ai,:])) / area
vals[i,:] = bary_a * target_vertices[ai,:] + bary_b * target_vertices[bi,:] + (1 - bary_a - bary_b) * target_vertices[ci,:]
count += 1
return vals
J, K = 40, 40
thetas, phis = evenly_spaced_angles(J, K)
vals = generate_samples(thetas.flatten(), phis.flatten(), target, sphere)
vals = vals.reshape(J, K, 3)
if VERBOSE:
print("Displaying samples from target mesh")
display_point_cloud(vals.reshape(-1, 3))
Displaying samples from target mesh
Renderer(camera=PerspectiveCamera(children=(DirectionalLight(color='white', intensity=0.6, position=(0.0779113…
Given point samples $z(\theta_i, \phi_i)$, we seek to construct function $$Z(\theta, \phi) = \sum_{\Gamma} r^{\Gamma} R^{\Gamma} F^{\Gamma}(\theta, \phi)$$ such that $Z(\theta_i, \phi_i) = z(\theta_i, \phi_i)$ and $r^{\Gamma} \in \mathbb{R}_{>0}$, $R^{\Gamma} \in SO(3)$, $F^{\Gamma}: S^2 \to S^2$.
The plan is to choose a family of functions $F^{\Gamma}$ and then solve for $r^{\Gamma}$ and $R^{\Gamma}$ to fit the equation for these known arguments. Note that this is why we separated $F^{\Gamma}$ and $R^{\Gamma}$, even though $R^{\Gamma}F^{\Gamma}$ was already on the sphere. It allowed us to have a larger degree of freedom when solving for parameters.
It may be convenient to first take a look at a slight reduction of our problem by replacing the orthogonal matrices with general transformations $T^{\Gamma} \in \mathbb{R}^{3\times 3}$. Say we have chosen samples $z$ and functions $F^{\Gamma}$ and would like to solve the sum $$Z(\theta, \phi) = \sum_{\Gamma} T^{\Gamma} F^{\Gamma}(\theta, \phi).$$
This gives a linear system in which we can directly solve for the entries of $T$. In particular, for the x-coordinates:
$$\begin{pmatrix}z(\theta_1, \phi_1)_x \\ z(\theta_2, \phi_2)_x \\ \vdots \end{pmatrix} = \begin{pmatrix}F^{\Gamma_1}(\theta_1, \phi_1)_x & F^{\Gamma_1}(\theta_1, \phi_1)_y & F^{\Gamma_1}(\theta_1, \phi_1)_z & F^{\Gamma_2}(\theta_1, \phi_1)_x & \ldots \\ F^{\Gamma_1}(\theta_2, \phi_2)_x & F^{\Gamma_1}(\theta_2, \phi_2)_y & F^{\Gamma_1}(\theta_2, \phi_2)_z & F^{\Gamma_2}(\theta_2, \phi_2)_x & \ldots \\ \vdots & \vdots & \vdots & \vdots & \end{pmatrix} \begin{pmatrix}T^{\Gamma_1}_{11} \\ T^{\Gamma_1}_{12} \\ T^{\Gamma_1}_{13} \\ \vdots\end{pmatrix}.$$To avoid solving this dense system, we choose $F^{\Gamma}$ as $$F^{0, j, k}(\theta, \phi) = \begin{pmatrix} \sin(j\theta)\cos(k\phi) \\ \sin(j\theta)\sin(k\phi) \\ \cos(j\theta) \end{pmatrix}, \, F^{1, j, k}(\theta, \phi) = \begin{pmatrix} \cos(j\theta)\cos(k\phi) \\ \cos(j\theta)\sin(k\phi) \\ \sin(j\theta) \end{pmatrix}$$ for integers $j \in [0, \frac{J}{2})$ and $k\in [0, \frac{K}{2})$ for some even numbers $J, K$. Then the functions $F^{\Gamma}_x, F^{\Gamma}_y :S^2 \to \mathbb{R}$ consisting of the top two entries of $F$ are pairwise orthogonal for the inner product $$\langle f, g \rangle = \sum_{j=0}^J \sum_{k=0}^K f\left(\frac{2\pi j}{J}, \frac{2\pi k}{K}\right)g\left(\frac{2\pi j}{J}, \frac{2\pi k}{K}\right).$$ In fact, this is equivalent to the two-dimensional DFT.
We can now solve for the entries of $T^{\Gamma}$. For example, $$T^{0,j,k}_{11} = \frac{\langle z_x, F^{0,j,k}_x\rangle}{\langle F^{0,j,k}_x, F^{0,j,k}_x\rangle}.$$ The code below does this.
def compute_general_transforms(vals):
"""
Computes the general 3x3 transformations
Does this by solving the linear system explicitly with inner products
"""
J, K = vals.shape[:2]
thetas, phis = evenly_spaced_angles(J, K)
# scs = sin * cos 's, etc.
scs = np.zeros((J//2, K//2, 3))
sss = np.zeros((J//2, K//2, 3))
ccs = np.zeros((J//2, K//2, 3))
css = np.zeros((J//2, K//2, 3))
ERR = 1e-12
for j in tqdm(range(J//2)):
for k in range(K//2):
scs_denom = np.sum((np.sin(j * thetas) * np.cos(k * phis))**2)
if np.abs(scs_denom) > ERR:
scs[j,k,:] = np.sum(vals * (np.sin(j * thetas) * np.cos(k * phis))[:,:,None], axis=(0,1)) / scs_denom
sss_denom = np.sum((np.sin(j * thetas) * np.sin(k * phis))**2)
if np.abs(sss_denom) > ERR:
sss[j,k,:] = np.sum(vals * (np.sin(j * thetas) * np.sin(k * phis))[:,:,None], axis=(0,1)) / sss_denom
ccs_denom = np.sum((np.cos(j * thetas) * np.cos(k * phis))**2)
if np.abs(ccs_denom) > ERR:
ccs[j,k,:] = np.sum(vals * (np.cos(j * thetas) * np.cos(k * phis))[:,:,None], axis=(0,1)) / ccs_denom
css_denom = np.sum((np.cos(j * thetas) * np.sin(k * phis))**2)
if np.abs(css_denom) > ERR:
css[j,k,:] = np.sum(vals * (np.cos(j * thetas) * np.sin(k * phis))[:,:,None], axis=(0,1)) / css_denom
transforms = np.zeros((2, J//2, K//2, 3, 3))
transforms[0,:,:,:,0] = scs
transforms[0,:,:,:,1] = sss
transforms[1,:,:,:,0] = ccs
transforms[1,:,:,:,1] = css
return transforms
transforms = compute_general_transforms(vals)
if VERBOSE:
check_vals = np.zeros((J, K, 3))
for j in tqdm(range(J)):
for k in range(K):
for ii in range(2):
for jj in range(J//2):
for kk in range(K//2):
ray = spherical_to_cartesian(jj * thetas[j,k] if ii == 0 else np.pi/2 - jj * thetas[j,k], kk * phis[j,k])
check_vals[j,k] += transforms[ii,jj,kk] @ ray
print("Mean abs error comparing target values with sum of spheres under general transformations:", np.abs(vals - check_vals).mean())
print("Displaying samples generated as a sum of spheres under general transformations")
display_point_cloud(check_vals.reshape(-1, 3))
0%| | 0/20 [00:00<?, ?it/s]
0%| | 0/40 [00:00<?, ?it/s]
Mean abs error comparing target values with sum of spheres under general transformations: 0.001194455735189159 Displaying samples generated as a sum of spheres under general transformations
Renderer(camera=PerspectiveCamera(children=(DirectionalLight(color='white', intensity=0.6, position=(0.0766123…
Now we must convert our general transforms into orthogonal transforms. We can prove that any 3 x 3 matrix can be written as the sum of three orthogonal 3 x 3 matrices, possibly with real coefficients. (Note that two is not generally possible since that would have eight degrees of freedom while 3 x 3 matrices have nine.) Given a matrix $T$, the singular value decomposition gives $T = U \Sigma V^T$ for orthogonal $T, U$. Then note that $$R_1 = U\begin{pmatrix}1&&\\&-1&\\&&-1\end{pmatrix} V^T,\, R_2 = U\begin{pmatrix}-1&&\\&1&\\&&-1\end{pmatrix} V^T,\, R_3 = U\begin{pmatrix}-1&&\\&-1&\\&&1\end{pmatrix} V^T$$ are orthogonal as they're each the product of orthogonal matrices. Furthermore, $$T = -\left(\left(\frac{\sigma_2 + \sigma_3}{2}\right) R_1 + \left(\frac{\sigma_3 + \sigma_1}{2}\right) R_2 + \left(\frac{\sigma_1 + \sigma_2}{2}\right) R_3 \right).$$ So the proof is done and the formula is implemented below.
def general_to_orthogonal(transform):
"""
Converts general 3x3 matrix into sum of three 3x3 orthogonal matrices with coefficients
"""
u, s, vh = np.linalg.svd(transform)
s_tot = np.sum(s)
coeffs = -0.5 * (s_tot - s)
matrices = np.array([
u @ np.diag([1,-1,-1]) @ vh,
u @ np.diag([-1,1,-1]) @ vh,
u @ np.diag([-1,-1,1]) @ vh
])
return coeffs, matrices
We can finally combine our results to get our desired formula.
orthogonal_transforms = np.zeros((2, J//2, K//2, 3, 3, 3))
radii = np.zeros((2, J//2, K//2, 3))
flips = np.zeros((2, J//2, K//2, 3), dtype=np.int32)
js = np.zeros((2, J//2, K//2, 3), dtype=np.int32)
ks = np.zeros((2, J//2, K//2, 3), dtype=np.int32)
for ii in range(2):
for j in range(J//2):
for k in range(K//2):
flips[ii,j,k,:] = ii
js[ii,j,k,:] = j
ks[ii,j,k,:] = k
coeffs, matrices = general_to_orthogonal(transforms[ii,j,k])
radii[ii,j,k] = np.abs(coeffs)
orthogonal_transforms[ii,j,k] = coeffs[:,None,None] * matrices
if VERBOSE:
print("Mean abs error from converting general transforms to orthogonal ones:", np.abs(transforms - np.sum(orthogonal_transforms, axis=3)).sum())
Mean abs error from converting general transforms to orthogonal ones: 1.127887155586323e-15
"""
Filter out small spheres and sort by radii for nicer visualization
"""
flips = flips.reshape(-1,)
js = js.reshape(-1,)
ks = ks.reshape(-1,)
radii = radii.reshape(-1)
orthogonal_transforms = orthogonal_transforms.reshape(-1, 3, 3)
mask = radii > 1e-9
flips = flips[mask]
js = js[mask]
ks = ks[mask]
radii = radii[mask]
orthogonal_transforms = orthogonal_transforms[mask,:,:]
order = np.argsort(-np.abs(radii))
flips = flips[order]
js = js[order]
ks = ks[order]
radii = radii[order]
orthogonal_transforms = orthogonal_transforms[order,:,:]
save = False
if save:
np.savez(f"data/parameters_{J}_{K}.npz", flips=flips, js=js, ks=ks, radii=radii, orthogonal_transforms=orthogonal_transforms)
if VERBOSE:
check_vals = np.zeros((J, K, 3))
for j in tqdm(range(J)):
for k in range(K):
for i in range(radii.shape[0]):
ray = spherical_to_cartesian(js[i] * thetas[j,k] if flips[i] == 0 else np.pi/2 - js[i] * thetas[j,k], ks[i] * phis[j,k])
check_vals[j,k] += orthogonal_transforms[i] @ ray
print("Overall mean abs error:", np.abs(vals - check_vals).mean())
print("Displaying samples generated as a sum of spheres")
display_point_cloud(check_vals.reshape(-1, 3))
0%| | 0/40 [00:00<?, ?it/s]
Overall mean abs error: 0.0011944557351891622 Displaying samples generated as a sum of spheres
Renderer(camera=PerspectiveCamera(children=(DirectionalLight(color='white', intensity=0.6, position=(0.0766123…
We use these parameters for our web visualization.
"""
Output obtained parameters to javascript file
"""
np.set_printoptions(threshold = 2*np.prod(orthogonal_transforms.shape))
with open('data.js', 'w') as f:
f.write('export const radii_data = ' + np.array2string(radii, separator=',', precision=4))
f.write(';\n')
f.write('export const flips_data = ' + np.array2string(flips, separator=',', precision=4))
f.write(';\n')
f.write('export const js_data = ' + np.array2string(js, separator=',', precision=4))
f.write(';\n')
f.write('export const ks_data = ' + np.array2string(ks, separator=',', precision=4))
f.write(';\n')
f.write('export const orthogonal_transforms_data = ' + np.array2string(orthogonal_transforms.swapaxes(1, 2).reshape(-1, 9), separator=',', precision=4))
f.write(';\n')
We render our results using three.js and host it here: https://www.lessvrong.com/cs/spheres.
We compute the error of our reconstruction against the intended samples and against random points. Unsurprisingly, the error goes down as the number of terms increases.
def apply(parameters, thetas, phis):
output = np.zeros((thetas.shape[0], 3))
js = parameters["js"]
ks = parameters["ks"]
flips = parameters["flips"]
radii = parameters["radii"]
orthogonal_transforms = parameters["orthogonal_transforms"]
for j in tqdm(range(thetas.shape[0])):
for i in range(radii.shape[0]):
ray = spherical_to_cartesian(js[i] * thetas[j] if flips[i] == 0 else np.pi/2 - js[i] * thetas[j], ks[i] * phis[j])
output[j,:] += orthogonal_transforms[i] @ ray
return output
JKs = [
(10*i, 10*i) for i in range(1, 11)
]
parameters = []
for j, k in JKs:
parameters.append(np.load(f"data/parameters_{j}_{k}.npz"))
structured_errors = []
for i in range(len(parameters)):
thetas, phis = evenly_spaced_angles(JKs[i][0], JKs[i][1])
thetas = thetas.flatten()
phis = phis.flatten()
if thetas.shape[0] > 1000:
rand_idx = np.random.choice(thetas.shape[0], (1000,))
thetas = thetas[rand_idx]
phis = phis[rand_idx]
true_vals = generate_samples(thetas, phis, target, sphere)
pred_vals = apply(parameters[i], thetas, phis)
error = np.abs(true_vals - pred_vals).mean()
structured_errors.append(error)
fig = plt.figure()
ax = fig.add_subplot(111)
ax.plot([6*JKs[i][0]**2 for i in range(len(JKs))], structured_errors, marker="o")
ax.set_title("Error at Samples vs Number of Parameters")
ax.set_xlabel("Number of Terms")
ax.set_ylabel("Mean Abs Error")
ax.set_xscale("log")
ax.set_yscale("log")
ax.set_aspect("equal")
plt.show()
0%| | 0/100 [00:00<?, ?it/s]
0%| | 0/400 [00:00<?, ?it/s]
0%| | 0/900 [00:00<?, ?it/s]
0%| | 0/1000 [00:00<?, ?it/s]
0%| | 0/1000 [00:00<?, ?it/s]
0%| | 0/1000 [00:00<?, ?it/s]
0%| | 0/1000 [00:00<?, ?it/s]
0%| | 0/1000 [00:00<?, ?it/s]
0%| | 0/1000 [00:00<?, ?it/s]
0%| | 0/1000 [00:00<?, ?it/s]
random_errors = []
for i in range(len(parameters)):
num = 1000
thetas = np.random.uniform(0, np.pi, (num,))
phis = np.random.uniform(0, 2*np.pi, (num,))
true_vals = generate_samples(thetas, phis, target, sphere)
pred_vals = apply(parameters[i], thetas, phis)
error = np.abs(true_vals - pred_vals).mean()
random_errors.append(error)
fig = plt.figure()
ax = fig.add_subplot(111)
ax.plot([6*JKs[i][0]**2 for i in range(len(JKs))], random_errors, marker="o")
ax.set_title("Error at Random Points vs Number of Terms")
ax.set_xlabel("Number of Terms")
ax.set_ylabel("Mean Abs Error")
ax.set_xscale("log")
ax.set_yscale("log")
ax.set_aspect("equal")
plt.show()
0%| | 0/1000 [00:00<?, ?it/s]
0%| | 0/1000 [00:00<?, ?it/s]
0%| | 0/1000 [00:00<?, ?it/s]
0%| | 0/1000 [00:00<?, ?it/s]
0%| | 0/1000 [00:00<?, ?it/s]
0%| | 0/1000 [00:00<?, ?it/s]
0%| | 0/1000 [00:00<?, ?it/s]
0%| | 0/1000 [00:00<?, ?it/s]
0%| | 0/1000 [00:00<?, ?it/s]
0%| | 0/1000 [00:00<?, ?it/s]
In this project, we generalized the epicycle animation to 3D. This was possible by formulating the problem as a linear algebra problem and applying standard techniques. Ultimately, we obtained a novel surface representation, although it appears that it is not very useful without an infeasible number of terms. Finally, this was made into a web demo, which can be accessed here. Though the results are not ideal, seeing the spheres spin through space is pretty cool.
There were many different ways we could have solved our problem formulation. For example, the maps $F^{\Gamma}$ could have simply mapped to great circles on the sphere, removing a dimension of the problem. Then each coordinate could independently be solved. That would have been boring though, which is why I proposed a different set of functions. Even in 2D, the DFT approach does not have to be the only one. So what is special about it? In fact, the circular harmonics are 1D irreducible representations of $SO(2)$ and the DFT is essentially the Peter-Weyl theorem on the cyclic subgroups of $SO(2)$. On the other hand, $SO(3)$ does not readily have an infinite set of 3D irreducible representations. Furthermore, its subgroups are not convenient for sampling a surface. The elegance of the 2D experiment had to be replaced with careful consideration in 3D.