A statistical shape model is a concise depiction of a cluster of shapes. This representation is based on Principal Component Analysis (PCA) performed on this cluster. Principal Component Analysis (PCA) aims to reduce the size of the model by selectively retaining eigenvectors that capture most of the original variability.
The statistical model generates shapes that are geometrically consistent with the shapes assimilated by the model, achieved by varying a single parameter influencing the expression of the eigenvalues encoding the variability.
A 3D shape is represented as a vector: $$x = (x_1, y_1, z_1,…,x_n, y_n, z_n)$$
To obtain the statistical model, we must first align and normalise the shapes in the dataset using the Procrustes analysis algorithm. This registration algorithm estimates the similarity transformation (excluding reflection and with isotropic scaling). Its purpose is to align the shapes two by two by optimizing the distance between them. All the shapes are generally aligned with a randomly chosen reference shape.
The shapes aligned by means of Procustes analysis are assembled in the form of a matrix $k \times n$, where $k$ is the number of shapes in the model:
$$X = \begin{pmatrix} x1_1 & y1_1 & z1_1 & … & x1_n & y1_n & z1_n \ & & & … & & &\ xk_1 & yk_1 & zk_1 & … & xk_n & yk_n & zk_n \end{pmatrix}$$
From $X$ it is possible to find the mean ($x_{mean}$).
import scipy.spatial as spatial
import numpy as np
import vtk
def loadVTKMesh(filename):
reader = vtk.vtkPolyDataReader()
reader.SetFileName(filename)
reader.Update()
return reader.GetOutput()
def arrayToVTKMesh(a, mesh):
a = a.reshape(-1, 3)
return [mesh.GetPoints().SetPoint(idx, a[idx,:]) for idx in range(mesh.GetNumberOfPoints())]
def procrusteAnalysis(meshes):
alignedMeshes = list()
for idx in range(1,len(meshes)):
mtx1, mtx2, disparity = spatial.procrustes(meshes[0], meshes[idx])
alignedMeshes.append((mtx1, mtx2))
return alignedMeshes
if __name__ == '__main__':
meshes = []
for filename in glob.glob("folder-with-meshes/*.vtk"):
meshes.append(vtkMeshPointsToArray(loadVTKMesh(filename)))
alignedMeshes = procrusteAnalysis(meshes)
X = np.stack([alignedMeshes[0][0].flatten()] + [alignedMeshes[idx][1].flatten() for idx in range(len(alignedMeshes))], axis=-1)
x_mean = np.mean(X, axis=-1)
Principal component analysis is used to extract the primary axes from the $X$ matrix. The eigen-decomposition of the covariance matrix $X$ allows only the most important eigenvectors to be retained. The decomposition is based on a criterion that guarantees that almost all the variability of the shapes chosen to construct the shape base can be expressed. From this decomposition, we derive the $P$ matrix which encodes the modes of variation of the model (see here, here and here for more details).
$$P = \begin{pmatrix} \sqrt{eigval[0]} * eigvec[0] \ … \ \sqrt{eigval[i_{max}]} * eigvec[i_{max}] \end{pmatrix}^T$$
Where the eigenvalues and eigenvectors have been sorted in decreasing order and truncated at $i_{max}$.
# inverse cov matrix
covMatrix = np.cov(X)
eigVals, eigVecs = np.linalg.eig(covMatrix)
# sort eigenvectors
idx = eigVals.argsort()[::-1]
eigVals = eigVals[idx]
eigVecs = eigVecs[:, idx]
# retain only significant eigenvectors
variance_explained = np.cumsum(eigVals/np.sum(eigVals))
for index,item in enumerate(variance_explained > 0.99):
if item:
npcs = index
break
P = np.real(np.array([np.sqrt(eigVals[i]) * eigVecs[:, i] for i in range(0, npcs)]).squeeze().T)
By varying the values of a vector $b$ of size $i_{max}$ we can then obtain a new form:
$$y = x_{mean} + Pb$$
b = np.zeros(P.shape[1])
y = (x_mean + P.dot(b)).reshape(-1,3)