Question 1

Pre-requisites

• An Signed Distance Function (SDF) is a mathematical function that gives the shortest distance from any point in space to the surface of a shape, with a sign indicating whether the point is inside or outside the shape.

• A Signed Distance Function (SDF) represents a surface implicitly and provides numerous benefits. It’s particularly useful for boolean operations, collision detection, ray marching, modeling, and physics simulations.

## Run this shell for setting up env
!pip install trimesh pyglet networkx shapely rtree
import numpy as np
import trimesh
import plotly.graph_objects as go
from skimage import measure


# Utility funcs for visualization
def get_pcl_trace(surface_points,size,color,opacity=1):
    points_trace = go.Scatter3d(
        x=surface_points[:, 0],
        y=surface_points[:, 1],
        z=surface_points[:, 2],
        mode='markers',
        marker=dict(
            size=size,
            color=color,
            opacity=0.8
        )
    )
    return points_trace

def get_mesh_trace(mesh,color,opacity=1):
    vertices = mesh.vertices
    faces = mesh.faces

    mesh_trace = go.Mesh3d(
        x=vertices[:, 0],
        y=vertices[:, 1],
        z=vertices[:, 2],
        i=faces[:, 0],
        j=faces[:, 1],
        k=faces[:, 2],
        color=color,
        opacity=opacity,
    )
    return mesh_trace

def get_cone_trace(surface_points,surface_normals,size=1):
    sn = surface_normals / np.linalg.norm(surface_normals, axis=1, keepdims=True)
    cone_trace = go.Cone(
        x=surface_points[:, 0],
        y=surface_points[:, 1],
        z=surface_points[:, 2],
        u=sn[:, 0],
        v=sn[:, 1],
        w=sn[:, 2],
        sizemode="absolute",
        sizeref=size,
        anchor='tail',
        colorscale=[[0, 'blue'], [1, 'blue']],
        showscale=False,
    )
    return cone_trace

def view_trace(data):
    fig = go.Figure(data=data)
    fig.update_layout(
        scene=dict(
            xaxis_title='X',
            yaxis_title='Y',
            zaxis_title='Z'
        ),
        margin=dict(r=0, l=0, b=0, t=0)
    )

    fig.show()

Requirement already satisfied: trimesh in /usr/local/lib/python3.11/dist-packages (4.6.4)
Requirement already satisfied: pyglet in /usr/local/lib/python3.11/dist-packages (2.1.3)
Requirement already satisfied: networkx in /usr/local/lib/python3.11/dist-packages (3.4.2)
Requirement already satisfied: shapely in /usr/local/lib/python3.11/dist-packages (2.0.7)
Requirement already satisfied: rtree in /usr/local/lib/python3.11/dist-packages (1.4.0)
Requirement already satisfied: numpy>=1.20 in /usr/local/lib/python3.11/dist-packages (from trimesh) (1.26.4)

Q 1.a

Define the implicit Signed Distance Function (SDF) for a sphere centered at the origin with a radius of 0.5. The function should accept a point as input and return its signed distance to the sphere’s surface — positive if the point is outside, negative if it’s inside.

import numpy as np
def sphere_sdf(point):
    center =  np.array([0, 0, 0])
    radius = 0.5
    ##
    # Write your code here
    x = np.linalg.norm(point-center)
    if x > radius:
        sdf_at_point = -1 * (radius - x)
    else:
        sdf_at_point = (x - radius)
    ## Could have had directly written
    # sdf_at_point = np.linalg.norm(point-center) - radius
    # Just for explanation and readability wrote in expanded form
    ## Remove the following
    return sdf_at_point

Visualization

Let’s create a 64x64x64 voxel grid for a unit cube centered at the origin, with a side length of 2, ranging from -1 to 1 and evaluate sdf at all voxels.

We can determine if a voxel lies on the surface by checking if its SDF value is close to zero (within a given threshold). The provided code can help visualize all voxels-- in red, and the points on (near) the surface-- in blue.

Learn about marching cubes algorithm-- it allows us to extract the surface implicitly represented by SDF.

grid_size = 64
x = np.linspace(-1, 1, grid_size)
y = np.linspace(-1, 1, grid_size)
z = np.linspace(-1, 1, grid_size)

X, Y, Z = np.meshgrid(x, y, z, indexing='ij')
points = np.stack((X, Y, Z), axis=-1)
sdf_values = np.apply_along_axis(sphere_sdf, -1, points)

# Visualize Grid (Red) and Points near sphere surface (Blue).

pt = get_pcl_trace(points.reshape(-1, 3),size=1,color='red',opacity=0.3)

surface_threshold = 0.02
surface_points = points[np.abs(sdf_values) < surface_threshold]
spt = get_pcl_trace(surface_points.reshape(-1, 3),size=5,color='blue',opacity=1)

data = [pt,spt]
view_trace(data)

