Yi Gu

Yoshito Otake
email: otake@is.naist.jp

Keisuke Uemura

Mazen Soufi

Masaki Takao

Hugues Talbot

Seiji Okada

Nobuhiko Sugano

Yoshinobu Sato

Bone mineral density estimation from a plain X-ray image by learning decomposition into projections of bone-segmented computed tomography

Keywords: 41A05, 41A10, 65D05, 65D17 Bone mineral density (BMD), Computed tomography (CT), Generative adversarial network (GAN), Radiography, Regression network, Representation learning

Osteoporosis is a prevalent bone disease that causes fractures in fragile bones, leading to a decline in daily living activities. Dual-energy X-ray absorptiometry (DXA) and quantitative computed tomography (QCT) are highly accurate for diagnosing osteoporosis; however, these modalities require special equipment and scan protocols. To frequently monitor bone health, low-cost, low-dose, and ubiquitously available diagnostic methods are highly anticipated. In this study, we aim to perform bone mineral density (BMD) estimation from a plain X-ray image for opportunistic screening, which is potentially useful for early diagnosis. Existing methods have used multi-stage approaches consisting of extraction of the region of interest and simple regression to estimate BMD, which require a large amount of training data. Therefore, we propose an efficient method that learns decomposition into projections of bone-segmented QCT for BMD estimation under limited datasets. The proposed method achieved high accuracy in BMD estimation, where Pearson correlation coefficients of 0.880 and 0.920 were observed for DXA-measured BMD and QCT-measured BMD estimation tasks, respectively, and the root mean square of the coefficient of variation values were 3.27 to 3.79% for four measurements with different poses. Furthermore, we conducted extensive validation experiments, including multi-pose, uncalibrated-CT, and compression experiments toward actual application in routine clinical practice.

Introduction

In the last decades, a significant increase in the prevalence of osteoporosis has been observed. It has become the most com-clinical practice, dual-energy X-ray absorptiometry (DXA) [START_REF] Iki | Bone mineral density of the spine, hip and distal forearm in representative samples of the japanese female population: Japanese population-based osteoporosis (JPOS) study[END_REF][START_REF] Mazess | Dual-energy xray absorptiometry for total-body and regional bone-mineral and soft-tissue composition[END_REF][START_REF] Blake | Role of dual-energy X-ray absorptiometry in the diagnosis and treatment of osteoporosis[END_REF][START_REF] Kröger | Bone densitometry of the spine and femur in children by dual-energy x-ray absorptiometry[END_REF][START_REF] O'malley | Trends in dual-energy x-ray absorptiometry in the United States, 2000-2009[END_REF][START_REF] Pisani | Screening and early diagnosis of osteoporosis through X-ray and ultrasound based techniques[END_REF] is considered the gold standard for BMD measurement. In addition, quantitative computed tomography (QCT) [START_REF] Engelke | Clinical use of quantitative computed tomography and peripheral quantitative computed tomography in the management of osteoporosis in adults: The 2007 iscd official positions[END_REF][START_REF] Mueller | Phantom-less QCT BMD system as screening tool for osteoporosis without additional radiation[END_REF][START_REF] Aggarwal | Opportunistic diagnosis of osteoporosis, fragile bone strength and vertebral fractures from routine CT scans; a review of approved technology systems and pathways to implementation[END_REF][START_REF] Ziemlewicz | Opportunistic Osteoporosis Screening: Addition of Quantitative CT Bone Mineral Density Evaluation to CT Colonography[END_REF] has been explored with a view to provide opportunistic screening. Despite the prevalence of osteoporosis worldwide, the diagnosis rates for adults remain unacceptably low in many regions [START_REF] Snodgrass | Osteoporosis diagnosis, management, and referral practice after fragility fractures[END_REF][START_REF] Choi | The prevalence of osteoporosis in Korean adults aged 50 years or older and the higher diagnosis rates in women who were beneficiaries of a national screening program: The Korea National Health and Nutrition Examination Survey 2008-2009[END_REF][START_REF] Papaioannou | The osteoporosis care gap in men with fragility fractures: The Canadian multicentre osteoporosis study[END_REF][START_REF] Mccloskey | The osteoporosis treatment gap in patients at risk of fracture in European primary care: a multi-country cross-sectional observational study[END_REF]. Therefore, there is a strong demand for developing more straightforward methods of measuring BMD, providing opportunistic screening for the early detection of osteoporosis, an asymptomatic disease unless a fracture happens.

Recent studies have focused on the use of deep learning to estimate BMD or diagnose osteoporosis from plain X-ray images [START_REF] Hsieh | Automated bone mineral density prediction and fracture risk assessment using plain radiographs via deep learning[END_REF][START_REF] Yamamoto | Deep learning for osteoporosis classification using hip radiographs and patient clinical covariates[END_REF][START_REF] Jang | Prediction of osteoporosis from simple hip radiography using deep learning algorithm[END_REF][START_REF] Ho | Application of deep learning neural network in predicting bone mineral density from plain X-ray radiography[END_REF][START_REF] Gu | BMD-GAN: Bone mineral density estimation using x-ray image decomposition into projections of bone-segmented quantitative computed tomography using hierarchical learning[END_REF] [START_REF] Wang | Lumbar Bone Mineral Density Estimation From Chest X-Ray Images: Anatomy-Aware Attentive Multi-ROI Modeling[END_REF], which are more widespread modalities than DXA and QCT. Most existing methods directly regress BMD or classify osteoporosis from regions of interest (ROI) of X-ray images, some of which achieved high correlations with DXA-measured BMD [START_REF] Hsieh | Automated bone mineral density prediction and fracture risk assessment using plain radiographs via deep learning[END_REF][START_REF] Wang | Lumbar Bone Mineral Density Estimation From Chest X-Ray Images: Anatomy-Aware Attentive Multi-ROI Modeling[END_REF]. However, these methods required large-scale training datasets and did not provide the spatial density distribution of the target bone. Furthermore, they did not leverage information from CT. The critical requirement of a large-scale training database would limit the application of these methods when the target bone of BMD measurement extends to other anatomical areas (e.g., different positions of the vertebrae, pelvis, sacrum, etc.).

Different from end-to-end regression, [START_REF] Whitmarsh | Reconstructing the 3D shape and bone mineral density distribution of the proximal femur from dual-energy X-ray absorptiometry[END_REF] proposed a method of estimating 3D distributions of QCT-measured BMD from DXA by reconstructing 3D shapes and BMD distributions with statistical shape modeling, improving the diagnosis of osteoporosis and fracture risk assessment. From the viewpoint of X-ray image processing, bone suppression is one of the main topics [START_REF] Suzuki | Image-processing technique for suppressing ribs in chest radiographs by means of massive training artificial neural network (MTANN)[END_REF][START_REF] Yang | Cascade of multi-scale convolutional neural networks for bone suppression of chest radiographs in gradient domain[END_REF][START_REF] Liu | Automatic delineation of ribs and clavicles in chest radiographs using fully convolutional DenseNets[END_REF][START_REF] Eslami | Image-to-images translation for multi-task organ segmentation and bone suppression in chest x-ray radiography[END_REF], enhancing the visibility of other soft tissues to increase the diagnosis rate of softtissue diseases by machines and clinicians. These studies used convolutional neural networks to suppress bones [START_REF] Suzuki | Image-processing technique for suppressing ribs in chest radiographs by means of massive training artificial neural network (MTANN)[END_REF][START_REF] Yang | Cascade of multi-scale convolutional neural networks for bone suppression of chest radiographs in gradient domain[END_REF] or decompose into target soft tissues directly [START_REF] Eslami | Image-to-images translation for multi-task organ segmentation and bone suppression in chest x-ray radiography[END_REF], some of which utilized generative adversarial networks (GANs) [START_REF] Isola | Image-to-image translation with conditional adversarial networks[END_REF] for their realistic synthesis abilities. Because those studies focused on the diagnosis of soft tissues, they did not address the quantitative evaluation of bone decomposition for BMD estimation. While those studies inspired this one, we propose to learn Xray image decomposition into projections of bone-segmented QCT for BMD estimation, targeting both DXA-measured BMD (hereafter "DXA-BMD") and QCT-measured BMD (hereafter "QCT-BMD").

