source
/
nenuscanner
3
0
Fork 0
Code Issues Pull Requests Packages Projects Releases Wiki Activity

Cleaning

This commit is contained in:
Thomas Forgione 2024-06-18 10:00:54 +02:00
parent 242595e119
commit 3d5d0a6d32
14 changed files with 589 additions and 745 deletions
Show Stats Download Patch File Download Diff File Expand all files Collapse all files

90
calibration.py
View File

@ -4,25 +4,99 @@ import functools
import numpy as np
import os
import sys
from PIL import Image
import utils
def print_error(msg: str) -> None:
# To extract a few images and resize them at 20% of their size:
# for file in ALL/led_00[0-9]0*; do; magick $file -resize 20% SMALL/$(basename $file); done
def print_error(msg: str):
print('\x1b[1;31m[ERR]' + msg + '\x1b[0m', file=sys.stderr)
def calibrate(input_dir: str):
images = [np.asarray(Image.open(os.path.join(input_dir, x))) for x in os.listdir(input_dir)]
# Camera parameters
nu, nv, nc = images[0].shape
nspheres = 5
focal_mm = 35
matrix_size = 24
focal_pix = nu * focal_mm / matrix_size
K = utils.build_K_matrix(focal_pix, nu/2, nv/2)
# Max image: image of brightest pixels, helps spheres segmentation
max_image = functools.reduce(np.maximum, images)
# Normalize and reshape the image pixels for Gaussian Mixture Model (GMM) estimation
pixels = np.reshape(max_image / 255.0, (-1, 3))
# Initialize parameters for GMM
init_params = np.ones(2), np.broadcast_to(np.eye(3) * 0.1, (2, 3, 3)), np.asarray([[0, 0, 0], [1, 1, 1]])
# Estimate GMM parameters and classify pixels
estimated_params = utils.gaussian_mixture_estimation(pixels, init_params, it=10)
classif = np.asarray(utils.maximum_likelihood(pixels, estimated_params), dtype=bool)
# Refine classification to select the appropriate binary mask
rectified_classif = utils.select_binary_mask(classif, lambda mask: np.mean(pixels[mask]))
# Identify the largest connected components (spheres) and extract their borders
sphere_masks = utils.get_greatest_components(np.reshape(rectified_classif, (nu, nv)), nspheres)
border_masks = np.vectorize(utils.get_mask_border, signature='(u,v)->(u,v)')(sphere_masks)
# Fit quadratic forms (ellipses) to the borders
def fit_on_mask(border):
return utils.fit_quadratic_form(utils.to_homogeneous(np.argwhere(border)))
ellipse_quadratics = np.vectorize(fit_on_mask, signature='(u,v)->(t,t)')(border_masks)
# Calibrate the ellipses using the camera intrinsic matrix
calibrated_quadratics = np.swapaxes(K, -1, -2) @ ellipse_quadratics @ K
# Deproject the ellipse quadratics to sphere centers
sphere_centers = utils.deproject_ellipse_to_sphere(calibrated_quadratics, 1)
# Create coordinates and calculate camera rays
coordinates = np.stack(np.meshgrid(range(nu), range(nv), indexing='ij'), axis=-1)
rays = utils.get_camera_rays(coordinates, K)
# Find the intersections between the camera rays and the spheres
sphere_points_map, sphere_geometric_masks = \
utils.line_sphere_intersection(sphere_centers[:, np.newaxis, np.newaxis, :], 1, rays[np.newaxis, :, :, :])
sphere_points = np.asarray([sphere_points_map[i, sphere_geometric_masks[i]] for i in range(nspheres)], dtype=object)
sphere_normals = np.vectorize(utils.sphere_intersection_normal, signature='(v),()->()', otypes=[object])(sphere_centers, sphere_points)
# Load grey values from images for the identified sphere regions
def to_grayscale(image):
return [np.mean(image, axis=-1)[sphere_geometric_masks[i]] / 255.0 for i in range(nspheres)]
grey_values = np.asarray(list(map(to_grayscale, images)), dtype=object)
# Estimate lighting conditions from sphere normals and grey values
estimated_lights = np.vectorize(
utils.estimate_light,
excluded=(2,),
signature='(),()->(k)',
otypes=[float]
)(sphere_normals, grey_values, (0.1, 0.9))
# Calculate the positions of the light sources
light_positions = utils.lines_intersections(sphere_centers, estimated_lights)
return light_positions
def main():
if len(sys.argv) < 2:
print_error('Expected path to images as argument')
sys.exit(1)
# Load images
input_dir = sys.argv[1]
images = [np.asarray(Image.open(os.path.join(input_dir, x))) for x in os.listdir(input_dir)]
# Max image
max_image = functools.reduce(np.maximum, images)
print(calibrate(sys.argv[1]))
if __name__ == '__main__':

507
utils.py Normal file
View File

