# Detection of Orbital Floor Fractures by Principal Component Analysis

Daniel Krpelik, Milan Jaros, Marta Jarosova, Petr Strakos, Tereza Buresova,
Alena Vasatova, Tomas Karasek

► To cite this version:

Daniel Krpelik, Milan Jaros, Marta Jarosova, Petr Strakos, Tereza Buresova, et al. Detection of Orbital Floor Fractures by Principal Component Analysis. 15th IFIP International Conference on Computer Information Systems and Industrial Management (CISIM), Sep 2016, Vilnius, Lithuania. pp.129-138, 10.1007/978-3-319-45378-1_12. hal-01637516

HAL Id: hal-01637516

https://inria.hal.science/hal-01637516

Submitted on 17 Nov 2017

HAL is a multi-disciplinary open access
archive for the deposit and dissemination of sci-
entific research documents, whether they are pub-
lished or not. The documents may come from
teaching and research institutions in France or
abroad, or from public or private research centers.

L'archive ouverte pluridisciplinaire HAL, est
destinée au dépôt et à la diffusion de documents
scientifiques de niveau recherche, publiés ou non,
émanant des établissements d'enseignement et de
recherche français ou étrangers, des laboratoires
publics ou privés.

<figure><img src="image_1.png" /></figure>

Distributed under a Creative Commons Attribution 4.0 International License
# Detection of Orbital Floor Fractures by Principal Component Analysis

Daniel Krpelik¹, Milan Jaros¹, Marta Jarosova¹, Petr Strakos¹, Tereza Buresova², Alena Vasatova¹ and Tomas Karasek¹

¹IT4Innovations National Supercomputing Center, VSB-Technical University of Ostrava, Czech Republic

²The University Hospital in Ostrava, Czech Republic

tomas.karasek@vsb.cz

**Abstract.** Principal component analysis (PCA) is a statistical method based on orthogonal transformation, which is used to convert possibly correlated datasets into linearly uncorrelated variables called principal components. PCA is one of the simplest methods based on the eigenvector analysis. This method is widely used in many fields, such as signal processing, quality control or mechanical engineering. In this paper, we present the use of PCA in area of medical image processing. In the medical image processing with subsequent reconstruction of 3D models, data from sources such as Computed Tomography (CT) or Magnetic Resonance Imagining (MRI) are used. Series of images representing axial slices of human body are stored in Digital Imaging and Communications in Medicine (DICOM) format. Physical properties of different body tissues are characterized by different shades of grey of each pixel correlated to the tissue density. Properties of each pixel are then used in image segmentation and subsequent creation of 3D model of human organs. Image segmentation splits digital image into regions with similar properties which are later used to create 3D model. In many cases accurate detections of edges of such objects are necessary. This could be for example the case of a tumour or orbital fracture identification. In this paper, identification of the orbital fracture using PCA method is presented as an example of application of the method in the area of medical image processing.

## 1 Introduction

Despite CT and MRI are widely used by medical doctors for diagnostic purposes their potential is still not fully utilized. Those modalities provide 2D images of 3D objects and although medical doctors learn how to read and interpret them, they provide only partial information. Extraction of additional information such as size, area or volume of desired object has to be done manually and it is labour-intensive. Information about volume is crucial for example for surgeons, who are planning liver resections where remaining volume of liver is matter of life
or death. Volume of upper respiratory tract serves in diagnoses of sleep apnea. Size of orbital floor fracture is important when doctors are making decision whether patient should undergo surgery or conservative treatment. Those are just examples of additional information requested by doctors which cannot be obtained automatically from CT or MRI images.

Image processing could be an option how to reduce labour intensity, shorten the time of diagnostic processes and increase fidelity of the results. Above mentioned quantities are measured on 3D models which are created usually in several following steps. First, image noise is reduced by image filtering. Than image segmentation consisting of localization of objects and its boundaries is performed. Last step is 3D surface reconstruction.

If the 3D model is created by the mentioned steps, further detection of significant features like edges can be highly important. These features can localize major changes in overall shape of the model. While working with medical data, localization of such changes can lead, e.g., to direct detection of fractures in bones or tumours in soft tissues.

PCA method is used here to detect the points on the edge of the fracture. The edges are identified by calculating eigenvalues and eigenvectors of covariance matrix and establishing ratios between the calculated eigenvalues. In this paper PCA method is applied for detection of orbital floor fracture.

This paper is organized as follows. In Section 2 we describe example of identification of orbital floor fracture. In Section 3 PCA method is described followed by chapter 4 where object identification is explained. In chapter 5 details about fracture area reconstruction can be found. Chapter 6 concludes whole paper.

