Karyna Isaieva
email: karyna.isaieva@univ-lorraine.fr

Justine Leclère

Jacques Felblinger

Romain Gillet

Xavier Dubernard

Pierre-André Vuissoz

Methodology for quantitative

Keywords: temporomandibular joint, mandible, condyles, articulation, segmentation, real-time MRI

INTRODUCTION

The temporomandibular joint (TMJ) is a synovial joint composed of an articular disc, mostly fibrous, interposed between two bone surfaces (mandibular condyle and glenoid fossa). TMJ kinematics is complex [START_REF] Woodford | Measurement of normal and pathological mandibular and temporomandibular joint kinematics: A systematic review[END_REF] because it involves numerous muscular structures interacting with the bone and cartilaginous tissues and has six degrees of freedom. The progressive deterioration of the articular surfaces of the TMJ or the muscle imbalance leads to the appearance of temporomandibular dysfunction (TMD) [START_REF] Schiffman | Diagnostic Criteria for Temporomandibular Disorders (DC/TMD) for Clinical and Research Applications: Recommendations of the International RDC/TMD Consortium Network* and Orofacial Pain Special Interest Group[END_REF]. These disorders are common and have reported a prevalence of 25% in the adult population [START_REF] Valesan | Prevalence of temporomandibular joint disorders: a systematic review and meta-analysis[END_REF]. TMD could lead to joint pain, clicking or grinding sounds when performing mandibular movements. In the long term and in the absence of appropriate treatment, these disorders can impact the quality of life of patients by limiting the amplitude of mandibular movements, altering masticatory muscle function, prematurely wearing down the articular disc, or even generating osteoarthritis. So, dynamical visualization of the temporomandibular joint is important [START_REF] Bennett | A Contribution to the Study of the Movements of the Mandible[END_REF] to improve the knowledge about the specific movements of its structures [START_REF] Orthlieb | Cinématique mandibulaire[END_REF].

Multiple approaches exist to study TMJ kinematics [START_REF] Woodford | Measurement of normal and pathological mandibular and temporomandibular joint kinematics: A systematic review[END_REF]. They include commonly used mechanical axiographs [START_REF] Kucukkeles | Compatibility of mechanical and computerized axiographs: A pilot study[END_REF], electromagnetic articulography [START_REF] Lezcano | Symmetry of mandibular movements: A 3D electromagnetic articulography technique applied on asymptomatic participants[END_REF], and various dedicated devices [START_REF] Santosa | The influence of the leaf gauge and anterior jig on jaw muscle electromyography and condylar head displacement: a pilot study[END_REF][START_REF] Agbaje | Assessment of occlusion with the T-Scan system in patients undergoing orthognathic surgery[END_REF]. These devices are affordable and could be directly accessible in the clinician's office. However, the principal drawback of these techniques is that extra dermic and dental reference points are used to record the global motion of the mandibular bone but do not give any information regarding the real TMJ soft and hard tissue variations [START_REF] Palla | Dynamic stereometry of the temporomandibular joint: Dynamic stereometry of the temporomandibular joint[END_REF], and the acquired data is insufficiently reliable [START_REF] Zimmer | The reliability of axiographic tracings[END_REF]. Furthermore, these devices are cumbersome and uncomfortable. Other devices such as infrared [START_REF] Furtado | A specialized motion capture system for real-time analysis of mandibular movements using infrared cameras[END_REF] and optical tracking systems [START_REF] Pinheiro | Measurement of jaw motion: the proposal of a simple and accurate method[END_REF] are considerably more comfortable but still rely on external tracking markers. Currently, only axiographs are widely used by clinicians; electromagnetic articulography is commercialized but used mostly for research purposes, and other devices are in the experimental prototype stage.

In parallel with these devices, TMJ's anatomical structures can be observed due to radiological acquisitions such as Computed Tomography scan (CT-Scan) [START_REF] Hilgenberg-Sydney | Diagnostic validity of CT to assess degenerative temporomandibular joint disease: a systematic review[END_REF] and Magnetic Resonance Imaging (MRI) [START_REF] Abdalla-Aslan | Diagnostic correlation between clinical protocols and magnetic resonance findings in temporomandibular disorders: A systematic review and meta-analysis[END_REF]. MRI has the advantage to be a non-invasive and non-ionizing technology. However, conventional acquisition protocols imply static positions and thus do not allow assessment of the entire TMJ's movement [START_REF] Stehling | High-Resolution Magnetic Resonance Imaging of the Temporomandibular Joint: Image Quality at 1.5 and 3.0 Tesla in Volunteers[END_REF]. Standard static MR images, generally made in two key positions (closed and open jaw), have been shown to frequently not detect symptomatic cases (e.g. temporomandibular disk displacement [START_REF] Ren | Magnetic resonance imaging of the temporomandibular joint Value of pseudodynamic images[END_REF]). Moreover, keeping the mouth in the open position for several minutes is uncomfortable for patients, but is required to avoid motion artifacts. An attempt at TMJ's motion MR imaging has also been made by pseudo-dynamic MRI that enables imaging of intermediate positions [START_REF] Ren | Magnetic resonance imaging of the temporomandibular joint Value of pseudodynamic images[END_REF]. However, it was shown that acquired positions were different from those of natural jaw opening [START_REF] Seemann | Pseudodynamic MRI differs from natural opening of the temporomandibular joint[END_REF].

