pp

of relaxation of the tetrahedral mesh of the program of tetrahedral
mesh generation

.gitignore

@@ -0,0 +1,16 @@
+.vs
+.vscode
+out
+build*
+CmakeSettings.json
+*.vtk
+*.aux
+*.bbl
+*.blg
+*.log
+*.lbl
+rapport/*.pdf
+*.toc
+*.snm
+*.out
+*.nav

CMakeLists.txt

@@ -2,7 +2,41 @@ cmake_minimum_required(VERSION 3.15)
 
 project(pfe)
 
-add_subdirectory(external/vtk)
+set(VTK_COMPONENTS
+  VTK::CommonCore
+  VTK::IOCore
+  VTK::FiltersCore
+  VTK::CommonColor
+  VTK::CommonDataModel
+  VTK::IOLegacy
+  VTK::IOGeometry
+  VTK::IOXML
+  VTK::FiltersModeling
+  VTK::FiltersGeometry
+  VTK::vtksys)
+set(ENABLE_VIEWER OFF CACHE BOOL "Enable the 3D viewer, depends on Qt.")
+if(ENABLE_VIEWER)
+  list(APPEND VTK_COMPONENTS
+	VTK::RenderingCore
+	VTK::ViewsCore
+	VTK::GUISupportQt
+	VTK::RenderingQt
+	VTK::ViewsQt
+	VTK::RenderingVolume
+	VTK::RenderingVolumeOpenGL2)
+endif()
+
+set(USE_SYSTEM_VTK NO CACHE BOOL
+  "Use the version of vtk installed in the system instead of downloading and compiling it ourselves.")
+if(USE_SYSTEM_VTK)
+  list(TRANSFORM VTK_COMPONENTS REPLACE "VTK::" ""
+	OUTPUT_VARIABLE VTK_PACKAGE_COMPONENTS)
+  message("VTK_COMPONENTS: ${VTK_COMPONENTS}")
+  message("VTK_PACKAGE_COMPONENTS: ${VTK_PACKAGE_COMPONENTS}")
+  find_package(VTK COMPONENTS ${VTK_PACKAGE_COMPONENTS})
+else()
+  add_subdirectory(external/vtk)
+endif()
 
 add_executable(pfe)
 
@@ -10,23 +44,33 @@ target_compile_features(pfe PRIVATE cxx_std_17)
 
 target_sources(pfe PRIVATE
   src/main.cc
-  src/angles_filter.cc
-  src/angles_filter.h
-  src/aspect_ratio_filter.cc
-  src/aspect_ratio_filter.h)
+  src/analysis/angles_filter.cc
+  src/analysis/angles_filter.h
+  src/analysis/aspect_ratio_filter.cc
+  src/analysis/aspect_ratio_filter.h
+  src/analysis/dihedral_angles_filter.cc
+  src/analysis/dihedral_angles_filter.h
+  src/surface_points_filter.cc
+  src/surface_points_filter.h
+  src/kd_tree.cc
+  src/kd_tree.h
+  src/fitting/mesh_fit_filter.cc
+  src/fitting/mesh_fit_filter.h
+  src/point_tris_dist.cc
+  src/point_tris_dist.h
+  src/fitting/project_surface_points_on_poly.cc
+  src/fitting/project_surface_points_on_poly.h
+  src/fitting/remove_external_cells_filter.cc
+  src/fitting/remove_external_cells_filter.h
+  src/fitting/relaxation_filter.cc
+  src/fitting/relaxation_filter.h
+  src/analysis/max_distance_filter.cc
+  src/analysis/max_distance_filter.h
+  src/closest_polymesh_point.cc
+  src/closest_polymesh_point.h)
 
-target_link_libraries(pfe PRIVATE
-  VTK::CommonCore
-  VTK::ViewsCore
-  VTK::RenderingCore
-  VTK::IOCore
-  VTK::FiltersCore
-  VTK::CommonColor
-  VTK::GUISupportQt
-  VTK::RenderingQt
-  VTK::ViewsQt
-  VTK::RenderingVolume
-  VTK::CommonDataModel
-  VTK::IOLegacy
-  VTK::IOXML
-  VTK::RenderingVolumeOpenGL2)
+target_link_libraries(pfe PRIVATE ${VTK_COMPONENTS})
+
+if(ENABLE_VIEWER)
+  target_compile_definitions(pfe PRIVATE ENABLE_VIEWER)
+endif()

LISEZMOI

@@ -4,9 +4,24 @@ Projet de fin d'étude de M2 Info Géométrie et Informatique Graphique.
 Compilation :
 	cmake -Bbuild -DCMAKE_BUILD_TYPE=Release && cmake --build build
 
-
 Compilation (pour développeurs) :
 	cmake -Bbuild \
 		-DCMAKE_BUILD_TYPE=RelWithDebInfo \
 		-DCMAKE_CXX_FLAGS="-Wall -Wextra" \
 	&& cmake --build build
+
+Pour utiliser une installation système de VTK :
+	-DUSE_SYSTEM_VTK=ON
+
+Pour désactiver la visualisation :
+	-DENABLE_VIEWER=OFF
+
+Exécution :
+	build/pfe fit|analyze tetmesh polydata [radiusScale relaxationIterCount]
+
+	 Produit un fichier nommé out.vtu. En mode `analyze`, ce fichier
+	 contient le même modèle que tetmesh mais avec les informations de
+	 qualité ajoutée, et affiche des statistiques dans la console.
+
+	 En mode `fit`, adapte le tetmesh au polydata, puis effectue
+	 l'analyse comme `analyze`.

cmds.txt

@@ -0,0 +1,2 @@
+cmake -Bbuild -DCMAKE_BUILD_TYPE=RelWithDebInfo -DCMAKE_CXX_FLAGS="-Wall -Wextra" -DUSE_SYSTEM_VTK=ON -DENABLE_VIEWER=OFF
+cmake --build build

data/TetMesh2.vtk

@@ -1,23 +0,0 @@
-# vtk DataFile Version 3.0
-vtk output
-ASCII
-DATASET UNSTRUCTURED_GRID
-POINTS 4 float
-0 0 0
-1 0 0
-0 1 0
-0 0 1
-CELLS 1 5
-4 0 1 2 3
-CELL_TYPES 1
-10
-POINT_DATA 4
-SCALARS Indication float
-LOOKUP_TABLE default
-1
-1
-1
-1
-FIELD law 1
-cell_law 1 1 int
-0

data/tetrahedron.obj

@@ -0,0 +1,10 @@
+o Solid
+v 0.000000 0.769688 0.000000
+v 0.000000 0.000000 0.769688
+v 0.769688 0.000000 0.000000
+v 0.000000 0.000000 0.000000
+s off
+f 1 2 3
+f 1 3 4
+f 1 4 2
+f 2 4 3

rapport/présentation.txt

@@ -0,0 +1,117 @@
++------------+
+| Partie 1/3 |
++------------+
+
+1. Présentation du sujet
+
+[4/37] Maillages volumiques
+ - Structure composée de cellules polyédriques (tétraèdres, hexaèdres…)
+   reliées entre elles comme pour un maillage polygonal
+   classique, mais les cellules sont en 3D.
+
+[5/37] Anisotropie
+Notion de dépendance à l'orientation.
+Image : plasma dans un réacteur à fusion nucléaire. On distingue clairement
+une direction principale qui est la « boucle ». Cet exemple présente
+donc de l'anisotropie.
+
+[6/37] Anisotropie
+Plus formellement, un champ de métriques.
+Une métrique définit la déformation locale.
+
+À gauche : isotropie.
+À droite : espace anisotropique circulaire.
+
+[7/37] Applications
+Principalement dans le domaine de la simulation physique
+Notamment la mécanique des fluides, où ils sont utilisés comme base de
+discrétisation pour porter les calculs.
+
+Par rapport à une simple voxélisation : résolution adaptative, plus
+haute résolution dans les zones « intéressantes », pas de perte de temps
+sur les calculs dans les zones moins intéressantes.
+
+Par rapport à une voxélisation multi-niveaux comme un octree ou un
+maillage isotropique, la forme des cellules peut être adaptée au
+problème pour encore améliorer la précision des calculs. Le problème
+reste maintenant de générer ces maillages anisotropiques.
+
+[8/37]
+Présente des approches pour la génération de maillages surfaciques (et
+pas volumiques) anisotropes.
+D'après leur conclusion : « Bien qu'elles soient toutes prouvablement
+correctes, aucune des méthodes présentées n'est universellement
+pratique ».
+De plus, celle avec les résultats les « moins pires » est coûteuse en
+temps de calcul, ce qui ne ferait qu'empirer après adaptation sur
+maillages volumiques.
+
++------------+
+| Partie 2/3 |
++------------+
+
+[11/37]
+Difficile d'imaginer en trois mois qu'on puisse produire une solution
+qui fonctionne quand une thèse entière qui ne se consacre qu'au
+surfacique y peine.
+
+[17/37]
+Présentation de la chaîne de traitement « pipeline »
+Deux modes : analyse et adaptation
+Image : en mode adaptation. En mode analyse on enlève simplement les
+étapes qui modifient la géométrie du maillage.
+
+[19/37]
+point-in-polygon : permet de détecter si un point est dans un
+polygone. On tire un rayon depuis le point vers l'infini et compte le
+nombre d'intersections, pair : le points est hors du polygone,
+impair : le point est dans le polygone.
+
+[20/37]
+Les cellules dont tous les points sont à l'extérieur du maillage
+surfaciques sont supprimées.
+
+[21/37]
+Les points de surface sont ensuite marqués. Pour cela on détecte
+quelles faces sont mitoyennes, ce qu'on peut déduire car on sait à
+quelle cellules appartiennent chaque point.
+Si les trois points d'une face appartiennent tous à une autre cellule,
+c'est que la face qu'ils forment et mitoyenne. Sinon, elle est en
+surface.
+
+[22/37]
+Les points de surface sont ensuite plaqués sur le point du maillage
+surfacique le plus proche.
+
+[23/37]
+En déplaçant les points, on accumule l'information du déplacement de
+chaque point dans le voisinage dans un certain rayon, pondéré par la
+distance.
+
+Le rayon est proportionel à la distance de déplacement multipliée par
+un facteur, ce qui assure qu'on ne produit pas de cellules inversées
+dites « chaussettes ».
+
+Ceci est similaire à la « modification proportionelle » de blender.
+
+[25/37] et [26/37]
+Influence du facteur de multiplication du rayon d'action sur la
+qualité. On voit qu'un rayon de 3 semble maximiser les angles les plus
+faibles, sans trop réduire les angles les plus élevés pour autant.
+
+[27/37]
+Relaxation : on déplace chaque point dans le barycentre de son
+voisinage. Une itération améliore la qualité du maillage, plus d'une
+itération la dégrade.
+
++------------+
+| Partie 3/3 |
++------------+
+
+implémentation, moi
+
++------------+
+| Partie 4/3 |
++------------+
+
+résultats et conclusion

rapport/rapport.bib

@@ -0,0 +1,14 @@
+@phdthesis{rouxellabbe:tel-01419457,
+  TITLE = {{Anisotropic mesh generation}},
+  AUTHOR = {Rouxel-Labbe, Mael},
+  URL = {https://tel.archives-ouvertes.fr/tel-01419457},
+  NUMBER = {2016AZUR4150},
+  SCHOOL = {{Universit{\'e} C{\^o}te d'Azur}},
+  YEAR = {2016},
+  MONTH = Dec,
+  KEYWORDS = {Riemannian Geometry ; Voronoi diagram ; Delaunay triangulation ; Mesh generation ; Anisotropy ; G{\'e}n{\'e}ration de maillages ; G{\'e}ometrie Riemannienne ; Diagramme de Voronoi ; Triangulation de Delaunay ; Anisotropie},
+  TYPE = {Theses},
+  PDF = {https://tel.archives-ouvertes.fr/tel-01419457v2/file/2016AZUR4150.pdf},
+  HAL_ID = {tel-01419457},
+  HAL_VERSION = {v2},
+}

rapport/rapport.tex

@@ -0,0 +1,166 @@
+\documentclass{gig}
+
+\usepackage[utf8]{inputenc}
+\usepackage[french]{babel}
+\usepackage[T1]{fontenc}
+\usepackage{url}
+\usepackage{caption}
+\usepackage{listings}
+\usepackage{graphicx}
+
+\title{Génération de maillages anisotropes}
+\author[Jérémy André, Cyril Colin et Christine Cozzolino]
+       {Jérémy André$^1$, Cyril Colin$^1$ et Christine Cozzolino$^1$
+         \\
+         $^1$Aix-Marseille Université}
+
+\begin{document}
+
+\maketitle
+
+\begin{abstract}
+  Nous étudions dans ce rapport la génération de maillages volumiques
+  anisotropes à partir de maillages surfaciques. Nous commençons par
+  une analyse des travaux de la thèse de Mael
+  Rouxel-Labbé\cite{rouxellabbe:tel-01419457}. Nous redéfinissons
+  ensuite les objectifs du projet, et présentons notre méthode pour y
+  répondre. Enfin, nous analysons les résultats obtenus.
+\end{abstract}
+
+\keywords{Génération de maillages, Maillages volumiques, Anisotropie}
+
+\section{Présentation du sujet}
+Les maillages volumiques sont des structures qui représentent la
+surface mais aussi l'intérieur d'un objet, contrairement aux maillages
+surfaciques classiques. Ils sont composés de cellules qui sont des
+polyèdres liés entre eux par leurs sommets, et peuvent être agencés de
+manière structurée ou non structurée. Un maillage volumique peut par
+exemple être un assemblage de voxels, auquel cas le maillage est
+structuré et les cellules sont des cubes.
+
+Une des applications de ces maillage est dans le domaine de la
+simulation, où l'on peut discrétiser un problème continu en se servant
+des cellules du maillage comme base de la discrétisation.
+
+Dans certains domaines, comme par exemple la mécanique des fluides,
+les conditions de la simulation rendent innéficace un maillage
+structuré : leur résolution fixe limite la précision dans les zones
+« intéressantes » comme celles ou un fluide intéragit avec une
+surface, tout en étant inutilement coûteuse dans les zones moins
+intéressantes.
+
+En plus d'augmenter la résolution dans les zones intéressantes, la
+précision des calculs peut aussi être augmentée en adaptant la forme
+des cellules aux conditions locales : il faut donc que celles-ci
+puissent être irrégulières pour pouvoir les transformer et étirer pour
+les adapter.
+
+Ainsi, l'utilisation de maillages volumiques non structurés s'est
+imposée comme solution pour tous les problèmes de ce type, posant
+cependant l'autre problème qu'est leur génération.
+
+\section{Présentation de la thèse}
+\section{Justification de la non-utilisation des travaux de la thèse}
+La thèse de Rouxel-Labbé offre une excellente introduction du domaine
+et des méthodes utilisées, mais les résultats des approches qu'elle
+développe sont peu satisfaisant.
+\begin{quote}
+  We have sought approaches that combine theoretical soundness, that
+  is the capacity to prove that the desired result is obtained under
+  certain conditions, and practicality. Our results are mixed: while
+  we have proposed a number of methods that were all provably correct,
+  none of these methods were universally practical.
+\end{quote}
+De plus, parmis les approches proposées, celle qui présente les
+meilleurs résultats à le défaut d'être lourde en temps de calcul sur
+des maillages surfaciques. L'adapter à des maillages volumiques la
+rendrait donc encore moins pratiquable.
+
+\section{Redéfinition des objectifs}
+Au vu de la complexité du problème, notre encadrant à décidé de
+redéfinir nos objectifs.
+
+Nous partons désormais d'un maillage volumique tétrahédrique existant
+mais qui approxime que grossièrement le maillage surfacique, et devons
+tenter de le rendre plus fidèle au maillage surfacique en évitant de
+dégrader la qualité.
+
+%-----------------------------------------------------------------------
+\section{Redéfinition de nos objectifs}
+La thèse de Rouxel-Labbé présente une excellente introduction du
+domaine en plus d'un état de l'art remarquable, mais les approches
+qu'elle propose restent assez limitée. L’auteur lui-même indique
+que bien que prouvablement correctes, elles sont limitées dans leurs
+applications. De plus, une des méthodes présentées qui fournit les
+meilleurs résultats à le défaut d'être coûteuse en temps de calcul,
+défaut qui ne ferait que s'aggraver si on la traduisait dans le
+domaine des maillages volumiques.
+
+Pour cette raison, et au vu de la complexité de la tâche, notre
+encadrant a redéfini nos objectifs : il s'agit désormais de prendre un
+maillage tétraédrique généré à partir d'un maillage surfacique de
+manière simpliste, et de l'adapter pour réduire l'écart entre le
+maillage surfacique et le maillage tétraédrique sans trop dégrader la
+qualité de celui-ci.
+%-----------------------------------------------------------------------
+\section{Résultats et validation}
+Notre méthode présente des résultats en demi-teinte. Elle réussi
+effectivement à grandement diminuer la distance entre le maillage
+polygonal source et le maillage tétraédrique, mais pas sans générer
+par endroit des zones de très faible qualité. De plus, elle ne peut
+pas toujours corriger les défauts causés par la tétraédrisation
+initiale.
+
+Les points les plus éloignés du maillage surfacique après application
+de notre méthode sont en effet ceux où trop de détails et de géométrie
+ont été perdus dans la tétrahédrisation (voir figure
+\ref{fig:cow-head}). On peut également y voir un défaut de notre
+méthode de projection qui choisit les points les plus proches.
+
+\begin{figure*}
+  \centering
+  \begin{tabular}{ccc}
+    \includegraphics[width=4.7cm]{img/cow-head-poly.png} &
+    \includegraphics[width=4.7cm]{img/cow-head-tet.png} &
+    \includegraphics[width=4.7cm]{img/cow-head-tet-out.png} \\
+  \end{tabular}
+  \caption{Tête de la vache dans le maillage polygonal source, le
+    maillage tétraédrique initial, et le maillage tétraédrique après
+    application de notre méthode, de gauche à droite respectivement.}
+  \label{fig:cow-head}
+\end{figure*}
+
+L'étape de relaxation améliore la qualité du maillage en minimisant le
+nombre de cellules qui ont des angles faibles. Cependant en faire plus
+d'une itération dégrade le maillage plus que de ne pas faire de
+relaxation du tout. Les angles moins faibles ne sont eux que peu
+affectés (Voir figure \ref{fig:courbe-relaxation}).
+
+\begin{figure*}
+  \centering
+  \includegraphics[width=14cm]{img/influence-relaxation.png} \\
+  \caption{Influence de la relaxation sur les angles dièdres}
+  \label{fig:courbe-relaxation}
+\end{figure*}
+
+Pour ce qui est du rayon d'action de la modification proportionelle,
+un multiplicateur de rayon de 3 réduit grandement la proportion de
+cellules avec des angles très faible, mais au-delà de ça l'effet est
+négligeable (Voir figure \ref{fig:courbe-rayon}).
+
+\begin{figure*}
+  \centering
+  \begin{tabular}{cc}
+    \includegraphics[width=8cm]{img/influence-rayon-0-3.png} &
+    \includegraphics[width=8cm]{img/influence-rayon-5-20.png} \\
+  \end{tabular}
+  \caption{Influence du rayon de modification proportionelle sur les
+    angles dièdres}
+  \label{fig:courbe-rayon}
+\end{figure*}
+
+\newpage
+\bibliography{rapport}
+\bibliographystyle{refig-alpha}
+
+\end{document}

rapport/slides.bib

@@ -0,0 +1,15 @@
+@article{plasma,
+author = {Helander, P and Beidler, C and Bird, T and Drevlak, M and Feng, Y. and Hatzky, R and Jenko, Frank and Kleiber, R. and Proll, Josefine and Turkin, Yu and Xanthopoulos, Panagiotis},
+year = {2012},
+month = {11},
+pages = {124009},
+title = {Stellarator and tokamak plasmas: A comparison},
+volume = {54},
+journal = {Plasma Physics and Controlled Fusion},
+doi = {10.1088/0741-3335/54/12/124009}
+}
+
+@misc{gamma,
+author = {Frédéric Alauzet et al.},
+howpublished = {https://team.inria.fr/gamma/}
+}

rapport/slides.tex

@@ -0,0 +1,272 @@
+\documentclass{beamer}
+
+\usepackage[french]{babel}
+\usepackage[utf8]{inputenc}
+\usepackage[T1]{fontenc}
+\usepackage{csquotes}
+\usepackage[backend=biber]{biblatex}
+
+\addbibresource{slides.bib}
+
+\usetheme{JuanLesPins}
+\beamertemplatenavigationsymbolsempty
+\setbeamertemplate{footline}[frame number]
+
+\makeatletter
+\newlength\graphwidth
+\setlength\graphwidth{10.5cm}
+\newlength\graphheight
+\setlength\graphheight{5.5cm}
+\makeatother
+
+\title{Génération de maillages anisotropes}
+\author{Jérémy André, Cyril Colin, Christine Cozzolino}
+\institute{Aix-Marseille Université} \date{\today}
+\titlegraphic{\includegraphics[width=2.5cm]{img/logo-gig.png}}
+
+\begin{document}
+
+\begin{frame}
+  \titlepage
+\end{frame}
+
+\begin{frame}{Plan}
+  \tableofcontents[hideallsubsections]
+\end{frame}
+
+\AtBeginSection[] {
+  \begin{frame}{Plan}
+    \tableofcontents[currentsection]
+  \end{frame}
+}
+
+\section{Présentation du sujet}
+
+\begin{frame}{Maillages volumiques}
+  \centering
+  \begin{figure}
+    \begin{tabular}{cc}
+      \includegraphics[width=4.5cm]{img/cow-00.png} &
+      \includegraphics[width=4.5cm]{img/cow-03.png} \\
+      \includegraphics[width=4.5cm]{img/cow-06.png} &
+      \includegraphics[width=4.5cm]{img/cow-09.png} \\
+    \end{tabular}
+    \caption{Différentes coupes du maillage volumique d'une vache}
+  \end{figure}
+\end{frame}
+
+\begin{frame}{Anisotropie}
+  \centering
+  \begin{figure}
+    \includegraphics[height=\graphheight]{img/plasma.png}
+    \caption{Plasma dans un réacteur à fusion nucléaire \cite{plasma}}
+  \end{figure}
+\end{frame}
+
+\begin{frame}{Anisotropie}
+  \centering
+  \begin{figure}
+    \begin{tabular}{cc}
+      \includegraphics[width=4.5cm]{img/isotropique.png} &
+      \includegraphics[width=4.5cm]{img/anisotropique.png}
+    \end{tabular}
+    \caption{À gauche : plan isotropique. À droite : plan anisotropique circulaire.}
+  \end{figure}
+\end{frame}
+
+\begin{frame}{Applications}
+  \begin{figure}[htb]
+    \includegraphics[width=9cm]{img/symmetry-plane.png}
+    \caption{\label{fig:inria-gamma}
+      Image issue du projet GAMMA d'Inria \cite{gamma}
+    }
+  \end{figure}
+\end{frame}
+
+
+\section{Thèse de Rouxel-Labbé}
+
+\begin{frame}
+  \begin{itemize}
+  \item Excellente introduction et état de l'art
+  \item Plusieurs approches proposées
+  \item Des résultats en demi-teinte
+  \end{itemize}
+\end{frame}
+
+
+\section{Redéfinition du sujet}
+
+\begin{frame}
+  Le sujet étant jugé trop complexe, notre encadrant nous guide vers
+  un sujet un peu simplifié.
+\end{frame}
+
+\begin{frame}
+  \titlepage
+\end{frame}
+
+\begin{frame}[plain]{}
+  \hspace*{-12mm}
+  \includegraphics[width=\paperwidth]{img/explosion.png}
+\end{frame}
+
+\title{Approximation et amélioration de la qualité des maillages tétrahédriques}
+\begin{frame}
+  \titlepage
+\end{frame}
+
+\begin{frame}
+  \begin{block}{Nouveau sujet}
+    Adapter un maillage tétrahédrique simple existant au maillage
+    surfacique utilisé pour le générer, tout en préservant sa qualité.
+  \end{block}
+\end{frame}
+
+\section{Présentation de notre méthode}
+
+\begin{frame}
+  \centering
+  \begin{figure}
+    \includegraphics[height=\graphheight]{img/pipeline2.png}
+    \caption{Chaîne de traitement de notre approche}
+  \end{figure}
+\end{frame}
+
+\subsection{Pré-traitement}
+
+\begin{frame}
+  \begin{itemize}
+  \item Suppression des cellules externes
+  \item Détection des points de surface
+  \end{itemize}
+\end{frame}
+
+\begin{frame}
+  \centering
+  \begin{figure}
+    \includegraphics[height=\graphheight]{img/point_in_polygon.png}
+    \caption{Point-dans-polygone}
+  \end{figure}
+\end{frame}
+
+\begin{frame}{Suppression des cellules externes}
+  \centering
+  \begin{figure}
+    \includegraphics[width=\graphwidth]{img/remove_external.png}
+    \caption{Cellules supprimées}
+  \end{figure}
+\end{frame}
+
+\begin{frame}
+  \centering
+  \begin{figure}
+    \includegraphics[width=\graphwidth]{img/surface_points.png}
+    \caption{Détection des points de surface}
+  \end{figure}
+\end{frame}
+
+\subsection{Projection des points}
+
+\begin{frame}{Recherche des points les plus proches}
+  \centering
+  \begin{figure}
+    \includegraphics[height=\graphheight]{img/project_triangle_schema.png}
+    \caption{Méthode de recherche du point le plus proche}
+  \end{figure}
+\end{frame}
+
+\begin{frame}
+  \centering
+  \begin{figure}
+    \includegraphics[width=\graphwidth]{img/project.png}
+    \caption{Vache après projection des points}
+  \end{figure}
+\end{frame}
+
+\subsection{« Modification proportionelle »}
+
+\begin{frame}{Sans}
+  \centering
+  \includegraphics[height=\graphheight]{img/project_clip_no_displace.png}
+\end{frame}
+
+\begin{frame}{Avec}
+  \centering
+  \includegraphics[height=\graphheight]{img/project_clip_displace.png}
+\end{frame}
+
+\begin{frame}{Rayon d'influence}
+  \centering
+  \includegraphics[width=\graphwidth]{img/influence-rayon-0-3.png}
+\end{frame}
+
+\begin{frame}{Rayon d'influence}
+  \centering
+  \includegraphics[width=\graphwidth]{img/influence-rayon-5-20.png}
+\end{frame}
+
+\subsection{Relaxation}
+
+\begin{frame}
+  \centering
+  \includegraphics[width=\graphwidth]{img/influence-relaxation.png}
+\end{frame}
+
+
+\section{Implémentation}
+
+\begin{frame}
+  \centering
+  Bibliothèque de Kitware sous license libre pour la manipulation de
+  données scientifiques 3D.
+\end{frame}
+
+\begin{frame}{VTK}
+  \centering
+  \includegraphics[width=\graphwidth]{img/vtk-unstructuredgrid-0.png}
+\end{frame}
+
+\begin{frame}{VTK}
+  \centering
+  \includegraphics[width=\graphwidth]{img/vtk-unstructuredgrid-1.png}
+\end{frame}
+
+\begin{frame}{VTK}
+  \centering
+  \includegraphics[width=\graphwidth]{img/vtk-unstructuredgrid-2.png}
+\end{frame}
+
+\begin{frame}{VTK}
+  \centering
+  \includegraphics[width=\graphwidth]{img/vtk-unstructuredgrid-3.png}
+\end{frame}
+
+\begin{frame}{Paraview}
+  \centering
+  \begin{figure}
+    \includegraphics[height=5cm]{img/paraview.png}
+    \caption{Capture d'écran de Paraview}
+  \end{figure}
+\end{frame}
+
+
+\section{Résultats}
+
+\subsection{Cas dégénéré}
+
+\begin{frame}
+  \centering
+  \includegraphics[width=\graphwidth]{img/cas-échec.png}
+\end{frame}
+
+\begin{frame}
+  \centering
+  \includegraphics[width=\graphwidth]{img/cas-échec-2.png}
+\end{frame}
+
+\begin{frame}
+  \printbibliography
+\end{frame}
+
+\end{document}

src/analysis/aspect_ratio_filter.cc

@@ -1,18 +1,30 @@
 #include "aspect_ratio_filter.h"
+#include "vtkDataObject.h"
 
 #include <vtkUnstructuredGrid.h>
 #include <vtkPointData.h>
 #include <vtkCellData.h>
 #include <vtkDoubleArray.h>
 #include <vtkCellIterator.h>
+#include <vtkInformation.h>
+#include <vtkInformationVector.h>
 
 
 vtkStandardNewMacro(AspectRatioFilter);
 
 
+// Ensure there is an "angles", 12 floats tuple array in the DataCell.
 int AspectRatioFilter::FillInputPortInformation(int port, vtkInformation *info) {
 	vtkUnstructuredGridAlgorithm::FillInputPortInformation(port, info);
-	// Ensure there is an "angles" array in the DataCell.
+	vtkNew<vtkInformation> anglesField;
+	anglesField->Set(vtkDataObject::FIELD_ASSOCIATION(), vtkDataObject::FIELD_ASSOCIATION_CELLS);
+	anglesField->Set(vtkDataObject::FIELD_NAME(), "angles");
+	anglesField->Set(vtkDataObject::FIELD_NUMBER_OF_COMPONENTS(), 12);
+	anglesField->Set(vtkDataObject::FIELD_ARRAY_TYPE(), VTK_DOUBLE);
+
+	vtkNew<vtkInformationVector> fields;
+	fields->Append(anglesField);
+	info->Set(vtkAlgorithm::INPUT_REQUIRED_FIELDS(), fields);
 	return 1;
 }
 

src/analysis/aspect_ratio_filter.h

@@ -6,6 +6,12 @@
 #include <vtkIdList.h>
 
 
+/*
+  Computes the worst triangle aspect ratio for each cell and stores it
+  in a "aspect_ratio", 1 float cell data array. Requires the input to
+  have an "angles", 12 float cell data array containing the angles of
+  every triangle for each cell, grouped by triangle.
+*/
 class AspectRatioFilter : public vtkUnstructuredGridAlgorithm {
 public:
 	static AspectRatioFilter *New();

src/analysis/dihedral_angles_filter.cc

@@ -0,0 +1,82 @@
+#include "dihedral_angles_filter.h"
+
+#include <algorithm>
+#include <limits>
+#include <vtkUnstructuredGrid.h>
+#include <vtkPointData.h>
+#include <vtkCellData.h>
+#include <vtkDoubleArray.h>
+#include <vtkCellIterator.h>
+#include <vtkTriangle.h>
+
+
+vtkStandardNewMacro(DihedralAnglesFilter);
+
+
+vtkTypeBool DihedralAnglesFilter::RequestData(
+	vtkInformation *request,
+	vtkInformationVector **inputVector,
+	vtkInformationVector *outputVector) {
+	(void) request;
+	vtkUnstructuredGrid* input =
+		vtkUnstructuredGrid::GetData(inputVector[0]);
+	vtkUnstructuredGrid* output =
+		vtkUnstructuredGrid::GetData(outputVector);
+	output->CopyStructure(input);
+	output->GetPointData()->PassData(input->GetPointData());
+	vtkCellData *cellData = output->GetCellData();
+	cellData->PassData(input->GetCellData());
+
+	vtkNew<vtkDoubleArray> dihedralAnglesArray;
+	dihedralAnglesArray->SetName("dihedral_angles");
+	dihedralAnglesArray->SetNumberOfComponents(6);
+	dihedralAnglesArray->SetNumberOfTuples(input->GetNumberOfCells());
+	double *dihedralAnglesBase = dihedralAnglesArray->GetPointer(0);
+
+	vtkNew<vtkDoubleArray> minDihedralAngleArray;
+	minDihedralAngleArray->SetName("min_dihedral_angle");
+	minDihedralAngleArray->SetNumberOfTuples(input->GetNumberOfCells());
+	double *minDihedralAngleBase = minDihedralAngleArray->GetPointer(0);
+
+	AverageMinDegrees = 0;
+	MinMinDegrees = std::numeric_limits<double>::infinity();
+
+	size_t i = 0;
+	vtkCellIterator *it = input->NewCellIterator();
+	for (it->InitTraversal(); !it->IsDoneWithTraversal(); it->GoToNextCell()) {
+		double a[3], b[3], c[3], d[3];
+		vtkPoints *points = it->GetPoints();
+		points->GetPoint(0, a);
+		points->GetPoint(1, b);
+		points->GetPoint(2, c);
+		points->GetPoint(3, d);
+
+		double nacb[3], nadc[3], nabd[3], nbdc[3];
+		vtkTriangle::ComputeNormal(a, c, b, nacb);
+		vtkTriangle::ComputeNormal(a, d, c, nadc);
+		vtkTriangle::ComputeNormal(a, b, d, nabd);
+		vtkTriangle::ComputeNormal(b, d, c, nbdc);
+
+		dihedralAnglesBase[i+0] = vtkMath::AngleBetweenVectors(nacb, nadc);
+		dihedralAnglesBase[i+1] = vtkMath::AngleBetweenVectors(nadc, nabd);
+		dihedralAnglesBase[i+2] = vtkMath::AngleBetweenVectors(nabd, nacb);
+		dihedralAnglesBase[i+3] = vtkMath::AngleBetweenVectors(nbdc, nacb);
+		dihedralAnglesBase[i+4] = vtkMath::AngleBetweenVectors(nbdc, nadc);
+		dihedralAnglesBase[i+5] = vtkMath::AngleBetweenVectors(nbdc, nabd);
+
+		minDihedralAngleBase[i / 6] =
+			*std::min_element(dihedralAnglesBase + i,
+			                  dihedralAnglesBase + i + 6);
+		AverageMinDegrees += minDihedralAngleBase[i / 6];
+		MinMinDegrees = std::min(MinMinDegrees, minDihedralAngleBase[i / 6]);
+
+		i += 6;
+	}
+	it->Delete();
+	AverageMinDegrees /= input->GetNumberOfCells();
+	AverageMinDegrees = AverageMinDegrees * 180. / 3.141592653589793;
+	MinMinDegrees = MinMinDegrees * 180. / 3.141592653589793;
+	cellData->AddArray((vtkAbstractArray *) dihedralAnglesArray);
+	cellData->AddArray((vtkAbstractArray *) minDihedralAngleArray);
+	return true;
+}

src/analysis/dihedral_angles_filter.h

@@ -0,0 +1,25 @@
+#ifndef DIHEDRAL_ANGLES_FILTER_H
+#define DIHEDRAL_ANGLES_FILTER_H
+
+#include <vtkUnstructuredGridAlgorithm.h>
+#include <vtkIdList.h>
+
+
+class DihedralAnglesFilter : public vtkUnstructuredGridAlgorithm {
+public:
+	static DihedralAnglesFilter *New();
+	vtkTypeMacro(DihedralAnglesFilter, vtkUnstructuredGridAlgorithm);
+	vtkTypeBool RequestData(vtkInformation *request,
+							vtkInformationVector **inputVector,
+							vtkInformationVector *outputVector) override;
+	vtkGetMacro(AverageMinDegrees, double);
+	vtkGetMacro(MinMinDegrees, double);
+
+protected:
+	~DihedralAnglesFilter() override = default;
+	double AverageMinDegrees;
+	double MinMinDegrees;
+};
+
+
+#endif

src/analysis/max_distance_filter.cc

@@ -0,0 +1,75 @@
+#include "max_distance_filter.h"
+#include "../closest_polymesh_point.h"
+
+#include <limits>
+#include <vtkPolyData.h>
+#include <vtkInformation.h>
+#include <vtkKdTree.h>
+
+
+vtkStandardNewMacro(MaxDistanceFilter);
+
+
+MaxDistanceFilter::MaxDistanceFilter() {
+	SetNumberOfInputPorts(2);
+	SetNumberOfOutputPorts(0);
+}
+
+
+int MaxDistanceFilter::FillInputPortInformation(int port,
+                                                vtkInformation *info) {
+	if (port == 0) {
+		vtkPolyDataAlgorithm::FillInputPortInformation(port, info);
+	} else if (port == 1) {
+		info->Set(vtkDataObject::DATA_TYPE_NAME(), "vtkPolyData");
+	}
+	return 1;
+}
+
+
+vtkTypeBool MaxDistanceFilter::RequestData(vtkInformation *request,
+                                           vtkInformationVector **inputVector,
+                                           vtkInformationVector *outputVector) {
+	(void) request;
+	(void) outputVector;
+	vtkPolyData* input1 = vtkPolyData::GetData(inputVector[0]);
+	vtkPolyData* input2 = vtkPolyData::GetData(inputVector[1]);
+	vtkNew<vtkKdTree> tree1;
+	tree1->BuildLocatorFromPoints(input1->GetPoints());
+	vtkNew<vtkKdTree> tree2;
+	tree2->BuildLocatorFromPoints(input2->GetPoints());
+	vtkNew<vtkStaticCellLinks> links1;
+	links1->BuildLinks(input1);
+	vtkNew<vtkStaticCellLinks> links2;
+	links2->BuildLinks(input2);
+	MaxDist = 0;
+	for (vtkIdType i = 0; i < input1->GetNumberOfPoints(); i++) {
+		vtkIdType nCells = links1->GetNcells(i);
+		if (nCells == 0) continue;
+		double point[3];
+		input1->GetPoint(i, point);
+		double vec[3];
+		double dist;
+		closestPolyMeshPoint(input2, point, tree2, links2, vec, &dist);
+		if (dist > MaxDist) {
+			MaxDist = dist;
+			MaxId = i;
+			MaxInput = 0;
+		}
+	}
+	for (vtkIdType i = 0; i < input2->GetNumberOfPoints(); i++) {
+		vtkIdType nCells = links2->GetNcells(i);
+		if (nCells == 0) continue;
+		double point[3];
+		input2->GetPoint(i, point);
+		double vec[3];
+		double dist;
+		closestPolyMeshPoint(input1, point, tree1, links1, vec, &dist);
+		if (dist > MaxDist) {
+			MaxDist = dist;
+			MaxId = i;
+			MaxInput = 1;
+		}
+	}
+	return true;
+}

src/analysis/max_distance_filter.h

@@ -0,0 +1,27 @@
+#ifndef MAX_DISTANCE_FILTER_H
+#define MAX_DISTANCE_FILTER_H
+
+#include <vtkPolyDataAlgorithm.h>
+
+
+class MaxDistanceFilter : public vtkPolyDataAlgorithm {
+public:
+	static MaxDistanceFilter *New();
+	vtkTypeMacro(MaxDistanceFilter, vtkPolyDataAlgorithm);
+	MaxDistanceFilter();
+	int FillInputPortInformation(int, vtkInformation *info) override;
+	vtkTypeBool RequestData(vtkInformation *request,
+							vtkInformationVector **inputVector,
+							vtkInformationVector *outputVector) override;
+	vtkGetMacro(MaxDist, double);
+	vtkGetMacro(MaxId, vtkIdType);
+	vtkGetMacro(MaxInput, vtkIdType);
+
+protected:
+	double MaxDist;
+	vtkIdType MaxId;
+	vtkIdType MaxInput;
+};
+
+
+#endif

src/closest_polymesh_point.cc

@@ -0,0 +1,60 @@
+#include "closest_polymesh_point.h"
+#include "point_tris_dist.h"
+#include <limits>
+
+
+vtkIdType findClosestPoint(const double *point, vtkPointSet *pointSet) {
+	vtkIdType minPoint = 0;
+	double minDist2 = std::numeric_limits<double>::infinity();
+	for (vtkIdType i = 0; i < pointSet->GetNumberOfPoints(); i++) {
+		double point2[3];
+		pointSet->GetPoint(i, point2);
+		double dist2 = vtkMath::Distance2BetweenPoints(point, point2);
+		if (dist2 < minDist2) {
+			minDist2 = dist2;
+			minPoint = i;
+		}
+	}
+	return minPoint;
+}
+
+
+void closestPolyMeshPoint(vtkPolyData *polyMesh,
+                          const double *point,
+                          vtkKdTree *kdTree,
+                          vtkStaticCellLinks *links,
+                          double *toClosestPoint,
+                          double *distanceToClosestPoint) {
+	vtkIdType closestPolyMeshPoint;
+	{
+		double dummy;
+		// closestPolyMeshPoint = kdTree->FindClosestPoint((double *) point, dummy);
+		closestPolyMeshPoint = findClosestPoint(point, polyMesh);
+	}
+	vtkIdType *cellIds = links->GetCells(closestPolyMeshPoint);
+	vtkIdType nCells = links->GetNcells(closestPolyMeshPoint);
+	if (nCells == 0) {
+		*distanceToClosestPoint = -1;
+		return;
+	}
+	double minDist = std::numeric_limits<double>::infinity();
+
+	/* For each tri the closest point belongs to. */
+	for (vtkIdType i = 0; i < nCells; i++) {
+		vtkIdType cellId = cellIds[i];
+		double direction[3];
+		double dist = pointTriangleDistance(
+			(double *) point, polyMesh->GetCell(cellId), direction);
+		/* Find the closest one. */
+		if (dist < minDist) {
+			minDist = dist;
+			toClosestPoint[0] = direction[0];
+			toClosestPoint[1] = direction[1];
+			toClosestPoint[2] = direction[2];
+		}
+	}
+	vtkMath::MultiplyScalar(toClosestPoint, minDist);
+	*distanceToClosestPoint = minDist;
+}
+
+

src/closest_polymesh_point.h

@@ -0,0 +1,17 @@
+#ifndef CLOSEST_POLYMESH_POINT_H
+#define CLOSEST_POLYMESH_POINT_H
+
+#include <vtkPolyData.h>
+#include <vtkKdTree.h>
+#include <vtkStaticCellLinks.h>
+
+
+void closestPolyMeshPoint(vtkPolyData *polyMesh,
+                          const double *point,
+                          vtkKdTree *kdTree,
+                          vtkStaticCellLinks *links,
+                          double *toClosestPoint,
+                          double *distanceToClosestPoint);
+
+
+#endif

src/fitting/mesh_fit_filter.cc

@@ -0,0 +1,55 @@
+#include "mesh_fit_filter.h"
+
+#include <vtkUnstructuredGrid.h>
+#include <vtkPointData.h>
+#include <vtkCellData.h>
+#include <vtkDoubleArray.h>
+#include <vtkCellIterator.h>
+
+#include "../kd_tree.h"
+
+vtkStandardNewMacro(MeshFitFilter);
+
+vtkTypeBool MeshFitFilter::RequestData(
+	vtkInformation *request,
+	vtkInformationVector **inputVector,
+	vtkInformationVector *outputVector) {
+	(void) request;
+
+	vtkUnstructuredGrid* input =
+		vtkUnstructuredGrid::GetData(inputVector[0]);
+	vtkUnstructuredGrid* output =
+		vtkUnstructuredGrid::GetData(outputVector);
+    output->CopyStructure(input);
+	output->GetPointData()->PassData(input->GetPointData());
+
+
+    std::vector<KdTree::Tuple> points;
+
+	auto it = input->NewCellIterator();
+	for(it->InitTraversal(); !it->IsDoneWithTraversal(); it->GoToNextCell()) {
+		if(it->GetCellType() != VTK_TETRA) continue;
+		
+        vtkIdList *idList = it->GetPointIds();
+        
+        for(int i = 0; i < 4; ++i) {
+            vtkIdType id = idList->GetId(i);
+
+            double point[3];
+            input->GetPoint(id, point);
+
+            points.push_back({point, id});
+
+            //std::cout << "[" << point[0] << ", " << point[1] << ", " << point[2] << "] (" << id << ")\n";
+        }
+	}
+
+    KdTree kdTree;
+    kdTree.fill(points);
+
+    std::cout << "[0.3, 1.5, -0.1] => " << kdTree.query({0.3, 1.5, -0.1}) << std::endl;
+    std::cout << "[5.1, 1.1, 0.5] => " << kdTree.query({5.1, 1.1, 0.5}) << std::endl;
+    std::cout << "[0, 0, 1] => " << kdTree.query({0, 0, 1}) << std::endl;
+
+	return true;
+}

src/fitting/mesh_fit_filter.h

@@ -0,0 +1,19 @@
+#ifndef MESH_FIT_FILTER_H
+#define MESH_FIT_FILTER_H
+
+#include <vtkUnstructuredGridAlgorithm.h>
+#include <vtkIdList.h>
+
+class MeshFitFilter : public vtkUnstructuredGridAlgorithm {
+public:
+	static MeshFitFilter *New();
+	vtkTypeMacro(MeshFitFilter, vtkUnstructuredGridAlgorithm);
+	vtkTypeBool RequestData(vtkInformation *request,
+							vtkInformationVector **inputVector,
+							vtkInformationVector *outputVector) override;
+
+protected:
+	~MeshFitFilter() override = default;
+};
+
+#endif

src/fitting/project_surface_points_on_poly.cc

@@ -0,0 +1,168 @@
+#include "project_surface_points_on_poly.h"
+#include "../closest_polymesh_point.h"
+
+#include <vtkUnstructuredGrid.h>
+#include <vtkPointData.h>
+#include <vtkCellData.h>
+#include <vtkDoubleArray.h>
+#include <vtkCellIterator.h>
+#include <vtkInformation.h>
+#include <vtkInformationVector.h>
+#include <vtkDoubleArray.h>
+
+#include <limits>
+
+
+vtkStandardNewMacro(ProjectSurfacePointsOnPoly);
+
+
+ProjectSurfacePointsOnPoly::ProjectSurfacePointsOnPoly() {
+	SetNumberOfInputPorts(2);
+
+	RadiusScale = 2.;
+}
+
+
+int ProjectSurfacePointsOnPoly::FillInputPortInformation(int port, vtkInformation *info) {
+	if (port == 0) {
+		vtkUnstructuredGridAlgorithm::FillInputPortInformation(port, info);
+	} else if (port == 1) {
+		vtkNew<vtkInformation> input2;
+		input2->Set(vtkDataObject::DATA_TYPE_NAME(), "vtkPolyData");
+		info->Append(input2);
+	}
+	return 1;
+}
+
+
+void findPointsWithinRadius(double radius, vtkPointSet *pointSet,
+                            const double *point, vtkIdList *points) {
+	for (vtkIdType i = 0; i < pointSet->GetNumberOfPoints(); i++) {
+		double vec[3];
+		vtkMath::Subtract(point, pointSet->GetPoint(i), vec);
+		double dist = vtkMath::Norm(vec);
+		if (dist < radius) {
+			points->InsertNextId(i);
+		}
+	}
+}
+
+
+vtkSmartPointer<vtkIdList> removeSurfacePoints(vtkIdList *list,
+                                               vtkIntArray *isSurface) {
+	vtkNew<vtkIdList> affectedPoints;
+	affectedPoints->Allocate(list->GetNumberOfIds());
+	for (vtkIdType i = 0; i < list->GetNumberOfIds(); i++) {
+		vtkIdType id = list->GetId(i);
+		if (!isSurface->GetValue(id)) {
+			affectedPoints->InsertNextId(id);
+		}
+	}
+	return affectedPoints;
+}
+
+
+void ProjectSurfacePointsOnPoly::moveSurfacePoint(
+		vtkUnstructuredGrid *tetMesh,
+		vtkPolyData *polyMesh,
+		vtkIdType pointId,
+		std::vector<double> &motionVectors,
+		std::vector<unsigned> &numberOfAffectors,
+		vtkKdTree *polyMeshKdTree,
+		vtkKdTree *tetMeshKdTree,
+		vtkStaticCellLinks *links) {
+	double point[3];
+	tetMesh->GetPoint(pointId, point);
+	double direction[3];
+	double distance;
+	closestPolyMeshPoint(
+		polyMesh, point, polyMeshKdTree, links, direction, &distance);
+	vtkNew<vtkIdList> pointsInRadius;
+	findPointsWithinRadius(distance * RadiusScale, tetMesh, point, pointsInRadius);
+	auto affectedPoints = removeSurfacePoints(pointsInRadius, isSurface);
+	// tetMeshKdTree->FindPointsWithinRadius(distance * RadiusScale, point, affectedPoints);
+
+	for (vtkIdType j = 0; j < affectedPoints->GetNumberOfIds(); j++) {
+		vtkIdType affectedPointId = affectedPoints->GetId(j);
+		double affectedPoint[3];
+		tetMesh->GetPoint(affectedPointId, affectedPoint);
+		double dist2 = vtkMath::Distance2BetweenPoints(affectedPoint, point);
+		double expected = (distance * RadiusScale) * (distance * RadiusScale);
+		if (dist2 > 2 * expected) {
+			std::cerr << dist2 << " " << expected << "\n";
+			// throw std::runtime_error("");
+			continue;
+		}
+
+		double motion[3] = {direction[0], direction[1], direction[2]};
+		double factor = 1. - std::sqrt(dist2) / (distance * RadiusScale);
+		vtkMath::MultiplyScalar(motion, factor);
+		motionVectors[affectedPointId * 3 + 0] += motion[0];
+		motionVectors[affectedPointId * 3 + 1] += motion[1];
+		motionVectors[affectedPointId * 3 + 2] += motion[2];
+		numberOfAffectors[affectedPointId]++;
+	}
+	vtkMath::Subtract(point, direction, point);
+	tetMesh->GetPoints()->SetPoint(pointId, point);
+	tetMesh->Modified();
+	affectedPoints->Reset();
+}
+
+
+static void applyMotionVectors(vtkUnstructuredGrid *tetMesh,
+                               std::vector<double> &motionVectors,
+                               const std::vector<unsigned> &numberOfAffectors) {
+	for (vtkIdType i = 0; i < tetMesh->GetNumberOfPoints(); i++) {
+		if(numberOfAffectors[i] == 0) continue;
+		double point[3];
+		tetMesh->GetPoint(i, point);
+		double *motion = motionVectors.data() + (i * 3);
+		vtkMath::MultiplyScalar(motion, 1. / numberOfAffectors[i]);
+		vtkMath::Subtract(point, motion, point);
+		tetMesh->GetPoints()->SetPoint(i, point);
+		tetMesh->Modified();
+	}
+}
+
+
+vtkTypeBool ProjectSurfacePointsOnPoly::RequestData(
+	vtkInformation *request,
+	vtkInformationVector **inputVector,
+	vtkInformationVector *outputVector) {
+	(void) request;
+	vtkUnstructuredGrid* tetMesh =
+		vtkUnstructuredGrid::GetData(inputVector[0]);
+	vtkPolyData* polyMesh = vtkPolyData::GetData(inputVector[1]);
+	vtkUnstructuredGrid* output = vtkUnstructuredGrid::GetData(outputVector);
+
+	output->CopyStructure(tetMesh);
+	output->GetPointData()->PassData(tetMesh->GetPointData());
+	output->GetCellData()->PassData(tetMesh->GetCellData());
+
+	/* Generate cell links, these tell us which cell a point belongs to. */
+	vtkNew<vtkStaticCellLinks> links;
+	links->BuildLinks(polyMesh);
+
+	vtkNew<vtkKdTree> polyMeshKdTree;
+	polyMeshKdTree->BuildLocatorFromPoints(polyMesh->GetPoints());
+	vtkNew<vtkKdTree> tetMeshKdTree;
+	tetMeshKdTree->BuildLocatorFromPoints(vtkPointSet::SafeDownCast(tetMesh));
+
+	vtkIdTypeArray *surfacePoints =
+		vtkIdTypeArray::SafeDownCast(tetMesh->GetFieldData()->GetArray("surfacePoints"));
+	isSurface = vtkIntArray::SafeDownCast(
+		output->GetPointData()->GetArray("isSurface"));
+	std::vector<double> motionVectors(3 * tetMesh->GetNumberOfPoints(), 0);
+	std::vector<unsigned> numberOfAffectors(tetMesh->GetNumberOfPoints(), 0);
+	for (vtkIdType i = 0; i < surfacePoints->GetNumberOfTuples(); i++) {
+		moveSurfacePoint(tetMesh, polyMesh,
+		                 surfacePoints->GetValue(i),
+		                 motionVectors,
+		                 numberOfAffectors,
+		                 polyMeshKdTree, tetMeshKdTree,
+		                 links);
+	}
+	applyMotionVectors(tetMesh, motionVectors, numberOfAffectors);
+
+	return true;
+}

src/fitting/project_surface_points_on_poly.h

@@ -0,0 +1,43 @@
+#ifndef PROJECT_SURFACE_POINTS_ON_POLY_H
+#define PROJECT_SURFACE_POINTS_ON_POLY_H
+
+#include <vtkUnstructuredGridAlgorithm.h>
+#include <vtkIdList.h>
+#include <vector>
+#include <vtkPolyData.h>
+#include <vtkKdTree.h>
+#include <vtkStaticCellLinks.h>
+
+class ProjectSurfacePointsOnPoly : public vtkUnstructuredGridAlgorithm {
+public:
+	static ProjectSurfacePointsOnPoly *New();
+	vtkTypeMacro(ProjectSurfacePointsOnPoly, vtkUnstructuredGridAlgorithm);
+	ProjectSurfacePointsOnPoly();
+	int FillInputPortInformation(int, vtkInformation *info) override;
+	vtkTypeBool RequestData(vtkInformation *request,
+							vtkInformationVector **inputVector,
+							vtkInformationVector *outputVector) override;
+	vtkSetMacro(affectedNeighborsCount, int);
+	vtkGetMacro(affectedNeighborsCount, int);
+
+	vtkSetMacro(RadiusScale, double);
+
+protected:
+	int affectedNeighborsCount = 0;
+	double RadiusScale;
+	vtkIntArray *isSurface = nullptr;
+
+	~ProjectSurfacePointsOnPoly() override = default;
+
+	void moveSurfacePoint(vtkUnstructuredGrid *tetMesh,
+                             vtkPolyData *polyMesh,
+                             vtkIdType pointId,
+                             std::vector<double> &motionVectors,
+                             std::vector<unsigned> &numberOfAffectors,
+                             vtkKdTree *polyMeshKdTree,
+                             vtkKdTree *tetMeshKdTree,
+                             vtkStaticCellLinks *links);
+};
+
+
+#endif

src/fitting/relaxation_filter.cc

@@ -0,0 +1,81 @@
+#include "relaxation_filter.h"
+
+#include <vtkUnstructuredGrid.h>
+#include <vtkPointData.h>
+#include <vtkCellData.h>
+#include <vtkDoubleArray.h>
+#include <vtkCellIterator.h>
+
+#include <set>
+
+vtkStandardNewMacro(RelaxationFilter);
+
+RelaxationFilter::RelaxationFilter()
+	:IterCount(1) {
+}
+
+vtkTypeBool RelaxationFilter::RequestData(
+	vtkInformation *request,
+	vtkInformationVector **inputVector,
+	vtkInformationVector *outputVector) {
+	(void) request;
+	vtkUnstructuredGrid* input =
+		vtkUnstructuredGrid::GetData(inputVector[0]);
+	vtkUnstructuredGrid* output = vtkUnstructuredGrid::GetData(outputVector);
+	output->CopyStructure(input);
+	output->GetPointData()->PassData(input->GetPointData());
+	output->GetCellData()->PassData(input->GetCellData());
+	vtkNew<vtkPoints> newPoints;
+	newPoints->DeepCopy(output->GetPoints());
+	vtkNew<vtkIdList> neighborCells;
+	vtkNew<vtkIdList> cellPoints;
+	vtkIdTypeArray *surfacePoints = vtkIdTypeArray::SafeDownCast(
+		output->GetFieldData()->GetArray("surfacePoints"));
+	vtkIntArray *isSurface = vtkIntArray::SafeDownCast(
+		output->GetPointData()->GetArray("isSurface"));
+	std::set<vtkIdType> neighbors;
+	output->BuildLinks();
+	double avg[3];
+
+	for (int iter = 0; iter < IterCount; iter++) {
+		for (vtkIdType i = 0; i < surfacePoints->GetNumberOfValues(); i++) {
+			vtkIdType id = surfacePoints->GetValue(i);
+
+			output->GetPointCells(id, neighborCells);
+
+			if (neighborCells->GetNumberOfIds() != 0) {
+
+				for (vtkIdType j = 0; j < neighborCells->GetNumberOfIds(); j++) {
+
+					vtkIdType cellId = neighborCells->GetId(j);
+					output->GetCellPoints(cellId, cellPoints);
+
+					for (vtkIdType k = 0; k < cellPoints->GetNumberOfIds(); k++) {
+
+						vtkIdType neighborId = cellPoints->GetId(k);
+						if (isSurface->GetValue(neighborId)) {
+							neighbors.insert(neighborId);
+						}
+					}
+					cellPoints->Reset();
+				}
+				neighbors.erase(neighbors.find(id));
+				if (neighbors.size() != 0) {
+					avg[0] = 0; avg[1] = 0; avg[2] = 0;
+					for (const vtkIdType &j : neighbors) {
+						double point[3];
+						output->GetPoints()->GetPoint(j, point);
+						vtkMath::Add(point, avg, avg);
+					}
+					vtkMath::MultiplyScalar(avg, 1. / neighbors.size());
+					newPoints->SetPoint(id, avg);
+				}
+			}
+			neighborCells->Reset();
+			neighbors.clear();
+		}
+
+		output->SetPoints(newPoints);
+	}
+	return true;
+}

src/fitting/relaxation_filter.h

@@ -0,0 +1,25 @@
+#ifndef RELAXATION_FILTER_H
+#define RELAXATION_FILTER_H
+
+
+#include <vtkUnstructuredGridAlgorithm.h>
+#include <vtkIdList.h>
+
+class RelaxationFilter : public vtkUnstructuredGridAlgorithm {
+public:
+	static RelaxationFilter *New();
+	vtkTypeMacro(RelaxationFilter, vtkUnstructuredGridAlgorithm);
+	RelaxationFilter();
+	vtkTypeBool RequestData(vtkInformation *request,
+							vtkInformationVector **inputVector,
+							vtkInformationVector *outputVector) override;
+	vtkSetMacro(IterCount, int);
+	vtkGetMacro(IterCount, int);
+
+protected:
+	~RelaxationFilter() override = default;
+	int IterCount;
+};
+
+
+#endif

src/fitting/remove_external_cells_filter.cc

@@ -0,0 +1,68 @@
+#include "remove_external_cells_filter.h"
+#include "../point_tris_dist.h"
+
+#include <vtkCellType.h>
+#include <vtkUnstructuredGrid.h>
+#include <vtkDataObject.h>
+#include <vtkPolyData.h>
+#include <vtkPointData.h>
+#include <vtkCellData.h>
+#include <vtkCellIterator.h>
+#include <vtkSelectEnclosedPoints.h>
+#include <vtkUnsignedCharArray.h>
+
+
+vtkStandardNewMacro(RemoveExternalCellsFilter);
+
+
+RemoveExternalCellsFilter::RemoveExternalCellsFilter() {
+	SetNumberOfInputPorts(2);
+}
+
+
+int RemoveExternalCellsFilter::FillInputPortInformation(
+	int port, vtkInformation *info) {
+	if (port == 0) {
+		info->Set(vtkDataObject::DATA_TYPE_NAME(), "vtkUnstructuredGrid");
+	} else if (port == 1) {
+		info->Set(vtkDataObject::DATA_TYPE_NAME(), "vtkPolyData");
+	}
+	return 1;
+}
+
+
+vtkTypeBool RemoveExternalCellsFilter::RequestData(
+	vtkInformation *request,
+	vtkInformationVector **inputVector,
+	vtkInformationVector *outputVector) {
+	(void) request;
+	vtkUnstructuredGrid* tetMesh = vtkUnstructuredGrid::GetData(inputVector[0]);
+	vtkPolyData *polyMesh = vtkPolyData::GetData(inputVector[1]);
+	vtkUnstructuredGrid* output = vtkUnstructuredGrid::GetData(outputVector);
+
+	vtkNew<vtkSelectEnclosedPoints> selectEnclosedPoints;
+	selectEnclosedPoints->SetInputDataObject(0, tetMesh);
+	selectEnclosedPoints->SetInputDataObject(1, polyMesh);
+	selectEnclosedPoints->Update();
+	vtkUnsignedCharArray *enclosedPoints =
+		vtkUnsignedCharArray::SafeDownCast(
+			selectEnclosedPoints->GetOutput()
+			->GetPointData()->GetArray("SelectedPoints"));
+
+	vtkNew<vtkCellArray> outCells;
+	auto it = tetMesh->NewCellIterator();
+	for (it->InitTraversal(); !it->IsDoneWithTraversal(); it->GoToNextCell()) {
+		vtkIdList *pointIds = it->GetPointIds();
+		if (enclosedPoints->GetValue(pointIds->GetId(0))
+		    || enclosedPoints->GetValue(pointIds->GetId(1))
+		    || enclosedPoints->GetValue(pointIds->GetId(2))
+		    || enclosedPoints->GetValue(pointIds->GetId(3))) {
+			outCells->InsertNextCell(it->GetPointIds());
+		}
+	}
+	it->Delete();
+	output->SetCells(VTK_TETRA, outCells);
+	output->SetPoints(tetMesh->GetPoints());
+
+	return true;
+}

src/fitting/remove_external_cells_filter.h

@@ -0,0 +1,29 @@
+#ifndef REMOVE_EXTERNAL_CELLS_FILTER_H
+#define REMOVE_EXTERNAL_CELLS_FILTER_H
+
+
+#include <vtkUnstructuredGridAlgorithm.h>
+#include <vtkIdList.h>
+#include <vtkInformation.h>
+
+
+/*
+   Removes cells from the UnstructuredGrid passed as first input that
+   lay outside the PolyData passed as second input.
+*/
+class RemoveExternalCellsFilter : public vtkUnstructuredGridAlgorithm {
+public:
+	static RemoveExternalCellsFilter *New();
+	vtkTypeMacro(RemoveExternalCellsFilter, vtkUnstructuredGridAlgorithm);
+	RemoveExternalCellsFilter();
+	int FillInputPortInformation(int port, vtkInformation *info) override;
+	vtkTypeBool RequestData(vtkInformation *request,
+							vtkInformationVector **inputVector,
+							vtkInformationVector *outputVector) override;
+
+protected:
+	~RemoveExternalCellsFilter() override = default;
+};
+
+
+#endif

src/kd_tree.cc

@@ -0,0 +1,83 @@
+#include "kd_tree.h"
+#include <algorithm>
+#include <limits>
+
+KdTree::Point::Point(): x{0}, y{0}, z{0} {}
+KdTree::Point::Point(double x, double y, double z): x{x}, y{y}, z{z} {}
+KdTree::Point::Point(double *p): x{p[0]}, y{p[1]}, z{p[2]} {}
+
+double KdTree::Point::operator[](std::size_t i) const {
+    if(i == 0) return x;
+    if(i == 1) return y;
+    return z;
+}
+
+double KdTree::Point::dist2(Point const &other) const {
+    double x_ = x - other.x;
+    double y_ = y - other.y;
+    double z_ = z - other.z;
+
+    return x_ * x_ + y_ * y_ + z_ * z_;
+}
+
+ostream& operator<<(ostream &os, KdTree::Point const &point) {
+    os << "[" << point.x << ", " << point.y << ", " << point.z << "]";
+    return os;
+}
+
+KdTree::Node::Node(Point const &position, vtkIdType id): position{position}, id{id} {}
+
+/////////////////////////////////////////////////////////////////////
+
+void KdTree::fill(std::vector<Tuple> &points) {
+    nodes.reserve(points.size());
+    for(std::size_t i = 0; i < points.size(); ++i)
+        nodes.push_back({points[i].first, points[i].second});
+
+    root = fillRec(0, points.size(), 0);
+}
+
+KdTree::Node *KdTree::fillRec(std::size_t begin, std::size_t end, int axis) {
+    if(end <= begin) return nullptr;
+
+    std::size_t n = begin + (end - begin) / 2;
+    auto i = nodes.begin();
+    std::nth_element(i + begin, i + n, i + end, [&](Node &p1, Node &p2) {
+        return p1.position[axis] < p2.position[axis];
+    });
+
+    axis = (axis + 1) % 3;
+    nodes[n].leftChild = fillRec(begin, n, axis);
+    nodes[n].rightChild = fillRec(n + 1, end, axis);
+
+    return &nodes[n];
+}
+
+KdTree::Point KdTree::query(Point const &point) {
+    bestDist = std::numeric_limits<double>::max();
+    bestNode = nullptr;
+    queryRec(root, point, 0);
+    return bestNode->position;
+}
+
+KdTree::Point KdTree::query(double *position) {
+    return query({position[0], position[1], position[2]});
+}
+
+void KdTree::queryRec(Node *node, Point const &point, int axis) {
+    if(node == nullptr) return;
+
+    double d = point.dist2(node->position);
+    if(d < bestDist) {
+        bestDist = d;
+        bestNode = node;
+    }
+
+    double dx = node->position[axis] - point[axis];
+
+    axis = (axis + 1) % 3;
+
+    queryRec(dx > 0 ? node->leftChild : node->rightChild, point, axis);
+    if(dx * dx >= bestDist) return;
+    queryRec(dx > 0 ? node->rightChild : node->leftChild, point, axis);
+}

src/kd_tree.h

@@ -0,0 +1,52 @@
+#ifndef KD_TREE_H
+#define KD_TREE_H
+
+#include <vector>
+#include <utility>
+
+#include <vtkIdList.h>
+
+// https://rosettacode.org/wiki/K-d_tree#C.2B.2B
+
+class KdTree {
+public:
+    KdTree() = default;
+
+    struct Point {
+        Point();
+        Point(double x, double y, double z);
+        Point(double *position);
+        double x, y, z;
+        double operator[] (std::size_t i) const;
+        double dist2(Point const &other) const;
+
+        friend ostream& operator<<(ostream &os, Point const &point);
+    };
+
+    using Tuple = std::pair<Point, vtkIdType>;
+    void fill(std::vector<Tuple> &points);
+    Point query(Point const &point);
+    Point query(double *point);
+
+private:
+    struct Node {
+        Node(Point const &position, vtkIdType id);
+
+        Point position;
+        vtkIdType id;
+
+        Node *leftChild;
+        Node *rightChild;
+    };
+
+    std::vector<Node> nodes;
+    Node *root;
+
+    double bestDist;
+    Node *bestNode;
+
+    Node *fillRec(std::size_t begin, std::size_t end, int axis);
+    void queryRec(Node *node, Point const &point, int axis);
+};
+
+#endif

src/main.cc

@@ -1,90 +1,152 @@
-#include "angles_filter.h"
-#include "aspect_ratio_filter.h"
+#include "analysis/dihedral_angles_filter.h"
+#include "analysis/max_distance_filter.h"
+#include "fitting/remove_external_cells_filter.h"
+#include "fitting/project_surface_points_on_poly.h"
+#include "fitting/relaxation_filter.h"
+#include "surface_points_filter.h"
 
+#include <vtkGeometryFilter.h>
+#include <vtkCellArrayIterator.h>
 #include <vtkCellData.h>
+#include <vtkDataSetReader.h>
+#include <vtkOBJReader.h>
+#include <vtkPolyData.h>
+#include <vtkPolyDataReader.h>
+#include <vtkSetGet.h>
 #include <vtkUnstructuredGrid.h>
+#include <vtkUnstructuredGridReader.h>
+#include <vtkXMLPolyDataReader.h>
+#include <vtkXMLUnstructuredGridReader.h>
+#include <vtkXMLUnstructuredGridWriter.h>
+#include <vtksys/SystemTools.hxx>
+#include <string>
+
+#ifdef USE_VIEWER
 #include <vtkCamera.h>
 #include <vtkNamedColors.h>
-#include <vtkVolumeProperty.h>
+#include <vtkOpenGLProjectedTetrahedraMapper.h>
+#include <vtkPiecewiseFunction.h>
 #include <vtkRenderWindow.h>
 #include <vtkRenderWindowInteractor.h>
 #include <vtkRenderer.h>
-#include <vtkVolumeMapper.h>
 #include <vtkVolume.h>
-#include <vtkOpenGLProjectedTetrahedraMapper.h>
-#include <vtkUnstructuredGridReader.h>
-#include <vtkPiecewiseFunction.h>
-#include <vtkDoubleArray.h>
+#include <vtkVolumeMapper.h>
+#include <vtkVolumeProperty.h>
+#endif // USE_VIEWER
 
 #include <array>
-#include <vector>
-#include <algorithm>
 
 
-// template <typename T>
-// T average(const std::vector<T> &data) {
-// 	T avg = 0;
-// 	for (const T &t : data) {
-// 		avg += t;
-// 	}
-// 	return avg / data.size();
-// }
+vtkAlgorithm *readerFromFileName(const char *fileName) {
+	std::string extension =
+		vtksys::SystemTools::GetFilenameLastExtension(fileName);
+	if (extension == ".vtk") {
+		auto reader = vtkDataSetReader::New();
+		reader->SetFileName(fileName);
+		return {reader};
+	} else if (extension == ".vtu") {
+		auto reader = vtkXMLUnstructuredGridReader::New();
+		reader->SetFileName(fileName);
+		return {reader};
+	} else if (extension == ".vtp") {
+		auto reader = vtkXMLPolyDataReader::New();
+		reader->SetFileName(fileName);
+		return {reader};
+	} else if (extension == ".obj") {
+		auto reader = vtkOBJReader::New();
+		reader->SetFileName(fileName);
+		return {reader};
+	}
+	throw std::runtime_error("Invalid file extension.");
+}
 
 
-// template <typename T>
-// T standard_deviation(const std::vector<T> &data) {
-// 	T avg = average(data);
-// 	T stddev = 0;
-// 	for (const T &t : data) {
-// 		stddev += (t - avg) * (t - avg);
-// 	}
-// 	return std::sqrt(stddev / data.size());
-// }
+static void usage(char **argv) {
+	std::cerr << "Usage: " << argv[0] << " analyze|fit tetmesh polydata [radiusScale relaxationIterCount]" << std::endl;
+}
 
 
 int main(int argc, char **argv) {
-	if (argc != 2) {
-		std::cerr << "Usage: " << argv[0] << " FILE.vtk" << std::endl;
+	if (argc < 4) {
+		usage(argv);
+		return EXIT_FAILURE;
+	}
+	bool fit = false;
+	if (std::string(argv[1]) == "fit") {
+		fit = true;
+	} else if (std::string(argv[1]) != "analyze") {
+		usage(argv);
+		return EXIT_FAILURE;
+	}
+	if ((fit && argc > 6) || (!fit && argc != 4)) {
+		usage(argv);
 		return EXIT_FAILURE;
 	}
 
-	/* Reader */
-	vtkNew<vtkUnstructuredGridReader> reader;
-	reader->SetFileName(argv[1]);
-
-
-	/* Angle filter */
-	vtkNew<AnglesFilter> anglesFilter;
-	anglesFilter->SetInputConnection(reader->GetOutputPort());
-
-	/* Display angles in console. */
-	anglesFilter->Update();
-	vtkUnstructuredGrid *grid = anglesFilter->GetOutput();
-	auto *angles = vtkDoubleArray::SafeDownCast(grid->GetCellData()->GetArray("angles"));
-	for (ssize_t i = 0; i < angles->GetNumberOfTuples(); i++) {
-		std::cout << "cell " << i << ": ";
-		for (ssize_t j = 0; j < angles->GetNumberOfComponents(); j++) {
-			std::cout << angles->GetTuple(i)[j] << ", ";
-		}
-		std::cout << "\b\b \n";
+	double radiusScale = 2.;
+	int iterCount = 1;
+	if(argc > 4) {
+		radiusScale = std::stod(argv[4]);
+	}
+	if (argc > 5) {
+		iterCount = std::stoi(argv[5]);
 	}
 
-	/* Aspect ratio */
-	vtkNew<AspectRatioFilter> aspectRatioFilter;
-	aspectRatioFilter->SetInputConnection(anglesFilter->GetOutputPort());
+	auto tetMeshReader = readerFromFileName(argv[2]);
+	auto polyMeshReader = readerFromFileName(argv[3]);
 
-	/* Display aspect ratios in console. */
-	aspectRatioFilter->Update();
-	grid = aspectRatioFilter->GetOutput();
-	auto *aspectRatios = vtkDoubleArray::SafeDownCast(grid->GetCellData()->GetArray("aspect_ratio"));
-	for (ssize_t i = 0; i < aspectRatios->GetNumberOfTuples(); i++) {
-		std::cout << "cell " << i << ": "
-				  << aspectRatios->GetTuple1(i) << "\n";
+	vtkNew<RemoveExternalCellsFilter> removeExternalCellsFilter;
+	removeExternalCellsFilter->SetInputConnection(0,
+		tetMeshReader->GetOutputPort());
+	removeExternalCellsFilter->SetInputConnection(1,
+		polyMeshReader->GetOutputPort());
+
+	vtkNew<SurfacePointsFilter> surfacePointsFilter;
+	surfacePointsFilter->SetInputConnection(
+		removeExternalCellsFilter->GetOutputPort());
+
+	vtkNew<ProjectSurfacePointsOnPoly> project;
+	project->SetRadiusScale(radiusScale);
+	project->SetInputConnection(0, surfacePointsFilter->GetOutputPort());
+	project->SetInputConnection(1, polyMeshReader->GetOutputPort());
+
+	vtkNew<RelaxationFilter> relaxationFilter;
+	relaxationFilter->SetIterCount(iterCount);
+	relaxationFilter->SetInputConnection(project->GetOutputPort());
+
+	vtkNew<DihedralAnglesFilter> dihedralAnglesFilter;
+	if (fit) {
+		dihedralAnglesFilter->SetInputConnection(
+			relaxationFilter->GetOutputPort());
+	} else {
+		dihedralAnglesFilter->SetInputConnection(
+			tetMeshReader->GetOutputPort());
 	}
 
+	vtkNew<vtkGeometryFilter> geometryFilter;
+	geometryFilter->SetInputConnection(dihedralAnglesFilter->GetOutputPort());
+
+	vtkNew<MaxDistanceFilter> maxDistanceFilter;
+	maxDistanceFilter->SetInputConnection(0, geometryFilter->GetOutputPort());
+	maxDistanceFilter->SetInputConnection(1, polyMeshReader->GetOutputPort());
+
+	vtkNew<vtkXMLUnstructuredGridWriter> writer;
+	writer->SetInputConnection(dihedralAnglesFilter->GetOutputPort());
+	writer->SetDataModeToAscii();
+	writer->SetFileName("out.vtu");
+	writer->Write();
+
+	maxDistanceFilter->Update();
+	std::cerr << "Max distance: " << maxDistanceFilter->GetMaxDist() << "\n"
+	          << "Max dist point id: " << maxDistanceFilter->GetMaxId()
+	          << " in " << (maxDistanceFilter->GetMaxInput() == 0 ? "tetMesh\n" : "polyMesh\n")
+	          << "Average min angle: " << dihedralAnglesFilter->GetAverageMinDegrees() << "\n"
+	          << "Min min angle: " << dihedralAnglesFilter->GetMinMinDegrees() << "\n";
+
+#ifdef USE_VIEWER
 	/* Volume rendering properties */
 	vtkNew<vtkOpenGLProjectedTetrahedraMapper> volumeMapper;
-	volumeMapper->SetInputConnection(aspectRatioFilter->GetOutputPort());
+	volumeMapper->SetInputConnection(externalPointsFilter->GetOutputPort());
 	vtkNew<vtkVolume> volume;
 	volume->SetMapper(volumeMapper);
 	vtkNew<vtkPiecewiseFunction> transferFunction;
@@ -110,6 +172,8 @@ int main(int argc, char **argv) {
 	renderWindowInteractor->SetRenderWindow(renderWindow);
 	renderWindow->Render();
 	renderWindowInteractor->Start();
-
+#endif
+	tetMeshReader->Delete();
+	polyMeshReader->Delete();
 	return EXIT_SUCCESS;
 }

src/point_tris_dist.cc

@@ -0,0 +1,95 @@
+#include "point_tris_dist.h"
+
+#include <vtkPoints.h>
+#include <vtkTriangle.h>
+#include <vtkMath.h>
+
+double pointSegmentDistance(double *p, double *a, double *b, double* direction) {
+    double segSqrLength = vtkMath::Distance2BetweenPoints(a, b);
+
+    double ap[3]; vtkMath::Subtract(p, a, ap);
+    double ab[3]; vtkMath::Subtract(b, a, ab);
+
+	if(segSqrLength <= 0) {
+        direction[0] = ap[0];
+        direction[1] = ap[1];
+        direction[2] = ap[2];
+        return vtkMath::Normalize(direction);
+    }
+
+    double h = vtkMath::ClampValue(vtkMath::Dot(ap, ab) / segSqrLength, 0., 1.);
+
+    vtkMath::MultiplyScalar(ab, h);
+    double v[3]; vtkMath::Subtract(ap, ab, v);
+
+    direction[0] = v[0];
+    direction[1] = v[1];
+    direction[2] = v[2];
+    return vtkMath::Normalize(direction);
+}
+
+double pointTriangleDistance(double *p, vtkCell *triangle, double *direction) {
+    double a[3]; triangle->GetPoints()->GetPoint(0, a);
+    double b[3]; triangle->GetPoints()->GetPoint(1, b);
+    double c[3]; triangle->GetPoints()->GetPoint(2, c);
+
+    double ab[3]; vtkMath::Subtract(b, a, ab);
+    double bc[3]; vtkMath::Subtract(c, b, bc);
+    double ca[3]; vtkMath::Subtract(a, c, ca);
+
+    double n[3]; vtkTriangle::ComputeNormal(a, b, c, n);
+    double vecTA[3]; vtkMath::Cross(ab, n, vecTA);
+    double vecTB[3]; vtkMath::Cross(bc, n, vecTB);
+    double vecTC[3]; vtkMath::Cross(ca, n, vecTC);
+
+    double ap[3]; vtkMath::Subtract(p, a, ap);
+    double bp[3]; vtkMath::Subtract(p, b, bp);
+    double cp[3]; vtkMath::Subtract(p, c, cp);
+
+    double da = vtkMath::Dot(vecTA, ap);
+	double db = vtkMath::Dot(vecTB, bp);
+	double dc = vtkMath::Dot(vecTC, cp);
+
+	if(da <= 0 && db <= 0 && dc <= 0) {
+        double d = vtkMath::Dot(n, ap);
+        double n2 = vtkMath::Dot(n, n);
+
+        if(d >= 0) {
+            direction[0] = n[0];
+            direction[1] = n[1];
+            direction[2] = n[2];
+        } else {
+            direction[0] = -n[0];
+            direction[1] = -n[1];
+            direction[2] = -n[2];
+        }
+        
+        return sqrt(d * d / n2);
+	} else {
+        double normalA[3];
+        double normalB[3];
+        double normalC[3];
+
+		double lab = pointSegmentDistance(p, a, b, normalA);
+		double lbc = pointSegmentDistance(p, b, c, normalB);
+		double lca = pointSegmentDistance(p, c, a, normalC);
+
+		double min = vtkMath::Min(lab, vtkMath::Min(lbc, lca));
+
+        if(min == lab) {
+            direction[0] = normalA[0];
+            direction[1] = normalA[1];
+            direction[2] = normalA[2];
+        } else if(min == lbc) {
+            direction[0] = normalB[0];
+            direction[1] = normalB[1];
+            direction[2] = normalB[2];
+        } else {
+            direction[0] = normalC[0];
+            direction[1] = normalC[1];
+            direction[2] = normalC[2];
+        }
+
+        return min;
+	}
+}

src/point_tris_dist.h

@@ -0,0 +1,10 @@
+#ifndef POINT_TRIS_DIST_H
+#define POINT_TRIS_DIST_H
+
+/* #include <vtkTriangle.h> */
+#include <vtkCell.h>
+#include <utility>
+
+double pointTriangleDistance(double *point, vtkCell *triangle, double *direction);
+
+#endif

src/surface_points_filter.cc

@@ -0,0 +1,84 @@
+#include "surface_points_filter.h"
+
+#include <vtkIdList.h>
+#include <vtkUnstructuredGrid.h>
+#include <vtkPointData.h>
+#include <vtkCellData.h>
+#include <vtkIdTypeArray.h>
+#include <vtkCellIterator.h>
+#include <vtkTriangle.h>
+#include <vtkStaticCellLinks.h>
+#include <vtkIntArray.h>
+
+#include <set>
+
+
+
+vtkStandardNewMacro(SurfacePointsFilter);
+
+
+vtkTypeBool SurfacePointsFilter::RequestData(
+	vtkInformation *request,
+	vtkInformationVector **inputVector,
+	vtkInformationVector *outputVector) {
+	(void) request;
+	vtkUnstructuredGrid* input =
+		vtkUnstructuredGrid::GetData(inputVector[0]);
+	vtkUnstructuredGrid* output =
+		vtkUnstructuredGrid::GetData(outputVector);
+
+	/* Copy point and cell data unmodified. */
+	output->CopyStructure(input);
+	output->GetPointData()->PassData(input->GetPointData());
+	output->GetCellData()->PassData(input->GetCellData());
+
+	/* Create an array to store the IDs of surface points. */
+	vtkNew<vtkIdTypeArray> surfacePoints;
+	surfacePoints->SetName("surfacePoints");
+	surfacePoints->SetNumberOfComponents(1);
+	vtkNew<vtkIntArray> isSurface;
+	isSurface->SetName("isSurface");
+	isSurface->SetNumberOfComponents(1);
+	isSurface->SetNumberOfTuples(input->GetNumberOfPoints());
+	isSurface->FillValue(0);
+
+	/* Generate cell links, these tell us which cell a point belongs to. */
+	vtkNew<vtkStaticCellLinks> links;
+	links->BuildLinks(input);
+
+	vtkNew<vtkIdList> neighborCells;
+	vtkNew<vtkIdList> facePoints;
+	std::set<vtkIdType> surfacePointsSet;
+	auto *it = input->NewCellIterator();
+	for (it->InitTraversal(); !it->IsDoneWithTraversal(); it->GoToNextCell()) {
+		vtkIdList *cellPointIds = it->GetPointIds();
+		for (vtkIdType faceId = 0; faceId < 4; faceId++) {
+			vtkIdType idA = cellPointIds->GetId((faceId + 0) % 4);
+			vtkIdType idB = cellPointIds->GetId((faceId + 1) % 4);
+			vtkIdType idC = cellPointIds->GetId((faceId + 2) % 4);
+			facePoints->InsertNextId(idA);
+			facePoints->InsertNextId(idB);
+			facePoints->InsertNextId(idC);
+			input->GetCellNeighbors(it->GetCellId(), facePoints,
+			                        neighborCells);
+			facePoints->Reset();
+			if (neighborCells->GetNumberOfIds() == 0) {
+				surfacePointsSet.insert(idA);
+				surfacePointsSet.insert(idB);
+				surfacePointsSet.insert(idC);
+				isSurface->SetValue(idA, 1);
+				isSurface->SetValue(idB, 1);
+				isSurface->SetValue(idC, 1);
+			}
+		}
+	}
+	it->Delete();
+	surfacePoints->Allocate(surfacePointsSet.size());
+	for (const vtkIdType &id : surfacePointsSet) {
+		surfacePoints->InsertNextValue(id);
+	}
+
+	output->GetPointData()->SetScalars(isSurface);
+	output->GetFieldData()->AddArray(surfacePoints);
+	return true;
+}

src/surface_points_filter.h

@@ -0,0 +1,22 @@
+#ifndef SURFACE_POINT_FILTER_H
+#define SURFACE_POINT_FILTER_H
+
+
+#include <vtkUnstructuredGridAlgorithm.h>
+#include <vtkIdList.h>
+
+
+class SurfacePointsFilter : public vtkUnstructuredGridAlgorithm {
+public:
+	static SurfacePointsFilter *New();
+	vtkTypeMacro(SurfacePointsFilter, vtkUnstructuredGridAlgorithm);
+	vtkTypeBool RequestData(vtkInformation *request,
+							vtkInformationVector **inputVector,
+							vtkInformationVector *outputVector) override;
+
+protected:
+	~SurfacePointsFilter() override = default;
+};
+
+
+#endif