## 2 Model example

In this paper PCA is used for identification of orbital floor fracture from 3D virtual model of the human orbit. There is a very thin bone on the bottom side of the human eye orbit. Since this bone is the most fragile part of the orbit it is predisposed to fractures (see Fig. 1, Fig. 2). Fracture of this bone can be determined from CT scans, but to determine the extension of the fracture, only rough approximate methods have been used so far. These methods rely on manual determination of most distant points on fracture boundary. The fracture size can then be assessed by calculating the ratio of distance between two points of fracture boundary, to the distance between another two points, representing orbital floor boundary. All point are selected in one CT image where the fracture diameter is maximal [Ploder02]. Another way is to select points denoting fracture boundary in two perpendicular planes and approximate the fracture by geometrical primitives (e.g. ellipse). Area of such primitives can be determined from known equations [Goggin15]. In this approach CT image with maximal fracture diameter is also used. Although both above mentioned methods were shown to be inaccurate, they are widely used as they require little human expert intervention and less time compared to purely manual segmentation. Manual segmentation requires manual selection of all points denoting the fracture on
each individual CT image followed by area calculation from selected voxels. Our
implemented method is based on semi-automatic segmentation and keeps the
human expert intervention minimal.

<figure><img src="image_2.png" /><figcaption>**Fig. 1.** Model example: Left image depicts eye area denoted by yellow rectangle on CT scan image of a skull. Right image shows boundary (pink colour) of inner eye region identified by segmentation algorithm.</figcaption></figure>

<figure><img src="image_3.png" /><figcaption>**Fig. 2.** 3D reconstruction of orbital floor. Hole on the bottom side represents sought fracture.</figcaption></figure>

3 Covariance matrix and principal component analysis

First, let us introduce the basic terms which will be used later in this paper.
To distinguish, whether particular object is 1D, 2D or 3D we need to know the
relationship between points representing this object. Let P ⊂ R³ be the set of
all data points. For each point $p \in P$ set of neighbouring nodes could be found as

$$N_p = \left\{ q \mid p, q \in P, |pq| = \sqrt{(p_x - q_x)^2 + (p_y - q_y)^2 + (p_z - q_z)^2} < d_{3D} \right\}. \quad (1)$$

The example of set of points for which covariance matrix will be created is depicted on Fig. 3.

<figure><img src="image_4.png" /><figcaption>Fig. 3. Point distribution: a) linear, b) planar, c) spherical</figcaption></figure>

Recognition of geometric patterns in an image by the invariant moment was first introduced by 1962 Hu [Hu62]. In this work two-dimensional moment was used. Later, the invariant moment was extended to the third dimension by [Maas99] and [Gross02].

Let us consider the set $N_p$ to be solid in such a way that $B \subset \mathbb{R}^3$. Then $(i+j+k)$-th moment is defined by the equation

$$m_{ijk} = \iiint_B x^i y^j z^k f(x, y, z) dxdydz, \quad (2)$$

where $i, j, k \in \mathbb{Z}_0^+$ and $f(x, y, z)$ is the mass density function of solid $B$. Let $f(x, y, z) = 1$ for $\forall p \in P$.

If we set $(i+j+k) \le 2$, coordinates of the centre of gravity $\bar{p}$ of the solid $B$ could be defined as

$$\bar{p} = \left( \bar{x} = \frac{m_{100}}{m_{000}}, \bar{y} = \frac{m_{010}}{m_{000}}, \bar{z} = \frac{m_{001}}{m_{000}} \right), \quad (3)$$

where $m_{000}$ is the mass of the solid $B$ and $m_{100}, m_{010}, m_{001}$ are the static moments with respect to the coordinate planes $yz, xz, xy$.

The $(i+j+k)$-th centralized moment of the solid $B$ is defined by

$$\mu_{ijk} = \iiint_B (x - \bar{x})^i (y - \bar{y})^j (z - \bar{z})^k f(x, y, z) dxdydz. \quad (4)$$

To obtain central moment invariant to its size we have to normalize it
$$ \widehat{\mu_{ijk}} = \frac{\mu_{ijk}}{\mu_{000}} = \frac{\iiint_B (x-\bar{x})^i (y-\bar{y})^j (z-\bar{z})^k f(x,y,z) \, dxdydz}{\iiint_B (x-\bar{x})^0 (y-\bar{y})^0 (z-\bar{z})^0 f(x,y,z) \, dxdydz}. \quad (5) $$

After the numerical approximation of the normalized centralized moment we get