Unlike existing BMD estimation methods, we virtually and automatically mimicked the decomposition procedure (similar to the DXA scanning) and the measurement procedure, which is physically and manually conducted in DXA-based BMD mea- surement, by fully utilizing rich information from QCT in its training phase by predicting bone-segmented digitally reconstructed radiography (DRR), which are projections from segmented QCT. Using GAN, we trained a model to decompose an X-ray image into the density distribution of a partial bone defined clinically for BMD measurement. Despite the difficulty in GAN training, recent studies [START_REF] Zhao | Differentiable augmentation for data-efficient GAN training[END_REF][START_REF] Zhang | Consistency regularization for generative adversarial networks[END_REF][START_REF] Wu | Gradient normalization for generative adversarial networks[END_REF][START_REF] Tseng | Regularizing generative adversarial networks under limited data[END_REF] were able to stabilize the training and reduce the demand for a large-scale database.

For similar purposes, we introduced the hierarchical learning (HL) method so that the small-area regions of the target bone required for BMD estimation were decomposed accurately and stably even without a large training dataset. Fig. 1 illustrates the relationship between X-ray images and DXA-BMDs in our patient dataset, demonstrating the challenge in the task of BMD prediction based on X-ray images, which showed no correlation with BMD. The proposed method, which disentangled soft tissues and bone from an X-ray image, accurately predicted the target bone, showing a high correlation with BMD. Furthermore, we showed improvement by regressing BMD from the learned representation generated by the encoder of the decomposition model, demonstrating the effectiveness of the proposed decomposition learning. We summarize our contributions as follows:

• Proposal of a BMD estimation method through X-ray image decomposition.

• Proposal of a hierarchical learning framework for decomposition into a small target.

• Proposal of a representation learning method for BMD estimation using decomposition learning.

• Extensive validation with real clinical datasets involving variations in patient poses, image compression ratios, and CT calibration. [START_REF] Gu | BMD-GAN: Bone mineral density estimation using x-ray image decomposition into projections of bone-segmented quantitative computed tomography using hierarchical learning[END_REF] and then learns mapping to BMD by regressing from the learned representation. (Right) The visualized procedure of shaping representation and model attention shifting with learning stages. The representations were visualized after being compressed to one-channel images using principal component analysis; the model attention maps were visualized using FullGrad [START_REF] Srinivas | Full-gradient representation for neural network visualization[END_REF] This paper is built upon our conference paper [START_REF] Gu | BMD-GAN: Bone mineral density estimation using x-ray image decomposition into projections of bone-segmented quantitative computed tomography using hierarchical learning[END_REF] with extensions into several aspects. From the methodological aspect, we explicitly targeted QCT-BMD in addition to DXA-BMD. We conducted an ablation study on thresholds for obtaining binarized masks when calculating BMDs from the predicted DRR. We validated the effectiveness of the proposed decomposition learning by regressing BMD from the learned representation with a transformer-based regression network. For the dataset aspect, we included more patients who had undergone X-ray scanning and QCT scanning and had taken BMD measurements, resulting in 305 patients (increased from 200 patients) in [START_REF] Gu | BMD-GAN: Bone mineral density estimation using x-ray image decomposition into projections of bone-segmented quantitative computed tomography using hierarchical learning[END_REF]. We included more Xray images imaged with supine, abduction, and adduction poses in addition to the standing pose. For the experimental aspect, we evaluated the performance of the proposed method using the extended dataset and conducted additional validation for robustness against pose variation, image compression, and uncalibrated CT. We also investigated error sources by identifying outliers.

Related Work

Grading osteoporosis from X-ray images

Jang et al. ( 2021) used a VGG16 model with integrated Nonlocal Neural Networks to classify osteoporosis from a manually cropped X-ray image. [START_REF] Hsieh | Automated bone mineral density prediction and fracture risk assessment using plain radiographs via deep learning[END_REF] and [START_REF] Ho | Application of deep learning neural network in predicting bone mineral density from plain X-ray radiography[END_REF] developed BMD estimation frameworks using localization and regression models, requiring large-scale databases to achieve high accuracy. [START_REF] Wang | Lumbar Bone Mineral Density Estimation From Chest X-Ray Images: Anatomy-Aware Attentive Multi-ROI Modeling[END_REF] proposed a Transformerbased model for BMD estimation from chest X-ray images with multiple ROIs, simultaneously capturing local and global information to improve accuracy. Although conventional methods strove to localize the target regions from an X-ray image (which lies in 2D space), they did not address the entanglement of bone and soft tissues (which align in 3D space). We argue the unsolved tissue entanglement, proposing the utilization of CT information during training in order to allow the estimation of the spatial density of bone from a plain X-ray image. [START_REF] Yamamoto | Deep learning for osteoporosis classification using hip radiographs and patient clinical covariates[END_REF] proposed an ensemble model combining multiple convolutional neural networks, trying to classify osteoporosis with patients' clinical covariates that are fused in the last layer of the model. Though most clinical covariates are easy to obtain in routine clinical practice and reasonably improve prediction accuracy, we do not use patient characteristics in our method for better clarification of the improvement.

Representation learning

Representation learning aims to learn efficient representations of data whose information is easier to be extracted for downstream tasks, improving effective learning [START_REF] Bengio | Representation Learning: A Review and New Perspectives[END_REF]. Many representation learning methods have been proposed for solving tasks in the fields of natural language processing and general image processing. [START_REF] Devlin | BERT: Pre-training of Deep Bidirectional Transformers for Language Understanding[END_REF][START_REF] Noroozi | Representation Learning by Learning to Count[END_REF][START_REF] Chen | Exploring Simple Siamese Representation Learning[END_REF][START_REF] Donahue | Large Scale Adversarial Representation Learning[END_REF]. However, those methods usually rely on largescale pretraining (which is hard to achieve using medical images) and they did not utilize the specialty such as penetrability in medical images. [START_REF] Zhou | Latent Correlation Representation Learning for Brain Tumor Segmentation With Missing MRI Modalities[END_REF] proposed a representation method for brain tumor segmentation with missing modes of magnetic resonance images, utilizing the similarity between multi-modal MR images. [START_REF] Liu | Incomplete multi-modal representation learning for Alzheimer's disease diagnosis[END_REF] proposed a multiview missing data completion framework to learn multi-modal representation for diagnosing Alzheimer's disease. In addition to providing the spatial density of bone for deriving BMD, our decomposition learning can serve as representation learning for grading osteoporosis from a plain X-ray image.

