\chapter{Implementation}\label{four}
In previous chapter I have presented several methods I have chosen for visualization of facial data to meet the needs of scientists. In this chapter I am going to describe the process of implementation and techniques used for it.

\section{Surface superimposition}
For improvement of shape understandability, when comparing two models, I introduced method based on superimposition principle and transparency. Unfortunately, achieving proper, so called Order Independent Transparency for intersecting geometry is problematic, because all fragments corresponding to one pixel need to be sorted and blended in correct order. As I stated in section \ref{OIT} of chapter 2, there are several common methods for solving this well know problem.

The Depth Peeling and approximative techniques are limited to the extensions and modification. However, for the designed visualization techniques I needed access to information such as neighboring fragments or depth distance between fragments, which is not attainable by these techniques. Therefore as a basis for the implementation I chose Per Pixel Linked Lists constructed on GPU \cite{yang2010real}. I implemented an algorithm for Order Independent Transparency presented in \cite{shreiner2013opengl} and modified it, creating a rendering pipeline presented in Figure \ref{pip}.

\begin{figure}[h]
\centering
\scalebox{0.29}{\includegraphics{pictures/pipeline.png}}
\caption{Rendering pipeline.} \label{pip}
\end{figure} 

The pipeline consists of four basic steps, two of which are optional:
\begin{itemize}
\item \textbf{Depth map creation (optional).} Creation of depth map necessary for shadow-casting glyphs, glyph placement.
\item \textbf{Shading.} Shading of fragments, creation of linked lists, glyph placement.
 \item \textbf{Rendering.} Ordering of the linked lists, color and opacity modulations for visualization, computation of final color per pixel.
 \item \textbf{Contouring (optional).} Detection and rendering of the intersection contours.
\end{itemize}

The following sections describe the whole process of rendering in greater detail. All of the described methods were implemented in Java using JOGL -- the OpenGL wrapper library for Java -- and GLSL -- the OpenGL Shading Language. 

\subsection{Data preprocessing}
Before the models can be passed to rendering pipeline they need to be preprocessed. For the visualization of shadow-casting curvature glyphs the glyph placement and orientation need to be computed.

There are two typical approaches to placement of glyphs on surfaces -- either points of some significance can be selected (e.g. points with highest curvature value) or the glyphs can be evenly distributed over surface. Since the visualization technique should primarily express the differences between models on entire model surfaces, even distribution of glyphs is more suitable solution. The even distribution is approximated by algorithm which constructs a valued graph over the model mesh with each vertex storing its neighbors and distances to them. The algorithm iteratively selects a vertex as glyph center and marks the vertices in neighborhood that don't satisfy given distance condition along the graph edges (are closer to the glyph than glyph radius) as used, so they cannot be selected in following iterations.

At the selected vertices the principal curvature directions need to be determined in order to rotate the glyphs to indicate these directions. For this I used the normal averaging method described in \cite{interrante1997conveying}.

Firstly, for each glyph center vertex an orthogonal coordinate system ($\vec{u},\vec{v},\vec{n}$) must be selected. Here, vector $\vec{n}$ is the normal of the given vertex. The other two vectors can be arbitrary as long as the orthogonality is preserved. All the vectors must be of unit length. From this system the second fundamental tensor \textbf{\textrm{II}} can be defined: 
\begin{equation}
A = 
\begin{pmatrix}
 \frac{\partial \vec{n}}{\partial \vec{u}} \cdot \vec{u} &\frac{\partial \vec{n}}{\partial \vec{v}} \cdot \vec{u} \\
 \frac{\partial \vec{n}}{\partial \vec{u}} \cdot \vec{v} &\frac{\partial \vec{n}}{\partial \vec{v}} \cdot \vec{v} 
 \end{pmatrix}
.
\end{equation}

The components $\frac{\partial \vec{n}}{\partial \vec{a}} \cdot \vec{b}$ indicate the rate at which the surface normal tips in the direction of $\vec{a}$ when we move in the direction of $\vec{b}$. This tensor is defined and averaged for the close neighborhood of the vertex at which the curvature is to be computed. The directions of principal curvature can be then found by rotation of the orthogonal system so that the terms $\frac{\partial \vec{n}}{\partial \vec{v}} \cdot \vec{u}$ and $ \frac{\partial \vec{n}}{\partial \vec{u}} \cdot \vec{v}$ in the tensor $A$ "disappear" -- see Figure \ref{curv}. This is done by diagonalizing $A$ to obtain:

 \begin{equation}
