Detection of Orbital Floor Fractures by Principal Component Analysis

. 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 eigenvec-tor analysis. This method is widely used in many ﬁelds, 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 Res-onance 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 diﬀerent body tis-sues are characterized by diﬀerent 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 identiﬁcation. In this paper, identiﬁcation of the orbital fracture using PCA method is presented as an example of application of the method in the area of medical image processing.


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.

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.

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 3 be the set of all data points. For each point p ∈ P set of neighbouring nodes could be found as (1) The example of set of points for which covariance matrix will be created is depicted on Fig. 3.
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 ⊂ R 3 . Then (i + j + k)-th moment is defined by the equation where i, j, k ∈ Z + 0 and f (x, y, z) is the mass density function of solid B. Let f (x, y, z) = 1 for ∀p ∈ P .
If we set (i + j + k) ≤ 2, coordinates of the centre of gravityp of the solid B could be defined asp 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 To obtain central moment invariant to its size we have to normalize it After the numerical approximation of the normalized centralized moment we get Finally, covariance matrix for each point p ∈ P and its neighbours could be assembled as The matrix C p could be simplified by substituting the moments where q ∈ N p andp 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 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 1 |Np| 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 The diagonal elements of the matrix D are the eigenvalues λ 1 ≥ λ 2 ≥ λ 3 and columns of the matrix Q are the orthonormal eigenvectors v 1 , v 2 , v 3 of the matrix C p

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 λ 1 ≥ λ 2 ≥ λ 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 λ 1 >> λ 2 , λ 3 ∼ = 0, the point p is marked as linear. If λ 1 , λ 2 >> λ 3 ∼ = 0, the point p is assigned to the planar group. If λ 1 ∼ = λ 2 ∼ = λ 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 where d λ d is the normalization coefficient and is the tolerance parameter. If a planar > , 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 .

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 manual 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 fracture 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 interior 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 example 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 postprocessing 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.

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.