A robust, semi-automated approach for counting cementum increments imaged with X-ray computed tomography

Cementum, the tissue attaching mammal tooth roots to the periodontal ligament, grows appositionally throughout life, displaying a series of circum-annual incremental features. These have been studied for decades as a direct record of chronological lifespan. The majority of previous studies on cementum have used traditional thin-section histological methods to image and analyse increments. However, several caveats have been raised in terms of studying cementum increments in thin-sections. Firstly, the limited number of thin-sections and the two-dimensional perspective they impart provide an incomplete interpretation of cementum structure, and studies often struggle or fail to overcome complications in increment patterns that complicate or inhibit increment counting. Increments have been repeatedly shown to both split and coalesce, creating accessory increments that can bias increment counts. Secondly, identification and counting of cementum increments using human vision is subjective, and it has led to inaccurate readings in several experiments studying individuals of known age. Here, we have attempted to optimise a recently introduced imaging modality for cementum imaging; X-ray propagation-based phase-contrast imaging (PPCI). X-ray PPCI was performed for a sample of rhesus macaque (Macaca mulatta) lower first molars (n=10) from a laboratory population of known age. A new method for semi-automatic incrementcounting was then integrated into a purpose-built software package for studying cementum increments. Comparison with data from conventional cementochronology, based on histological examination of tissue sections, confirmed that X-ray PPCI reliably records cementum increments. Validation of the increment counting algorithm suggests that it is robust and provides accurate estimates of increment counts. In summary, we show that our new increment counting method has the potential to overcome caveats of conventional cementochronology approaches, when used to analyse 3D images provided by X-ray PPCI.

9 variation of six months for m1 replacement necessitated the use of a minimum 200 expected increment count, and a maximum expected increment count (one increment 201 higher than the minimum expected count) (Table 1)  In a preliminary experiment, the cementum tissue of specimen l56 was imaged 211 using X-ray PPCI through SR CT for a range of different experimental settings. Four 212 key experimental settings (X-ray energy, exposure time, number of X-ray projections, 213 and sample-to-detector distance) were individually varied according to Table 2, while 214 all other experimental settings were fixed at an X-ray energy of 20 keV, a voxel size 215 Image quality measures for each experimental setting were calculated for 10 225 CT slices, representing the same regions of the tooth root of l56 in each scan (Fig. 3).
where and represent the standard deviations of cementum and background, 236 respectively. Mean values of SNR and CNR were calculated from the values for these 237 10 slices and compared between all experimental settings ( Fig. 4 and Table 2). 238 Following this preliminary study, an X-ray energy of 20 keV and a sample-to-239 detector distance of 14 mm was chosen. This provided sufficient image contrast 240 between cementum increments. For each scan, the exposure time was set to 150 ms, 241 and 1501 projections were taken per scan. These settings provided sufficient image 242 quality at scan times that allowed the entire sample to be imaged during our fixed-243 time experiment. The phase of each X-ray projection was retrieved through the 244 index =3.7 10 -8 and the decrement of the real part of the refractive index =1.7 10 -247 10 , and hence the ratio   ⁄ = 218, were fixed for all scans. The phase images were 248 reconstructed using an in-house implementation of the Gridrec algorithm (Marone 249 and Stampanoni, 2012) at TOMCAT. The resulting PPCI CT reconstructions were 250 saved as 16-bit tiff stacks. 251