# Visualize Mesh
verts, faces, _, _ = measure.marching_cubes(sdf_values, level=0, spacing=(2/grid_size, 2/grid_size, 2/grid_size))
verts = verts - 1
mesh = trimesh.Trimesh(vertices=verts, faces=faces)
mt = get_mesh_trace(mesh,color='red',opacity=0.8)
vts = get_pcl_trace(mesh.vertices,size=2,color='yellow',opacity=1)
data = [mt,vts]
view_trace(data)

Q 1.b

Calculate the normals for the surface represented by SDF. A normal is a vector that is perpendicular to the surface at a given point, typically pointing outward from the object’s surface.

Hint: Consider calculating the spatial gradient of the SDF using finite differences for voxel grid and filter normals belonging to surface_points

Why does computing the spatial gradient of an SDF using finite differences at surface points yield the surface normals? Describe elaboratively and provide code in following cell.

###
# Add explanation here as comment
# Gradient -----> direction of steepest increase

# Gradient of a scalar and Divergence of a vector (mathematically)

# For any surface defined by f(x, y, z) = c (with c being a constant),
# the gradient vector ∇f = [∇f/∇x, ∇f/∇y, ∇f/∇z]
# This same gradient is also orthogonal to the surfaces where f remains constant.

# To illustrate, take any curve r(t) that lies on the surface. r : R -----> R^3
# Because the curve is entirely on the surface, f(r(t)) stays constant at c for all t.

# Differentiating f(r(t)) with respect to t gives:
#   d/dt [f(r(t))] = 0

# By the chain rule:
#   d/dt [f(r(t))] = ∇f(r(t)) ⋅ dr/dt

# Since the derivative is zero, we get:
#   ∇f ⋅ dr/dt = 0

# This result shows that the gradient ∇f is perpendicular to the tangent vector dr/dt of the curve.(because dot prod =0)
# Since this applies to every possible curve on the surface,
# it follows that ∇f is perpendicular to the entire tangent plane at that point.


###


def sdf_gradient(sdf_values, surface_points=None):
    ##
    # Write your code here, change input args as required
    # Computing gradients of the signed distance field (SDF) using central differences.
    # Gradient approximation:
    # ∇f(x, y, z) ≈ [(f(x + h, y, z) - f(x - h, y, z)) / 2h,
    #                (f(x, y + h, z) - f(x, y - h, z)) / 2h,
    #                (f(x, y, z + h) - f(x, y, z - h)) / 2h]

    # grad_x, grad_y, grad_z store the partial derivatives of the SDF along x, y, and z directions.
    grad_x= np.zeros_like(sdf_values)
    grad_y= np.zeros_like(sdf_values)
    grad_z= np.zeros_like(sdf_values)

    spacing = 2/64

    grad_x[1:-1, :, :] = (sdf_values[2:, :, :] - sdf_values[:-2, :, :]) / (2 * spacing)
    grad_y[:, 1:-1, :] = (sdf_values[:, 2:, :] - sdf_values[:, :-2, :]) / (2 * spacing)
    grad_z[:, :, 1:-1] = (sdf_values[:, :, 2:] - sdf_values[:, :, :-2]) / (2 * spacing)

    grad_x[0,:,:] = (sdf_values[1,:,:] - sdf_values[-1,:,:]) / (2 * spacing)
    grad_x[-1, :, :] = (sdf_values[0, :, :] - sdf_values[-2, :, :]) / (2 * spacing)

    grad_y[:,0,:] = (sdf_values[:,1,:] - sdf_values[:,-1,:]) / (2 * spacing)
    grad_y[:, -1, :] = (sdf_values[:, 0, :] - sdf_values[:, -2, :]) / (2 * spacing)

    grad_z[:,:,0] = (sdf_values[:,:,1] - sdf_values[:,:,-1]) / (2 * spacing)
    grad_z[:, :, -1] = (sdf_values[:, :, 0] - sdf_values[:, :, -2]) / (2 * spacing)


            # Create a grid of indices
    indices = np.round(((surface_points + 1) * (grid_size - 1) / 2)).astype(int)
    # Extract gradients at surface points
    normals = np.zeros_like(surface_points)
    for i in range(len(surface_points)):
        x, y, z = indices[i]
        # Combine gradient components into a vector
        normal = np.array([grad_x[x, y, z], grad_y[x, y, z], grad_z[x, y, z]])
        # Normalize to get unit vector
        normals[i] = normal / (np.linalg.norm(normal) + 1e-10)  # Add small epsilon to avoid division by zero

    return normals


surface_normals = sdf_gradient(sdf_values,surface_points=surface_points)
ct = get_cone_trace(surface_points,surface_normals,size=1) # Try changing size here, if cones too large or invisible!
data = [ct]
view_trace(data)

Visualizing Normals