D = 
\begin{pmatrix}
 \kappa_1 &0 \\
 0 & \kappa_2  
 \end{pmatrix}\quad
\text{and}\quad P=
 \begin{pmatrix}
   w_{uu} & w_{vu}  \\
    w_{uv} & w_{vv}
  \end{pmatrix}
,
\end{equation}
where |$\kappa_1$| > |$\kappa_2$| and $A= PDP^{-1}$. The columns of $P$ are the eigenvectors of $A$, while $\kappa_1$, $\kappa_2$ are the eigenvalues.

\begin{figure}[h]
\centering
\scalebox{0.7}{\includegraphics{pictures/curv.pdf}}
\caption{Curvature. Orthogonal system before rotation (black) and after rotation (green) indicating principal curvature directions.} \label{curv}
\end{figure} 

The principal curvature directions are given as:
 \begin{equation}
\begin{split}
\vec{u'} = w_{uu} \vec{u} + w_{uv} \vec{v}, \\
\vec{v'} = w_{vu} \vec{u} + w_{vv} \vec{v}.
\end{split}
\end{equation}

The vector $\vec{u'}$ indicates the direction of maximal principal curvature, while vector $\vec{v'}$ indicates the direction of minimal principal curvature. The values of the curvature are defined by $\kappa_1$ and $\kappa_2$ respectively.


\subsection{Depth map creation}
As soon as the computation of glyph centers and orientation is completed, the models can be passed to rendering pipeline. If the option of shadow-casting glyphs is enabled, the shadow placement must be computed. For this I chose basic shadow mapping approach presented in \cite{williams1978casting}, based on testing the visibility of fragment from the light position.

The geometry is rendered from the position of light in order to get depth of fragments nearest to the light. These depth values are stored in depth map -- a texture which is passed along the rendering pipeline to the point where the geometry is render from normal camera position. Then the depth value in the light space is computed for each fragment and it is compared with corresponding value in depth map. If the value in depth map is lower, the fragment lies in shadow. See Figure \ref{shadow} for example.

\begin{figure}[h]
\centering
\scalebox{0.7}{\includegraphics{pictures/shadow.pdf}}
\caption{Shadow mapping. At the location of fragment f, the surface (b) is shadowed by the surface (a).} \label{shadow}
\end{figure} 

In order to create the shadows of the glyphs, it must be tested at the time of rendering from light position whether the rendered fragment is part of the glyph or not.

My first idea for determining this was to generate texture coordinates for the vertices around glyph center, based on the vertex normal, curvature directions and distance from the vertex and map the glyph texture on the model. However, the algorithm based on this idea proved to be inefficient -- the texture coordinates computation taking up to several minutes. Therefore, I used the idea presented in \cite{interrante1997conveying} -- constructing a 3D representation of glyph around glyph center and testing whether fragment lies inside it -- see Figure \ref{slab} (a).

\begin{figure}[h]
\centering
\scalebox{0.5}{\includegraphics{pictures/slab.pdf}}
\caption{(a) 3D representation of glyph placed on the surface. (b)~Cuboid defining planes.} \label{slab}
\end{figure} 

The 3D representation of glyph consists of two crossing cuboids. Each cuboid can be defined by 6 planes (see Figure \ref{slab} (b)) and each of these planes is defined by its distance from glyph center and either by normal of the glyph center, or principal curvature direction vector. Therefore for each rendered triangle the closest glyph centers for its vertices along with corresponding normals and curvature vectors are passed to the fragment shader. These vectors and values are computed in object space. In order to operate in the same space, untransformed positions of vertices in object space need to be passed from vertex to fragment shader to obtain the fragment position in object space. 

Once the cuboid planes and the fragment position are defined in the same space, the test whether fragment lies inside the 3D glyph is a simple matter of testing the distances from fragment position to the planes.

If the rendered fragment lies inside the glyph, it is processed and its depth value may be written to the depth map, otherwise the fragment is discarded.

\subsection{Shading}
Once the depth map is created, the geometry is rendered from the camera position. In this phase shading of fragments and creation of Per Pixel Linked List takes place.

The fragments are shaded using Phong shading model \cite{phong1975illumination}. It is a shading model based on the interpolation of vertex attributes in fragment shader, thus achieving visually smooth results. It is therefore commonly used model for per pixel lighting. Figure \ref{light} shows the lighting scheme.

For each fragment the color computation consists of three parts -- ambient, diffuse and specular -- which are summed to get the final color. The ambient component (equation \ref{amb}) serves as rough approximation of indirect light reflections. It therefore doesn't depend on the positions, only on colors of the material and light. The diffuse component (equation \ref{diff}) is computed according to Lambert's law, which says the amount of reflected light is proportional to the cosine of the angle (a dot product) between normal vector of the surface and light direction vector. The specular component (equation \ref{spec}) specifies the direct reflections of light. The value depends on the cosine of the angle between reflection direction and eye (camera) direction from the surface.

\begin{figure}[h]
\centering
\scalebox{1.2}{\includegraphics{pictures/light.pdf}}
\caption{Lighting scheme.} \label{light}
\end{figure} 

 \begin{equation}\label{amb}
Color_{amb} = Light_{amb}\cdot Material_{amb}
\end{equation}
 \begin{equation}\label{diff}
Color_{diff} = Light_{diff} \cdot Material_{diff} \cdot cos(\alpha) 
\end{equation}
 \begin{equation}\label{spec}
\begin{split}
Color_{spec} = &Light_{spec} \cdot Material_{spec} \cdot cos(\beta)^s \\
&s = shininesses
\end{split}
\end{equation}


When the fragment is shaded, it is added to a fragment list. To be able to use the list for Order Independent Transparency, several things need to be set up. 

Firstly, there needs to be a way to access the first fragment list entry for each pixel. For this purpose a texture with the size of maximal expected image size is created. This texture is called Head Pointer texture and at each position in stores index of the last processed fragment corresponding to the given pixel position.

The shading of fragments and list building is computed on GPU, therefore several fragments will be processed in parallel. To precede concurrent modifications of the list, an Atomic Counter needs to be created. This counter keeps track of number of entries in the list, as it is increased every time new entry is to be added.

Finally, the list itself is implemented as \textit{RGB32UI} image buffer. This means that each fragment entry needs to be of \textit{uvec4} format -- a four component vector of \textit{uint} -- unsigned integer. 

A first component of each fragment entry will be the pointer to the index of next fragment at given pixel. As index of \textit{uint} format, this is a straightforward operation.

The second component should store the color of the fragment. However, the color is of \textit{vec4} format -- four component float vector. Fortunately, the function \textit{packUnorm4x8} defined in GLSL allows to pack four floating point values into 8-bit integers resulting in 32-bit unsigned integer. In this format the color can be added as the second component of the fragment entry. 

Third component consists of fragment depth. The depth is the float value, so it needs to be transformed to \textit{uint} using \textit{floatBitsToUint} function.

The fourth component is unused in original Order Independent Transparency implementation. However, for the implementation of the visualization techniques such as fogging, I needed additional information about fragment. I therefore used this space to store index of model to which the fragment belongs (passed to fragment shader as uniform), information about whether the fragment is part of the glyph (computed as described in previous section) and information about whether fragment lies in shadow. These three components with one additional unused zero component are represented as float values and packed using the \textit{packUnorm4x8} function again.

The complete scheme of linked list can be found in Figure \ref{list}. 

\begin{figure}[t]
\centering
\scalebox{0.75}{\includegraphics{pictures/list.pdf}}
\caption{Scheme of the linked list and structure of the fragment. Next -- the pointer to the next fragment at given pixel, RGBA -- color of the fragment, Depth -- depth of the fragment, M -- index of model to which the fragment belongs, X -- fragment is/isn't part of the glyph, S -- fragment is shadowed, -- unused} \label{list}
\end{figure} 

Once everything is set up, the process of storing a fragment in the list looks as follows:
\begin{enumerate}
\item Atomic Counter storing the number of fragments in the list is incremented. The new value is stored in \textit{index} variable. This is an atomic operation.
\item The value in Head Pointer texture corresponding to current pixel position is exchanged, the new value being \textit{index}. The old value is stored as \textit{old\_index}. This is an atomic operation as well.
\item The fragment is shaded and prepared for storing in the list as \textit{uvec4} -- the first component will be \textit{old\_index}, second packed color, third depth of the fragment and fourth the additional information for rendering.
\item The fragment information is added to the list at the position of \textit{index}.
\end{enumerate}

\subsection{Rendering}
When the list of fragments is complete, it is passed to the next step in rendering pipeline, where the final color for each pixel is computed. 

Prior to the color computation, the values in the Depth buffer after rendering models are analyzed in order to get their span -- this aids the later computations. Then the Color and Depth buffer are cleared and a rectangle covering the viewport is rendered, so that the new color can be computed for each pixel of viewport.

In fragment shader corresponding to this part of rendering pipe\-line, the list of fragments adherent to the pixel processed at given time needs to be extracted from the list of all fragments -- a reversal process to the list building. The index in Head Pointer texture at given pixel is looked up and an entry from the fragment list at this index is copied to the temporary list created in the shader. The first component of the copied entry points to the next fragment entry in the list corresponding to same pixel, so the next entry is retrieved and copied to the temporary list as well. The process is repeated as long as the first component of the newly retrieved entry points to valid index.

Once the per pixel list is retrieved, it is sorted according to depth value of each fragment (the third component of each fragment entry), so the nearest fragment is the first in the list. After this, the computation of color, based on the selected visualization options (passed to shader as uniform variables), may take place. 

If the option of making the inner surface opaque (see section \ref{transp_m} of chapter 3) is turned on, the alpha values of fragments need to be adjusted based on the inner-outer classification. The alpha value of first fragment is not modified -- it is user adjustable -- and the fragment is automatically considered outer. For the second fragment the test must be performed to decide whether it belongs to the same model as the first one, or is part of another model. If it is part of another model, the second fragment is turned opaque (alpha value set to 1), labeled inner and the fragments behind it may be discarded, since they will not be visible. If the second fragment belongs to the same model as the first one, it is usually the case of protuberant surface such as nose. In this case it is desirable to keep the alpha values unchanged, consider the fragment outer and proceed to process third fragment as if it was the first one. The process is repeated until one fragment is declared opaque, or the list of fragments ends. It should be noted that it will always be the the even (second, fourth, etc.) fragment that will be declared opaque -- see Figure \ref{nose}.

\begin{figure}[h]
\centering
\scalebox{0.1}{\includegraphics{pictures/nose.pdf}}
\caption{Example of case when the viewing ray intersect the surface belonging to same model several times. (a) Schematic view. (b) Rendered result.} \label{nose}
\end{figure} 

When the surfaces are classified as inner or outer, the solution for shadow-casting glyphs option is trivial. With each fragment the information about whether it is part of the glyph and whether it is shadowed is stored. So for fragments that were declared outer in previous step it is tested whether they are part of the glyph and if true, they are turned opaque. For the fragments that were declared inner it is tested whether they are shadowed, and if true, their brightness and saturation is decreased.

As for the fog, it is always function of depth difference between two fragments -- outer and inner, or, in special cases such as in Figure \ref{nose}, first and second, third and fourth fragment, etc. The $z_{ratio}$ which is used as depth parameter in computations is defined as following:
 \begin{equation}
z_{ratio} =k \cdot \frac{\lvert z_{a} - z_{b} \rvert}{\lvert z_{max} - z_{min} \rvert} ,
\end{equation}
where $z_{max}$ and $z_{min}$ are the maximum and minimum depth values detected when rendering geometry, $z_{a}$ and $z_{b}$ are the depth values of two fragments between which the fog is to be simulated and $k$ is a parameter for fog density.

There are three different fog implementations (see section \ref{fog_s} of chapter 3 for description):

\begin{itemize}
\item \textbf{Color overlay.} The first method highlights the distances by modifying the color of outer fragment -- the fragment nearer to the camera. It interpolates between fog color and surface color, using $z_{ratio}$ as interpolation parameter:
\begin{equation}\label{co}
Color_{new} = Color_{fragment}\cdot (1- z_{ratio}) + Color_{fog}\cdot {z_{ratio}}
\end{equation}
\item \textbf{Transparency mapping on outer surface.} With this method each outer fragment is colored to the color of the fog by changing the hue component of the HSV fragment color representation to hue of the fog. The saturation and value are thus preserved. The distance is then mapped on the alpha values (opacity) of outer fragment:
 \begin{equation}
Alpha_{new} = 0,2 \cdot Alpha_{fragment} + z_{ratio},
\end{equation}

 \item \textbf{Color mapping on inner surface.} The last method modifies the color of inner fragments in the same way as the first method modifies the color of outer fragments -- see equation \ref{co}.
\end{itemize}

After the aforementioned modifications to fragments are done, the final per pixel color may be computed. The variable $Color_{final}$ is first initialized to the color of the background. Then the sorted per pixel list of fragments is traversed, starting with the fragment furthest from camera and with each fragment the final color is modified in following way:
\begin{equation}
\begin{split}
Color_{final} = &Color_{final}\cdot (1- Alpha_{fragment}) + \\
				&+Color_{fragment}\cdot { Alpha_{fragment}}.
\end{split}
\end{equation}

At the end of the computation the pointers in original list of fragments are modified to represent the sorted version of per pixel lists.

\subsection{Contouring}
After the Rendering phase, the final image is nearly complete. The only remaining part is highlighting the model intersections. If this option is selected, one additional layer -- again a rectangle covering viewport -- is rendered over the image result from previous steps. In this layer only the pixels detected at the intersection points will be modified, the rest of them will be discarded and thus not interfering with the already rendered image.

To compute the intersection, the sorted list of fragments acquired in previous step is passed to the fragment shader. There, for the pixel corresponding to processed fragment the neighborhood of $3 \times 3$ pixels is searched. For each pixel within the neighborhood the first two fragments closest to the camera are looked up in the sorted list of fragments. Two conditions must be satisfied in order for the fragment to be declared as intersection point:

\begin{enumerate}
\item At least one of the first fragments in the neighborhood must belong to different model then the first fragment at the current pixel position.
\item The depth difference between first two fragments at given pixels must be lower than experimentally chosen threshold, to prevent false detections of intersections.
\end{enumerate}

If the fragment is identified as intersection point, it will be rendered, otherwise it will be discarded.

\section{Cross sections}
The technique I have designed for visualization of local shape variability shows the variability in a set of models at the place of their intersection with 3D plane. It was again implemented using JOGL library.

The base of this technique lies in computing the intersection of facial model with 3D plane. Since the meshes used in FIDENTIS Analyst application are triangulated, this problem reduces to computing an intersection of triangle with plane and this in turn reduces to intersection of line segment (triangle edge) with the plane. For a plane given by point $P$ and normal $\vec{n}$ and the line segment given by two points $A$ and $B$ -- see Figure \ref{lps}, it can be solved using following equations: 
\begin{equation}
\vec{u} = B-A \quad \vec{v} = A-P
\end{equation}
\begin{equation}
D = cos(\alpha) = \vec{n} \cdot \vec{u} \quad N = - cos(\beta)= - \vec{n} \cdot \vec{v} 
\end{equation}
\begin{figure}[t!]
\centering
\scalebox{0.58}{\includegraphics{pictures/plane.pdf}}
\caption{Line-plane intersection scheme.} \label{lps}
\end{figure} 
If $ cos(\alpha)$ is equal to zero, the segment is parallel to the plane. In that case if $cos(\beta)$ is equal to zero as well, the segment lies in the plane, otherwise there is no intersection. If segment is not parallel to the plane, the intersection point of line (given by the segment) with the plane must be computed and tested whether it lies within segment. The line corresponding to the segment AB may be represented as:
\begin{equation}
l = A + t \cdot \vec{u},
\end{equation}
and the $t$ parameter for line-plane intersection point is given as:
\begin{equation}
t = \frac{N}{D}.
\end{equation}
If $t \in <0, 1>$, then the segment intersects with the plane.

The described test is done for each edge of each triangle. The detected plane-triangle intersections -- points and line segments are stored in a list. When each triangle is processed, the intersecting line segments are connected at the points with same or very close coordinates to create polyline paths -- see Figure \ref{f} (a). Depending on the model and plane orientation, there may be one or more intersection paths per model.
 
These paths are then rotated to XY plane to be displayed in 2D view. As all the paths lie in same plane, the rotation corresponds to rotation of the plane normal to z-axis. The angle between the normal and z-axis can be computed from their dot product, while the rotation axis is given by their cross product:

\begin{equation}
 cos(\theta) = \vec{n} \cdot \vec{z} \quad \vec{u} = \vec{n} \times \vec{z}
\end{equation}
The rotation of each point to the XY plane is then done by following matrix for rotation around arbitrary axis $u$:
\begin{equation}\label{matica}
\scriptstyle 
R = 
\begin{bmatrix}
\scriptstyle 
cos\theta + u_x^2 \left ( 1- cos \theta \right ) & \scriptstyle u_xu_y \left ( 1- cos \theta \right ) - u_z sin \theta & \scriptstyle u_xu_z \left ( 1- cos \theta \right ) + u_y sin \theta\\ 
\scriptstyle u_yu_x \left ( 1- cos \theta \right ) + u_z sin \theta & \scriptstyle cos\theta + u_y^2 \left ( 1- cos \theta \right ) & \scriptstyle u_yu_z \left ( 1- cos \theta \right ) - u_x sin\\ 
\scriptstyle u_zu_x \left ( 1- cos \theta \right ) - u_y sin & \scriptstyle u_zu_y \left ( 1- cos \theta \right ) + u_x sin & \scriptstyle cos\theta + u_z^2 \left ( 1- cos \theta \right )
\end{bmatrix}
\end{equation}

Once in 2D, the intersection paths belonging to the model marked as primary (further referred to as primary paths) are sampled at equal distances and at each sample point a normal to the path is computed -- a trivial task, since the path consist of line segments. For each sampling point a ray given by the point and the normal to the primary path at that point is defined -- Figure \ref{f} (b). Then the intersections between these rays and paths of models from dataset are computed -- if there is more than one intersection for the model, the closest one is taken -- Figure \ref{f} (c). 
\begin{figure}[h]
\centering
\scalebox{0.3}{\includegraphics{pictures/faces.pdf}}
\caption{(a) Example of detected intersection path on triangular mesh. (b) Example of sampled intersection path. The point I on secondary model illustrates its intersection with the ray given by sampling point and normal to the primary path at that point. (c) Intersections with entire dataset. (d) Vectors indicating variability (scaled for illustration).} \label{f}
\end{figure} 

Finally, between each sampling point and intersections adherent to it (intersections of the line passing through the sampling point) a vector is computed. The vectors adherent to one sampling point are summed, the final vector providing information about variability at given sampling point -- Figure \ref{f} (d).


\section{Plots}
The last visualization technique I chose, consists of 2D plots for numerical data. The implementation of this technique was done using Java 2D API. 

For displaying the numerical results of Batch processing, the table-like plot is drawn. Numerical value is mapped on each cell of the table in a form of color. For computing of colors I operate in HSV color space, mapping the values onto hue component.

To compute the color which should represent given value, firstly the minimal and maximal values in the numerical results must be found. These values are then mapped onto boundary color values, which are by default red (hue = 0) for minimum and blue (hue = 240) for maximum. The hue color component is then computed according to following equation:

\begin{equation}
h = h_{min} + \frac{(v -v_{min})}{(v_{max}- v_{min})} \cdot (h_{max}- h_{min}),
\end{equation}
where $h$ denotes hue and $v$ denotes numerical values. 

The same color picking principle is also applied when displaying auxiliary results. The only additional calculations that needs to be done concern clustering the neighboring values to fit the whole dataset into plot area -- as there may be thousands of entries per model -- and reversal process for zooming and displaying the original values.

Both plots are drawn using Java 2D functions for drawing lines and rectangles and all interactivity must be achieved by tracking the mouse motion and computing its position in relation to interactive objects.