D.2 Method 151 training ear canals and Φ is found by a principal component analysis of the Procrustes aligned shapes. It was found that t = 7 modes of variation was enough to sufficiently explain the variation seen in the training data [201]. A new shape exhibiting the variance seen in the training set is synthesised by adding a combination of eigenvectors to the average shape xnew = x+ Φb, wherebis a vector of weights controlling the modes of shape variation. In this article a size-and-shape model of the ear canal is used meaning that the size variation of ear canal is included in Φ [83].
D.2.3 Alignment of Faceplates
The rigid body transformation,Ti, that in a least squares sense aligns ear canal xi to the average ear canal is also applied to the faceplate plane of that shape, thus transforming the faceplate plane to the average ear canal. All the faceplate planes from the training set transformed to the average canal and the average ear canal can be seen in Fig. D.3b.
A plane can be defined as the set of points satisfying:
P ={x∈R3|F(x) =λAx+λBy+λCz+λD= 0, λ6= 0}, (D.1) where the normal vector of the plane is n = (A, B, C) and x = (x, y, z) is a point. This means that a plane can be represented as a parameter vector, p= (A, B, C, D), which is defined up to a normalising factorλ6= 0, and that can be normalised to have unit length, kpk = 1. After normalising, all planes in the training set are now lying on the 3 Dimensional hypersphere.
Using this representation, planes can be looked upon as points of P3R, the real projective space of dimension 3 [103]. A metric can now be established in this differential manifold and thus it is possible specify a distance between two infinite planes,pi andpj; namely the arc-length on the hypersphere,da = min(arccos(pi,pj),arccos(−pi,pj)). Note that it is necessary to handle that a plane p is equal to the plane −p. An approximate distance measure is de = min(kpi−pjk,kpi+pjk). An average planepcan be determined from:
p= arg min
p
X
i
d(p,pi), kpk= 1, (D.2) and is found by gradient descent optimisation withde as the distance function.
Planevectors not pointing into the same hemisphere as the current estimate ofp are negated prior to each iteration of the optimisation algorithm. The average ear canal and the estimated average faceplate plane can be seen in D.3b. For visualisation purposes the planes are represented as finite planes limited to the
proximity of the ear canal. From an anatomical viewpoint the placement of the average plane in the average ear canal looks sensible.
D.2.4 Shape Fitting and Recognition
The statistical shape model built previously can be fitted to new ear canals not included in the training set following the approach used in the active shape models (ASM) [61]. The following method is similar to the method described in [143]. Given a sampled surface of an ear canal not included in the training set, xtarget, the goal is to fit a model mesh, xmodel toxtarget. The fitting procedure consists of a combination of the iterative closest point transformation (ICP) [25, 262] and a shape deformation using, xdeform = ¯x+ Φbf. In each stage the model mesh is therefore represented by a global transform, Tf, and a set of deformation parameters,bf: xmodel=Tf(¯x+ Φbf). The size-and-shape model is used and therefore Tf is a rigid-body transform [83]. Initially the average ear canal (bf =0) is placed in the scene and then iteratively transformed and deformed to fit the target ear canal. In each iterationTf is found by ICP. The ICP algorithm initially aligns the centre of mass of the model mesh to the centre of mass of the target mesh and then the closest point on the target surface is found for all points in the model mesh. The surface of the target ear canal is typically acquired using a laser scanner and therefore consists of a large number of points and triangles. The closest point on the target surface is found by the vtkCellLocator class found in the Visualization Toolkit [219]. It uses a uniform-level octree subdivision [188], where each octant carries an indication of whether it is empty or not, and each leaf octant carries a list of the polygons inside of it. Therefore the closest point is found as a point lying on a triangle and not just as a vertex. Tf is found as the transformation that minimises the distances between the model mesh and the closest points on the target mesh in a least squares sense [133].
After the model mesh has been rigidly aligned to the target mesh the clos-est point on the target shape is found for all points in the model mesh. The collection of closest points yields a new shape, x0. This new shape is trans-formed back to alignment with the shape model using the inverse transforma-tion: x00 = T−1f (x0). The model shape is now deformed to look as much like x00 as possible, but still being in the allowed shape space. This is done by finding the parameter vector, by back-projecting x00 on the eigenvector space:
bf = ΦT(x00−x). The parameters found are then clamped to three standard¯ deviations and the new model mesh is synthesised usingxdeform= ¯x+ Φbf and used in the next iteration.
The quality of the fit is defined to be the square root of the summed squared
D.2 Method 153 distances from all the points in the model mesh to the closest point on the target mesh
F= ( XN
i
d2(pmi , ptclosest))1/2, (D.3)
where d2(pmi , pti) is the squared Euclidean distance from point i in the model shape to the closest point on the surface of the target shape. The iteration stops when the change inFis sufficiently small. The parameter vector describing the approximated target shape, bf, can be used to test if the found shape belongs to the same class of shapes as the training set [61].
The algorithm described previously does not handle local minima very well and tends to get stuck if the initial orientation of the model is very different from the orientation of the target shape. To overcome this problem a guess of the starting placement is made by starting ICP with the model pre-transformed in 24 different configurations and then selecting the ICP with the best initial fit, F. The 24 configurations are made by an initial rotation around the x-axis followed by a rotation around either the y- or the z-axis. The allowed rotations are 90,180 and 270 degrees. So it is equal to the object first being rotated around itself and then rotated to lie on one of the 6 axis of symmetry of a cube. This is a primitive and brute-force approach, but it works very well. Examples of ICP fits with different pre-rotations can be seen in Fig. D.4. It is seen that the pre-rotation used in the top row is ideal.
D.2.5 Prediction of Faceplate Placement
The fitting procedure described above can also be used to place a faceplate plane on a new ear canal. The average faceplate plane, p, associated with the average ear canalxcan be fitted to the new ear canal using the transform used to fit the model shape to the new ear canal. The fitting consists of a rigid body transformTf and a deformation expressed withbf. It would be ideal to use both the rigid transformation and the shape deformation in the prediction of the faceplate placement. In this work the prediction is limited to the rigid transformation meaning that the predicted plane is calculated by applying Tf
top. A faceplate plane placed by the algorithm in a new ear canal can be seen together with the fitted model mesh in Fig. D.5.
The fitting and faceplate plane prediction algorithm is outlined in algorithm 2.
Figure D.4: ICP fitting. Here the green model ear is fitted to the blue target ear canal. The initial placement of the model ear is seen to the left and after ICP to the right. The fit quality after the converged ICP shown in the top row isF= 0.722 andF = 2.50 for the one in the bottom row.
D.2 Method 155
Algorithm 2Fitting the Shape Model to a Target Surface,xtarget, and Predict Faceplate Placement
bf =0,Tf =I→xmodel=x
Determine pre-transformationTf by multiple ICP while Not convergencedo
xmodel =Tf(¯x+ Φbf)
DetermineTi using ICP soxmodel fitsxtarget. UpdateTf withTi and updatexmodel
For each point inxmodel find the closest point on the surface ofxtarget
Generate a new shapex0 of these closest points ApplyT−1f tox0 generating x00
Setbf = ΦT(x00−x)¯ Convergence if ∆F is small end while
Final fitted model: xfit =Tf(¯x+ Φbf)
ApplyTf topto get the predicted faceplate plane forxtarget
Figure D.5: The shape model fitted to a new ear canal and the predicted face-plate plane. The target shape is blue and the fitted model is green.