@ -0,0 +1,507 @@
import numpy as np
import scipy.ndimage as ndimage
def dot_product(v1, v2):
"""Computes the dot product between two arrays of vectors.
Args:
v1 (Array ..., ndim): First array of vectors.
v2 (Array ..., ndim): Second array of vectors.
Returns:
Array ...: Dot product between v1 and v2.
"""
result = np.einsum('...i,...i->...', v1, v2)
return result
def norm_vector(v):
"""computes the norm and direction of vectors
Args:
v (Array ..., dim): vectors to compute the norm and direction for
Returns:
Array ...: norms of the vectors
Array ..., dim: unit direction vectors
"""
norm = np.linalg.norm(v, axis=-1)
direction = v/norm[..., np.newaxis]
return norm, direction
def to_homogeneous(v):
"""converts vectors to homogeneous coordinates
Args:
v (Array ..., dim): input vectors
Returns:
Array ..., dim+1: homogeneous coordinates of the input vectors
"""
append_term = np.ones(np.shape(v)[:-1] + (1,))
homogeneous = np.append(v, append_term, axis=-1)
return homogeneous
def cross_to_skew_matrix(v):
"""converts a vector cross product to a skew-symmetric matrix multiplication
Args:
v (Array ..., 3): vectors to convert
Returns:
Array ..., 3, 3: matrices corresponding to the input vectors
"""
indices = np.asarray([[-1, 2, 1], [2, -1, 0], [1, 0, -1]])
signs = np.asarray([[0, -1, 1], [1, 0, -1], [-1, 1, 0]])
skew_matrix = v[..., indices] * signs
return skew_matrix
def build_K_matrix(focal_length, u0, v0):
"""
Build the camera intrinsic matrix.
Parameters:
focal_length (float): Focal length of the camera.
u0 (float): First coordinate of the principal point.
v0 (float): Seccond coordinate of the principal point.
Returns:
numpy.ndarray: Camera intrinsic matrix (3x3).
"""
K = np.asarray([[focal_length, 0, u0],
[0, focal_length, v0],
[0, 0, 1]])
return K
def get_camera_rays(points, K):
"""Computes the camera rays for a set of points given the camera matrix K.
Args:
points (Array ..., 2): Points in the image plane.
K (Array 3, 3): Camera intrinsic matrix.
Returns:
Array ..., 3: Camera rays corresponding to the input points.
"""
homogeneous = to_homogeneous(points)
inv_K = np.linalg.inv(K)
rays = np.einsum('ij,...j->...i', inv_K, homogeneous)
return rays
def matrix_kernel(A):
"""Computes the eigenvector corresponding to the smallest eigenvalue of the matrix A.
Args:
A (Array ..., n, n): Input square matrix.
Returns:
Array ..., n: Eigenvector corresponding to the smallest eigenvalue.
"""
eigval, eigvec = np.linalg.eig(A)
min_index = np.argmin(np.abs(eigval), axis=-1)
min_eigvec = np.take_along_axis(eigvec, min_index[..., None, None], -1)[..., 0]
normed_eigvec = norm_vector(min_eigvec)[1]
return normed_eigvec
def evaluate_bilinear_form(Q, left, right):
"""evaluates bilinear forms at several points
Args:
Q (Array ...,ldim,rdim): bilinear form to evaluate
left (Array ...,ldim): points where the bilinear form is evaluated to the left
right (Array ...,rdim): points where the bilinear form is evaluated to the right
Returns:
Array ... bilinear forms evaluated
"""
result = np.einsum('...ij,...i,...j->...', Q, left, right)
return result
def evaluate_quadratic_form(Q, points):
"""evaluates quadratic forms at several points
Args:
Q (Array ...,dim,dim): quadratic form to evaluate
points (Array ...,dim): points where the quadratic form is evaluated
Returns:
Array ... quadratic forms evaluated
"""
result = evaluate_bilinear_form(Q, points, points)
return result
def merge_quadratic_to_homogeneous(Q, b, c):
"""merges quadratic form, linear term, and constant term into a homogeneous matrix
Args:
Q (Array ..., dim, dim): quadratic form matrix
b (Array ..., dim): linear term vector
c (Array ...): constant term
Returns:
Array ..., dim+1, dim+1: homogeneous matrix representing the quadratic form
"""
dim_points = Q.shape[-1]
stack_shape = np.broadcast_shapes(np.shape(Q)[:-2], np.shape(b)[:-1], np.shape(c))
Q_b = np.broadcast_to(Q, stack_shape + (dim_points, dim_points))
b_b = np.broadcast_to(np.expand_dims(b, -1), stack_shape+(dim_points, 1))
c_b = np.broadcast_to(np.expand_dims(c, (-1, -2)), stack_shape + (1, 1))
H = np.block([[Q_b, 0.5 * b_b], [0.5 * np.swapaxes(b_b, -1, -2), c_b]])
return H
def quadratic_to_dot_product(points):
"""computes the matrix W such that
x.T@Ax = W(x).T*A[ui,uj]
Args:
points ( Array ...,ndim): points of dimension ndim
Returns:
Array ...,ni: dot product matrix (W)
Array ni: i indices of central matrix
Array ni: j indices of central matrix
"""
dim_points = points.shape[-1]
ui, uj = np.triu_indices(dim_points)
W = points[..., ui] * points[..., uj]
return W, ui, uj
def fit_quadratic_form(points):
"""Fits a quadratic form to the given zeroes.
Args:
points (Array ..., n, dim): Input points.
Returns:
Array ..., dim, dim: Fitted quadratic form matrix.
"""
dim_points = points.shape[-1]
normed_points = norm_vector(points)[1]
W, ui, uj = quadratic_to_dot_product(normed_points)
H = np.einsum('...ki,...kj->...ij', W, W)
V0 = matrix_kernel(H)
Q = np.zeros(V0.shape[:-1] + (dim_points, dim_points))
Q[..., ui, uj] = V0
return Q
def gaussian_pdf(mu, sigma, x):
"""Computes the PDF of a multivariate Gaussian distribution.
Args:
mu (Array ...,k): Mean vector.
sigma (Array ...,k,k): Covariance matrix.
x (Array ...,k): Input vector.
Returns:
Array ...: Value of the PDF.
"""
k = np.shape(x)[-1]
Q = np.linalg.inv(sigma)
normalization = np.reciprocal(np.sqrt(np.linalg.det(sigma) * np.power(2.0 * np.pi, k)))
quadratic = evaluate_quadratic_form(Q, x - mu)
result = np.exp(-0.5 * quadratic) * normalization
return result
def gaussian_estimation(x, weights):
"""Estimates the mean and covariance matrix of a Gaussian distribution.
Args:
x (Array ...,n,dim): Data points.
weights (Array ...,n): Weights for each data point.
Returns:
Array ...,dim: Estimated mean vector.
Array ...,dim,dim: Estimated covariance matrix.
"""
weights_sum = np.sum(weights, axis=-1)
mu = np.sum(x*np.expand_dims(weights, axis=-1), axis=-2) / np.expand_dims(weights_sum, axis=-1)
centered_x = x - np.expand_dims(mu, axis=-2)
sigma = np.einsum('...s, ...si, ...sj->...ij', weights, centered_x, centered_x)/np.expand_dims(weights_sum, axis=(-1, -2))
return mu, sigma
def gaussian_mixture_estimation(x, init_params, it=100):
"""Estimates the parameters of a k Gaussian mixture model using the EM algorithm.
Args:
x (Array ..., n, dim): Data points.
init_params (tuple): Initial parameters (pi, sigma, mu).
pi (Array ..., k): Initial mixture weights.
sigma (Array ..., k, dim, dim): Initial covariance matrices.
mu (Array ..., k, dim): Initial means.
it (int, optional): Number of iterations. Defaults to 100.
Returns:
Tuple[(Array ..., k), (Array ..., k, dim, dim), (Array ..., k, dim)]:
Estimated mixture weights,covariance matrices, means.
"""
pi, sigma, mu = init_params
for _ in range(it):
pdf = gaussian_pdf(
np.expand_dims(mu, axis=-2),
np.expand_dims(sigma, axis=-3),
np.expand_dims(x, axis=-3)
) * np.expand_dims(pi, axis=-1)
weights = pdf/np.sum(pdf, axis=-2, keepdims=True)
pi = np.mean(weights, axis=-1)
mu, sigma = gaussian_estimation(x, weights)
return pi, sigma, mu
def maximum_likelihood(x, params):
"""Selects the best gaussian model for a point
Args:
x (Array ..., dim): Data points.
params (tuple): Gaussians parameters (pi, sigma, mu).
pi (Array ..., k): Mixture weights.
sigma (Array ..., k, dim, dim): Covariance matrices.
mu (Array ..., k, dim): Means.
Returns:
Array ...: integer in [0,k-1] giving the maximum likelihood model
"""
pi, sigma, mu = params
pdf = gaussian_pdf(mu, sigma, np.expand_dims(x, axis=-2))*pi
result = np.argmax(pdf, axis=-1)
return result
def get_greatest_components(mask, n):
"""
Extract the n largest connected components from a binary mask.
Parameters:
mask (Array ...): The binary mask.
n (int): The number of largest connected components to extract.
Returns:
Array n,...: A boolean array of the n largest connected components
"""
labeled, _ = ndimage.label(mask)
unique, counts = np.unique(labeled, return_counts=True)
greatest_labels = unique[unique != 0][np.argsort(counts[unique != 0])[-n:]]
greatest_components = labeled[np.newaxis, ...] == np.expand_dims(greatest_labels, axis=tuple(range(1, 1 + mask.ndim)))
return greatest_components
def get_mask_border(mask):
"""
Extract the border from a binary mask.
Parameters:
mask (Array ...): The binary mask.
Returns:
Array ...: A boolean array mask of the border
"""
inverted_mask = np.logical_not(mask)
dilated = ndimage.binary_dilation(inverted_mask)
border = np.logical_and(mask, dilated)
return border
def select_binary_mask(mask, metric):
"""Selects the side of a binary mask that optimizes the given metric.
Args:
mask (Array bool ...): Initial binary mask.
metric (function): Function to evaluate the quality of the mask.
Returns:
Array bool ...: Selected binary mask that maximizes the metric.
"""
inverted = np.logical_not(mask)
result = mask if metric(mask) > metric(inverted) else inverted
return result
def deproject_ellipse_to_sphere(M, radius):
"""finds the deprojection of an ellipse to a sphere
Args:
M (Array 3,3): Ellipse quadratic form
radius (float): radius of the researched sphere
Returns:
Array 3: solution of sphere centre location
"""
H = 0.5 * (np.swapaxes(M, -1, -2) + M)
eigval, eigvec = np.linalg.eigh(H)
i_unique = np.argmax(np.abs(np.median(eigval, axis=-1, keepdims=True) - eigval), axis=-1)
unique_eigval = np.take_along_axis(eigval, i_unique[..., None], -1)[..., 0]
unique_eigvec = np.take_along_axis(eigvec, i_unique[..., None, None], -1)[..., 0]
double_eigval = 0.5 * (np.sum(eigval, axis=-1) - unique_eigval)
z_sign = np.sign(unique_eigvec[..., -1])
dist = np.sqrt(1 - double_eigval / unique_eigval)
C = np.real(radius * (dist * z_sign)[..., None] * norm_vector(unique_eigvec)[1])
return C
def weighted_least_squares(A, y, weights):
"""Computes the weighted least squares solution of Ax=y.
Args:
A (Array ...,u,v): Design matrix.
y (Array ...,u): Target values.
weights (Array ...,u): Weights for each equation.
Returns:
Array ...,v : Weighted least squares solution.
"""
pinv = np.linalg.pinv(A * weights[..., np.newaxis])
result = np.einsum('...uv,...v->...u', pinv, y * weights)
return result
def iteratively_reweighted_least_squares(A, y, epsilon=1e-5, it=20):
"""Computes the iteratively reweighted least squares solution. of Ax=y
Args:
A (Array ..., u, v): Design matrix.
y (Array ..., u): Target values.
epsilon (float, optional): Small value to avoid division by zero. Defaults to 1e-5.
it (int, optional): Number of iterations. Defaults to 20.
Returns:
Array ..., v: Iteratively reweighted least squares solution.
"""
weights = np.ones(y.shape)
for _ in range(it):
result = weighted_least_squares(A, y, weights)
ychap = np.einsum('...uv, ...v->...u', A, result)
delta = np.abs(ychap-y)
weights = np.reciprocal(np.maximum(epsilon, np.sqrt(delta)))
return result
def lines_intersections_system(points, directions):
"""computes the system of equations for intersections of lines, Ax=b
where x is the instersection
Args:
points (Array ..., npoints, ndim): points through which the lines pass
directions (Array ..., npoints, ndim): direction vectors of the lines
Returns:
Array ..., 3*npoints, ndim: coefficient matrix A for the system of equations
Array ..., 3*npoints: right-hand side vector b for the system of equations
"""
n = norm_vector(directions)[1]
skew = np.swapaxes(cross_to_skew_matrix(n), -1, -2)
root = np.einsum('...uij, ...uj->...ui', skew, points)
A = np.concatenate(np.moveaxis(skew, -3, 0), axis=-2)
b = np.concatenate(np.moveaxis(root, -2, 0), axis=-1)
return A, b
def lines_intersections(points, directions):
"""computes the intersections of lines
Args:
points (Array ..., npoints, ndim): points through which the lines pass
directions (Array ..., npoints, ndim): direction vectors of the lines
Returns:
Array ..., ndim: intersection
"""
A, b = lines_intersections_system(points, directions)
x = iteratively_reweighted_least_squares(A, b)
return x
def line_sphere_intersection_determinant(center, radius, directions):
"""computes the determinant for the intersection of a line and a sphere,
Args:
center (Array ..., dim): center of the sphere
radius (Array ...): radius of the sphere
directions (Array ..., dim): direction of the line
Returns:
Array ...:intersection determinant
"""
directions_norm_2 = np.square(norm_vector(directions)[0])
center_norm_2 = np.square(norm_vector(center)[0])
dot_product_2 = np.square(dot_product(center, directions))
delta = dot_product_2 - directions_norm_2 * (center_norm_2 - np.square(radius))
return delta
def line_plane_intersection(normal, alpha, directions):
"""Computes the intersection points between a line and a plane.
Args:
normal (Array ..., ndim): Normal vector to the plane.
alpha (Array ...): Plane constant alpha.
directions (Array ..., dim): direction of the line
Returns:
Array ..., ndim: Intersection points between the line and the sphere.
"""
t = -alpha*np.reciprocal(dot_product(directions, normal))
intersection = directions*t[..., np.newaxis]
return intersection
def line_sphere_intersection(center, radius, directions):
"""Computes the intersection points between a line and a sphere.
Args:
center (Array ..., ndim): Center of the sphere.
radius (Array ...): Radius of the sphere.
directions (Array ..., ndim): Direction vectors of the line.
Returns:
Array ..., ndim: Intersection points between the line and the sphere.
Array bool ...: Mask of intersection points
"""
delta = line_sphere_intersection_determinant(center, radius, directions)
mask = delta > 0
directions_norm_2 = np.square(norm_vector(directions)[0])
distances = (dot_product(center, directions) - np.sqrt(np.maximum(0, delta))) * np.reciprocal(directions_norm_2)
intersection = np.expand_dims(distances, axis=-1) * directions
return intersection, mask
def sphere_intersection_normal(center, point):
"""Computes the normal vector at the intersection point on a sphere.
Args:
center (Array ..., dim): Coordinates of the sphere center.
point (Array ..., dim): Coordinates of the intersection point.
Returns:
Array ..., dim: Normal normal vector at the intersection point.
"""
vector = point - center
normal = norm_vector(vector)[1]
return normal
def estimate_light(normals, grey_levels, treshold=(0, 1)):
"""Estimates the light directions using the given normals, grey levels, and mask.
Args:
normals (Array ..., n, dim): Normal vectors.
grey_levels (Array ..., n): Grey levels corresponding to the normals.
threshold (tuple, optional): Intensity threshold for valid grey levels. Defaults to (0, 1).
Returns:
Array ..., dim: Estimated light directions.
"""
validity_mask = np.logical_and(grey_levels > treshold[0], grey_levels < treshold[1])
lights = weighted_least_squares(normals, grey_levels, validity_mask)
return lights