Image processing: straightening and filtering 252
Image processing of the raw data is often needed before digital image visualisation, 253 segmentation and quantification. This can involve a wide range of image processing 254 methods, the majority of which are based on the manipulation of 2D pixels and/or 3D 255 voxels using mathematical operations (Nixon and Aguado, 2012). Here, we applied 256 two principal image processing methods to individual CT slices: straightening and 257 isolation of cementum, and directional filtering of cementum increments. 258 For circumferential structures such as cementum increments, it is often 259 difficult to apply standard image processing tools and analyses without distorting 260 results, due to complexities in their patterns and boundaries. Hence, 2D straightening 261 algorithms are often applied in order to further analyse the data. We chose to use the 262 'Straighten' tool of the open source ImageJ/Fiji image analysis software (Schneider et 263 al., 2012). This tool applies cubic-spline interpolation across a segmented midline of 264 the feature of interest, which is defined by the user. Straightening is then performed 265 using a series of non-linear cubic splines for an arbitrary number of pixels on either 266 side of the midline that can be also be determined by the user. Here, we assigned this 267 number on an individual basis for each dataset, based on the (radial) thickness of 268 cementum being imaged, to ensure that all cementum, but no dentine, was included in 269 the processed image ( Fig. 5.a-b). This workflow was then repeated in a semi-270 automated fashion for all CT slices of each dataset following Newham et al. (2020a), wherein the same midline coordinates and thickness of the segmentation was applied 272 for as many slices as possible from the first slice of a dataset, before adjusting these 273 parameters once they became suboptimal at several points further through the 274 PPSRCT volume due to misalignment between the scanning and root axes. 275 Following straightening, cementum datasets were further processed using 276 directional filtering in order to enhance contrast between increments (Fig. 5.d). Gaussian filter oriented at 90° to the x-axis has been shown to substantially enhance 290 image contrast between straightened cementum increments (see Supplement, section 291 1. 'Image processing by directional filtering') ( Fig. 5). 292

Cementum increment counting algorithm 293
The cementum increment counting algorithm developed here was designed to count 294 cementum increments in a user-independent and semi-automated fashion, further 295 developing the rationale proposed by DCLA for distinguishing individual increments 296 by using a cut-off point that greyscale peaks/troughs must differ by, in order to be 297 counted as a genuine increment. Our algorithm employs population statistics (mean 298 and standard deviation of greyscale values) to count increments, based on the unique 299 distribution of greyscale values within each individual CT slice (see Supplement, 300 section 2 'Robustness testing for increment counting algorithm'). As in DCLA, this 301 method makes use of the average greyscale values along 10 pixel-thick transects through 302 cementum (Fig. 2a). Adapting methods used in tribological surface profiling 303  Therefore, the use of the mean greyscale value and its standard deviation for the entire 309 transect, as opposed to local values for individual sections of the transect, may 310 preclude the counting of genuine increments towards the outermost increment (as 311 greyscale peaks/troughs are below the mean value), and counting of increments 312 towards the cemento-dentine boundary (as greyscale peaks/troughs are above the 313 mean value). The mean and standard deviation of greyscale values in each section is 314 then calculated (Fig. 2b-d), and light-dark increment pairs are distinguished as peak-315 trough systems in greyscale that depart from the mean value beyond the local standard 316 deviation in each section (Fig. 2c). 317 Most importantly, this new method for increment counting can be operated in 318 a semi-automated fashion, following an algorithm implemented in the MATLAB 319 statistical environment (see Appendix for MATLAB script). In MATLAB, each of 1000 transects through the cementum chosen at random using a random number 322 generator (Fig. 2a). For each transect, the distance across the transect that the first 323 pixel above zero appears is saved, which gives the radial length of sampled 324 cementum along the transect. Any transect that is less than the lower standard toolbox' in MATLAB) is then used to identify peaks and troughs in their respective 342 datasets (following multiplication with -1 to convert troughs into peaks) and 343 calculate their difference from the mean greyscale value of that section. This allows 344 peaks and troughs that extend beyond the top and bottom cut-off values steps are taken within the algorithm to ensure that 'piggy-back' features (secondary 349 peaks/troughs along the ascending/descending limbs of genuine increment peaks 350 and troughs) do not affect increment counts (Fig. 2). No peaks are counted that 351 immediately proceed from the last respective peak; so only one peak is counted for 352 every trough (in Fig. 2c. peaks i and ii are not counted). Also, no peak/trough The robustness of the proposed algorithm was tested by applying it to a series 370 of digital sine wave patterns of known increment number between five and 30. 371 Random noise at different degrees was applied to these patterns in a controlled 372 manner by increasing their standard deviation along the y-axis (Fig. 6). Noise was 373 increased incrementally by SNR decrements of 0.1; starting from a SNR of 0.9 and 374 ending at an SNR of 0.1. For each SNR level, increments were counted for 30 sine 375 wave patterns for each count between five and 30. Increment estimates were 376 considered as accurate if the mean estimated count equalled the known increment 377 number to an accuracy of + 0.5. Estimates were considered robust for each count as 378 long as the standard deviation for the 30 counted sine wave patterns was < 1, as 379 values above this may produce estimated increment counts of over 1 year 380 above/below known/expected counts (see Supplement Section 2 'Robustness testing 381 for increment counting algorithm'). 382