$$ \widehat{\mu_{ijk}} = \frac{\sum_{n=1}^{|N_p|} (x_n - \bar{x})^i (y_n - \bar{y})^j (z_n - \bar{z})^k}{\sum_{n=1}^{|N_p|} (x_n - \bar{x})^0 (y_n - \bar{y})^0 (z_n - \bar{z})^0}, \quad (6) $$

$$ \widehat{\mu_{ijk}} = \frac{1}{|N_p|} \sum_{n=1}^{|N_p|} (x_n - \bar{x})^i (y_n - \bar{y})^j (z_n - \bar{z})^k. \quad (7) $$

Finally, covariance matrix for each point $p \in P$ and its neighbours could be assembled as

$$ C_p = \begin{pmatrix} \widehat{\mu}_{200} & \widehat{\mu}_{110} & \widehat{\mu}_{101} \\ \widehat{\mu}_{110} & \widehat{\mu}_{020} & \widehat{\mu}_{011} \\ \widehat{\mu}_{101} & \widehat{\mu}_{011} & \widehat{\mu}_{002} \end{pmatrix}. \quad (8) $$

The matrix $C_p$ could be simplified by substituting the moments

$$ C_p = \frac{1}{|N_p|} \sum_{n=1}^{|N_p|} \begin{pmatrix} (x_n - \bar{x})^2 & (x_n - \bar{x})(y_n - \bar{y}) & (x_n - \bar{x})(z_n - \bar{z}) \\ (x_n - \bar{x})(y_n - \bar{y}) & (y_n - \bar{y})^2 & (y_n - \bar{y})(z_n - \bar{z}) \\ (x_n - \bar{x})(z_n - \bar{z}) & (y_n - \bar{y})(z_n - \bar{z}) & (z_n - \bar{z})^2 \end{pmatrix}, \quad (9) $$

$$ C_p = \frac{1}{|N_p|} \sum_{n=1}^{|N_p|} \left( \begin{bmatrix} x_n - \bar{x} \\ y_n - \bar{y} \\ z_n - \bar{z} \end{bmatrix} \begin{bmatrix} x_n - \bar{x} & y_n - \bar{y} & z_n - \bar{z} \end{bmatrix} \right), \quad (10) $$

$$ C_p = \frac{1}{|N_p|} \sum_{n=1}^{|N_p|} (q_n - \bar{p})^T (q_n - \bar{p}), \quad (11) $$

where $q \in N_p$ and $\bar{p}$ is the center of gravity of the set $N_p$.

For purpose of implementation, equation (11) could be rewritten to the form of matrix multiplication

$$ Q_p = \begin{bmatrix} x_1 & y_1 & z_1 \\ x_2 & y_2 & z_2 \\ x_3 & y_3 & z_3 \\ \vdots & \vdots & \vdots \\ x_{|N_p|} & y_{|N_p|} & z_{|N_p|} \end{bmatrix}, \quad J = \begin{bmatrix} 1 & 1 & \dots & 1 \\ 1 & 1 & \dots & 1 \\ 1 & 1 & \dots & 1 \\ \vdots & \vdots & \ddots & \vdots \\ 1 & 1 & \dots & 1 \end{bmatrix}, \quad (12) $$

$$ S_p = Q_p - J^T Q_p \frac{1}{|N_p|}, \quad (13) $$
$$C_p = S_p^T S_p, \qquad (14)$$

where $Q_p$ is the matrix composed of coordinates of a set $N_p$, matrix J is square matrix of all ones and the matrix multiplication $J^T Q_p \frac{1}{|N_p|}$ represents the center of the gravity of the set $N_p$.

The covariance matrix $C_p$ contains only real numbers, it is symmetric positive semi-definite and its dimension is $n = 3$. It means that all its eigenvalues are real with a complete set of the orthonormal eigenvectors. If this is a case we could say that there exists an orthogonal matrix $Q$ and a diagonal matrix $D$ such that

$$C_p = QDQ^T. \qquad (15)$$

The diagonal elements of the matrix $D$ are the eigenvalues $\lambda_1 \ge \lambda_2 \ge \lambda_3$ and columns of the matrix $Q$ are the orthonormal eigenvectors $v_1, v_2, v_3$ of the matrix $C_p$

$$C_p = (v_1^T \ v_2^T \ v_3^T) \begin{pmatrix} \lambda_1 & 0 & 0 \\ 0 & \lambda_2 & 0 \\ 0 & 0 & \lambda_3 \end{pmatrix} \begin{pmatrix} v_1 \\ v_2 \\ v_3 \end{pmatrix}. \qquad (16)$$