34
utils/camera.py
View File

@ -1,34 +0,0 @@
import numpy as np
import utils.vector_utils as vector_utils
def build_K_matrix(focal_length, u0, v0):
"""
Build the camera intrinsic matrix.
Parameters:
focal_length (float): Focal length of the camera.
u0 (float): First coordinate of the principal point.
v0 (float): Seccond coordinate of the principal point.
Returns:
numpy.ndarray: Camera intrinsic matrix (3x3).
"""
K = np.asarray([[focal_length, 0, u0],
[0, focal_length, v0],
[0, 0, 1]])
return K
def get_camera_rays(points,K):
"""Computes the camera rays for a set of points given the camera matrix K.
Args:
points (Array ..., 2): Points in the image plane.
K (Array 3, 3): Camera intrinsic matrix.
Returns:
Array ..., 3: Camera rays corresponding to the input points.
"""
homogeneous = vector_utils.to_homogeneous(points)
inv_K = np.linalg.inv(K)
rays = np.einsum('ij,...j->...i',inv_K,homogeneous)
return rays

82
utils/gaussian.py
View File

@ -1,82 +0,0 @@
import numpy as np
import utils.quadratic_forms as quadratic_forms
def gaussian_pdf(mu,sigma,x):
"""Computes the PDF of a multivariate Gaussian distribution.
Args:
mu (Array ...,k): Mean vector.
sigma (Array ...,k,k): Covariance matrix.
x (Array ...,k): Input vector.
Returns:
Array ...: Value of the PDF.
"""
k = np.shape(x)[-1]
Q = np.linalg.inv(sigma)
normalization = np.reciprocal(np.sqrt(np.linalg.det(sigma)*np.power(2.0*np.pi,k)))
quadratic = quadratic_forms.evaluate_quadratic_form(Q,x-mu)
result = np.exp(-0.5*quadratic)*normalization
return result
def gaussian_estimation(x,weights):
"""Estimates the mean and covariance matrix of a Gaussian distribution.
Args:
x (Array ...,n,dim): Data points.
weights (Array ...,n): Weights for each data point.
Returns:
Array ...,dim: Estimated mean vector.
Array ...,dim,dim: Estimated covariance matrix.
"""
weights_sum = np.sum(weights,axis=-1)
mu = np.sum(x*np.expand_dims(weights,axis=-1),axis=-2)/np.expand_dims(weights_sum,axis=-1)
centered_x = x-np.expand_dims(mu,axis=-2)
sigma = np.einsum('...s,...si,...sj->...ij',weights,centered_x,centered_x)/np.expand_dims(weights_sum,axis=(-1,-2))
return mu,sigma
def gaussian_mixture_estimation(x,init_params,it=100):
"""Estimates the parameters of a k Gaussian mixture model using the EM algorithm.
Args:
x (Array ..., n, dim): Data points.
init_params (tuple): Initial parameters (pi, sigma, mu).
pi (Array ..., k): Initial mixture weights.
sigma (Array ..., k, dim, dim): Initial covariance matrices.
mu (Array ..., k, dim): Initial means.
it (int, optional): Number of iterations. Defaults to 100.
Returns:
Tuple[(Array ..., k), (Array ..., k, dim, dim), (Array ..., k, dim)]:
Estimated mixture weights,covariance matrices, means.
"""
pi,sigma,mu = init_params
for _ in range(it):
pdf = gaussian_pdf(np.expand_dims(mu,axis=-2),
np.expand_dims(sigma,axis=-3),
np.expand_dims(x,axis=-3))*np.expand_dims(pi,axis=-1)
weights = pdf/np.sum(pdf,axis=-2,keepdims=True)
pi=np.mean(weights,axis=-1)
mu,sigma = gaussian_estimation(x,weights)
return pi,sigma,mu
def maximum_likelihood(x,params):
"""Selects the best gaussian model for a point
Args:
x (Array ..., dim): Data points.
params (tuple): Gaussians parameters (pi, sigma, mu).
pi (Array ..., k): Mixture weights.
sigma (Array ..., k, dim, dim): Covariance matrices.
mu (Array ..., k, dim): Means.
Returns:
Array ...: integer in [0,k-1] giving the maximum likelihood model
"""
pi,sigma,mu = params
pdf = gaussian_pdf(mu,sigma,np.expand_dims(x,axis=-2))*pi
result = np.argmax(pdf,axis=-1)
return result

