Automatic Localization of Condylar Process Centre Based on Generalized Hough Transform from CBCT images

Cone Beam Computed Tomography (CBCT) is a technique for capturing maxillofacial images that provide tissue views with low radiation doses. In this paper, we propose an automatic method to locate the centre of the mandibular condylar process in CBCT images. The bone edges of the condylar process are extracted to detect the mandibular ramus in each slice based on the generalized Hough Transform. We propose an upward region-growing method to segment the condylar and coracoid processes. The centre of the condylar process is then determined as the centre of the condylar region in the axial slice with the maximum condyle area. The centre located by the proposed method was compared with the centre manually located by two experts. Our results indicate that the proposed method can accurately locate the condylar process centre.


Introduction
The mandibular condylar process is on the mandible that articulates with the disk of the temporomandibular joint (TMJ). Its position should be considered in the diagnosis and treatment of condyle fracture, orthognathic surgery, occlusive reconstruction and temporomandibular disorders [1]. The diseases connected to TMJ is the well-known temporomandibular joint disorder (TMD) which is a kind of common diseases that happen to middle and old age people. Segmentation of TMJ and study on TMJ movement can help a doctor to understand normal movement and abnormal TMJ movement which causes people suffer TMD. Then the doctor can make correct physical treatment on TMJ [1]. The mandibular condylar process is the axis of the mandibular jaw rotation dynamically with translation to open and close mouth. Segmentation of TMJ and study TMJ movement play important roles in TMD treatment.
With the development of clinical applications of cone-beam computed tomography (CBCT) in dentistry, it can also be used to measure condyle position for three-dimensional display and high accuracy [2,3].
The condylar processes have diversified and complex morphologies such as the low condylar bone density, proximity of the discus articulation and the overshadowing glenoid fossa. It is a challenging task to locate the condylar process centre [4]. Segmentation of the condylar process is an essential step for the localization of the condylar process centre. A region-growing method was proposed to segment the condylar process from CBCT images but manual interactivity was needed [5].
In order to study TMJ movement, segmentation of TMJ and location of the condylar process centres are the first things to do. This paper proposes an automatic algorithm to locate the centre of the mandibular condylar process from CBCT images based on the generalized Hough transform. The idea is to use a typical mandibular ramus template on a two-dimensional slice to detect the ramus in each  [9]. Since the mandibular condylar process is roughly a sphere, the axial slice containing the maximum area will be selected and its centre will be also the centra of the mandibular condylar process.

Head Segmentation and Edge Detection
The head is segmented by thresholding and the threshold is determined automatically by the method proposed by Mantas etc. [8]. Then, the Sobel operator is employed to detect the edges of the head. Here (Hounsfield Unit) HU values of different tissues in CBCT images are same as in [6,7]. These detected edges will be used to detect the shape of the mandibular ramus.

Shape Detection of Mandibular Ramus
The shape detection of mandibular ramus is the essential step of the proposed method. The mandibular ramus presents a specific shape in the axial plane as marked in red rectangular boxes shown in figure 1. The typical Hough transform is a feature extraction method to detect simple geometric shapes such as line, circle, ellipse etc. A generalized Hough transform [9] can detect arbitrary shape. A typical shape template containing the mandibular ramus is selected as shown in figure 2. It is used to detect the mandibular ramus from the edges of the condylar process as in the last step. Details of shape detection are shown as follows.

Establish R-table according to the shape template.
It is necessary to arbitrarily select a reference point on the template image such as the image centre as the reference point in our experiment. The R-table describes the distance between the edge point and the reference point and the direction from the edge point to the reference point. Compute the gradient direction ∅ for each edge point and displacement between edge point and reference point: , , is the intensity of the template image at , ; and denote the horizontal and vertical gradients; and denote the horizontal and vertical displacements between the edge point and the reference point. Set ∅ as the index and , as values to establish R-table as shown in Table 1. In addition, an index ∅ may map multiple , . Note that the gradient direction ∅ was expressed as degree rounded down to integer in our experiment. For example, 6.2 and 6.7 are all rounded down to 6 degrees.

Compute a matching score of each slice according to R-table.
This procedure aims to determine the matching slice with the highest matching score to the shape template in the stack of CBCT images. For each edge point p in a slice of the stack, the gradient direction ∅ is computed. If there exists ∅ ∅ in the R-table, the corresponding displacements , can be found in the R-table and then , will be considered as one point to update the reference point. These points generated in this way from all edge points will vote for a new reference point. The point with the most votes is selected as the new reference point and the number of votes stands for the matching score of the slice. The matching score is computed in each slice and the slice with the highest matching score is selected. With such automatically detected slice containing the ramus, region growing process to accurately find the condylar process can be performed and then to find its centre.