# 4 Object recognition

An auxiliary 3D model is acquired by prior segmentation of CT images. Surface of this model is then converted into cloud of points. Following task is to select points which lie on the fracture boundary. This is solved by the classification of neighbourhood of each single point in the cloud. Auxiliary model surface is a 2D object in 3D space. Points of the fracture boundary are determined based on violation of planarity in their neighbourhoods.

Object recognition, respectively classification whether object is 1D, 2D or 3D is based on properties of eigenvalues $\lambda_1 \ge \lambda_2 \ge \lambda_3$. Eigenvalues are obtained by spectral decomposition of matrix $C_p$ in equation (14) and they will represent the shape of the set $N_p$ by the ellipsoid [Gross02,Demantke11]. Ellipsoid representing linear, planar or spheric behaviour of neighbours is depicted in Fig. 3. If $\lambda_1 >> \lambda_2, \lambda_3 \cong 0$, the point $p$ is marked as linear. If $\lambda_1, \lambda_2 >> \lambda_3 \cong 0$, the point $p$ is assigned to the planar group. If $\lambda_1 \cong \lambda_2 \cong \lambda_3$, the point $p$ belongs to spheric group.

For the detection of a fracture we divide the points into two groups, $P_{planar}$ and $P_{boundary}$ using the criteria

$$a_{planar} = \frac{\lambda_3}{\sum_d \lambda_d}, \qquad (17)$$

$$P_{planar} = \{p | p \in P : a_{planar} \le \epsilon\}, \qquad (18)$$

$$P_{boundary} = \{p | p \in P : a_{planar} > \epsilon\}, \qquad (19)$$

where $\sum_d \lambda_d$ is the normalization coefficient and $\epsilon$ is the tolerance parameter. If $a_{planar} > \epsilon$, the point is on the edge and it is assigned to the group $P_{boundary}$. Otherwise the point is on the plane and it is assigned to the group $P_{planar}$.
<figure><img src="image_5.png" /><figcaption>**Fig. 4.** Conversion of voxels to point cloud (left) selection of cloud subset highlighted in green color (middle) and result of automatic boundary points extraction (right)</figcaption></figure>

Points on the fracture boundary are not the only ones, which are differen-
tiated by the feature calculated in equation (17). A subset of model surface
is chosen prior to the segmentation. Resulting segmented points are shown on
Fig. 4.

5 Reconstruction of fracture area

Fracture is in fact the absence of bone and as such it cannot be localized based
solely on pixel/voxel intensity values but rather on its relative position to the
other parts of human body. The course of work involves prior segmentation and
reconstruction of either undamaged bone area or eye orbit interior as sought
fracture boundary lies on the intersection of their surfaces. The eye orbit interior
is composed of soft tissue and in many cases leaks through the fracture into air
cavity bellow, creating a droplet-like structure. The fracture boundary itself is
then searched as a distinctive subsection of selected auxiliary object surface,
mainly determined by relatively extreme curvature in its neighbourhood.

This curvature is then exploited in two ways. Firstly, in case of both man-
ual and automated point selection, this curvature is used for reconstruction of
fracture boundary. Working hypothesis is that the shortest Jordan curve on the
object surface passing through selected points is a good approximation of frac-
ture boundary. Secondly, since the curvature is a distinctive feature, it may be
used to extract boundary points via automated segmentation of auxiliary model
vertices.

In applied method, we are reconstructing the fracture boundary. Fracture in-
terior itself is impossible to determine, because we cannot evaluate, which specific
surface shape is the correct one. Fracture interior is reconstructed heuristically
as a piecewise planar surface bounded by found fracture boundary. Alternative
methods could prefer smooth surfaces and reconstruct fracture interior for ex-
ample via solving Laplace problem with Dirichlet boundary.

Proposed method for Jordan curve localization assume that selected points
are vertices of auxiliary object mesh and are distributed sparsely around sought
boundary. The points in the case of manual point selection are selected on the
CT scan images and might not correspond with any vertex in the auxiliary object mesh. This is solved simply by selecting the closest vertex in the mesh for each manually selected point. In case of automated point selection we usually acquire a dense cloud of points, which denies requirement on sparsity. This might be solved automatically by, e.g., the k-means clustering, in which case we select output k-means centroids positions and treat those as a compression of initially dense cloud of points.

The following boundary reconstruction consists of two steps. First we construct a rough approximation of boundary consisting only of selected points and edges, and directly connecting them across Euclidean space. In terms of creating boundary, these points have to be properly ordered. This is a special case of travelling salesman problem for Euclidean metric space, which is generally NP-hard. Although a good polynomial approximation is known, we, for the sake of simplicity, use a greedy algorithm.

