Internal:Adaptive Meshing For APBNRR

Problem Description

Reducing the error for Non-Rigid Registration (NRR) applications depends among other factors on the Equi-distribution of the registration points within the mesh. The goal of this project is to explore metric adaptive methods from CFD and apply them if possible in the area of NRR.

Background

The NRR model is given by

As per section 2.2 of

" the matrix H^TSH has a non-zero 3 × 3 entry for each mesh vertex and edge with the cell complexes containing registration points. The corresponding submatrices can be expressed as the sum over the registration points in a cell. "

Also note that " We define the cell complex as a set of mesh tetrahedra incident on a mesh k-cell, e.g., mesh vertex is a 0-cell, and mesh edge 1-cell"

Mesh Requirements

In a set of application specific mesh requirements for NRR were presented.

• (R1) Equi-distribution of the registration points w.r.t. mesh vertex cells  :

  small number of registration points (e.g., less than 10, but greater than 0) within the vertex cell complex makes the formulation more sensitive to outliers and introduces additional displacement error

• (R2) Minimization of the approximation error at registration points :

  error of the displacements recovered at registration points can be reduced locally by using smaller mesh elements.

• (R3) Prevention of tetrahedron inversion during mesh deformation:

  while the interpolation error shown does not depend on the tetrahedron shape, inversion or collapse of a tetrahedron will result in an unrealistic deformation field, e.g., points inside the different tetrahedra can map to the same image location. We can attempt to remedy this problem by adjusting element size locally according to the expected deformation.

Alongside the application-specific requirements we would also like to adhere to the following generic requirements which are related with FEM calculations and image-based mesh generatio in general

• (R4) The mesh boundary (triangulation) should be close to the surface of the image.
• (R5) Mesh size should be minimized to reduce the computational costs.
• (R6) Mesh elements should not have small angles

some of these rules work against each other. For example we require small elements (R2) to reduce the error but we also need a small size (R5) of a mesh.

Methods

Method 1 :Isotropic Sizing Functions based on knn

Already in two attempts were made.

1. Creating new splitting rules in a BCC-based mesh generator and
2. Creating a sizing function for TetGen and NETGEN.

The sizing function of the second approach is based on a distance criterion. In particular:

The sizing value at each mesh vertex was initialized with the distance to the kth registration point closest to it, to reflect the density of the registration point distribution. The idea here is that the shape of a perfect mesh vertex cell complex is close to a ball with the radius prescribed by the background mesh. Where vertex cell complex is defined as the union of all the elements attached to a vertex.

Method 2 : Metric-based Anisotropic adaptation based on knn

The method above creates adaptive meshes but is does not capture the density of the landmarks efficiently due to the reason than only the k-th point is used and the relative locations of the rest of the k-1 landmarks is ignored. Going back to the paper

The idea here is that the shape of a perfect mesh vertex cell complex is close to a ball with the radius prescribed by the background mesh.

Building upon this observation one could create bounding volumes around each vertex that meet the density requirements of the problem. For linear tetrahedral meshes which is the focus in this work bounding volumes can be replaced with bounding ellipsoids and specifically the smallest bounding ellipsoid that is centered at a vertex and encloses the k-closest registration points. See for example Figure 1.

Figure 1 Comparison of the two methods: Left isotropic metric based on distance of k-th neighbour. Right Smallest volume enclosing ellipse centred at the mesh point and containing the k registration points

The smallest bounding ellipsoid is a geometrical construction that used a lot in optimisation and specifically in integer programming [CITE]. As such there are existing algorithms like the Khachiyan algorithm which we use here in order to construct such ellipsoids. The choice of ellipsoid has another significant advantage: It has a natural mapping to the set of positive-definite matrices. This means that for every ellipsoid we can construct a 3x3 positive-definite matrix that corresponds to it.

Another field where positive-definite matrices are used is in metric-adaptive mesh adaptation where an a posteriori error estimate of the interpolation error guides the mesh adaptation . One of the popular methods is to construct a metric tensor expressed as positive definite matrices based on the Hessian reconstructed from a scalar variable.

In this work however we have a purely geometrical approach where we use the metric derived from the ellipsoids directly as a way to control the density of the mesh and ultimately the error between the registered and the intra-operative image. Section XX presents the results.

Experimental Setup

Augmenting mesh adaptation to PAPBNRR:

1. Generate mesh (A) from registered segmented image using an image-to-mesh tool, PODM was used in this work.
2. Using the identified landmarks from the block-matching step of PABPNRR find the k-nearrest landmarks at each point p of A (VTK).
3. Evaluate an approximate minimum volume bounding ellispoid E using the Khachiyan algorithm (implemented in Eigen [CITE].)
4. Construct a metric tensor based on ellipsoid E and assign it to point p.
5. adapt the mesh using MMG3D [1]
6. return the mesh to PAPBNRR

Results

k : number of landmarks in each vertex-cell

a : inflation parameter