Image-to-image translation

An image-to-image translation (I2I) model tries to map between source and target image domains. Typically, GANs were often used to train the I2I model [START_REF] Isola | Image-to-image translation with conditional adversarial networks[END_REF][START_REF] Zhu | Unpaired Image-to-Image Translation Using Cycle-Consistent Adversarial Networks[END_REF][START_REF] Wang | High-resolution image synthesis and semantic manipulation with conditional GANs[END_REF]. We regard the decomposition task as an I2I task that aims to translate a plain X-ray image into the decomposed DRR. Depending on the availability of data, the I2I model can be trained in a paired [START_REF] Isola | Image-to-image translation with conditional adversarial networks[END_REF] and unpaired [START_REF] Zhu | Unpaired Image-to-Image Translation Using Cycle-Consistent Adversarial Networks[END_REF] manner while the paired-trained model usually learns faster and better. Our method appreciated the paired [START_REF] Uemura | Automated segmentation of an intensity calibration phantom in clinical CT images using a convolutional neural network[END_REF], bone segmentation [START_REF] Hiasa | Automated muscle segmentation from clinical CT using Bayesian U-Net for personalized musculoskeletal modeling[END_REF], 2D-3D registration to the X-ray image [START_REF] Otake | Intraoperative image-based multiview 2D/3D registration for image-guided orthopaedic surgery: incorporation of fiducialbased C-arm tracking and GPU-acceleration[END_REF] and DRR generation by projecting QCT.

training by utilizing 2D-3D registration [START_REF] Otake | Intraoperative image-based multiview 2D/3D registration for image-guided orthopaedic surgery: incorporation of fiducialbased C-arm tracking and GPU-acceleration[END_REF] to align X-ray images and decomposed DRRs, which are generated from 3D CT images.

Method

Fig. 2 shows an overview of the proposed method that uses decomposition-based representation learning for BMD estimation. In decomposition training with the proposed BMD-GAN, an image synthesis model G = {M E , M D }, which consists of an encoder backbone M E and a decoder M D , decomposes the X-ray image into the DRR of the proximal femur region [hereafter "proximal femur region DRR (PF-DRR)"], whose average intensity provides the predicted BMD. We adopted HRFormer [START_REF] Yuan | HRFormer: High-resolution vision Transformer for dense predict[END_REF] as the encoder backbone. Our BMD-GAN applies a hierarchical framework during the training in which the synthesis model is first trained to extract the pelvis and femur bones in stage-one and then the proximal femur region in stage-two. During regression (stage-three) training, the encoder backbone M E is inherited and connected to a Transformer-based regression head M R to become a regression network R = {M E , M R } for end-to-end training of regressing BMD. The regression head tries to regress BMD from a learned representation instead of the decomposed PF-DRR for efficient forwarding. Note that our method is network-agnostic. Other backbone models and regression heads can be adopted to achieve decomposition learning and regression learning, respectively.

Dataset construction

Fig. 3 illustrates the overall dataset construction procedure. In this study, we constructed two datasets for training different stages: 1) dataset A, which is used for stage-one training containing X-ray images; QCT; the 3D segmentation masks of the pelvis and femur, which were automatically obtained by applying the Bayesian U-net [START_REF] Hiasa | Automated muscle segmentation from clinical CT using Bayesian U-Net for personalized musculoskeletal modeling[END_REF]; and GT bone DRR, which were created from QCT using 2D-3D registration [START_REF] Otake | Intraoperative image-based multiview 2D/3D registration for image-guided orthopaedic surgery: incorporation of fiducialbased C-arm tracking and GPU-acceleration[END_REF] followed by projection with the 3D mask, and 2) dataset B, which is used for other stages' training containing X-ray images; QCT; the 3D mask of the proximal femur region, which is obtained via manually labeled bony landmarks defined in [START_REF] Uemura | Development of an open-source measurement system to assess the areal bone mineral density of the proximal femur from clinical CT images[END_REF]) by an expert clinician; and the GT PF-DRR. The construction procedure of dataset B followed [START_REF] Uemura | Development of an open-source measurement system to assess the areal bone mineral density of the proximal femur from clinical CT images[END_REF]. The intensity-based 2D-3D registration using gradient correlation similarity metrics and the CMA-ES optimizer [START_REF] Otake | Intraoperative image-based multiview 2D/3D registration for image-guided orthopaedic surgery: incorporation of fiducialbased C-arm tracking and GPU-acceleration[END_REF] was performed on each patient's X-ray image and QCT. GT DXA-BMD and GT QCT-BMD values were associated with each patient's data in dataset B. All Xray images and DRRs were normalized into the 256 × 512 size via central cropping and resizing. The aspect ratio of the original X-ray images varied from 0.880 to 1.228 (width/height). We first split them horizontally in half at the center. Then the side with the target hip was reshaped to a predefined image size (256 × 512 in this experiment) by aligning the center of the image and cropping the region outside the image after resizing to fit the shorter edge of the width and height.

X-ray image Decomposition

The GAN with conditional discriminators was used to train the decomposition model. We followed most settings used in Pix2PixHD [START_REF] Wang | High-resolution image synthesis and semantic manipulation with conditional GANs[END_REF], including the multi-scale discriminators and the Feature Matching loss, among others. Instead of the ResNet Generator used in Pix2PixHD [START_REF] Wang | High-resolution image synthesis and semantic manipulation with conditional GANs[END_REF] and CycleGAN [START_REF] Zhu | Unpaired Image-to-Image Translation Using Cycle-Consistent Adversarial Networks[END_REF], we adopted the state-of-the-art model HRFormer [START_REF] Yuan | HRFormer: High-resolution vision Transformer for dense predict[END_REF], a transformer-based model for segmentation, to be the backbone of the generator G = {M E , M D }, namely, the HRFormer Generator. Instead of using the hierarchical structure of the generator used in Pix2PixHD [START_REF] Wang | High-resolution image synthesis and semantic manipulation with conditional GANs[END_REF], we applied the HL framework in which a two-stage training is used. In stage one, the model is trained to decompose an X-ray image into pelvis and femur bones; in stage two, the target is transferred to the proximal femur region. During training, we used the adversar- ial loss L GAN , which is defined as

L GAN (G, D) = E (I Xp ,I DRR ) log D(I Xp , I DRR ) + E I Xp log 1 -D(I Xp , G(I Xp )) , (1)
where D, I Xp , and I DRR are the discriminator, X-ray image, and decomposed DRR, respectively. We denote E I Xp ≜ E I Xp ∼p data (I Xp ) and E (I Xp ,I DRR ) ≜ E (I Xp ,I DRR )∼p data (I Xp ,I DRR ) . Furthermore, we used the Feature Matching loss L FM proposed in [START_REF] Wang | High-resolution image synthesis and semantic manipulation with conditional GANs[END_REF], given by

L FM (G, D) = E (I Xp ,I DRR ) T i=1 1 N i FM(G, D (i) ), (2)
where D (i) , T , and N i denote the ith-layer feature extractor of the discriminator D, the total number of layers, and the number of elements in each layer, respectively. The layer feature matching criterion FM is defined as