40
utils/geometric_fit.py
View File

@ -1,40 +0,0 @@
import numpy as np
import utils.kernels as kernels
import utils.vector_utils as vector_utils
import utils.quadratic_forms as quadratic_forms
import utils.kernels as kernels
def sphere_parameters_from_points(points):
"""evaluates sphere parameters from a set of points
Args:
points (Array ... npoints ndim): points used to fit the sphere, homogeneous coordinates
Returns:
Array ... ndim: coordinates of the center of the sphere
Array ...: values of radius of the sphere
"""
homogeneous = vector_utils.to_homogeneous(points)
Q = quadratic_forms.fit_quadratic_form(homogeneous)
scale = np.mean(np.diagonal(Q[...,:-1,:-1],axis1=-2,axis2=-1))
scaled_Q = Q*np.expand_dims(np.reciprocal(scale),axis=(-1,-2))
center = -(scaled_Q[...,-1,:-1]+scaled_Q[...,:-1,-1])/2
centered_norm = vector_utils.norm_vector(center)[0]
radius = np.sqrt(np.square(centered_norm)-scaled_Q[...,-1,-1])
return center,radius
def plane_parameters_from_points(points):
"""Computes the parameters of a plane from a set of points.
Args:
points (Array ..., dim): Coordinates of the points used to define the plane.
Returns:
Array ..., dim: Normal vector to the plane.
Array ...: Plane constant alpha.
"""
homogeneous = vector_utils.to_homogeneous(points)
E = np.einsum('...ki,...kj->...ij',homogeneous,homogeneous)
L = kernels.matrix_kernel(E)
n,alpha = L[...,:-1],L[...,-1]
return n, alpha