iso : use of isotropic metric instead (diagonal matrix)

case 18 - (original casename CHF)

method	HD	Landmark error			# vertices	#cells
		min	max	avg
basecase	4.06	2.06	5.37	3.65	2,833	11,040
k100	3.96	2.15	6.77	3.74	13,114	79,822
k100iso	4.05	2.46	5.98	3.74	13,114	67,198
k200	3.94	1.97	5.64	3.59	7,960	41,131
k200a1.2	5.28	2.15	5.91	3.90	6,471	32,809
k200a1.5	4.78	1.80	5.89	3.55	5,095	25,465
k200a2	4.39	2.04	5.46	3.60	3,876	18,017
k200a2iso	4.31	1.89	5.19	3.48	3,309	15,346
k200a1.2iso	4.55	2.01	5.48	3.91	5,831	29,874
k200a1.5iso	3.85	1.74	5.48	3.46	3,343	22,054
k200iso	3.94	2.15	6.77	3.74	7,199	37,703
k300	3.59	1.67	5.81	3.86	5,699	28,667
k300a1.2	4.53	2.01	5.82	3.65	4,797	23,836
k300a1.5	3.24	2.26	5.32	3.59	3,784	18,003
k300iso	3.28	1.66	5.54	3.82	6,012	26,022
k300a1.2iso	4.05	2.61	5.50	3.94	4,269	20,624
k300a1.5iso	4.51	1.76	5.39	3.74	3,243	15,277
k400	3.71	2.12	5.50	3.96	4,460	22,342
k400a1.2	4.05	2.06	5.05	3.61	3,766	18,077
k400a1.5	4.05	1.92	5.17	3.65	2,983	13,812
k400a2	4.42	1.90	5.59	3.69	2,211	9,759
k500	3.85	1.90	5.21	3.82	3,762	18,508
k600	4.39	2.09	5.25	3.58	3,321	16,060
k400iso	3.42	2.29	5.76	3.92	4,008	19,466
k500iso	3.41	2.05	5.11	3.85	3,546	16,870
k600iso	4.29	2.17	4.97	3.58	3,334	15,529

case 26 - (original casename DDL)

method	HD	Landmark error			# vertices	#cells
		min	max	avg
basecase	3.56	0.72	3.57	2.50	2,504	9,650
200k	3.24	0.68	3.87	2.46	6,412	32,990
200ka1.2	3.24	1.26	3.36	2.57	5,025	25,565
200kiso	4.00	1.18	4.06	2.61	5,864	29,520
300k	3.58	0.90	3.37	2.55	4,347	21,904
300ka1.2	4.03	0.43	3.42	2.42	3,484	17,322
300kiso	3.56	0.18	3.16	2.31	4,224	20,973
400k	3.24	1.01	3.30	2.41	3,421	17,050
500k	3.58	0.79	3.7	2.55	2,785	13,703
500kiso	3.71	0.65	3.46	2.45	2,716	13,185

case12

	HD	Landmark error			# vertices	#cells
method		min	max	avg
basecase	2.24	1.07	5.9	3.51	3,264	13,210
k100	2.44141	1.74	7.79	3.79	18,856	97,772
k100a1.5	2.013	1.18	6.27	3.58	11,110	56,357
k100iso	2	1.16	5.74	3.24	15,415	76,803
k200	2	0.23	7.04	3.28	9,760	50,085
k200a1.2	2	1.03	7.2	3.53	7,816	39,519
k200a1.5	2.194	1.73	7.83	4.12	5,995	29,831
k200iso	2	1.1	6.88	3.44	8,763	43,449
k200a1.2iso	2.013	1.32	6.15	3.62	6,832	33,353
k300	2.024	1.17	6.76	3.72	6,854	34,812
k300a1.5	2.01	1.86	6.76	3.95	4,387	21,170
k300iso	2.024	0.95	6.38	3.39	6,245	30,601
k30	3.19	1.77	5.33	3.77	58,382	305,602
k400	2.02	2.12	7.59	3.91	5,364	26,874
k400iso	2.19	1.10	7.28	3.83	4,990	24,125
k400a1.2	2.24	1.52	6.46	3.45	4,456	21,708
k400a1.5	2.01	1.82	6.07	3.67	3,424	16,356
k500iso	1.95	1.22	7.53	3.71	4,177	19,893
k500	2.22	0.55	7.85	3.99	4,520	22,383
k500a1.2	2	1.01	7.10	3.7	3,629	17,593
k500a1.5	2.64	0.93	6.15	3.25	2,838	13,291

Future Work

• As part of the mesh adaptation , the quality of the mesh is improved significantly.

This is possible since the mesh adaptation software is not constrained by the Delaunay

criterion and moreorev is has mesh improvment cababilities not yet implemented in the image-to-mesh tool.

It would be of interested to study the effect of adaptation on the condition number of stiffness matrix

of the solver as well as to the running time of the solver.

• The effects of applying gradation: to the metric were not studied in order to keep this study short.