Application of cementum increment counting algorithm 383
The increment counting algorithm presented here was used to generate estimates of 384 increment counts for straightened and filtered CT slices for each lower first molar 385 specimen of the 10 Macaca mulatta individuals. We applied the cementum 386 increment counting algorithm to 30 CT slices for each individual, representative of 387 highest cementum increment contrast and quality for each individual (Naji et al., 388 2016). Each straightened and filtered CT dataset was examined by eye, in order to 389 find the regions of highest increment contrast and minimum amounts of complexity 390 in increment patterns (i.e. splitting and coalescence of increments). 1000 transects 391 were plotted through the cementum in each CT slice, and increment pair counts 392 were generated for each transect. The mean increment pair count for all 1000 393 transects was then used as the estimated increment pair count for the slice, and the 394 mean count of the 30 slices rounded to the nearest integer was defined as the 395 estimated increment pair count for that Macaca mulatta individual.

Optimisation of cementum imaging 399
When the X-ray energy was changed in isolation, SNR became consistently lower 400 with increasing X-ray energy, whereas CNR peaked at 20 keV, before steadily falling 401 with increasing X-ray energy beyond this point (Fig. 4a). SNR and CNR steadily 402 improved with increasing exposure time (Fig. 4b). SNR and CNR also improved with 403 increased number of projections, although the relative increase in CNR was marginal 404 between 3001 and 4501 projections (Fig. 4c). SNR steadily increased with larger 405 sample-to-detector distances up to 60 mm (Fig. 4d). Whereas CNR steadily rose from 406 14 mm sample-to-detector distance to a peak at 28 mm, it fell between a sample-to-407 detector distance of 28 mm and 100 mm (Fig. 4d). 408 The image quality of the dataset imaged at 28 mm sample-to-detector distance 409 (14 mm -20 mm), high image contrast resulted between increments, but their 416 boundaries were smoothed and less well defined (case Fig. 3b). For sample-to-417 detector distances above this, the increasing amounts of smoothing diminished the 418 differences in greyscale values between light and dark increments such that by 100 419 mm sample-to-detector distance, they were difficult to distinguish by eye (case Fig. through the same region of cementum in each dataset (Fig. 3e) acquired at different 422 sample-to-detector distances. 423

Cementum imaging results 424
Cementum was clearly visible in each CT dataset as an incremental tissue wrapping 425 around the dentine of tooth roots and comprising a series of radial increments (Figs. 5 426 and 7-8). The cementum could be distinguished from the dentine due to its 427 significantly lower mean grey values, and the cemento-dentine boundary was marked 428 by the characteristic tissues of the granular layer of Tomes and the high-density 429 hyaline layer of Hopewell Smith (Fig. 7). Individual increments were clearly visible 430 within the cementum and could be followed through the entire dataset, both 431 transversely and longitudinally (Fig. 8). Imaging') suggests that both imaging techniques represent the same cementum 436 increments (Fig. 7). Optical differences between increments in histological data were 437 reflected as grey value differences in CT data. Thick, light increments in histological 438 data corresponded to thick, light increments in CT data, and so absorbed a higher 439 proportion of X-rays relative to thin, dark increments (Fig. 7). Volumetric CT data 440 could further be used to help elucidate primary increments from accessory increments 441 in several specimens. Complexities in increment patterns were witnessed 442 intermittently in every Macaca mulatta individual, with individual increments 443 splitting and coalescing to create apparent accessory increments. Following Newham 444 et al. (2020a,b), individual increments could be mapped through the cementum tissue, 445 and the same primary increments could be plotted through the entire scanned tissue volume (Fig. 8) across these complexities, and distinguished from the accessory 447 increments created. Therefore, regions that were confounded by splitting and 448 coalescing of these increments could be distinguished and excluded for analysis of 449 increment counts (Fig. 8). Also, cellular cementum, the tissue with the least 450 chronological precision in its increment periodicity (Naji et al., 2016), could be 451 distinguished from acellular cementum by the presence of cellular voids, and so could 452 be avoided when identifying high-contrast regions of increments with a circum-453 annual periodicity (Fig. 7). These two factors, possible due to the entire coronal 454 (crownward) third of the cementum tissue being imaged, led to the identification of 455 the highest quality regions of circum-annual cementum increments for each specimen 456 (Fig. 7). 457

Validation of cementum increment counting algorithm 458
Robustness tests for the proposed increment counting algorithm suggest that it is 459 reliable for SNRs down to 0.2 (Fig. 5). For each simulated pattern of known 460 increment number, the average value of 30 automated counts was identical to the 461 known count for SNRs between 0.9-0.5. The upper/lower standard deviations of these 462 samples did not exceed one integer above/below the known count (Fig. 6). Between 463 SNR levels of 0.5-0.2, average automated counts only differed from known increment 464 number by a value of one in a single sample (with a known increment number of 465 eight). The standard deviations of automated counts exceed one integer above/below 466 the known count for known increment counts of 22 and 28 (Fig. 6). SNRs of 0.1 467 introduced more errors of between one and two in increment count when compared to 468 the known increment number, and the automated count was outside the region of one 469 standard deviations around the known increment number (Fig. 6). compared to expected counts for our sample based on known age, a Spearman's r of 472 0.77 (p<0.009) and Kendall's τ of 0.71 (p = 0.004) suggest significant correlation 473 between semi-automated increment pair counts and expected numbers of cementum 474 increment pairs. The mean of the semi-automated increment pair counts for each 475 Macaca mulatta individual either met the minimum or maximum expected count 476 based on their known age or fell in-between the two for every sample, apart from the 477 juvenile individual t46 whose mean estimated count was 0.5 years more than the 478 maximum expected count (Table 1 and Fig. 9). Juvenile cementum has been 479 previously shown to contain more complex incrementation and greater amounts of opposite relationship seen between SNR and X-ray energy can also be explained by a 492 diminished X-ray absorption with increased X-ray energy due to an exponentially 493 decreased probability of photoelectric interactions between X-rays and the tissue. 494 CNR has a more complex relationship with each experimental setting, with an 495 optimum setting at a different level compared to SNR for each experimental setting. 496 For instance, we located the optimal energy at TOMCAT for cementum increments in 497 terms of CNR at around 20-21 keV, while SNR was highest at 19 KeV and 498 continuously decreased with higher X-ray energies. 499 The steady increase in SNR with increasing sample-to-detector distance is in with increasing sample-to-detector distance versus the peak in mean greyscale value 507 of cementum (̅̅̅) between 28-60 mm (Fig. 4e). The Paganin phase retrieval 508 algorithm acts as a low pass filter, reducing the image noise in the resultant CT 509 reconstructions. This filtering has been enhanced here with increased sample-to- The first objective of this study was to optimise PPCI through SR CT for studying 522 cementum increments in 3D tomographic data, as an alternative strategy to 523 destructive thin-sectioning and light microscopy for imaging and counting cementum 524 increments. We have shown here that optimised PPCI strategies can overcome the 525 principal caveats of thin-section imaging: namely the destructive sample preparation 526 process and the limited 2D view of tissue that is actually 3D in nature, so lacking 527 context for interpreting complexities in increment patterns; and also limited control 528 over which cementum tissue type is imaged (AEFC versus CIFC). The high image 529 quality offered by SR CT, including phase retrieval offered in PPCI, has provided 530 comparable fidelity for counting individual cementum increments to thin-section 531 histological images of the same regions. The volumetric nature of CT datasets allows 532 navigation through the entire cementum tissue at an isotropic and sub-micrometre 533 nominal spatial resolution. Individual cementum increments can be followed across or tissue distortion through the mechanical cutting process, which can obscure or alter 544 However, our PPCI of cementum through SR CT has also highlighted the 546 sensitivity of cementum image quality to experimental settings. This suggests that 547 optimisation of experimental settings should be conducted preliminary to every 548 cementochronological PPCI experiment using SR CT, in order to ensure optimised 549 image quality for identifying and counting cementum increments. Optimal 550 experimental settings are specific to the optics of the synchrotron beamline and the 551 material properties, size and morphology of the specimen, and so should be 552 investigated when any of these factors are changed. Also, although CT is generally 553 considered to be non-destructive it became apparent during scanning that micrometre-554 scale cracks, which are not visible macroscopically, have formed within the 555 cementum tissue (Supplementary Fig. S3) due the interaction of the hard X-rays with 556 the teeth and/or related effects due to this interaction. Although this damage could not 557 be seen macroscopically, it may indicate that further preparation of teeth and/or 558 adaptation of experimental conditions for SR CT imaging is needed, including tissue 559 dehydration and/or cooling (Peña Fernández et al., 2019). 560

Cementum increment counting algorithm 561
The second objective of this study was to provide a validated, semi-automated and 562 robust algorithm for user-independent counting of cementum increments. The manual 563 counting of cementum increments amongst a restricted number of thin-sections per 564 tooth, plays a central role in the current user-dependent approach for counting 565 cementum increments. This subjectivity has led to a wide range of different 566 accuracies and precisions reported for increment counts and their correlation with 567 known age in animal and human samples. Both accuracy and precision in estimated 568 increment counts correlate with the experience of the researcher when compared to 569 known age (i.e. expected increment count) in validation studies (Naji et al., 2016).

Our algorithm offers a new method for objectively counting cementum increments 571
in a user-independent and semi-automated fashion. This substantially decreases the 572 subjectivity and propensity for human error involved in increment counting. Within 573 the same selected sample of straightened, isolated and filtered PPCI slices, our 574 algorithm requires no further human input for counting cementum increments and will 575 estimate the same increment count regardless of the experience of the researcher. The 576 accuracy and precision of this algorithm has been validated here for both simulated 577 data and our experimental sample of Macaca mulatta cementum. It could also be 578 further assessed in the same quantitative manner with other PPCI cementum data 579 from animals of known age. Such assessment will allow for further optimisation of 580 our algorithm and tailoring for a wide range of PPCI cementum data. 581 Finally, although we state the advantages of PPCI imaging over traditional 582 thin-section histological imaging here, the validation of our application on thin-583 section data of cementum from animals of known age may afford its application for 584 thin-section images. If found to be an accurate method for counting thin-section 585 increments, implementation of our algorithm for thin-sections has potential as an 586 important tool for validating the accuracy of counts estimated by-eye, or even 587 discounting the need for counting increments by-eye completely. 588 589

Conclusion 590
In conclusion, we have undertaken a first systematic experimental study on cementum 591 increment counting for non-fossilised dental tissue, based on a comparison between 592 optimised PPCI through SR CT and thin-section histological imaging. Comparison 593 between these two imaging techniques has shown that PPCI SR CT data can provide increments in the cementum tissue. CT reconstructions are of sufficient quality to 596 count increments semi-automatically using image processing, by defining them as 597 peaks and troughs in greyscale values along transects through the cementum. We have 598 implemented this semi-automated method of increment counting as part of a novel 599 workflow of image processing (cementum isolation, straightening and filtering) and    shown are from m1 teeth of ten Macaca mulatta individuals of known age (see Table  899 1). Cyan boxes indicate the interquartile range around the mean estimated increment