Real-time MRI is an emerging technology that is largely used for speech studies [START_REF] Isaieva | Multimodal dataset of real-time 2D and static 3D MRI of healthy French speakers[END_REF] and numerous medical applications [START_REF] Zhang | Real-time magnetic resonance imaging of normal swallowing[END_REF][START_REF] Dreha-Kulaczewski | Inspiration Is the Major Regulator of Human CSF Flow[END_REF][START_REF] Nayak | Real-Time Magnetic Resonance Imaging[END_REF]. Real-time MRI allows imaging of all intermediate dynamic jaw positions during natural movements. It could replace existing methods for TMJ's investigation by combining advantages of the existing methods: accessibility of motion information (similarly to kinematicsmeasuring devices), and precise position and general anatomical structure (similarly to conventional MRI).

The quality and the efficacity of real-time MRI in oblique sagittal planes were already assessed on cohorts of healthy subjects [START_REF] Krohn | Real-time MRI of the temporomandibular joint at 15 frames per second-A feasibility study[END_REF] and patients [START_REF] Krohn | Diagnosis of disk displacement using real-time MRI: Clinical report of two patients[END_REF]. It was shown that good-quality quasi-simultaneous multislice imaging was also possible [START_REF] Krohn | Multi-slice real-time MRI of temporomandibular joint dynamics[END_REF]. These pioneers' works offered a real-time MRI of TMJ feasibility evidence and a demonstration of the advantages of the real-time technique over conventional static imaging. They also made an exhaustive qualitative description of condylar motion and introduced the possibility of biomechanical analysis of the TMJ's motion [START_REF] Krohn | Biomechanical analysis of temporomandibular joint dynamics based on real-time magnetic resonance imaging[END_REF]. However, they mostly reported on the characterization of the condyles' motion based on visual observation or manual measurements.

The purpose of the current work was to present an acquisition protocol and a processing pipeline enabling automatizable quantitative evaluation of condylar motion from real-time MRI in the axial plane (which allows strictly simultaneous imaging of both condyles). Oppositely to existing works, we are focusing not on qualitative motion description, but on designing quantitative parameters aiming to provide a clinical physician with quick and comprehensible information describing mandibula's motion symmetricity. Such complement to acquired real-time MRI films would help with patient follow-up and estimation of rehabilitation efficacity based on objective criteria and without the employment of additional motion-tracking devices. This work also reports on the quality evaluation of each step of the proposed pipeline: actual slice orientation for the cohort of volunteers, segmentation, mid-sagittal axis determination, and robustness of the selected parameters. And, finally, we demonstrate the motion parameters variability and their correlation with symmetricity scores based on visual observation.

MATERIALS AND METHODS

Participants

The participants were 1+10 (S0 and S1-S10) healthy volunteers with no known temporomandibular disorders or orofacial pain. The subjects S1-S10 were involved in the inter-subject variability investigation, and S0 was involved in the estimation of the dependence of the proposed parameters on the slice position. The age of S1-S10 varied from 25 to 62 (with a mean of 42 years) and the age of S0 was 52 years. All participants provided written informed consent. The data was acquired under the approved ethical protocol "METHODO" (ClinicalTrials.gov Identifier: NCT02887053). The study was approved by the institutional ethics review board (CPP EST-III, 08.10.01).

MR images acquisition and visual evaluation

The MR images were acquired on a Siemens Prisma 3T. The imaging sequence was a radial radiofrequency (RF)-spoiled fast low angle shot (FLASH) sequence [START_REF] Uecker | Real-time MRI at a resolution of 20 ms[END_REF] with TR=2.34 ms, TE=1.47 ms, the flip angle was 5°, pixel bandwidth was 1670 Hz, and slice thickness was 8 mm. Image size was 136×136 pixels, and the in-plane resolution was 1.41 mm (field of view 192x192 mm). The number of radial spokes was 9 and two slices were imaged quasi-simultaneously for all sequences, leading to a temporal resolution of 42 ms (24 frames per second). The participants were asked to repetitively open and close the mouth with the maximal possible amplitude during the sequence noise. A reference opening/closing (with a period of 3+3 seconds) video was shown to the volunteers before the acquisition start, however, the actual movement speed varied depending on the volunteers. Total acquisition time varied from 12 seconds to 55 seconds covering at least 2 full opening/closing cycles. The acquisition protocol was different for subjects S1-S10 and for S0.

Note that the images were acquired in anatomical position so that the left condyle on the image represents the right condyle of the subject.

Inter-subject variability analysis

For S1-S10, the following protocol was applied. Firstly, oblique sagittal slices going through each of the condyles were imaged quasi-simultaneously. Then, image series for two parallel coronal slices were acquired. Finally, the following instructions were given to the MRI operator for two parallel axial slices' positioning:

1) Parallel to the horizontal Frankfurt plane (localized in sagittal images) which links the lowermost point of the orbital floor to the uppermost point of the external auditory canal.

2) Parallel to the uppermost point of each condyle, pointed in coronal plane images.

3) Parallel to the most posterior point of each condyle pointed in axial plane images.