56
utils/image_loader.py
View File

@ -1,56 +0,0 @@
from PIL import Image
import functools
from tqdm import tqdm
import numpy as np
def loader(file_list):
"""
Load images from the file list, convert them to numpy arrays, and show a progress bar.
Parameters:
file_list (list of str): List of file paths to the images.
Returns:
map: A map object containing numpy arrays of the images.
"""
return map(np.asarray, map(Image.open, tqdm(file_list)))
def load_reduce(file_list, reducer):
"""Loads and reduces a list of image files using a reduction function.
Args:
file_list (list of str): List of paths to image files.
reducer (function): Function to reduce the images.
Returns:
array: Reduced image.
"""
reduced_image = functools.reduce(reducer, loader(file_list))
return reduced_image
def load_map(file_list, mapper):
"""Loads and maps a list of image files using a mapping function.
Args:
file_list (list of str): List of paths to image files.
mapper (function): Function to map the images.
Returns:
List of array: Mapped images.
"""
mapped_images = map(mapper, loader(file_list))
return mapped_images
def load_max_image(file_list):
"""Loads a list of image files and computes the element-wise maximum.
Args:
file_list (list of str): List of paths to image files.
Returns:
array: Image with the element-wise maximum values.
"""
max_image = load_reduce(file_list, np.maximum)
return max_image

103
utils/intersections.py
View File