FM(G, D (i) ) = D (i) (I Xp , I DRR ) -D (i) (I Xp , G(I Xp )) 1 . (3)
We did not use a perceptual loss because a well-pretrained perceptual model is difficult to obtain using a limited-size dataset.

We instead used a simple L1 loss L L1 defined as

L L1 (G) = E (I Xp ,I DRR ) I DRR -G(I Xp ) 1 . (4)
To maintain the consistency of the structure between the fake DRR G(I Xp ) and the true DRR I DRR , we regularized the generator with the gradient matching constraints proposed in [START_REF] Penney | A comparison of similarity measures for use in 2-D-3-D medical image registration[END_REF], using the gradient correlation loss [START_REF] Hiasa | Cross-Modality Image Synthesis from Unpaired Data Using CycleGAN[END_REF] L GC defined as

L GC (G) = E (I Xp ,I DRR ) GC(G).
(5)

The criterion GC is defined as

GC(G) = NCC(∇ x I DRR , ∇ x G(I Xp )) + NCC(∇ y I DRR , ∇ y G(I Xp )), ( 6
)
where NCC(A, B) is the normalized cross-correlation of A and B, and ∇ x and ∇ y are the x and y components of the gradient vector, respectively. We introduced hyper-parameters λ L1 , λ GC , and λ F M to balance the importance of the loss terms summarized as

L dec = λ L1 L L1 (G) + λ GC L GC (G) + λ FM k=1,2,3 L FM (G, D k ), ( 7
)
where the multi-scale discriminators D 1 , D 2 , and D 3 were used under three resolutions as in [START_REF] Wang | High-resolution image synthesis and semantic manipulation with conditional GANs[END_REF]. Thus, the full objective for decomposition training is defined as

min G L dec + max D 1 ,D 2 ,D 3 k=1,2,3 L GAN (G, D k ) . ( 8
)
Both stage one and stage two training use the same loss functions.

Once the decomposition was completed, the average intensity of the predicted PF-DRR was calculated. Note that pixels whose intensities were equal to or greater than the threshold t were averaged for deriving BMD. The PF-DRR-average of all training datasets was linearly fitted to DXA-BMD and QCT-BMD to obtain the slope and intercept, which were used to convert the PF-DRR-average to the BMDs of the test dataset. The threshold t tuning procedure with GT-PF-DRR is shown in Fig. 4, where we empirically chose values for different situations.

Regression from learned representation

For the regression head M R , we used several Transformer blocks before the final fully-connected layer to extract BMDrelated features from a representation learned by the encoder backbone from an X-ray image. We borrowed the architecture of the CoAtNet [START_REF] Dai | CoAtNet: Marrying convolution and attention for all data sizes[END_REF] as our regression head for its high performance on image classification tasks by combining convolution blocks and Transformer blocks. However, we only used the Transformer blocks from it, which applies inputindependent relative attention. Unlike in CoAtNet whose first Transformer block takes a representation (i.e., feature maps) with 1 8 size of the input image, our regression head takes a higher-resolution representation with 1 4 size learned by the encoder M E . The inherited encoder and the regression head form our regression network R.

During training, a simple L1 loss was used. In addition, we introduced a predefined sample weight into the loss to improve unbalanced learning since BMD values are not uniformly distributed in the real world. The sample weight is calculated according to the distance d between the sample's BMD y and the average BMD ŷ of the dataset. The distance is defined as

d = |y -ŷ|. ( 9
)
For each X-ray image, its sample weight w is defined as where the d min and d max are the minimum and maximum distances in the dataset. Thus, the objective of regression training is defined as min

w = 1.5 - d -d min d max -d min , ( 10
)
R E (I Xp ,y,w) w y -R(I Xp ) , (11)
where we denote E (I Xp ,y,w) ≜ E (I Xp ,y,w)∼p data (I Xp ,y,w) . Because we treat the DXA-BMD and QCT-BMD estimations as independent tasks, the sample weight for an X-ray image can be different when regressing different BMDs.

Experiments and Results

We validated our method under the constraint of limited datasets. Based on the constructed datasets (A and B), we conducted two experiments for 1) evaluating the performance of the BMD estimation of our method and 2) conducting an ablation study on the proposed hierarchical learning under different backbone models.

Dataset

The two datasets (A and B) used in the experiments are summarized in Table 1. Ethical approval was obtained from the Institutional Review Boards (IRBs) of the institutions participating in this study (IRB approval numbers: 21115 for Osaka University Hospital and 2021-M-11 for the Nara Institute of Science and Technology). The constructed dataset A contained 275 cases. Each case had an X-ray image, and its paired bone DRRs of the left and right sides were split by the vertical middle line, resulting in 525 image pairs after excluding the images with hip implants. The constructed dataset B contained 305 cases obtained retrospectively from 305 patients (250 females) who underwent primary total hip arthroplasty between May 2011 and December 2015. Each case has three or four Xray images with different poses and their paired PF-DRR of one side with their ground-truth DXA-BMD and QCT-BMD. Fig. 5 visualized the T-SNE of the collected multi-pose X-ray images of dataset B, which includes poses of standing, adduction, abduction, and supine, suggesting high potential for conducting multiple clinical validations. The calibration phantom (B-MAS200, Kyoto Kagaku, Kyoto, Japan) [START_REF] Uemura | Automated segmentation of an intensity calibration phantom in clinical CT images using a convolutional neural network[END_REF], which is used to convert radiodensity [in Hounsfield units] to the ground truth QCT-BMD (in mg/cm 3 ), contains known densities of hydroxyapatite Ca 10 (PO 4 ) 6 (OH) 2 . All CT images used in this study were obtained using the OptimaCT660 scanner (GE Healthcare Japan, Tokyo, Japan), and all DXA images of the proximal femur were acquired for the operative side (Discovery A, Hologic Japan, Tokyo, Japan) to obtain the ground truth DXA-BMD.

BMD estimation experiment

Experimental settings

In this experiment, we tried to estimate DXA-BMD and QCT-BMD as independent tasks. Dataset A was used for training stage one of the proposed method, containing 492 pairs of X-ray images and GT bone DRR. Five-fold cross-validation was performed on stages two and three using dataset B. We compared our method with the conventional one proposed by [START_REF] Hsieh | Automated bone mineral density prediction and fracture risk assessment using plain radiographs via deep learning[END_REF], which directly regresses BMD from an ROI of an X-ray image obtained using a localization model that tries to find the proximal femur. When applying the conventional method, we used pre-cropped ROI of the proximal femur of X-ray images to simulate the most ideal situation. The proposed and conventional methods were trained using the same number of patient data to evaluate the effects of utilizing CT data in training by the proposed method. We evaluated BMD estimation performance by Pearson correlation coefficient (PCC), intraclass correlation coefficient (ICC), mean absolute error (MAE), and standard error of estimate (SEE). Statistical significance was evaluated using the single-factor repeated measures analysis of the variance model. P-values were used to denote statistical significance. We also evaluated the reproducibility under the four poses (standing, supine, abduction, and adduction) using root mean square of the coefficient of variation (RMS-CV). The results will be reported in Sec. 4.2.2. In addition to using QCT and X-ray images of four poses, we also conducted three more experiments by 1) training with the limited-pose dataset, 2) training with uncalibrated GT-DXA-BMD (g/cm2) GT-DXA-BMD (g/cm2) GT-DXA-BMD (g/cm2) Predicted DXA-BMD (g/cm2) 1) Convenional reg.

