View source for 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 =

see section 2.2 in <ref name=FedorovIMR>  A. Fedorov and N. Chrisochoides, “Tetrahedral Mesh Generation for Non-rigid Registration of Brain MRI: Analysis of the Requirements and Evaluation of Solutions,” in Proceedings of the 17th International Meshing Roundtable, Berlin, Heidelberg, 2008, pp. 55–72, doi: 10.1007/978-3-540-87921-3_4.   </ref> about why equidistribution might help

= Mesh Requirements =
In<ref name=FedorovIMR /> 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 <ref name=FedorovIMR /> two attempts were made. 
# Creating new splitting rules in a BCC-based mesh generator and 
# Creating a sizing function for TetGen and NETGEN. 
The sizing function of the second approach is based on a distance criterion. In particular:
<blockquote>
'' 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. <ref name=FedorovIMR />'' 
Where vertex cell complex is defined as the union of all the elements attached to a vertex.
</blockquote>

== 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 efficienty 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 arround each vertex that meet the density requirements of the problem. For linear tetrahedral meshes which is the interset in this work boundaing 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.
[[File:APBNRR-adaptivity-method-comparison.png|500px|thumb|none|Figure 1 Comparison of two methods]]

The smallest bounding ellipsoid is a geometrical construction that used a lot in optimization and specifically in interger programing [CITE]. As such there are existing algorithms like the Khachiyan algorithm which we use here in order to consutrct such ellipsoids. <ref> Khachiyan, L. G. 1979. "A Polynomial Algorithm in Linear Programming". Doklady Akademii Nauk SSSR 244, 1093-1096 (translated in Soviet Mathematics Doklady 20, 191-194, 1979) </ref> The choise 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 <ref> P.J. Frey and F. Alauzet. Anisotropic mesh adaptation for CFD computations. Computer methods in applied mechanics and engineering, 194(48):5068{5082, 2005 </ref>. 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 uise the metric derived from the ellipsoids directly as a way to control the density of the mesh and ultimately the error between the registred and the intra-operative image. Section XX presents the results.


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

= Results =
== DDL case ==
k : number of landmarks in each vertex-cell

a : inflation parameter

iso : use of isotropic metric instead (diagonal matrix)
{| class="wikitable" style="text-align:right;"
|- style="text-align:center;"
! rowspan="2" style="text-align:left;" | HD
! colspan="3" style="vertical-align:bottom;" | Landmark error
! rowspan="2" style="text-align:left;" | # vertices
! rowspan="2" style="text-align:left;" | #cells
|- style="text-align:left;"
| style="vertical-align:bottom;" | min
| style="vertical-align:bottom;" | max
| style="vertical-align:bottom;" | avg
|- style="vertical-align:bottom;"
| 3.56
| 0.72
| 3.57
| 2.50
| 2,504
| 9,650
|- style="vertical-align:bottom;"
| style="background-color:#B7E1CD;" | 3.24
| style="background-color:#B7E1CD;" | 0.68
| 3.87
| style="background-color:#B7E1CD;" | 2.46
| 6,412
| 32,990
|- style="vertical-align:bottom;"
| style="background-color:#B7E1CD;" | 3.24
| 1.26
| style="background-color:#B7E1CD;" | 3.36
| 2.57
| 5,025
| 25,565
|- style="vertical-align:bottom;"
| 4.00
| 1.18
| 4.06
| 2.61
| 5,864
| 29,520
|- style="vertical-align:bottom;"
| 3.58
| 0.90
| style="background-color:#B7E1CD;" | 3.37
| 2.55
| 4,347
| 21,904
|- style="vertical-align:bottom;"
| 4.03
| style="background-color:#B7E1CD;" | 0.43
| style="background-color:#B7E1CD;" | 3.42
| style="background-color:#B7E1CD;" | 2.42
| 3,484
| 17,322
|- style="vertical-align:bottom;"
| 3.56
| style="background-color:#B7E1CD;" | 0.18
| style="background-color:#B7E1CD;" | 3.16
| style="background-color:#B7E1CD;" | 2.31
| 4,224
| 20,973
|}

== CHF == 
{| class="wikitable" style="text-align:right;"
|- style="text-align:left;"
! rowspan="2" | method
! rowspan="2" | HD
! colspan="3" style="text-align:center; vertical-align:bottom;" | Landmark error
! rowspan="2" | # vertices
! rowspan="2" | #cells
|- style="text-align:left;"
| style="vertical-align:bottom;" | min
| style="vertical-align:bottom;" | max
| style="vertical-align:bottom;" | avg
|- style="vertical-align:bottom;"
| style="text-align:left;" | basecase
| 4.06
| 2.06
| 5.37
| 3.65
| 2,833
| 11,040
|- style="vertical-align:bottom;"
| style="text-align:left;" | k100
| style="background-color:#B7E1CD;" | 3.96
| 2.15
| 6.77
| 3.74
| 13,114
| 79,822
|- style="vertical-align:bottom;"
| style="text-align:left;" | k100iso
| style="background-color:#B7E1CD;" | 4.05
| 2.46
| 5.98
| 3.74
| 13,114
| 67,198
|- style="vertical-align:bottom;"
| style="text-align:left;" | k200
| style="background-color:#B7E1CD;" | 3.94
| style="background-color:#B7E1CD;" | 1.97
| 5.64
| style="background-color:#B7E1CD;" | 3.59
| 7,960
| 41,131
|- style="vertical-align:bottom;"
| style="text-align:left;" | k200a1.2
| 5.28
| 2.15
| 5.91
| 3.90
| 6,471
| 32,809
|- style="vertical-align:bottom;"
| style="text-align:left;" | k200a1.5
| 4.78
| style="background-color:#B7E1CD;" | 1.80
| 5.89
| style="background-color:#B7E1CD;" | 3.55
| 5,095
| 25,465
|- style="vertical-align:bottom;"
| style="text-align:left;" | k200a2
| 4.39
| style="background-color:#B7E1CD;" | 2.04
| 5.46
| style="background-color:#B7E1CD;" | 3.60
| 3,876
| 18,017
|- style="vertical-align:bottom;"
| style="text-align:left;" | k200a2iso
| 4.31
| style="background-color:#B7E1CD;" | 1.89
| style="background-color:#B7E1CD;" | 5.19
| style="background-color:#B7E1CD;" | 3.48
| 3,309
| 15,346
|- style="vertical-align:bottom;"
| style="text-align:left;" | k200a1.2iso
| 4.55
| style="background-color:#B7E1CD;" | 2.01
| 5.48
| 3.91
| 5,831
| 29,874
|- style="vertical-align:bottom;"
| style="text-align:left;" | k200a1.5iso
| style="background-color:#B7E1CD;" | 3.85
| style="background-color:#B7E1CD;" | 1.74
| 5.48
| style="background-color:#B7E1CD;" | 3.46
| 3,343
| 22,054
|- style="vertical-align:bottom;"
| style="text-align:left;" | k200iso
| style="background-color:#B7E1CD;" | 3.94
| 2.15
| 6.77
| 3.74
| 7,199
| 37,703
|- style="vertical-align:bottom;"
| style="text-align:left;" | k300
| style="background-color:#B7E1CD;" | 3.59
| style="background-color:#B7E1CD;" | 1.67
| 5.81
| 3.86
| 5,699
| 28,667
|- style="vertical-align:bottom;"
| style="text-align:left;" | k300a1.2
| 4.53
| style="background-color:#B7E1CD;" | 2.01
| 5.82
| 3.65
| 4,797
| 23,836
|- style="vertical-align:bottom;"
| style="text-align:left;" | k300a1.5
| style="background-color:#B7E1CD;" | 3.24
| 2.26
| style="background-color:#B7E1CD;" | 5.32
| style="background-color:#B7E1CD;" | 3.59
| 3,784
| 18,003
|- style="vertical-align:bottom;"
| style="text-align:left;" | k300iso
| style="background-color:#B7E1CD;" | 3.28
| style="background-color:#B7E1CD;" | 1.66
| 5.54
| 3.82
| 6,012
| 26,022
|- style="vertical-align:bottom;"
| style="text-align:left;" | k300a1.2iso
| style="background-color:#B7E1CD;" | 4.05
| 2.61
| 5.50
| 3.94
| 4,269
| 20,624
|- style="vertical-align:bottom;"
| style="text-align:left;" | k300a1.5iso
| 4.51
| style="background-color:#B7E1CD;" | 1.76
| 5.39
| 3.74
| 3,243
| 15,277
|}

== case12 ==
{| class="wikitable" 
|-
! style="vertical-align:bottom;" | 
! rowspan="2" | HD
! colspan="3" style="vertical-align:bottom;" | Landmark error
! rowspan="2" style="vertical-align:bottom;" | # vertices<br />
! rowspan="2" style="vertical-align:bottom;" | #cells<br />
|-
| style="vertical-align:bottom;" | method
| style="vertical-align:bottom;" | min
| style="vertical-align:bottom;" | max
| style="vertical-align:bottom;" | avg
|- style="vertical-align:bottom;"
| basecase
| 2.24
| 1.07
| 5.9
| 3.51
| 3,264
| 13,210
|- style="vertical-align:bottom;"
| k100
| 2.44141
| 1.74
| 7.79
| 3.79
| 18,856
| 97,772
|- style="vertical-align:bottom;"
| k100a1.5
| style="background-color:#B7E1CD;" | 2.013
| 1.18
| 6.27
| 3.58
| 11,110
| 56,357
|- style="vertical-align:bottom;"
| k100iso
| style="background-color:#B7E1CD;" | 2
| 1.16
| style="background-color:#B7E1CD;" | 5.74
| style="background-color:#B7E1CD;" | 3.24
| 15,415
| 76,803
|- style="vertical-align:bottom;"
| k200
| style="background-color:#B7E1CD;" | 2
| style="background-color:#B7E1CD;" | 0.23
| 7.04
| style="background-color:#B7E1CD;" | 3.28
| 9,760
| 50,085
|- style="vertical-align:bottom;"
| k200a1.2
| style="background-color:#B7E1CD;" | 2
| style="background-color:#B7E1CD;" | 1.03
| 7.2
| 3.53
| 7,816
| 39,519
|- style="vertical-align:bottom;"
| k200a1.5
| style="background-color:#B7E1CD;" | 2.194
| 1.73
| 7.83
| 4.12
| 5,995
| 29,831
|- style="vertical-align:bottom;"
| k200iso
| style="background-color:#B7E1CD;" | 2
| 1.1
| 6.88
| style="background-color:#B7E1CD;" | 3.44
| 8,763
| 43,449
|- style="vertical-align:bottom;"
| k200a1.2iso
| style="background-color:#B7E1CD;" | 2.013
| 1.32
| 6.15
| 3.62
| 6,832
| 33,353
|- style="vertical-align:bottom;"
| k300
| style="background-color:#B7E1CD;" | 2.024
| 1.17
| 6.76
| 3.72
| 6,854
| 34,812
|- style="vertical-align:bottom;"
| k300a1.5
| style="background-color:#B7E1CD;" | 2.01
| 1.86
| 6.76
| 3.95
| 4,387
| 21,170
|- style="vertical-align:bottom;"
| k300iso
| style="background-color:#B7E1CD;" | 2.024
| style="background-color:#B7E1CD;" | 0.95
| 6.38
| style="background-color:#B7E1CD;" | 3.39
| 6,245
| 30,601
|- style="vertical-align:bottom;"
| k400
| style="background-color:#B7E1CD;" | 2.02
| 2.12
| 7.59
| 3.91
| 5,364
| 26,874
|- style="vertical-align:bottom;"
| k30
| 3.19
| 1.77
| 5.33
| 3.77
| 58,382
| 305,602
|}

= 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:<ref> F. Alauzet, ‘Size gradation control of anisotropic meshes’, Finite Elements in Analysis and Design, vol. 46, no. 1–2, pp. 181–202, Jan. 2010, doi: 10.1016/j.finel.2009.06.028. </ref> to the metric were not studied in order to keep this study short.

==References==
{{Reflist}}