Participants and real-time MRIThis study was approved by the Institutional Ethics Committee (no. 6/7/21) following the Declaration of Helsinki. Written informed consent was obtained from all participants.Subjects were ten adult orthodontic patients (mean age: 31.5 ± 5.9 years, female: n = 5) with full dentition, skeletal class I (mean Wits appraisal: 0.17(± 2.06) mm), and no symptoms of temporomandibular disorder (TMD). Exclusion criteria were age below 18, craniofacial anomaly, large tattoos, and intraoral or intracorporal metal components such as orthodontic braces, pacemakers, cochlear or joint implants due to the MRI constraints.In each participant the right TMJ was scanned using a Siemens Magneton Prisma fit with T2/T1 contrast (refocused FLASH) at ten frames per second according to a previously published protocol20. The in-plane resolution was 0.75 × 0.75 mm for a field of view of 128 × 128 mm. The slice thickness was set at 6.0 mm, echo time (TE) at 1.56 ms, repetition time (TR) at 2.56 ms, and the number of excitations (NEX) at 1. Three slices (6 mm of inter-slices) were positioned on the center of the condyle, localized through a static calibration scan in sagittal, coronal, and transversal planes. The acquisition plane was rotated to include the mandibular ramus. The patients were instructed to do four cycles of mouth opening and closing (intercuspal position – maximal mouth opening (MMO) – intercuspal position), each cycle lasting ten seconds. Standardized instructions on the monitor for opening and closing the mouth were displayed to the patient. Due to the time required by the manual tracking, only the first cycle of each patient was analyzed in this study.Pre-processing for automatic and manual tracking of the mandibular movementThe pre-processing of all data has been done by a the same dentist (F.S.) on a custom Matlab R2019a applicationrunning on an ASUS ROG G751JT with 16Go RAM and an Intel i7 processor. Through this application, a sequence of tasks was applied.First, the contours of the condyle and the fossa articularis were marked by the same experienced user, as illustrated in Fig. 5. To ensure the precision of the position of the condyle’s center, it was not placed by the rater but defined automatically as the center of the circle fitting the contour of the condyle’s head. This circle was computed using LMS algorithm(Pratt, 1987).Fig. 5Screenshot of the GUI once the rater is finished. Manual landmarks are not displayed for visibility. They share positions 1 and 3 with automatic tracking. The areas are drawn by the rater and represent the pixels used for the registration. The contours are also drawn by the user but are only used to determine the center of the condyle. Due to the inclination of the acquisition plane, the Processus coronoideus mandibulae is not visible.Then, from the initial 100 frames, seven keyframes were selected from each MRI scan as follows:
first and last frame of each scan.
first and last frame of each movement cycle.
every 25th frame in between.
On these keyframes, the operator placed two landmarks on the temporal bone (1. the most inferior point of the eminence crest and 2. the most superior point of the fossa articularis) and three landmarks on the mandible (1. the most superior point of the condyle, 2. The most inferior point of the incisura mandibulae and 3. the gonion).Method I—Automatic tracking of the movementFor an automatic tracking of the mandible movement in the rtMRI data a three steps procedure has to be performed based on a mathematically algorithm: (1) the rough transformation based on landmarks, (2) the pixel extraction, and (3) the LMS registration.Rough transformation based on landmarksFirst, the rough transformation was extracted from the manually placed landmarks at each of the frames using SVD algorithm. At each frame, the position of the landmarks was interpolated from the two closest keyframes. The draft transformation described the transformation of the landmarks from their position in the first frame to their position in the current frame using the formula:$$\:{Y}_{rough}={T}_{rough}+{R}_{rough}*Coor{d}_{2D}$$With Coord2D the coordinates of the pixel in the scan, Yrough its rough transformed position, and Trough and Rrough the rough translation and rough rotation respectively.Extraction of the coordinates and tangent of the pixel belonging to the mandibulaIn order to filter the pixels of interest for the registration, the rough transformation (Trough, Rrough) is applied to the mandibula at each frame to act as a mask and filter out every pixel that is not covered. The greyvalues of the masked area are then normalized to achieve two goals:First, the grevalues are normalized to be in the same range of value as the spatial dimensions:$$\:Scal{e}_{FirstFrame}=\frac{interdecile\left(dot\left(\overrightarrow{tangente},\overrightarrow{Coor{d}_{2D}}\right)\right)}{8*interdecile\left(G{V}_{FirstFrame}\right)}$$With ScaleFirstFrame the normalization coefficient, GVFirstFrame the grey values of the pixels of the mandibula area in the first frame, \(\:\overrightarrow{tangente}\) the last component of the principal component analysis applied to the pixels’ coordinates, and \(\:interdecile\) the difference between the first and the last deciles (from 10 to 90%).The second goal of the normalization is to keep the grey-values of every frame in the same range. The other frames are therefore normalized using the formula:$$\:{3}^{rd}Coor{d}_{3D}=\left(GV-D{1}_{currFrame}\right)*\frac{interdecile\left(G{V}_{FirstFrame}\right)}{interdecile\left(G{V}_{currFrame}\right)}+D{1}_{FirstFrame}$$With 3rd Coord3D the third coordinate of the 3D representation, GV the grey value of the pixel and D1FirstFrame the first deciles of the grey values of the first and current frames respectively. The result is displayed by Figs. 6 and 7 Once the 3D representation of the mandibula is computed, the tangent of each of the voxels is computed as being the last component of the principal component analysis applied to the voxel and its neighbors. Finally, pixels whose tangent is too close to the vertical (< 8°) are removed to reduce the needed computational power.Fig. 6Grey-values of the pixels within the area drawn by the rater for the mandibula for the current frame in black and for the first frame in red. Horizontal pixels (with vertical tangents) are discarded to reduce computing power. [Matlab, Version R2019a, https://www.mathworks.com/].Fig. 73D shape extracted from the grey-values for the current frame in black and for the first frame in red This 3D shape of the first frame will be used as a model for the automatic tracking to compute the optimal superposition. [Matlab, Version R2019a, https://www.mathworks.com/].2D LMS-registration based on the 3D representationOnce the 3D representations of both slices are computed for every frame, the LMS registration is applied. The algorithm is designed to minimize the distance between the transformed model (mandible observed at the first frame) and the mandible observed on the current frame. This distance (the error of superimposition) is computed as the mean of the distance between each pixel of the model and its projection, following its tangent, on the observed frame. This is achieved through a repeated cycle of three parts: (1) the points pairing, (2) the rigid body transformation, and (3) the scaling. During this process, every pixel of the model does not have the same weight. In order to focus the precision of the superposition on the area of interest, the pixels of the condyle have their weight doubled, and pixels from the middle slice, centered on the TMJ, have a supplementary weight increase of 150%.At each cycle, each pixel of the first frame is paired with the closest point from the observed scan. The average distance between each pair is called the error of superimposition Then, the rigid body transformation aims to find the translation TLMS and the rotation RLMS using an LMS approach to minimize the error of superimposition. Eventually, scaling aims to find the translation TLMS and the homothety HLMS, minimizing the mean of the squared error. Since this approach searches for a local optimum, the initial transformation is given by the rough translation previously computed. Eventually, the final 2D transformation (TLMS, RRMS, HRMS) can be represented as$$\:{Y}_{LMS}={T}_{LMS}+{R}_{LMS}*{H}_{LMS}*X \quad {\text{with}}\quad {H}_{LMS}=\:\left(\begin{array}{ll}{H}_{11}&\:0\\\:0& \:{H}_{22}\end{array}\right) \;\;\; {\text{and }} \;\; {{\text{R}}_{{\text{LMS}}}} \;\;\;\; {\text{the rotation matrix}}$$Method II—Manual trackingTo compare our automatic tracking method to the actual state of the art, we implemented manual tracking according to Krohn et al. (2020). The same operator placed two landmarks, condylar superior and gonion, at the first frame and each frame of the movement during the first movement cycle, resulting in 200 landmarks placed on the 100 frames of the first opening – closing cycle for each participant. For each frame, the translation was computed based on condyle superior, while the rotation was adjusted using gonion.AnalysisFor each subject and for both methods, three metrics were recorded and compared: (1) the time taken by the rater to track the mandible, (2) the error of superimposition between the two methods, and (3) the distance between the pathways of the superior condyle, the condyle’s center, the gonion and the instantaneous center of rotation (ICR) (see6 in the most distant third of the movement as well as the angle between the rotation of the mandibula given by both tracking methods.The normal distribution of the data has been assessed using the Shapiro-Wilk Test. With the exception of the manual tracking method’s superimposition error, every group we compared did not follow a normal distribution (p < 0.05). Therefore, the Wilcoxon sign rank test has been used, with p < 0.05, to assess the difference between the durations and the superimposition error and if the pathways and the inclination were more distant than their thresholds (right tail). For all metrics, median and the interquartile range were reported as durations in seconds, distances in millimeters, and angles in degrees. The error of superposition is a constructed value based on distances and greyvalue and, therefore, has no dimension. The statistical assessments and the figures were computed using Matlab R2019a.