@ -1,103 +0,0 @@
import numpy as np
import utils.vector_utils as vector_utils
import utils.kernels as kernels
def lines_intersections_system(points,directions):
"""computes the system of equations for intersections of lines, Ax=b
where x is the instersection
Args:
points (Array ..., npoints, ndim): points through which the lines pass
directions (Array ..., npoints, ndim): direction vectors of the lines
Returns:
Array ..., 3*npoints, ndim: coefficient matrix A for the system of equations
Array ..., 3*npoints: right-hand side vector b for the system of equations
"""
n = vector_utils.norm_vector(directions)[1]
skew = np.swapaxes(vector_utils.cross_to_skew_matrix(n),-1,-2)
root = np.einsum('...uij,...uj->...ui',skew,points)
A = np.concatenate(np.moveaxis(skew,-3,0),axis=-2)
b = np.concatenate(np.moveaxis(root,-2,0),axis=-1)
return A,b
def lines_intersections(points,directions):
"""computes the intersections of lines
Args:
points (Array ..., npoints, ndim): points through which the lines pass
directions (Array ..., npoints, ndim): direction vectors of the lines
Returns:
Array ..., ndim: intersection
"""
A,b = lines_intersections_system(points,directions)
x = kernels.iteratively_reweighted_least_squares(A,b)
return x
def line_sphere_intersection_determinant(center,radius,directions):
"""computes the determinant for the intersection of a line and a sphere,
Args:
center (Array ..., dim): center of the sphere
radius (Array ...): radius of the sphere
directions (Array ..., dim): direction of the line
Returns:
Array ...:intersection determinant
"""
directions_norm_2 = np.square(vector_utils.norm_vector(directions)[0])
center_norm_2 = np.square(vector_utils.norm_vector(center)[0])
dot_product_2 = np.square(vector_utils.dot_product(center,directions))
delta = dot_product_2-directions_norm_2*(center_norm_2-np.square(radius))
return delta
def line_plane_intersection(normal,alpha,directions):
"""Computes the intersection points between a line and a plane.
Args:
normal (Array ..., ndim): Normal vector to the plane.
alpha (Array ...): Plane constant alpha.
directions (Array ..., dim): direction of the line
Returns:
Array ..., ndim: Intersection points between the line and the sphere.
"""
t = -alpha*np.reciprocal(vector_utils.dot_product(directions,normal))
intersection = directions*t[...,np.newaxis]
return intersection
def line_sphere_intersection(center,radius,directions):
"""Computes the intersection points between a line and a sphere.
Args:
center (Array ..., ndim): Center of the sphere.
radius (Array ...): Radius of the sphere.
directions (Array ..., ndim): Direction vectors of the line.
Returns:
Array ..., ndim: Intersection points between the line and the sphere.
Array bool ...: Mask of intersection points
"""
delta = line_sphere_intersection_determinant(center,radius,directions)
mask = delta>0
dot_product = vector_utils.dot_product(center,directions)
directions_norm_2 = np.square(vector_utils.norm_vector(directions)[0])
distances = (dot_product-np.sqrt(np.maximum(0,delta)))*np.reciprocal(directions_norm_2)
intersection = np.expand_dims(distances,axis=-1)*directions
return intersection,mask
def sphere_intersection_normal(center,point):
"""Computes the normal vector at the intersection point on a sphere.
Args:
center (Array ..., dim): Coordinates of the sphere center.
point (Array ..., dim): Coordinates of the intersection point.
Returns:
Array ..., dim: Normal normal vector at the intersection point.
"""
vector = point-center
normal = vector_utils.norm_vector(vector)[1]
return normal

83
utils/kernels.py
View File

@ -1,83 +0,0 @@
import numpy as np
import utils.vector_utils as vector_utils
def weighted_least_squares(A,y,weights):
"""Computes the weighted least squares solution of Ax=y.
Args:
A (Array ...,u,v): Design matrix.
y (Array ...,u): Target values.
weights (Array ...,u): Weights for each equation.
Returns:
Array ...,v : Weighted least squares solution.
"""
pinv = np.linalg.pinv(A*weights[...,np.newaxis])
result = np.einsum('...uv,...v->...u',pinv,y*weights)
return result
def least_squares(A,y):
"""Computes the least squares solution of Ax=y.
Args:
A (Array ...,u,v): Design matrix.
y (Array ...,u): Target values.
Returns:
Array ...,v : Least squares solution.
"""
result = weighted_least_squares(A,y,np.ones(A.shape[0]))
return result
def iteratively_reweighted_least_squares(A,y, epsilon=1e-5, it=20):
"""Computes the iteratively reweighted least squares solution. of Ax=y
Args:
A (Array ..., u, v): Design matrix.
y (Array ..., u): Target values.
epsilon (float, optional): Small value to avoid division by zero. Defaults to 1e-5.
it (int, optional): Number of iterations. Defaults to 20.
Returns:
Array ..., v: Iteratively reweighted least squares solution.
"""
weights = np.ones(y.shape)
for _ in range(it):
result = weighted_least_squares(A,y,weights)
ychap = np.einsum('...uv,...v->...u',A,result)
delta = np.abs(ychap-y)
weights = np.reciprocal(np.maximum(epsilon,np.sqrt(delta)))
return result
def matrix_kernel(A):
"""Computes the eigenvector corresponding to the smallest eigenvalue of the matrix A.
Args:
A (Array ..., n, n): Input square matrix.
Returns:
Array ..., n: Eigenvector corresponding to the smallest eigenvalue.
"""
eigval, eigvec = np.linalg.eig(A)
min_index = np.argmin(np.abs(eigval),axis=-1)
min_eigvec = np.take_along_axis(eigvec,min_index[...,None,None],-1)[...,0]
normed_eigvec = vector_utils.norm_vector(min_eigvec)[1]
return normed_eigvec
def masked_least_squares(A,y,mask):
"""Computes the least squares solution of Ax = y for masked data.
Args:
A (Array ..., n, p): Design matrix.
y (Array ..., n): Target values.
mask (Array ..., n, bool): Mask to select valid data points.
Returns:
Array ..., p: Least squares solution for the masked data.
"""
masked_solver = lambda A,y,mask : least_squares(A[mask,:],y[mask])
vectorized = np.vectorize(masked_solver,signature='(n,p),(n),(n)->(p)')
result = vectorized(A,y,mask)
return result

35
utils/mashal.py
View File

@ -1,35 +0,0 @@
import numpy as np
def marshal_arrays(arrays):
"""
Flatten a list of numpy arrays and store their shapes.
Parameters:
arrays (list of np.ndarray): List of numpy arrays to be marshalled.
Returns:
tuple: A tuple containing:
- flat (np.ndarray): A single concatenated numpy array of all elements.
- shapes (list of tuple): A list of shapes of the original arrays.
"""
flattened = list(map(lambda a : np.reshape(a,-1),arrays))
shapes = list(map(np.shape,arrays))
flat = np.concatenate(flattened)
return flat, shapes
def unmarshal_arrays(flat,shapes):
"""
Rebuild the original list of numpy arrays from the flattened array and shapes.
Parameters:
flat (np.ndarray): The single concatenated numpy array of all elements.
shapes (list of tuple): A list of shapes of the original arrays.
Returns:
list of np.ndarray: The list of original numpy arrays.
"""
sizes = list(map(np.prod,shapes))
splits = np.cumsum(np.asarray(sizes,dtype=int))[:-1]
flattened = np.split(flat,splits)
arrays = list(map(lambda t : np.reshape(t[0],t[1]),zip(flattened,shapes)))
return arrays