3) Proposed decomp. + reg. 2) Proposed decomp.

1) Convenional reg.

3) Proposed decomp. + reg. 2) Proposed decomp. is the direct regression method trained and tested using the pre-cropped ROI of the proximal femur of X-ray images. 2) Proposed decomp. is the proposed decomposition method that predicts 2D BMD distribution. 3) Proposed decomp. + reg. is the proposed regression method based on method 2). CT, and 3) testing with compressed X-ray images for robustness validation against 1) pose variation, 2) un-calibration, and 3) image compression reported in Sec. 4.2.3, Sec. 4.2.5, and Sec. 4.2.4, respectively. These additional experiments were designed to validate the generalizability of the proposed method considering the situations that are more likely to happen in realworld applications. For the training of all stages, we used the learning rate policy of stochastic gradient descent with warm restarts [START_REF] Loshchilov | SGDR: Stochastic gradient descent with warm restarts[END_REF] and AdamW optimizer [START_REF] Loshchilov | Decoupled weight decay regularization[END_REF]. In stage three training, we fixed the weights of encoder M E by several epochs at the beginning and then trained R by the remaining epochs.

Overall performance of BMD estimation

Comparisons of BMD estimation performance between methods are shown in Fig. 6. The proposed decomposition method ("proposed decomp." in Fig. 6), which is the G model trained with hierarchical learning, improved DXA-BMD estimation PCC and ICC from 0.667 and 0.609 to 0.864 and 0.841, respectively. A similar trend was also observed for the QCT-BMD estimation task, where the PCC and ICC were improved from 0.764 and 0.730 to 0.861 and 0.841 by the proposed decomposition method, respectively. The estimation performance was further improved by learning regression from a representa-tion learned by encoder M E via decomposition training ("proposed decomp. + reg." in Fig. 6), achieving the highest PCC and ICC of 0.880 and 0.848s, respectively, for DXA-BMD estimation, and 0.920 and 0.909, respectively, for QCT-BMD estimation. An evaluation summary for full experiments and metrics is shown in Table 2. Under the guidance of decomposition training, regressing BMD already achieved high scores while freezing model weights of encoder M E , demonstrating successful representation learning by decomposition. Furthermore, final tuning for the encoder and regression head improved performance, achieving the highest scores on all metrics we used. Notice that the decomposition-guided regression even outperformed the derivation of QCT-BMD from GT-DRR on most metrics. RMS-CV values for our methods were 3.27 to 3.79%, while the acceptable RMS-CV in the DXA-BMD measurement is 1.8% [START_REF] Lewiecki | Best Practices for Dual-Energy X-ray Absorptiometry Measurement and Reporting: International Society for Clinical Densitometry Guidance[END_REF]. However, this criterion assumes that the patient pose should be reproduced as much as possible (whereas in our conditions, the patient poses vary considerably). A video showing prediction samples can be found in the supplemental materials.

Robustness against pose variation

Unlike DXA scanning, which requires patients to be in a supine pose, X-ray imaging does not restrict the scanning pose, Sta., Sup., Abd., Add.

Sta., Sup.

Sta., Sup. Abd., Add. Sta., Sup.

Abd., Add. resulting in pose variations, which imply challenges when predicting the BMD from different poses. To validate robustness against pose variation, we conducted additional experiments using training data with poses limited to standing and supine, which are more commonly accessible in routine clinical practice, reducing the number of training data from 963 to 481 per fold.

Poses in test dataset

Fig. 7 compares the BMD estimation performance of the proposed method between training with full poses (standing, supine, abduction, adduction) and limited poses (standing, supine). For the seen poses (standing and supine), limitedpose training showed a small reduction of PCCs from 0.882 to 0.873 and from 0.919 to 0.916 for DXA-BMD and QCT-BMD, respectively. For the unseen poses, abduction and adduction, limited-pose training showed a reduction of PCCs from 0.879 to 0.862 and from 0.922 to 0.912 for DXA-BMD and QCT-BMD, respectively. Although such slight degradations were observed, statistical tests via the Tukey HSD pose-wise showed p-values all higher than 0.05, failing to indicate significant differences between training with full and limited-pose data. Furthermore, we checked the prediction consistency between poses. The proposed method showed all PCCs higher than 0.9 for the densely compared poses on not only seen poses but also unseen poses. Detailed pose prediction consistency can be found in the supplementary materials.

Robustness against image compression

We validated the proposed method for robustness against image compression since the compressed images are much easier to obtain in routine clinical practice. This compression validation is based on full-pose training experiments reported in Sec 4.2.2. We evaluated the performance of the trained models using X-ray images compressed by the JPEG lossy method.

Fig. 8 shows the relations between the average file sizes of compressed images and the performance of BMD estimation. As the file size has been significantly reduced by lower precision and the JPEG algorithm with 100 quality, the proposed decomposition-based regression method's performance remains almost unchanged with those with no compression. Furthermore, we explored the tolerance limitation by inputting super-compressed images after reducing the JPEG quality. We found that using the half quality, JPEG 50, slightly hurt overall performance, reducing PCCs from 0.880 and 0.920 to 0.877 and 0.918 for DXA-BMD and QCT-BMD estimations, respectively, with only 1.3% average file size of the original, demonstrating high robustness against significant loss of available information in X-ray images. Severe performance degradation was observed from JPEG 10, where the PCCs were reduced to 0.790 and 0.783 for DXA-BMD and QCT-BMD, respectively. Visualizations of the compression impact on our method can be found in the supplementary materials.

Robustness uncalibrated CT

Our method utilized CT information to improve accuracy and efficiency; however, quantitative CT with phantom calibration is required, while large-scale CT datasets, which do not contain phantoms for calibration, were not utilized. On the other hand, uncalibrated CT introduced the domain-shifting problem, which is a long-standing challenge in image processing. We validated the proposed method by using uncalibrated CT.

Without using QCT, the proposed method of Decomp. + Reg. achieved a PCC of 0.872 and 0.922 and ICC of 0.838 and 0.910 for the DXA-BMD and QCT-BMD estimation tasks with full-pose training, respectively. For limited-pose training, the proposed method achieved a PCC of 0.847 and 0.903 and ICC of 0.823 and 0.893 for the DXA-BMD and QCT-BMD estimation tasks, respectively. The Tukey HSD test was performed on a pose-by-pose basis on the predictions given by the proposed method with QCT and uncalibrated CT for full-pose training and limited-pose training separately, and no significant differences were found between them. The results were expected as the CT data obtained in our dataset were all from the same type of scanner. However, the results could be considered as an upper bound of what can be achieved by training using uncalibrated CT.

Source for large error