The last two verifications ensured the symmetricity of our measurements (by taking into account the subject's asymmetry). The slices were placed at the bottom of the temporal tubercule (eminence) to see the condyle's displacement all along its path from its posterior position in the glenoid fossa to its most anterior position ahead of the temporal eminence. The axial planes were imaged for 30-60 seconds. Only the best one of the two (where the condyles were well distinguishable from the surrounding tissues) was used for further processing.

The dynamic MRI series in the axial plane were visually evaluated independently by two clinicians: Expert 1 (JL, dental surgeon with 7 years of experience in TMD treatment) and Expert 2 (RG, radiologist with 8 years of experience in head and neck radiology), and a symmetricity score was given for each subject for the opening and closing movement separately. In order to evaluate the intra-observer agreement, Expert 1 scored the symmetricity twice, with a several-month interval between the scorings.

The most posterior point of each condyle has been marked, in the closed position (condyle in the most posterior position). A line has been made between these two points to determine a "bi-condylar plane" taking into account the natural asymmetry of volunteers. Then, the clinician eye-followed opening and closing movements making particular attention to the distance of each condyle in comparison to the traced line:

 Grade 0: corresponds to an absence of visual difference between right and left points compared to the coronal line.

 Grade 1: corresponds to a relative asymmetry, i.e. not present on each movement or slight asymmetry which can be identified visually only after slowing down the video sequence  Grade 2: corresponds to a visible difference between right or left distance regarding the coronal line, i.e. one side distance compared to the line was most important than the other one.

The experts were blinded to the subject number and to the information extracted from the postprocessing.

Intra-subject variability analysis

For S0, the sagittal and the coronal slices were acquired in the same manner, and the acquisition of the axial slices was repeated several times in order to evaluate slice positioning impact. The MRI operator was instructed to turn the slice in the sagittal plane (using the oblique sagittal image of the left condyle for positioning) and in the coronal plane, as well as translate the plane. The rotation in the sagittal plane was denoted by angle α and was 0 for the strict axial plane and positive for clockwise rotation. The rotation in the coronal plane was denoted by angle β and was 0 for the perfectly symmetric anatomical axial plane defined retrospectively, and positive for the clockwise rotation. The distance between the slice and the reference point (lying in the referent symmetrized axial plane intersecting the uppermost point of each condyle) was denoted as d. The parameters definition is illustrated in Figures 1A and1B, and their precise values were calculated retrospectively from DICOM-files headers. For some of the slices, the condyles were very poorly visible, and they were thus excluded. The following slices were kept for further analysis ( For the analysis, the following series were used.

Automatic data processing 2.3.1 Segmentation

For training and evaluation of the segmentation, we used a single axial slice of subjects S1-S10.

Considering a relatively small number of subjects, we have chosen an intra-subject training strategy. 20 images were manually segmented for each subject, and a U-Net-based neural network was then trained for each subject separately in order to automatically segment the remaining 280-1300 images.

Additionaly, a separate test was performed to test the generalizability of the proposed segmentation approach. The model was trained on the S1-S10 subject set and applied to the subject S0.

The mean square difference (MSD) between the projections calculated from the automatically and the manually segmented condyles has been chosen as the quality metrics. The segmentation has been considered reliable if MSD is inferior to half-pixel size. Additionally, mean displacements of centers of mass of automatically segmented condyles (with respect to those of the manually segmented ones), were calculated. The details on the segmentation method are provided in Appendix.

Mass center projections

The mass center of each condyle was projected on the mid-sagittal axis. This axis was supposed to be perpendicular to the line connecting the condyles' centers and intersecting their common mass center in the closed jaw position. To decrease the impact of position variation and segmentation imperfections on midsagittal axis determination, instead of taking one image per subject to calculate mean condyles centers in closed jaw position, we took the 10 first images. The automatically determined mid-sagittal axes were visually checked by Expert 1 and corrected if needed.

To remove the outliers appearing in the case of erroneous automatic segmentation and manifesting in form of narrow steep peaks, the projection curves were filtered with a median filter (with a width of 5).

Subsequent post-processing included low-pass filtering (with a 1-st order Butterworth filter) which helped to avoid oscillations related to discrete grid and imperfect segmentation. Cutoff frequency was calculated automatically for each subject assuming the MSD between filtered and input data should not exceed 1% of maximal amplitude (in the interval from 0.1 to 10 Hz with a step of 0.1 Hz). Velocity was calculated as a 5-point derivative of the projection. Elementwise difference between projections of condyles centers was also calculated (displacements difference).

All cycle phases were found automatically. Opening and closing intervals were defined as those where velocity was greater than 10% of the velocity maximum. In order to keep only one opening and one closing interval per cycle, all short intervals (less than 300 ms) were removed. The remaining intervals were concatenated if they were closer than 1 sec and had the same velocity direction. All necessary additional checks were performed to exclude the misclassification.

Calculation of the motion parameters

The motion parameters were calculated separately for each cycle and can be divided into three groups:

1. The first group included the motion amplitude (the difference between condyles' open position and close position projections) for the left and the right condyle.

2. The second group contained the mean and standard deviation of the velocity projections, calculated independently for opening and closing.

3. Finally, the third group aimed to describe motion symmetricity. The parameters were normalized to not integrate the subject's morphology and motion speed.

 Spatial asymmetry was represented as the maximal displacement difference (MDD) between the condyles during opening or closing movements. This parameter was normalized by maximal displacement from the whole image series. Additionally, temporal positions of each displacement difference maximum were expressed as a number that linearly varies from 0 (posterior position) to 1 (anterior position).

 Temporal asymmetry parameters were: latency (L) defined as the delay between onsets of two events (opening or closing), and velocity peak delay (VPD) (as illustrated in Figure 1C). The parameters were normalized on opening/closing duration.

To finally combine all three parameters describing the motion symmetricity (MDD, L, and VPD), we applied principal component analysis (PCA), where the variables were absolute values of normalized symmetricity-describing parameters and the observations were the values averaged across all the repetitions of each subject (concatenated for every opening and closing). Only the first principal component was taken as the unique parameter describing motion symmetricity.

To analyze the impact of slice position on the method reliability, we selected three parameters: a geometrical parameter (motion amplitude), a temporal motion asymmetry parameter (VPD), and a spatial motion asymmetry parameter (MDD). Velocity dependence on slice position has not been analyzed due to its significant intra-variability during repetitions. For the parameter α variability, only one of the condyles (the left one, which was chosen for the slice positioning) was taken for the analysis.

Statistical analysis

We performed a weighted Cohen kappa test to determine if the experts' scores are consistent and interpreted them as it is suggested in [START_REF] Landis | The Measurement of Observer Agreement for Categorical Data[END_REF]. The first score set of Expert 1 was taken for comparison with Expert 2 and for further analysis.

To evaluate the correlation between the symmetricity parameters (SP) and the slice position parameters, the Spearman correlation coefficient (denoted as r) was calculated. The correlation was considered negligible for r<0.1, weak for 0.1<r<0.39, moderate for 0.4<r<0.69, strong for 0.7<r<0.89, and very strong for r>0.9 [START_REF] Schober | Correlation Coefficients: Appropriate Use and Interpretation[END_REF]. We considered the quantities not significantly correlated if p > 0.2 (beta level) [START_REF] Mascha | Significance, Errors, Power, and Sample Size: The Blocking and Tackling of Statistics[END_REF]. The same test was used to find the pairwise strength of correlation between the SP.

To find the correlation between the experts' scores and the SP, the scores were averaged between two experts, and the Spearman correlation coefficient (r) was calculated. We considered the quantities significantly correlated if p<0.05 (alpha level) [START_REF] Mascha | Significance, Errors, Power, and Sample Size: The Blocking and Tackling of Statistics[END_REF].

RESULTS

Automatic data processing 3.1.1 Semi-automatic segmentation time efficiency and reliability

With the chosen semi-automatic segmentation strategy, the processing took 12 minutes per volunteer on average for all stages (from DICOM load and delineation to prediction and parameter calculation) and 4 minutes for delineation. Mean MSD was 0.51 mm (with a pixel size of 1.41 mm) and thus can be considered reliable. The mean displacement was +0.26 mm (the segmentation masks are, on average, higher than the manually created ones). From the generalizability test, mean MSD value was 0.55 mm and mean displacement was +0.42 mm which are both inferior to the half of the pixel size. The detailed results on the segmentation are presented in Appendix.

Mass center projections

Mid-sagittal axis determined automatically was estimated as acceptable for 8 subjects and corrected for two subjects: S7 and S9. The resulting motion curves and their difference were visually analyzed and compared with the generated videos (Video 1-10). Examples of motion curves are given in Figure 2: the majority of the volunteers show almost the same right and left displacement difference (Figure 2A). However, some volunteers demonstrated an asymmetry during the motion (for example, as can be seen in Video 1). It manifests in form of a steep peak on the displacement difference curve (black curve) (e.g. in Figure 2B). Such peaks at the movement beginning are related to the latency between the right and left condyles to initiate the displacement. They differ from the peaks of most of the volunteers by their shape, amplitude, and temporal position (for comparison, see Video 2-8). Steep peaks can be distinguished from small peaks related to minor segmentation errors, and from wide peaks that can be observed at the ending, which can be explained by an incomplete performing opening or closing (e.g. first repetitions of S3, see Video 3).

Similarly, in the case of asynchronous motion, latency, and velocity peak delay were observed on anterior-posterior velocity projections (e.g. right condyle delay (blue curve) for S1 in Figure 2D). However, a velocity peak delay could also be observed in the case of almost perfectly symmetrical motion (Figure 2C, the third repetition of S2, and Video 2).

Intra-subject variability analysis

All processed videos from real-time MRI of the subject S0 acquired in different planes are provided in Analysis of the dependence of the three selected parameters (motion amplitude, VPD, and MDD) on distance from the reference point has shown that motion amplitude is smaller for the slices placed lower (Figure 3A) with r right = -0.94 (p<0.001) and r left = -0.97 (p<0.001). This was an expected result considering the mandibular bone is a rigid structure whose rotation center is located above the imaging plane.

The MDD demonstrated a not significant correlation (r open = 0.28, p open =0.27; r clos =-0.31, p clos =0.30) to the slice displacement for the slices whose center was placed from 4 to 14 mm down from the reference point (Figure 3C). However, the slice displaced on more than 14 mm has shown completely different results. Indeed, the curves and the video analysis (Video 18) have shown a completely different motion trajectory of this part of the mandible.

Figure 3D shows that the motion amplitude was greater for positive values of α, and very small for a large negative value of α (around 3.5 mm for -39° compared to 6.5 mm for 23°) which decreases the method sensitivity. Furthermore, it has been noted that for too large positive values of α, another issue (partial volume effects related to an extended visualization of the mandibular bone parts in the slice) manifests (Video 15). This phenomenon could complicate the automatic segmentation task to differentiate the condyle from the ramus and crucially affect the precision of the mass center determination. Thus, in the axial plane, we defined anatomical markers which have to be seen to ensure an optimal placement while slice positioning (parallel to the Horizontal Frankfurt Plane (HFP)): the nasal septum, mandibular condyles (with an oval shape), the external auditory canals and the clivus between sphenoid and occipital bones (appearing as a horseshoe shape). Other validations can be done in other orientations: in the coronal plane, axial slices have to be parallel to the line connecting the uppermost point of each condyle and in the sagittal plane, axial slices have to be parallel to the HFP connecting the lowermost point of the orbital with the uppermost point of the internal auditory canal. As expected, the variations of β had a strong impact on the motion parameters. According to our previous explanations, in Figure 3G, it can be seen that in the case of greater β, a decrease in the motion amplitude of the left condyle (r=-0.94, p < 0.001) and an increase for the right one (r=0.94, p < 0.001) can be observed. VPD was not significantly correlated (r open =-0.33, p open =0.30; r clos = -0.34, p clos =0.31) to β changes (Figure 3H). Finally, MDD was strongly correlated with β (r opening =0.94, r closing =0.94, p < 0.001 for both). It means that slice placement greatly affects this last parameter and produced measurement errors in the case of asymmetric positioning.

In order to estimate possible error values of the MDD, actual β angles of the axial slices selected by MRI operators have been retrospectively analyzed for S1-S10. This analysis has shown that 8 of 10 axial slices were almost symmetrically positioned with a β between -0.68° and 0.77° compared to the retrospectively determined ground truth position, and one slice was inclined at 1.35° (S10). One more slice was inclined at 2.87° (S5) and was removed from further analysis. For this last slice, we noticed that the slice was strictly axial placed instead of anatomically axial placed due to a human mistake during the performed exam.

The β values for axial slices for S1-S10 are depicted as vertical lines in Figure 3F and the accepted values are denoted as a green area in Figures 3F-4I. Considering the slope of the MDD (on β) was of the same order for all volunteers, one can estimate the maximal possible error as ±0.5.

Inter-subject variability analysis 3.3.1 Visual symmetricity evaluation

The scores are presented in Table 2. The Weighted Cohen kappa value between two different score sets of the Expert 1 was κ=0.80, demonstrating a strong intra-observer agreement. The kappa-value between Expert 1 and Expert 2 was κ=0.58, therefore the scores of the two experts thus demonstrated a moderate agreement. One can note that the intra-observer agreement was much stronger than the inter-observer agreement, and Expert 1 (a dental surgeon) tended to choose higher scores than Expert 2 (a radiologist).

Subject Opening Closing Expert 1 Expert 2 Expert 1 bis Expert 1 Expert 2 Expert 1 bis S1 2 2 2 1 1 1 S2 0 0 0 0 0 0 S3 1 0 1 0 0 0 S4 0 0 1 0 0 0 S5 - - - - - - S6 2 1 2 1 1 1 S7 1 0 1 0 0 0 S8 1 0 0 0 0 0 S9 2 2 2 1 0 1 S10 2 2 2 1 1 2

Inter-subject variability of the motion parameters

The displacement between anterior and posterior plateaus (Figure 4A) varied from 8 mm (S10) to 17 mm (S1, S2, S7) and depended both on the subjects' anatomy and slice position. A large difference between right and left condyle posterior plateaus was observed mostly in the case of the incomplete jaw closing from one of the sides (S2, S9). Such difference for anterior plateaus had greater amplitude and was related to both incomplete opening from one side (S2) and global motion asymmetricity (S1). 4B and4C that the mean velocity largely varied between subjects and between repetitions for the same subject (even though a pre-recorded video showing the desirable motion rate was shown to the subjects at the beginning of each sequence). While some subjects (S4, S10) showed stable opening-closing rhythm, others (S3) demonstrated high dispersion of mean velocity values.

It has been observed in Figures

The motion symmetricity parameters for the subjects S1-S10 are depicted in Figures 4D-4F. Despite all these parameters characterizing the motion symmetricity, only some of them were correlated with each other:

-The two parameters describing temporal displacement (L and the VPD) were not significantly correlated (r=0.21 p = 0.41).

-L and the MDD have demonstrated a weak correlation (r =-0.39, p=0.09).

-A significant correlation was observed between VPD and MDD (r=-0.64, p = 0.004).

DISCUSSION

This work presents an acquisition and processing protocol that allows reliable quantification of directly visualized condylar motion symmetricity. The precision of the slice placement and its impact on the measurements have been assessed. The first estimation of inter-subject variability in a healthy population is established. The proposed parameters could serve as markers helping to follow up with patients, in particular, to quantitatively evaluate the efficacity of the therapy.

Despite the axial plane orientation reducing the volume of available motion information compared to the oblique sagittal slices [START_REF] Krohn | Real-time MRI of the temporomandibular joint at 15 frames per second-A feasibility study[END_REF][START_REF] Krohn | Diagnosis of disk displacement using real-time MRI: Clinical report of two patients[END_REF][START_REF] Krohn | Multi-slice real-time MRI of temporomandibular joint dynamics[END_REF][START_REF] Krohn | Biomechanical analysis of temporomandibular joint dynamics based on real-time magnetic resonance imaging[END_REF], it represents a simpler way to detect motion symmetry. In this orientation, the condyles can be images simultaneously as the mandibular bones' crosssections (contrary to the oblique sagittal orientation where they represent the extremities of mandibular bones). Therefore, the definition of a condyle center is more straightforward and objective which allows processing automation.

We have shown that the proposed symmetricity parameters can be extracted in an semi-automated manner, by manually segmenting only 20 images (from several hundreds). It was also demonstrated that this approach is potentially generalizable in terms of inter-subject segmentation prediction. It can be noted, that the automatic masks were usually higher than the manually created ones. This phenomenon is more pronounced in the case of large negative values of α and attributed to the partial volume effect caused by the presence of mandibular ramus.

In an example of one subject, we have shown that the temporal symmetricity parameters (latency, velocity peak delay) do not highly depend on the slice position. We showed that slice positioning allows a suitable tolerable margin of placement error regarding the anteroposterior inclination α. However, large values of right-left inclination β may lead to false asymmetry occurrence. Maximal displacement difference was shown to significantly vary within a 1-degree interval which was in practice respected by the involved MRI technicians. The variation extracted from the linear fit for intra-subject data had the same order as the inter-subject variation. Thus, the results indicate that maximal displacement difference is the least reliable from the proposed symmetricity-describing parameters.

The visual scores were highly reproducible by the same reader but in a moderate agreement between the physicians. This imperfect association could probably be explained by the nature of the professional activity of the experts. Nevertheless, the average scores were significantly correlated with the automatically calculated parameters. Latency demonstrated a better overall correlation with the visual scores than other automatically calculated parameters. The high correlation between latency and the introduced parameters could probably be explained by peculiarities of human perception and does not necessarily mean that this parameter is the most reliable for the objective motion symmetricity assessment.

Limitations

One of the substantial limitations of this study is the relatively small number of volunteers. To demonstrate the segmentation quality for different morphologies and slice placements, the neural network was re-trained for each new subject (intra-subject prediction) which is not a generalized approach. Only one subject was involved to inter-subject segmentation prediction. The results are encouraging, however, the sample size in our case was not sufficient to consider these results highly reliable. Currently, the segmentation is relatively time-consuming and the manual part takes 4 minutes on average which is hardly acceptable in a real clinical routine. Nevertheless, in the case of a larger cohort (at least several dozen of subjects), one could envision a completely generalized approach with an appropriate generalizability validation, as it was successfully done for datasets of MR images of different organs [START_REF] Radke | Deep Learning-Based Post-Processing of Real-Time MRI to Assess and Quantify Dynamic Wrist Movement in Health and Disease[END_REF][START_REF] Fu | A novel MRI segmentation method using CNN-based correction network for MRI-guided adaptive radiotherapy[END_REF][START_REF] Wu | Left ventricle automatic segmentation in cardiac MRI using a combined CNN and U-net approach[END_REF][START_REF] Kijowski | Deep Learning for Lesion Detection, Progression, and Prediction of Musculoskeletal Disease[END_REF]. The processing in this case would be also much faster, since the manual segmentation and the network training were the most time-consuming parts of the processing.

Similarly, only one subject was involved in the validation of the method's robustness with respect to different plane orientations. Nevertheless, these preliminary results allow a highly reliable conclusion (due to a substantial difference in correlation coefficients) that distance-based parameters (such as MDD) are more sensitive to slice orientation than time-based parameters (such as latency and velocity peak delay).

Another limitation of this study is that it was focused on the evaluation of motion symmetricity of asymptomatic volunteers with no known pathologies. The proposed parameters were linked to the symmetricity, but not to the clinically important signs (such as the intensity of orofacial pain, or the presence of clicking sounds). Therefore, a larger population (including symptomatic patients), and evaluation of clinically valuable parameters will be added to the future data collection protocol.

Analogously, the mid-sagittal axis determination based on close jaw position was shown to be reliable for most healthy subjects that participated in the study. However, this approach might fail for patients presenting pathological jaw anatomy. Thus, in the future, the mid-sagittal axis determination may need a more robust definition and could include extraoral anatomical markers.

CONCLUSION

An imaging protocol and a data processing pipeline were proposed for the investigation of condylar motion, and three normalized parameters describing motion symmetricity (latency, velocity peak delay, and maximal displacement difference) were introduced. It was shown that this approach potentially enables automatic segmentation and extraction of the introduced motion parameters.

Two of the three proposed parameters (latency and velocity peak delay) did not significantly change with a slight (observed in practice) change in the slice position. However, the maximal displacement difference was found to be dependent on the slice inclination in the coronal plane and varied considerably.

The physicians' symmetricity scores based on the videos acquired in the axial plane were in moderate agreement. All three motion symmetricity parameters, separately, have demonstrated a significant correlation with the average experts' scores, and a PCA-based parameters combination has shown the best correlation. Therefore, the presented results prove the reliability of the method and can serve for future clinical studies.

APPENDIX. SEMI-AUTOMATIC SEGMENTATION

Accurate image segmentation is a key challenge in quantitative analysis, however, recent developments in machine learning have made this task more manageable. Nowadays, the U-Net convolutional neural network [START_REF] Ronneberger | U-Net: Convolutional Networks for Biomedical Image Segmentation[END_REF] has demonstrated its reliability for the segmentation of real-time MRI of various organs:

heart [START_REF] Hennemuth | Automatic analysis of multicycle realtime MRI for the assessment of variable cardiac function based on multi-orientation U-net segmentation[END_REF], wrist bones [START_REF] Radke | Deep Learning-Based Post-Processing of Real-Time MRI to Assess and Quantify Dynamic Wrist Movement in Health and Disease[END_REF], or tongue [START_REF] Isaieva | Automatic Tongue Delineation from MRI Images with a Convolutional Neural Network Approach[END_REF]. The emergence of machine learning frameworks has also created new opportunities to extract motion characteristics from organs of interest. In our work, we employed a semi-automatic approach based on the U-Net which consisted in applying a model trained on a small number of images, to segmentation of the remaining several hundreds of images of the same subject.

Additionally, we tested the possibility to extend this approach to a fully automatic segmentation.

Methods

Intra-subject semi-automatic segmentation

The development set consisted of 20 images per subject: 16 for training and 4 for validation. In order to ensure the best diversity in the development set, a k-means clustering of all images (with k=20) has been performed. The distance metric was the Euclidean norm in the intensity space of vectorized images.

Images having the least distance to each of the 20 cluster centers were chosen as 20 non-similar images.

The condyles on the selected 20 images were manually delineated by KI who was previously trained by Expert 1. Then, a square region of interest (ROI) of 32 pixels was manually selected for each condyle (the same for the whole series) whose center corresponds to the mean position of the center of the delineated condyle. The resulting ROIs of images and masks were resized to 128 pixels using functions from the

OpenCV library [START_REF] Bradski | The OpenCV Library[END_REF] (bicubic interpolation over 4x4 pixel neighborhood in case of image increase and resampling using pixel area relation in case of image decrease). To increase the development set size, we used data augmentation (translations on -4 and +4 pixels in both horizontal and vertical directions).

The resulting images and masks served as input and output for the convolutional neural network. The network was based on the classic U-Net [START_REF] Ronneberger | U-Net: Convolutional Networks for Biomedical Image Segmentation[END_REF] with 4 downsampling and 4 upsampling steps. Each downsample step consisted of two 3x3 convolutions, followed by a rectified linear unit (ReLU) and a 2x2 max pooling operation with stride 2 for downsampling, and the upsampling path was constructed symmetrically. A zero-padding was applied after each convolution in order to conserve the matrix size, so that the image sizes were always the powers of 2, and a dropout regularization was applied on each downsampling or upsampling step. The number of feature channels started from 32, doubled on each downsampling step, and halved on each upsampling step. The final step consisted of the application of a 1x1 convolution and a sigmoid activation function to map 32 features to the probability of belonging to a condyle. The loss function was selected to be binary cross-entropy which usually demonstrates the best performance when the classes are well-balanced [START_REF] Jadon | A survey of loss functions for semantic segmentation[END_REF]. The network was trained separately for each subject with an Adam optimizer with a batch of 8 samples and an automatic early stopping with patience 2. The code was implemented using a GPU version of the Keras framework [35] on a machine with a 28-cores Intel Xeon 2.5 GHz processor, 64 GB RAM, and an NVIDIA Quadro P5000 graphics card. The U-Net model was trained in two steps. First, a model was trained starting from random weights and a fully automatic k-means-based segmentation of the data from subject S1. In the second step, the U-Net model was trained on 20 images specifically for each subject, with the first step as initialization. After the training, the models of each subject were applied for condyles segmentation of the whole series of images.

For the quality test of the automatic segmentation, each 50th image from each axial series of S1-S10 was manually delineated by Expert 1 in ImageJ software [START_REF] Schindelin | Fiji: an open-source platform for biological-image analysis[END_REF]. All intersections with the training and validation sets were excluded from the testing set to avoid bias (but kept on the graphs for better visualization).

The mean square difference (MSD) has been chosen between the projections calculated from the automatically and the manually segmented condyles as the quality metrics considering all further measures are based on projection values. Our segmentation has been considered reliable if MSD is inferior to halfpixel size. Additionally, mean displacements of centers of mass of automatically segmented condyles (with respect to those of the manually segmented ones), were calculated.

Generalizability of the segmentation approach

In order to test the generalizability of the proposed segmentation approach, an inter-subject test was performed. Firstly, the neural network was trained with both manually-delineated condyles of S1 as the validation set, and both manually-delineated condyles of subjects S2-S10 as the training set. Then, a series with the best orientation of S0 (series 4), was selected as the testing set (20 images selected with the kmeans algorithm, both condyles).

In order to test the impact of the region of interest choice, the inter-subject segmentation evaluation explained above has been evaluated with different ROI positions. Firstly, the center of the region of interest was selected as the center of mass of the manually segmented condyles, and its size was fixed to 32 pixels. Then, translations of 4 pixels (approximately 6 mm) in horizontal and vertical directions were applied, resulting in 9 sets of the model test results.

Results

Intra-subject semi-automatic segmentation

Details on the metrics values are presented in Table A

  Fixed α and β, variable d: series 2, 3, 7, 8, 9; Planes intersecting the upper left part of the mandibular bone at different points (by a clockwise rotation), variable α: series 3, 4, 5, 6; Fixed d and α, variable β: series 3, 10, 12.

Figure 1 .

 1 Figure 1. (A) and (B) Definition of the parameters that describe relative slice position. The white slice denotes the reference axial slice, and the yellow slice corresponds to the imaged slice. (C) Definition of latency and peak delay illustrated on the 2-nd opening of S6. The red corresponds to the right condyle and the blue to the left condyle. The dashed vertical lines denote opening onsets defined as the time point when velocity exceeds 10% from its maximal value (denoted by horizontal dotted lines which are almost indistinguishable from each other). The solid vertical lines denote velocity peak positions.

Figure 2 .

 2 Figure 2. (A) and (B): post-processed anterior-posterior projections of condylar centers during the opening-closing motion for subjects S2 (A) and S1 (B). (C) and (D): anterior-posterior velocities of condylar centers for subjects S2 and S1, correspondingly. Opening and closing time intervals calculated automatically are denoted by color codes.

Figure 3 .

 3 Figure 3. (A)-(E), (G)-(I): Condylar motion parameters measured for S0 in various slice orientations. VPD abbreviates "velocity peak delay", and MDD abbreviates "maximal displacement difference". Dashed lines correspond to a linear fit of the points. The green area corresponds to the accepted range of the parameters. Horizontal dotted lines in subplot (I) denote values of normalized maximal displacement difference for extremities of the accepted β zone. (F): Vertical lines depict the actual β of axial slices selected by MRI technicians with respect to the retrospectively selected symmetric slice for S1-S10. Horizontal lines and circles show the correspondence between β and the subject ID. The red circle points to the subject S5 rejected from the further analysis. No significant correlation between α and peak delays was observed (r open =0.23, p open =0.34; r clos =0.08,

- - Figure 4 .

 4 Parameters calculated from posterior-anterior motion of condyles centers in the axial plane. Markers on plots (A)-(C) correspond to mean values and error bars correspond to standard deviations. (D)-(F): Symmetricity parameters. VPD abbreviates "velocity peak delay", and MDD abbreviates "maximal displacement difference". The color code on plot (F) denotes the temporal position of the peak, which varies linearly from 0 (posterior position) to 1 (anterior position). The dashed lines on plots (G)-(I) depict the linear fit of the data (red for opening, blue for closing). -Nevertheless, all these parameters were significantly correlated with the experts' scores with a correlation coefficient of r=0.78 (p<0.001) for latency, r=0.58 (p = 0.01) for velocity peak delay and r=0.53 (p = 0.02) for MDD. As can be seen from Figures 4G-4I, the outliers (where the automatically calculated quantities did not correspond to the visual scores), were different for different parameters: for example, the opening of S6 with an average visual score = 1.5, demonstrated MDD = 0.01 (too small for such visual score) and L = 0.14 (consistent with the visual score). Similarly, the opening of S3 (visual score = 0.5) demonstrated L = 0.4 (too high for this visual score) and VPD = 0.11 (consistent with the visual score). The PCA-based combination of the symmetricity parameters based on principal component analysis has shown the strongest correlation with the physicians' scores with an average correlation coefficient of r=0.80 (p<0.001).

  Examples of successful and average segmentations are depicted in Figures A.1A and A.1B. The occurred segmentation errors can be divided into two groups: those related to partial volume effects and complete failures.

Figure A. 1 .

 1 Figure A.1. Examples of (A) good (S1 with d right =0.09, d left =0.15), (B) average (S6 with d right = 0.5 mm, d left = 0.76), (C) poor (S10 with d=1.3 mm for right condyle) and (D) failed (S5 with d=2.2 for left condyle) segmentation of condyles. Red and blue curves denote automatic segmentation respectively of the right and the left condyle, and green and yellow (if present) denote manual segmentation respectively of the right and the left condyle. (E), (F) Red and blue circles show the mass centers for respectively manually delineated right and left condyles, and red and blue curves represent the automatically determined anterior-posterior projections of condylar centers during opening-closing for two subjects (before any filtering); (E) For S5, blue peak (left condyle) represents a complete failure of segmentation. (F) For S6, red and blue curves demonstrate a successful automatic segmentation.

Table 1 )

 1 :

	Series	α	β	d
	number
	1	-1,2	2,72	12
	2	-1,2	1,23	4
	3	-1,2	1,23	12
	4	-12	-2,42	11
	5	23	0,33	0
	6	-39	0,77	24
	7	-1,2	1,23	8
	8	-1,2	1,23	16
	9	-1,2	1,23	14
	10	-1,2	5,83	11
	11	-1,2	-3,57	4
	12	-1,2	-3,57	12

Table 1 .

 1 Positions of the slices used for intra-subject variability analysis.

Table 2 .

 2 Symmetricity scores provided by the experts from visual observation only. Expert 1 is a dental surgeon and Expert 2 is a radiologist.

  .1.Table A.1. Mean square difference and mean displacement values that show the discrepancy between the automatic and the manual determination of a condyle center for intra-subject prediction for each subject.

	Subject	MSD left condyle (mm)	MSD right condyle (mm)	Mean displacement left condyle (mm)	Mean displacement right condyle (mm)
	S1	0.37	0.40	0.01	0.14
	S2	0.35	0.44	0.19	-0.07
	S3	0.64	0.47	0.22	-0.30
	S4	0.54	0.61	0.41	0.37
	S5	0.77	0.35	0.67	0.04
	S6	0.51	0.53	0.45	0.25
	S7	0.66	0.48	0.61	0.38
	S8	0.43	0.47	0.19	-0.03
	S9	0.44	0.53	0.28	0.39
	S10	0.76	0.43	0.66	0.28

Video Content Video 1.avi -Processed real-time MRI in an axial plane acquired during jaw opening and closing of the subject S1. Video 2.avi -Processed real-time MRI in an axial plane acquired during jaw opening and closing of the subject S2. Video 3.avi -Processed real-time MRI in an axial plane acquired during jaw opening and closing of the subject S3. Video 4.avi -Processed real-time MRI in an axial plane acquired during jaw opening and closing of the subject S4. Video 5.avi -Processed real-time MRI in an axial plane acquired during jaw opening and closing of the subject S5. Video 6.avi -Processed real-time MRI in an axial plane acquired during jaw opening and closing of the subject S6. Video 7.avi -Processed real-time MRI in an axial plane acquired during jaw opening and closing of the subject S7. Video 8.avi -Processed real-time MRI in an axial plane acquired during jaw opening and closing of the subject S8. Video 9.avi -Processed real-time MRI in an axial plane acquired during jaw opening and closing of the subject S9. Video 10.avi -Processed real-time MRI in an axial plane acquired during jaw opening and closing of the subject S10. Video 11.avi -Processed real-time MRI in axial plane 1 acquired during jaw opening and closing of the subject S0. Video 12.avi -Processed real-time MRI in axial plane 2 acquired during jaw opening and closing of the subject S0. Video 13.avi -Processed real-time MRI in axial plane 3 acquired during jaw opening and closing of the subject S0. Video 14.avi -Processed real-time MRI in axial plane 4 acquired during jaw opening and closing of the subject S0. Video 15.avi -Processed real-time MRI in axial plane 5 acquired during jaw opening and closing of the subject S0. Video 16.avi -Processed real-time MRI in axial plane 6 acquired during jaw opening and closing of the subject S0. Video 17.avi -Processed real-time MRI in axial plane 7 acquired during jaw opening and closing of the subject S0. Video 18.avi -Processed real-time MRI in axial plane 8 acquired during jaw opening and closing of the subject S0. Video 19.avi -Processed real-time MRI in axial plane 9 acquired during jaw opening and closing of the subject S0. Video 20.avi -Processed real-time MRI in axial plane 10 acquired during jaw opening and closing of the subject S0. Video 21.avi -Processed real-time MRI in axial plane 11 acquired during jaw opening and closing of the subject S0.

FUNDING

This work was supported by Full3DTalkingHead (ANR-20-CE23-0008) and MECAMAX (University of Lorraine, Interdisciplinary call) projects, CPER IT2MP, Région Lorraine and FEDER.

DATA AVAILABILITY STATEMENT

The full processed dataset is available in the video content. The raw dataset generated and analyzed during the current study is not publicly available, but are available from the corresponding author on reasonable request.

Generalizability of the segmentation approach

The results of inter-subject segmentation prediction on S0 for 9 different ROI positions are presented in

ACKNOWLEDGMENT

The authors thank Gabriela Hossu for her advices.

DECLARATION OF COMPETING INTEREST

The authors declare that they have no known competing financial interests or personal relationships that could have appeared to influence the work reported in this paper.