48
utils/masks.py
View File

@ -1,48 +0,0 @@
import numpy as np
import scipy.ndimage as ndimage
def get_greatest_components(mask, n):
"""
Extract the n largest connected components from a binary mask.
Parameters:
mask (Array ...): The binary mask.
n (int): The number of largest connected components to extract.
Returns:
Array n,...: A boolean array of the n largest connected components
"""
labeled, _ = ndimage.label(mask)
unique, counts = np.unique(labeled, return_counts=True)
greatest_labels = unique[unique != 0][np.argsort(counts[unique != 0])[-n:]]
greatest_components = labeled[np.newaxis,...] == np.expand_dims(greatest_labels,axis=tuple(range(1,1+mask.ndim)))
return greatest_components
def get_mask_border(mask):
"""
Extract the border from a binary mask.
Parameters:
mask (Array ...): The binary mask.
Returns:
Array ...: A boolean array mask of the border
"""
inverted_mask = np.logical_not(mask)
dilated = ndimage.binary_dilation(inverted_mask)
border = np.logical_and(mask,dilated)
return border
def select_binary_mask(mask,metric):
"""Selects the side of a binary mask that optimizes the given metric.
Args:
mask (Array bool ...): Initial binary mask.
metric (function): Function to evaluate the quality of the mask.
Returns:
Array bool ...: Selected binary mask that maximizes the metric.
"""
inverted = np.logical_not(mask)
result = mask if metric(mask)>metric(inverted) else inverted
return result

66
utils/photometry.py
View File

@ -1,66 +0,0 @@
import numpy as np
import utils.kernels as kernels
import utils.vector_utils as vector_utils
def estimate_light(normals,grey_levels, treshold = (0,1)):
"""Estimates the light directions using the given normals, grey levels, and mask.
Args:
normals (Array ..., n, dim): Normal vectors.
grey_levels (Array ..., n): Grey levels corresponding to the normals.
threshold (tuple, optional): Intensity threshold for valid grey levels. Defaults to (0, 1).
Returns:
Array ..., dim: Estimated light directions.
"""
validity_mask = np.logical_and(grey_levels>treshold[0],grey_levels<treshold[1])
lights = kernels.weighted_least_squares(normals,grey_levels,validity_mask)
return lights
def geometric_shading_parameters(light_point, principal_directions, points):
"""Computes geometric parameters for shading based on light source and points.
Args:
light_point (Array ..., dim): Coordinates of the light source.
principal_directions (Array ..., dim): Principal directions of each light.
points (Array ..., dim): Coordinates of the points on the surface.
Returns:
Array ..., dim: Distances from each point to the light source.
Array ...: Computed light direction at each point.
Array ...: Angular factors based on the principal directions and light direction.
"""
distance, light_direction = vector_utils.norm_vector(light_point-points)
angular_factor = np.maximum(vector_utils.dot_product(principal_directions,-light_direction),0)
return distance, light_direction, angular_factor
def estimate_anisotropy(light_point, principal_directions, points,normals, grey_levels, min_grey_level = 0.1, min_dot_product = 0.2):
"""Estimates anisotropy parameters based on geometric shading and grey levels.
Args:
light_point (Array ..., dim): Coordinates of the light source.
principal_directions (Array ..., dim): Principal directions of each light.
points (Array ..., dim): Coordinates of the points on the surface.
normals (Array ..., dim): Normal vectors at each point.
grey_levels (Array ...): Observed grey levels at each point.
min_grey_level (float, optional): Minimum valid grey level. Defaults to 0.1.
min_dot_product (float, optional): Minimum valid dot product for shading and angular factors. Defaults to 0.2.
Returns:
Array ..., 1: Estimated anisotropy parameter mu.
Array ..., 1: Estimated flux parameter.
"""
distance, light_direction, angular_factor = geometric_shading_parameters(light_point, principal_directions, points)
computed_shading = np.maximum(vector_utils.dot_product(normals,light_direction),0)
validity_mask = np.logical_and(grey_levels>min_grey_level,np.logical_and(computed_shading>min_dot_product,angular_factor>min_dot_product))
log_flux = np.log(np.maximum(grey_levels,min_grey_level)*np.square(distance)*np.reciprocal(np.maximum(computed_shading,min_dot_product)))
log_factor = vector_utils.to_homogeneous(np.expand_dims(np.log(np.maximum(angular_factor,min_dot_product)),axis=-1))
eta = kernels.weighted_least_squares(log_factor,log_flux,validity_mask)
mu,log_phi = eta[...,0], eta[...,1]
estimated_flux = np.exp(log_phi)
return mu,estimated_flux
def light_conditions(light_point, principal_directions, points, mu, flux):
distance, light_direction, angular_factor = geometric_shading_parameters(light_point, principal_directions, points)
light_conditions = light_direction*(np.reciprocal(np.square(distance))*np.power(angular_factor,mu)*flux)[...,np.newaxis]
return light_conditions

98
utils/quadratic_forms.py
View File