We investigated the error source of the proposed method in the experiments that use full-pose training data and QCTderived DRR. Fig. 9 shows Bland-Altman plots of the DXA-BMD estimation task (left) and the QCT-BMD estimation task (right), respectively, by the proposed decomposition-guided regression method. Points beyond the 95% limits of agreement (i.e., mean ± 1.96SD) are considered sample outliers, where the SD is the standard deviation of the difference between prediction and ground truth. Each point represents an X-ray image in the corresponding task. The trend showed a bias underestimating and overestimating high-BMD and low-BMD patients, respectively. One possible bias source could be the lack of training data in those BMD areas, which introduced unbalanced learning, which is a long-standing problem in real-world applications [START_REF] Li | Targeted Supervised Contrastive Learning for Long-Tailed Recognition[END_REF][START_REF] Alshammari | Long-Tailed Recognition via Weight Balancing[END_REF]. Bland-Altman plots by the conventional and proposed decomposition methods can be found in supplemental materials.

For the DXA-BMD estimation task, there are 57 outlier samples in Fig. 9. We carefully checked their X-ray imaging conditions, including the tension peak of the X-ray tube in kV, exposure, and magnification; however, we found no particular trend or strong correlation with the estimation error. Additionally, we checked the correlation between the estimation error and the patient's femur position with respect to the pelvis (i.e., joint angle). Fig. 10 shows the correlation to the flexion angle, where the estimation outlier rates are 3.93% and 22.41% in the normal-and abnormal-degree patients, respectively, implying potential estimation error introduced by patient positioning [START_REF] Uemura | The effect of patient positioning on measurements of bone mineral density of the proximal femur: a simulation study using computed tomographic images[END_REF]. When correlating to other angles, the estimation outlier rates in the normal-and abnormal-degree pa- tients are 4.18% and 11.11%, respectively, for the abduction angle, and 5.37% and 8.96%, respectively, for the external rotation (ER) angle. Scatter plots representing the correlation between the estimation error and abduction and ER angles can be found in supplemental materials. Another possible error source could be hip diseases in patients. In our dataset, most patients had developed hip diseases (e.g., osteonecrosis and osteoarthritis), and some even had developed severe forms of these diseases, which may have caused the measurement to be inaccurate. We identified four outlier patients (outliers 1 to 4) in Fig. 9 (left), whose points of all poses were indicated as outlier samples. The DXA-BMD of the four outlier patients were in the high-BMD area. Fig. 11 (left) shows the X-ray images and their paired PF-DRRs of the outlier-1 patient. The patient's Kellgren-Lawrence grade [START_REF] Kohn | Classifications in Brief: Kellgren-Lawrence Classification of Osteoarthritis[END_REF] was 4 (the highest level), and the patient was classified as level II of Crowe [START_REF] Sugano | The morphology of the femur in developmental dysplasia of the hip[END_REF]; i.e., high-severity hip osteoarthritis. Referring to the QCT-BMD of outlier 1, we found that the QCT-BMD is below the distribution average, and our method accurately predicted it.

A similar analysis was performed for the QCT-BMD estimation task as well.

Implementation details

The encoder backbone we adopted is HRFormer-B, which is a middle-size model. Our decoder followed the upsampling parts of the GlobalGenerator in Pix2PixHD, consisting of several convolution layers and transpose-convolution layers. In the regression head, we used Transformer blocks of S3-TFM Rel and S4-TFM Rel from CoAtNet We replaced the Batch Normalization and Instance Normalization with Group Normalization, except for the Layer Normalization used in Transformer. We set λ L1 , λ GC , and λ F M to 100, 1, and 10, respectively. The initial learning rate and weight decay were 6 × 10 -5 and 1 × 10 -2 , respectively, except for tuning the end-to-end regression in stage three, which used 8 × 10 -6 as the initial learning rate. For all our methods (i.e., decomposition and decomposition + regression), we train 630 epochs in stage one. For the decomposition method, we trained 150 epochs in stage two. For decomposition + regression, we trained 630 and 310 epochs in stages two and three, respectively, and tuned end-to-end regression with 70 epochs. Data augmentation included random rotation(±25), shear(±9), translation(±0.3), scaling(±0.3), brightness(±0.5), contrast(±0.5), and horizontal and vertical flipping. When implementing the conventional method for BMD estimation, we followed most settings and training protocols; however, the total number of epochs was set to 400, and we did not perform validation to capture the best model during training because we observed worse performance with our limited-size dataset.

Ablation study on hierarchical learning 4.3.1. Experimental settings

The ablation study aims to compare decomposition performances between 1) model backbones and 2) models with and without hierarchical learning. In addition to the ResNet and HRFormer, we compared two more models that were used for semantic segmentation as the backbone of the generator-DAFormer [START_REF] Hoyer | DAFormer: Improving network architectures and training strategies for domain-adaptive semantic segmentation[END_REF]HRNetV2 Wang et al. (2021), namely, DAFormer Generator and HRNet Generator, respectively. Though these state-of-the-art models have been proven to have high performance in segmentation tasks, their ability to decompose images has not been thoroughly assessed. We set the ResNet Generator without HL as the baseline method. Dataset A was used to train stage-one models. For stage two, we only used 200 cases of standing pose from dataset B with five-fold cross-validation. We also compared our best decomposition method with the conventional method [START_REF] Hsieh | Automated bone mineral density prediction and fracture risk assessment using plain radiographs via deep learning[END_REF] on the DXA-BMD estimation task using our limited dataset. To evaluate the performance of image decomposition, we used the peak signal-to-noise ratio (PSNR), multi-threshold dice coefficient (DC), ICC, and PCC of the average intensity of PF-DRR. To evaluate BMD estimation, we used ICC, PCC, MAE, and SEE.

Decomposition performance

The decomposition accuracy, as evaluated by PSNR and DC is shown in Fig. 12 (a), where a significant improvement in HL was observed. The high performance of the HRFormer Generator with HL in DC indicated the ability to maintain the silhouette of the decomposed structure. The high PSNR suggested the superior capability of the quantitative decomposition compared with the same generator without HL and the baseline method. Fig. 12 (b) shows the training progress for each backbone with and without HL, in which the robust convergence was achieved consistently using HL, even with only a few epochs. One case was randomly chosen to track progress during training, which is shown in Fig. 12 (d) and (e). The qualitative comparison demonstrated that the target region was well-formed in the early epoch using HL, suggesting the effectiveness of HL. A summary of the experimental results for all backbones is shown in Table 3.

BMD estimation accuracy

A comparison of the BMD estimation performance between the conventional method and the proposed HRFormer Generator with HL is shown in Fig. 13 The proposed method achieved high PCC of 0.882 and 0.888 (compared to 0.361 and 0.447), respectively, demonstrating the effectiveness of the estimation strategy of the proposed method that extracts the density distribution of the target region of the bone. We evaluated 13 more cases for which repeated X-ray images (acquired in the standing and supine positions on the same day) were available. The average coefficient of variation was 3.06% ± 3.22% when the best model, HRFormer with HL, was used.

Discussion

Generalizability

We showed the robustness of the proposed method against pose variation, image compression, and uncalibrated CT toward clinical application. As shown in Fig. 7, training with fewer poses did not result in significant degradation even during evaluations with X-ray images of unseen poses. The dense comparison between predictions from different poses suggested the high consistency against pose variations for X-ray images of seen and unseen poses. In the image compression experiment, our method showed high robustness against image compression. The performance was stable unless the X-ray images were heavily compressed, as shown in Fig. 8. Training our method with uncalibrated CT achieved PCC of 0.872 and 0.922 for DXA-BMD and QCT-BMD estimation, respectively, which are slightly different from training with QCT as shown in Table 2. Our results suggested that the calibration may not be needed as long as CT images are taken using the same type of scanner.