Obtain the ramus shape on matched slice.
To determine whether an edge point , is located on the ramus shape in the matching slice, a metric expressed as IsOnRamus is defined as: Here , is the reference point determined from last step. ′, ′ is the point voted by the edge point , . is a threshold to determine if a point is on the ramus (here d = 8 is used). If , = 1, the edge point at , is located on the ramus shape. The minimum bounding rectangle (MBR) of the points on the ramus can be determined. The MBR would be the range for selection of seed points in subsequent steps. If the boundary shape of current slice is much different from that of the template, then the number of the points on ramus would be less but the MBR is basically not changed and then the seed points can be found.

Upward region growing
The datasets which we used were from an iCAT scanner (Hatfield, PA). The range of HU values for bones including trabecular bones with low HU values is roughly from 0 to 3000 from scanned images. Due to the lower HU values of the condylar process, the condylar process is not fully segmented in Section 2.1. In order to accurately segment the condylar process, the matched slice is the first slice for two-dimensional region-growing and slice-by-slice region growing is then performed continuously in each upper slice: HU values of the condylar process are roughly in the range [0, 800] and so and are set as 200 and 0 respectively in our upper region growing criteria for better inclusion; The pixels in the MBR of the matched slice with HU values greater than are set as initial seeds of the first slice and grow the pixels with HU values greater than the low threshold on the first slice; On the immediate upper slice, the pixels with HU values greater than and belonging to the 9-upper neighbors of the initial seeds from immediate lower slice are selected as the initial seeds and the pixels with HU value greater than are included to grow in the immediate upper slice; The process stops when there are no initial seeds in an upper slice. Figure 3 shows the 4-neighbors and 9-upper-neighbors and 9-lowerneighbors. Figure 3. Schematics of 4-neighbors for two-dimensional region-growing and 9-neighbors for threedimensional region-growing: seed point (yellow) and its 4-neighbors (green) and 9-upper-neighbors (red) and 9-lower-neighbors (blue) The upward region-growing algorithm can segment condylar and coracoid processes well in general. But the growing result may contain a part of maxilla when the condylar process is close to the maxilla. For such case, the number of bone pixels will increase abruptly compared to the previous slice due to the local bone structures. When the number of bone pixels in the current slice is more than 30% of the previous slice. The two-dimensional region-growing is carried out again under an additional condition that a pixel is included if one of its 9-lower-neighbors is a grown pixel.

Localization of Condylar Process Centre
The centre of the condylar region in the axial slice with the maximum condyle area of the left/right condylar process is considered as the centre of the left/right condylar process. The gap or distance between the left/right condylar and the coracoid process is computed slice-by-slice starting from the matched slice based on the upward region-growing result as shown in the axial view of figure 4. When the gap is firstly larger than 1 cm (an empirical value) at the lower part of the condylar process, the corresponding slice is assumed as the bottom slice where the condylar process lies in.   5 We calculate the distance between the first pixel and last bone pixel in each row of the left/right condylar process for each slice on or above the bottom slice. The left/right condylar process area on the slice is the sum of all these distances. The slice with the maximum area of left/right condylar process is selected and the average coordinate of all left/right condylar process in this slice is the centre of the left/right condylar process. Note that the pixel between the first and last bone pixels in each row of the left/right condylar process should be regarded as the condylar process pixel for calculating the average coordinate.

Results and Discussion
In order to verify the performance of the proposed algorithm, 20 CBCT datasets were used. These CBCT datasets were acquired using i-CAT system with X-ray tube current 5mA and voltage 120 KV. The voxel size of these datasets is 0.2×0.2×0.2 mm 3 and the resolution of them is 776×776×432. The proposed method was implemented using C++ on aLenovo desktop (Win 10, Intel Core i7-9700K 3.60 GHz, 16 GB RAM). It took 13.8 ± 1.15s to locate the centre of the condylar process for one single CBCT dataset by proposed method. The located centres of the condylar processes were labelled with red point as shown in figure 5. The maxilla was removed to observe the location results clearly from a jaw segmentation procedure [10]. We calculated the spatial distance between the centre located by the proposed method and ground truth centres manually located by two experts to evaluate the accuracy. The average of such spatial distances of the 40 centres (left and right) of 20 CBCT datasets was 0.646 ± 0.824 mm (3.23 ± 4.12 voxels).
The errors were mainly caused by the incomplete segmentation of the condylar process which led to the inaccurate calculation of the condylar process area. As mentioned above, the HU values of the condylar process are in a low level. Therefore, some edges of the condylar process might be missing after two-dimensional region growing.

Conclusion
In this paper, we proposed a method to automatically locate the centre of the condylar process from CBCT images. The proposed method used the generalized Hough Transform [9] to detect the specific shape of the mandibular ramus. The condylar and coracoid processes were segmented by our proposed upward region-growing method. The condylar process centre was determined in a slice where the maximum area of the condylar process was achieved. The proposed method provided accurate location of the condylar process centre and was verified by the manual location. Based on localization of the condylar process centres, we will study the movement model of the TMJ in the future work.