@ -1,98 +0,0 @@
import numpy as np
import utils.kernels as kernels
import utils.vector_utils as vector_utils
def evaluate_bilinear_form(Q,left,right):
"""evaluates bilinear forms at several points
Args:
Q (Array ...,ldim,rdim): bilinear form to evaluate
left (Array ...,ldim): points where the bilinear form is evaluated to the left
right (Array ...,rdim): points where the bilinear form is evaluated to the right
Returns:
Array ... bilinear forms evaluated
"""
result = np.einsum('...ij,...i,...j->...',Q,left,right)
return result
def evaluate_quadratic_form(Q,points):
"""evaluates quadratic forms at several points
Args:
Q (Array ...,dim,dim): quadratic form to evaluate
points (Array ...,dim): points where the quadratic form is evaluated
Returns:
Array ... quadratic forms evaluated
"""
result = evaluate_bilinear_form(Q,points,points)
return result
def merge_quadratic_to_homogeneous(Q,b,c):
"""merges quadratic form, linear term, and constant term into a homogeneous matrix
Args:
Q (Array ..., dim, dim): quadratic form matrix
b (Array ..., dim): linear term vector
c (Array ...): constant term
Returns:
Array ..., dim+1, dim+1: homogeneous matrix representing the quadratic form
"""
dim_points = Q.shape[-1]
stack_shape = np.broadcast_shapes(np.shape(Q)[:-2],np.shape(b)[:-1],np.shape(c))
Q_b = np.broadcast_to(Q,stack_shape+(dim_points,dim_points))
b_b = np.broadcast_to(np.expand_dims(b,-1),stack_shape+(dim_points,1))
c_b = np.broadcast_to(np.expand_dims(c,(-1,-2)),stack_shape+(1,1))
H = np.block([[Q_b,0.5*b_b],[0.5*np.swapaxes(b_b,-1,-2),c_b]])
return H
def quadratic_to_dot_product(points):
"""computes the matrix W such that
x.T@Ax = W(x).T*A[ui,uj]
Args:
points ( Array ...,ndim): points of dimension ndim
Returns:
Array ...,ni: dot product matrix (W)
Array ni: i indices of central matrix
Array ni: j indices of central matrix
"""
dim_points = points.shape[-1]
ui,uj = np.triu_indices(dim_points)
W = points[...,ui]*points[...,uj]
return W,ui,uj
def fit_quadratic_form(points):
"""Fits a quadratic form to the given zeroes.
Args:
points (Array ..., n, dim): Input points.
Returns:
Array ..., dim, dim: Fitted quadratic form matrix.
"""
dim_points = points.shape[-1]
normed_points = vector_utils.norm_vector(points)[1]
W,ui,uj = quadratic_to_dot_product(normed_points)
H = np.einsum('...ki,...kj->...ij',W,W)
V0 = kernels.matrix_kernel(H)
Q = np.zeros(V0.shape[:-1]+(dim_points,dim_points))
Q[...,ui,uj]=V0
return Q
# import matplotlib.pyplot as plt
# Q = np.random.randn(3,3)
# x0, y0 = np.linspace(-1,1,300),np.linspace(-1,1,300)
# x,y = np.meshgrid(x0,y0,indexing='ij')
# points = vector_utils.to_homogeneous(np.stack([x,y],axis=-1))
# f = evaluate_quadratic_form(Q,points)
# mask = np.abs(f)<0.01
# u,v = np.where(mask)
# zeros = vector_utils.to_homogeneous(np.stack([x0[u],y0[v]],axis=-1))+np.random.randn(5,u.shape[0],3)*0.1
# Qc = fit_quadratic_form(zeros)
# fchap = evaluate_quadratic_form(Qc,points[...,None,:])
# print()

23
utils/sphere_deprojection.py
View File

@ -1,23 +0,0 @@
import numpy as np
import utils.vector_utils as vector_utils
def deproject_ellipse_to_sphere(M, radius):
"""finds the deprojection of an ellipse to a sphere
Args:
M (Array 3,3): Ellipse quadratic form
radius (float): radius of the researched sphere
Returns:
Array 3: solution of sphere centre location
"""
H = 0.5*(np.swapaxes(M,-1,-2)+M)
eigval, eigvec = np.linalg.eigh(H)
i_unique = np.argmax(np.abs(np.median(eigval,axis=-1,keepdims=True)-eigval),axis=-1)
unique_eigval = np.take_along_axis(eigval,i_unique[...,None],-1)[...,0]
unique_eigvec = np.take_along_axis(eigvec,i_unique[...,None,None],-1)[...,0]
double_eigval = 0.5*(np.sum(eigval,axis=-1)-unique_eigval)
z_sign = np.sign(unique_eigvec[...,-1])
dist = np.sqrt(1-double_eigval/unique_eigval)
C = np.real(radius*(dist*z_sign)[...,None]*vector_utils.norm_vector(unique_eigvec)[1])
return C

69
utils/vector_utils.py
View File

@ -1,69 +0,0 @@
import numpy as np
def norm_vector(v):
"""computes the norm and direction of vectors
Args:
v (Array ..., dim): vectors to compute the norm and direction for
Returns:
Array ...: norms of the vectors
Array ..., dim: unit direction vectors
"""
norm = np.linalg.norm(v,axis=-1)
direction = v/norm[...,np.newaxis]
return norm,direction
def dot_product(v1,v2):
"""Computes the dot product between two arrays of vectors.
Args:
v1 (Array ..., ndim): First array of vectors.
v2 (Array ..., ndim): Second array of vectors.
Returns:
Array ...: Dot product between v1 and v2.
"""
result = np.einsum('...i,...i->...',v1,v2)
return result
def cross_to_skew_matrix(v):
"""converts a vector cross product to a skew-symmetric matrix multiplication
Args:
v (Array ..., 3): vectors to convert
Returns:
Array ..., 3, 3: matrices corresponding to the input vectors
"""
indices = np.asarray([[-1,2,1],[2,-1,0],[1,0,-1]])
signs = np.asarray([[0,-1,1],[1,0,-1],[-1,1,0]])
skew_matrix = v[...,indices]*signs
return skew_matrix
def to_homogeneous(v):
"""converts vectors to homogeneous coordinates
Args:
v (Array ..., dim): input vectors
Returns:
Array ..., dim+1: homogeneous coordinates of the input vectors
"""
append_term = np.ones(np.shape(v)[:-1]+(1,))
homogeneous = np.append(v,append_term,axis=-1)
return homogeneous
def one_hot(i,imax):
"""Converts indices to one-hot encoded vectors.
Args:
i (Array ...): Array of indices.
imax (int): Number of classes.
Returns:
Array ..., imax: One-hot encoded vectors.
"""
result = np.arange(imax)==np.expand_dims(i,axis=-1)
return result