Difference in DXA-BMD and QCT-BMD estimation

We noticed the performance gap between estimating DXA-BMD and QCT-BMD. When deriving BMD from GT-PF-DRR, the performance for DXA-BMD, which achieved a PCC of 0.930, is much higher than that for QCT-BMD, which achieved a PCC of 0.883. However, our decomposition-guided regression method showed the opposite by achieving a PCC of 0.880 and 0.920 for DXA-BMD and QCT-BMD, respectively in Fig. 6. The 2D modalities such as DXA and PF-DRR are projected from 3D spaces, losing partial information of posterior parts. The average intensity of PF-DRR may be better to fit DXA-BMD because PF-DRR and DXA are in the same 2D space. However, one may need to consider the information loss that occurs during the projection when using DXA-BMD as the ground truth. On the other hand, the QCT-BMD, which is derived via a 3D modality, appears to be easier to regress by models not only for our method but also the conventional method we compared with. QCT-BMD could serve as additional supervision for potential improvement in DXA-BMD estimation. Furthermore, this phenomenon may suggest the model's ability to infer 3D information from a 2D modality, encouraging us to recover the 3D spatial density of bone from a plain X-ray image in the future.

Learning and predicting 2D BMD distribution

A unique feature of our method is that it produces the pixelwise BMD estimation, representing the 2D BMD distribution of the PF bone in addition to the predicted DXA-BMD in stage 2, improving explainability. Indeed, one can tell that the predicted BMD is incorrect before knowing the predicted BMD value if the model fails to produce a reasonable PF-DRR, which inspires us to develop a failure detection mechanism for actual clinical applications. This ability was obtained by learning the mapping from an X-ray image to PF-DRR, suggesting better use of the texture pattern inside the bone, which is potentially helpful for the fracture risk prediction, as it has been shown to be factorized by the texture of the bone [START_REF] Dong | Random field assessment of inhomogeneous bone mineral density from DXA scans can enhance the differentiation between postmenopausal women with and without hip fractures[END_REF][START_REF] Farzi | Quantitating Age-Related BMD Textural Variation from DXA Region-Free-Analysis: A Study of Hip Fracture Prediction in Three Cohorts[END_REF].

Method limitations

Although the proposed method demonstrated high performance and robustness in limited datasets, the large-error source is unclear, and the performance on a large-scale dataset was not evaluated. We will collect more data and validate the method using larger datasets collected from different centers, involving more variation. The error source analysis showed that hip disease in patients resulted in erroneous segmentation of the PF region, which introduced estimation error. We will improve the segmentation processing pipeline to handle the diseased hip more robustly. Although our method can localize the proximal femur region from an X-ray image, detecting key points is unavailable. We will incorporate key points detection, which could be useful in actual applications and potentially improve the model's understanding of the patient's hip joint pose. Our method used 2D-3D registration, which may limit the application. We will investigate the performance with unpaired decomposition, which does not require registration. Furthermore, we will generalize our method to other target structures, such as the spine.

Summary

We developed a technique for estimating BMD from a plain X-ray image through decomposition learning using CT, which additionally provides the spatial density distribution of the target bone. Our method disentangles the soft tissue and bone from an X-ray image to improve the BMD estimation accuracy, efficiently leveraging information from a small number of datasets consisting of X-ray images, BMD values, and CT of the same patient. We first used the proposed HL framework to train a model to decompose an X-ray image into PF-DRR, then trained a regressor to estimate BMD from the representation learned via decomposition learning. Our method showed significant performance on DXA-BMD estimation and QCT-BMD estimation tasks with limited datasets. The proposed HL framework stabilized the decomposition training and improved decomposition accuracy as shown in the ablation study. Furthermore, we demonstrated the high robustness of the proposed method of BMD estimation with the validations of multipose X-ray images, image compression, and uncalibrated CT, demonstrating the high potential for opportunistic screening in clinical practice.

Fig. 1 .

 1 Fig. 1. Relationship between the intensities of the proximal-femur region of the DRR and BMD values in datasets of 305 patients. (Left) scatter plots showing the correlation of the average intensity of ground-truth (GT) DRR and predicted DRR with DXA-BMD and QCT-BMD evaluated by Pearson correlation coefficient. (Right) proximal femur ROIs of five representative cases. ROIs #1 and #2 have similar X-ray intensities but significantly different BMDs, whereas ROIs #3 and #4 have similar BMDs but significantly different X-ray intensities. The predicted DRRs correctly recovered the intensity of QCT DRR, regardless of the intensity of the input X-ray image.

  Fig. 2. (Left) Overview of the proposed method consisting of three training stages, where the model first learns X-ray image representation by decomposition training in stage-one and stage-two with HL applied (Gu et al., 2022) and then learns mapping to BMD by regressing from the learned representation.(Right) The visualized procedure of shaping representation and model attention shifting with learning stages. The representations were visualized after being compressed to one-channel images using principal component analysis; the model attention maps were visualized using FullGrad[START_REF] Srinivas | Full-gradient representation for neural network visualization[END_REF]

Fig. 3 .

 3 Fig. 3. Construction of the training dataset consisting of the intensity calibration of CT (Uemura et al., 2021), bone segmentation[START_REF] Hiasa | Automated muscle segmentation from clinical CT using Bayesian U-Net for personalized musculoskeletal modeling[END_REF], 2D-3D registration to the X-ray image[START_REF] Otake | Intraoperative image-based multiview 2D/3D registration for image-guided orthopaedic surgery: incorporation of fiducialbased C-arm tracking and GPU-acceleration[END_REF] and DRR generation by projecting QCT.

Fig. 4 .

 4 Fig. 4. Tuning process of the threshold used to derive BMD from PF-DRR in stage two. GT-BMD is linearly fitted to the average intensity of the PF region in PF-DRR. PF-DRR is binarized by thresholding to determine the pixel count inside the PF region for calculating the average intensity. The line plot shows the threshold against PCC between GT-BMD and the BMD derived from PF-DRR, where the black points indicate the chosen threshold that maximizes PCC. The label images show the binarized results from a representative case using the chosen thresholds.

Fig

  Fig. 5. T-SNE visualization of X-ray images in dataset B with data collected from 305 patients. Each data point represents an X-ray image. The horizontal and vertical axes are the first and second dimensions in T-SNE space, respectively. The X-ray images were collected with different poses (standing, adduction, abduction, and supine), suggesting a high variability in the dataset and a high potential for various validations.

Fig. 6 .

 6 Fig. 6. Correlation and absolute error between predicted and GT BMDs. Prediction of DXA-BMD. (b) Prediction of QCT-BMD. 1) Conventional reg.is the direct regression method trained and tested using the pre-cropped ROI of the proximal femur of X-ray images. 2) Proposed decomp. is the proposed decomposition method that predicts 2D BMD distribution. 3) Proposed decomp. + reg. is the proposed regression method based on method 2).

  -BMD (mg/cm3) Pose x-ray images from a patient Predicted QCT-BMD (mg/cm3) Poses in training dataset

