Mahalanobis Distance

In the following, we will visualize the Mahalanobis distance between two planes using a Python script.

We sample 20 points at random and organize them on two planes, centered at (0, 0.5, 0.5) and (1, 0.5, 0.5) respectively. We then rotate one of the two planes around z axis for 90 degrees, and see how the distance changes. Also, we set different noise levels on the normal direction of the planes, and see how the distance changes.

Before Rotation, Noise Level = 0.1, M-distance = 6-9

After Rotation, Noise Level = 0.1, M-distance = 2-4

Before Rotation, Noise Level = 0.001, M-distance = 600-900

After Rotation, Noise Level = 0.001, M-distance = 2-4

Furthermore, if we set noise level to 0.01, the results are still consistent. Figures omitted.

Before Rotation, Noise Level = 0.01, M-distance = 60-90
After Rotation, Noise Level = 0.01, M-distance = 2-4

If we rotate the plane slightly, such as even only 10 degrees, the distance will decrease drastically.

Before Rotation, Noise Level = 0.001, M-distance = 600-900
After Rotation, Noise Level = 0.001, M-distance = 20-40

Python Code

import numpy as np
import matplotlib.pyplot as plt
from mpl_toolkits.mplot3d import Axes3D

# Function to generate points on a plane
def generate_plane_points(mean, normal, n=20, noise_level=0.001):
    # Generate random points
    points = np.random.rand(n, 3)
    # Adjust points to lie on a given plane
    for i in range(n):
        points[i] = points[i] - np.dot(points[i] - mean, normal) * normal
        # Add some noise
        points[i] += np.random.normal(scale=noise_level, size=3)
    return points

# Function to calculate Mahalanobis distance
def mahalanobis_distance(mu1, mu2, sigma1, sigma2):
    cov_inv = np.linalg.inv(sigma1 + sigma2)
    delta = mu1 - mu2
    return np.sqrt(np.dot(np.dot(delta.T, cov_inv), delta))

# Function to plot planes
def plot_planes(plane1, plane2, ax):
    ax.scatter(plane1[:,0], plane1[:,1], plane1[:,2], color='blue', alpha=0.5)
    ax.scatter(plane2[:,0], plane2[:,1], plane2[:,2], color='red', alpha=0.5)

# Define two planes
mean1 = np.array([0, 0, 0])
normal1 = np.array([1, 0, 0])
plane1 = generate_plane_points(mean1, normal1)

mean2 = np.array([1, 0, 0])
normal2 = np.array([-1, 0, 0])
plane2 = generate_plane_points(mean2, normal2)

# Calculate means and covariances
mu1 = np.mean(plane1, axis=0)
sigma1 = np.cov(plane1, rowvar=False)

mu2 = np.mean(plane2, axis=0)
sigma2 = np.cov(plane2, rowvar=False)

# Calculate Mahalanobis distance before transformation
distance_before = mahalanobis_distance(mu1, mu2, sigma1, sigma2)
print("Mahalanobis Distance Before Transformation:", distance_before)

# Apply transformation to second plane (e.g., rotation)
theta = np.radians(90)  # Rotation angle around z axis
rotation_matrix = np.array([
    [np.cos(theta), -np.sin(theta), 0],
    [np.sin(theta), np.cos(theta), 0],
    [0, 0, 1]
])
plane2_transformed = np.dot(plane2 - mu2, rotation_matrix) + mu2

# Recalculate mean and covariance for the transformed plane
mu2_transformed = np.mean(plane2_transformed, axis=0)
sigma2_transformed = np.cov(plane2_transformed, rowvar=False)

# Calculate Mahalanobis distance after transformation
distance_after = mahalanobis_distance(mu1, mu2_transformed, sigma1, sigma2_transformed)
print("Mahalanobis Distance After Transformation:", distance_after)

# Visualization
fig = plt.figure(figsize=(10, 8))
ax = fig.add_subplot(111, projection='3d')
plot_planes(plane1, plane2_transformed, ax)
ax.set_title("Mahalanobis Distance Between Two Planes")
plt.show()

PreviousGeneralized ICP NextComputer Science

Last updated 1 year ago

Was this helpful?

import numpy as np import matplotlib.pyplot as plt from mpl_toolkits.mplot3d import Axes3D # Function to generate points on a plane def generate_plane_points(mean, normal, n=20, noise_level=0.001): # Generate random points points = np.random.rand(n, 3) # Adjust points to lie on a given plane for i in range(n): points[i] = points[i] - np.dot(points[i] - mean, normal) * normal # Add some noise points[i] += np.random.normal(scale=noise_level, size=3) return points # Function to calculate Mahalanobis distance def mahalanobis_distance(mu1, mu2, sigma1, sigma2): cov_inv = np.linalg.inv(sigma1 + sigma2) delta = mu1 - mu2 return np.sqrt(np.dot(np.dot(delta.T, cov_inv), delta)) # Function to plot planes def plot_planes(plane1, plane2, ax): ax.scatter(plane1[:,0], plane1[:,1], plane1[:,2], color='blue', alpha=0.5) ax.scatter(plane2[:,0], plane2[:,1], plane2[:,2], color='red', alpha=0.5) # Define two planes mean1 = np.array([0, 0, 0]) normal1 = np.array([1, 0, 0]) plane1 = generate_plane_points(mean1, normal1) mean2 = np.array([1, 0, 0]) normal2 = np.array([-1, 0, 0]) plane2 = generate_plane_points(mean2, normal2) # Calculate means and covariances mu1 = np.mean(plane1, axis=0) sigma1 = np.cov(plane1, rowvar=False) mu2 = np.mean(plane2, axis=0) sigma2 = np.cov(plane2, rowvar=False) # Calculate Mahalanobis distance before transformation distance_before = mahalanobis_distance(mu1, mu2, sigma1, sigma2) print("Mahalanobis Distance Before Transformation:", distance_before) # Apply transformation to second plane (e.g., rotation) theta = np.radians(90) # Rotation angle around z axis rotation_matrix = np.array([ [np.cos(theta), -np.sin(theta), 0], [np.sin(theta), np.cos(theta), 0], [0, 0, 1] ]) plane2_transformed = np.dot(plane2 - mu2, rotation_matrix) + mu2 # Recalculate mean and covariance for the transformed plane mu2_transformed = np.mean(plane2_transformed, axis=0) sigma2_transformed = np.cov(plane2_transformed, rowvar=False) # Calculate Mahalanobis distance after transformation distance_after = mahalanobis_distance(mu1, mu2_transformed, sigma1, sigma2_transformed) print("Mahalanobis Distance After Transformation:", distance_after) # Visualization fig = plt.figure(figsize=(10, 8)) ax = fig.add_subplot(111, projection='3d') plot_planes(plane1, plane2_transformed, ax) ax.set_title("Mahalanobis Distance Between Two Planes") plt.show()