## Illustration for plotting normals of extracted mesh by marching cubes
vertex_normals = mesh.vertex_normals
vertices = mesh.vertices
ct = get_cone_trace(vertices,vertex_normals,size=1)  # Try changing size here, if cones too large or invisible!
data = [ct,mt]
view_trace(data)

Question 2

How can Singular Value Decomposition (SVD) be used to compress and reconstruct images efficiently, and what role do singular values play in preserving image features?

• Reconstruct the image using only the top k singular values.
• Experiment with different values of k and visualize the results.
• Plot reconstruction error by measuring Mean squared error between reconstructed and original image for different values of k.
• Write a brief note on how this can be beneficial for image compression applications.

Singular Value Decomposition (SVD) for Image Compression

Singular Value Decomposition (SVD) is used for image compression by decomposing an image into U, S, and V matrices. The S matrix contains singular values in descending order, representing the importance of image features. By retaining only the top k singular values, you can compress the image while preserving its essential features.

Key Points:

• Decomposition: $$ A = U \Sigma V^T $$
• Compression: Retain top k singular values in Σ
• Role of Singular Values: Preserve important features by filtering out less significant ones.

import numpy as np
import matplotlib.pyplot as plt
import imageio.v2 as imageio
from skimage.color import rgb2gray

url = "https://upload.wikimedia.org/wikipedia/commons/a/af/Golden_retriever_eating_pigs_foot.jpg"  # If this image becomes unavailable (rare!) read some other from web
image = imageio.imread(url)

gray_image = rgb2gray(image)

plt.figure(figsize=(6, 6))
plt.imshow(gray_image, cmap='gray')
plt.title('Original Grayscale Image')
plt.axis('off')
plt.show()

def compress_image_using_topk_singular_values(image, k):
    ## Your code here
    U, Sigma, Vt = np.linalg.svd(image)
    U_k = U[:, :k]
    Sigma_k = np.diag(Sigma[:k])
    Vt_k = Vt[:k, :]
    image = U_k@Sigma_k@ Vt_k

    compressed_image= image
    return compressed_image

def reconstruction_error(original, reconstructed):
    mse = np.mean((original - reconstructed) ** 2)
    return mse
    ## Your code here

m, n = gray_image.shape
k_values = np.linspace(1, min(m, n), 15, dtype=int)  # Lets visualize 15 images
plt.figure(figsize=(12, 12))
for i, k in enumerate(k_values):
    compressed_image = compress_image_using_topk_singular_values(gray_image, k)
    re = reconstruction_error(gray_image, compressed_image)
    plt.subplot(3, 5, i + 1)
    plt.imshow(compressed_image, cmap='gray')
    plt.title(f'k = {k} , RE = {re:.4f}')
    plt.axis('off')

plt.tight_layout()
plt.show()

$$ MSE = \frac{1}{mn} \sum_{i=1}^{m} \sum_{j=1}^{n} \left( A_{i,j} - (A_k)_{i,j} \right)^2 $$

# Plot Reconstruction error w.r.t. 'k'
k_values = np.linspace(1, min(m, n), 100, dtype=int)  # Lets sample more 100 k values for analysing how k affects reconstruction quality.
## Your code here for plot

def reconstruction_error(original, reconstructed):
    mse = np.mean((original - reconstructed) ** 2)
    return mse

reconstruction_errors = []
for k in k_values:
    compressed_image = compress_image_using_topk_singular_values(gray_image, k)
    re = reconstruction_error(gray_image, compressed_image)
    reconstruction_errors.append(re)

plt.plot(k_values, reconstruction_errors, marker='o')
plt.xlabel('k')
plt.ylabel('Reconstruction Error')
plt.title('Reconstruction Error vs. k')

Text(0.5, 1.0, 'Reconstruction Error vs. k')

Write a brief note on

• (a) how this can be beneficial for image compression applications
• (b) how k affects reconstruction quality.

Advantages for Image Compression

Utilizing Singular Value Decomposition (SVD) for image compression offers several benefits:

• Efficient Storage: By retaining only a subset of singular values and their corresponding vectors, significant storage savings are achieved without compromising essential image features.

• Progressive Transmission: Transmitting singular values in descending order of significance allows for an initial quick display of a recognizable image, with finer details progressively refined.

• Noise Reduction: Eliminating smaller singular values, which often represent noise, enhances overall image clarity.

Influence of ( k ) on Reconstruction Quality

The parameter ( k ), representing the number of singular values retained, directly impacts the quality of the reconstructed image:

• Low ( k ) Values: Capturing the most critical image features with minimal storage, resulting in rapid quality improvements.

• High ( k ) Values: As ( k ) increases, each additional singular value contributes less to visual quality, leading to diminishing returns in reconstruction improvement.

• Optimal ( k ) Selection: Balancing ( k ) is crucial to achieve desired compression efficiency while maintaining acceptable image fidelity.