Fig. 7 .

 7 Fig. 7. Pose-wise correlations between predicted and ground-truth BMD by the proposed method in experiments of training with (upper row) four poses and (lower row) two poses. (a) Prediction of DXA-BMD. (b) Prediction of QCT-BMD. (c) Typical X-ray image of each pose.

Fig. 8 .

 8 Fig.8. Line plots showing the BMD estimation performance by our method affected by the compression of X-ray images with fewer bits and the JPEG compression algorithm when targeting DXA-BMD (left) and QCT-BMD (right). The horizontal axes indicate the average file sizes of images after compression under the size of 512 × 256. The vertical axes show the corresponding correlation scores obtained by our method.

Fig. 9 .Fig. 10 .

 910 Fig. 9. Bland-Altman plots for the DXA-BMD estimation task (left) and the QCT-BMD estimation task (right). Case-level outliers are indicated for whose points of all poses are either above +1.96SD or below -1.96SD.

Fig. 11 .

 11 Fig. 11. Outlier cases in the DXA-BMD estimation task (left) and the QCT-BMD estimation task (right). The PF-DRR and label images are cropped with the same bounding box to align the horizontal view. The red arrows indicate the edge of the femoral head included in the PF region due to hip deformation, potentially introducing errors in BMD estimation.

Fig. 12 .

 12 Fig. 12. Ablation study results of X-ray image decomposition. (a) Evaluation of image decomposition for ResNet (baseline) without HL, HRFormer without HL, and HRFormer with HL. (b) The Convergence analysis of each backbone without and with HL. (c) The ROI of the input X-ray image and GT-DRR, and the comparison of the training progress between (d) without HL and (e) with HL. Without HL in (d), predicted BMD became zero at 150 and 300 due to instability of decomposition in several cases

Fig. 13 .

 13 Fig. 13. Results of BMD estimation using a limited dataset. (a) Correlation of the predicted BMD with DXA-BMD and QCT-BMD. (b) Boxplot of the absolute error (AE) of the predicted BMD. The BMD predicted using the proposed method clearly shows a higher correlation with DXA-measured and QCT-measured BMDs and smaller absolute errors.

Table 1 .

 1 Datasets Cohort

	Item	Dataset A	Dataset B
	Number of X-ray images Number of patients Female, n (%) Mean age (sd), years Mean BMI (sd), kg/m 2 Mean DXA-BMD (sd), g/cm 2 Median T-score (IQR) Osteoporosis w.r.t. DXA-BMD, n (%) n/a 525 275 236 (85.8) 250 (82.0) 1204 305 58.6 (14.3) 57.6 (14.9) 23.1 (3.8) 23.5 (4.0) n/a 0.759 (0.151) n/a -1.2 (-2.1, 0.4) 50 (16.4) Mean CT-BMD (sd), g/cm 3 n/a 290 (57)

Table 2 .

 2 Summary of BMD estimation experiment result

	Pose in training dataset		Method	DXA-BMD Estimation Accuracy PCC↑ ICC↑ MAE↓ SEE↓ RMS-CV↓ PCC↑ ICC↑ MAE↓ SEE↓ RMS-CV↓ QCT-BMD Estimation Accuracy
		Conventional Regression	0.667 0.609	0.091	0.117	6.02%	0.764 0.730	30.3	38.6	5.99%
			Decomposition using QCT 0.864 0.841	0.058	0.077	3.27%	0.861 0.847	23.1	29.3	3.79%
	Standing, Supine, Abduction, Adduction	Proposed	+Regression (freeze M E ) +Regression (tune R) Decomposition using CT +Regression (freeze M E ) +Regression (tune R)	0.874 0.837 0.880 0.848 0.836 0.814 0.864 0.829 0.872 0.838	0.057 0.055 0.062 0.057 0.056	0.076 0.074 0.083 0.078 0.076	3.40% 3.35% 3.53% 3.54% 3.46%	0.914 0.904 0.920 0.909 0.871 0.852 0.915 0.905 0.922 0.910	18.4 17.9 22.2 18.1 17.6	23.6 22.9 28.4 23.4 22.7	3.61% 3.36% 3.78% 3.73% 3.54%
		GT-PF-DRR from QCT	0.930 0.928	0.040	0.055	1.39%	0.883 0.877	20.8	27.0	1.24%
		GT-PF-DRR from CT	0.921 0.918	0.044	0.059	1.49%	0.893 0.888	19.9	25.8	1.36%
			Decomposition using QCT 0.855 0.837	0.059	0.079	3.72%	0.851 0.833	23.8	30.3	3.77%
	Standing, Supine	Proposed	+Regression (freeze M E ) +Regression (tune R) Decomposition using CT	0.862 0.835 0.866 0.840 0.833 0.803	0.058 0.057 0.063	0.078 0.076 0.084	4.02% 3.93% 3.60%	0.909 0.897 0.913 0.900 0.848 0.832	18.9 18.6 24.0	24.3 23.8 30.5	3.99% 3.85% 4.49%
			+Regression (freeze M E ) +Regression (tune R)	0.833 0.808 0.847 0.823	0.063 0.061	0.084 0.081	4.18% 4.28%	0.898 0.888 0.903 0.893	19.9 19.4	25.7 25.0	4.19% 4.12%

Poses in test dataset

					Standing	Supine
					Sta.	Sup.
	PCC=0.882	PCC=0.879	PCC=0.919	PCC=0.922
	ICC=0.850	ICC=0.847	ICC=0.909	ICC=0.909
	MAE=0.054±0.050	MAE=0.055±0.049	MAE=17.7±14.6	MAE=18.0±13.9	Abduction	Adduction
	PCC=0.873	PCC=0.862	PCC=0.916	PCC=0.912	Abd.	Add.
	ICC=0.844	ICC=0.837	ICC=0.905	ICC=0.895
	MAE=0.055±0.051	MAE=0.059±0.050	MAE=18.3±14.6	MAE=18.9±15.0

Table 3 .

 3 Summary of the ablation study results.

	Method	Image Decomposition Accuracy Mean PSNR Mean DC ICC	PCC	Mean AE	DXA-BMD Estimation Accuracy SEE ICC PCC PCC w.r.t. QCT-BMD
	ResNet	30.688	0.658	-0.278 0.006	0.117	0.157 -0.024 -0.208	-0.130
	+ HL	39.105	0.952	0.866 0.894	0.057	0.074 0.872	0.879	0.818
	DAFormer	36.874	0.926	0.538 0.671	0.085	0.114 0.639	0.680	0.632
	+ HL	37.994	0.945	0.853 0.875	0.057	0.078 0.856	0.865	0.799
	HRNet	36.996	0.931	0.369 0.650	0.097	0.127 0.537	0.581	0.640
	+HL	39.971	0.958	0.883 0.920	0.057	0.074 0.870	0.878	0.843
	HRFormer	36.594	0.927	0.255 0.495	0.109	0.143 0.313	0.400	0.498
	+ HL	40.168	0.961	0.910 0.927	0.053	0.071 0.882	0.888	0.853

Acknowledgments

This work was funded by MEXT/JSPS KAKENHI (19H01176, 20H04550, 21K16655).

Code availability

The source code of our method is available from the authors (gu.yi.gu4@is.naist.jp, otake@is.naist.jp, yoshi@is.naist.jp) upon reasonable request for research activity.