Our greedy algorithm starts with arbitrary triangle (a cycle with 3 vertices) and iteratively lengthens the cycle by inserting new point in between a pair of consecutive points from previous iteration. The pair is selected in such a way that the resulting new cycle is the shortest one among those considered (length being measured as a sum of distances between all consecutive points in the cycle). This algorithm is polynomial and given our assumptions, it yields good approximation of sought shortest cycle.

As a next step we refine found boundary by replacing edges in between consecutive points by the shortest paths between them on the auxiliary object mesh. This can be done effectively via Dijkstra algorithm [Dijkstra59] or some of its approximating variants. To ensure properties of Jordan curve a simple post-processing consisting of removing the inner cycles, which may occur because of independent runs of Dijkstra algorithms, is carried out.

Resulting reconstructed fracture area is visualized on Fig. 5 and Fig. 6.

<figure><img src="image_6.png" /><figcaption>**Fig. 5.** Fracture extraction - backward projection: manual point selection (pink), automatic point selection (yellow)</figcaption></figure>
<figure><img src="image_7.png" /><figcaption>**Fig. 6.** 3D visualization of fracture extraction : manual point selection (pink), automatic point selection (yellow)</figcaption></figure>

# 6 Conclusion

A method for edge recognition applied to detect orbital floor fracture was examined in this paper. The method is based on utilization of auxiliary 3D models obtained by segmentation of CT scan images. The problem of fracture reconstruction was converted to a problem of segmentation of auxiliary model mesh vertices and edges.

Proposed method requires selection of specific vertices on auxiliary model mesh. This selection was carried out both manually and automatically, via discriminating features based on PCA, and the results have been compared visually. Once the points were selected, algorithm chose appropriate subset of edges from auxiliary model mesh to reconstruct fracture boundary. These edges were selected so that they correspond with the shortest paths between selected points. This was carried out by Dijkstra algorithm.

Although the throughout validation of the size of orbital floor fracture based on semi-automatic segmentation is still pending, preliminary results based on expert based validation of fracture areas projected backward into CT images show promising results. PCA based point selection provides satisfactory results and thus it deserves further development. For example, the used discriminative feature, the relative magnitude of the smallest eigenvalue, is unable to fully capture all the parts of the fracture boundary. Finding better discriminative feature which could overcome this problem will be our future work.

# Acknowledgement

This work was supported by The Ministry of Education, Youth and Sports from the National Programme of Sustainability (NPU II) project "IT4Innovations excellence in science - LQ1602" and from the Large Infrastructures for Research, Experimental Development and Innovations project "IT4Innovations National
Supercomputing Center - LM2015070". All CT images used in this paper were
provided by The University Hospital in Ostrava.

References

[Hu62] Ming-Kuei Hu: *Visual pattern recognition by moment invariants*. Information Theory, IRE Transactions, Volume:8, Issue:2. 1962, p. 179-187.

[Maas99] Hans-Gerd Maas, George Vosselman: *Two algorithms for extracting build-
ing models from raw laser altimetry data*. ISPRS Journal of Photogrammetry &
Remote Sensing, 54, 1999, p. 153-163.

[Gross02] Hermann Gross, Ulrich Thoennessen: *Extraction of lines from laser point clouds*. ISPRS Commission III Photogrammetric Computer Vision PCV, 2006.

[Demantke11] Jerome Demantke, Clement Mallet, Nicolas David, Bruno Vallet: *Dimensionality based scale selection in 3D lidar point clouds*. ISPRS Commission III - WG III/2, 2011.

[Goggin15] James Goggin, Daniel C. Jupiter, Marcin Czerwinski: *Simple Computed Tomography-Based Calculations of Orbital Floor Fracture Defect Size Are Not Sufficiently Accurate for Clinical Use*, Journal of Oral and Maxillofacial Surgery, Volume 73, Issue 1, January 2015, Pages 112-116, ISSN 0278-2391.

[Ploder02] Oliver Ploder, Clemens Klug, Martin Voracek, Georg Burggasser, Christian Czerny: *Evaluation of computer-based area and volume measurement from coronal computed tomography scans in isolated blowout fractures of the orbital floor*, Journal of Oral and Maxillofacial Surgery, Volume 60, Issue 11, November 2002, Pages 1267-1272, ISSN 0278-2391.

[Dijkstra59] Dijkstra, E. W. *A note on two problems in connexion with graphs*. Num-
merische Mathematik 1: 269-271.