Automatic facial axes standardization of 3D fetal ultrasound images
Antonia Alomar Ricardo Rubio Laura Salort Gerard Albaiges Antoni Payà Gemma Piella Federico Sukno


















A. Alomar et al.



Department of Information and Communications Technologies, Universitat Pompeu Fabra, 122-140 Tànger, Barcelona, Spain Department of Obstetrics and Gynecology, Hospital del Mar, 25-29 Passeig Marítim, Barcelona,Spain Department of Medicine and Life Sciences, Universitat Pompeu Fabra, 88 Doctor Aiguader,Barcelona, Spain Fetal Medicine Unit, Obstetrics Service, Department of Obstetrics, Gynecology and Reproductive Medicine, University Hospital Quirón Dexeus, Barcelona, Spain
Automatic facial axes standardization of 3D fetal ultrasound images
    Antonia Alomar1 Ricardo Rubio2,3Laura Salort 1 Gerard Albaiges4 Antoni Payà2,3Gemma Piella1Federico Sukno1
    September 9, 2024
==============================================================================================================





§ ABSTRACT 
Craniofacial anomalies indicate early developmental disturbances and are usually linked to many genetic syndromes. Early diagnosis is critical, yet ultrasound (US) examinations often fail to identify these features. This study presents an AI-driven tool to assist clinicians in standardizing fetal facial axes/planes in 3D US, reducing sonographer workload and facilitating the facial evaluation. Our network, structured into three blocks—feature extractor, rotation and translation regression, and spatial transformer—processes three orthogonal 2D slices to estimate the necessary transformations for standardizing the facial planes in the 3D US. These transformations are applied to the original 3D US using a differentiable module (the spatial transformer block), yielding a standardized 3D US and the corresponding 2D facial standard planes. The dataset used consists of 1180 fetal facial 3D US images acquired between weeks 20 and 35 of gestation. Results show that our network considerably reduces inter-observer rotation variability in the test set, with a mean geodesic angle difference of 14.12^∘ ± 18.27^∘ and an Euclidean angle error of 7.45^∘ ± 14.88^∘. These findings demonstrate the network's ability to effectively standardize facial axes, crucial for consistent fetal facial assessments. In conclusion, the proposed network demonstrates potential for improving the consistency and accuracy of fetal facial assessments in clinical settings, facilitating early evaluation of craniofacial anomalies.






§ INTRODUCTION

Craniofacial anomalies serve as indicators of developmental disturbances at early stages of life, encompassing a wide range of heterogeneous conditions associated with many genetic syndromes <cit.>. Estimates suggest that up to 40% of genetic syndromes produce alterations in the normal morphology of the face and the head. Although these associations have predominantly been identified in adult populations, there is increasing interest in early assessment <cit.>. Consequently, diagnostic efforts are moving towards prenatal and postnatal stages <cit.>. 

To evaluate the fetal development, 2D ultrasound (US) imaging is the standard procedure. Unfortunately, dysmorphology features are hard to identify in this way, due to the noisy nature of fetal US (low signal-to-noise ratio, fetal or probe movements, fetal position, and limbs in front of the face) <cit.>. Currently, 3D/4D US serves as a complement to 2D US. They prove to be particularly useful in diagnosing various fetal anomalies, especially those involving facial abnormalities, neural tube defects, and skeletal anomalies <cit.>. 

In this context, acquiring an US standard plane (SP) is crucial for performing an accurate fetal diagnosis, as the SP is used to measure and analyse biomarkers and abnormal features <cit.>. 3D US has the advantage of capturing multi-view planes allowing sonographers to manually select SPs from 3D US images or videos during prenatal exams. While this process is essential, it is also time-consuming and observer-dependent. This manual selection can be laborious and biased due to the extensive search space, the sonographer experience, and the variability of the fetus orientation <cit.>. 

In this study, we aim to reduce the sonographer's workload while enhancing the accuracy and interpretability of fetal facial SP detection. We propose an AI-driven tool designed to assist clinicians in standardizing the facial axes/planes in the 3D US. It aims to minimize variability across planes detection while mitigating the effects of clinician subjectivity in selecting accurate SPs for fetal facial assessment. Standardizing the fetal facial axes intends to facilitate the evaluation of facial biomarkers and biometric measurements to perform facial assessment. The proposed method consists in regressing the transformation necessary to standardize the sagittal, coronal and axial fetal facial axes, taking as input 3 orthogonal planes centered at the middle of the 3D US image. The novelty lies in that instead of combining the regression model with another task, such as the classification of the planes, we add a differentiable block that incorporates the image loss between estimated and ground truth (GT) planes as part of the minimization strategy. This helps the network learn the structures that should be present in the SPs. Additionally, the proposed algorithm offers the advantage of low computational cost and easy integration into in the echographer or in 3D Slicer as a built-in feature. 






§ RELATED WORK

Several methods have been proposed for automatically detecting 2D SPs in US images using deep learning. Usually, this task has been approached as an image classification problem, using convolutional neural networks (CNNs) or recurrent neural networks (RNNs) <cit.>. However, these methods only determine whether the acquired 2D slices are SPs, but do not inform on what correction shall be applied to them in case they are not SPs.

Another common strategy consists of regressing the plane parameters or transformation matrices to achieve the SP in the 3D US volume. For example, Feng et al.<cit.> proposed a constrained marginal space learning method that combines both 2D and 3D information for fetal face detection in 3D US. Nie et al. <cit.> introduced a deep belief network combined with a detection algorithm to provide a prior structural knowledge to the network. Li et al. <cit.> presented an iterative transformation network to detect SPs in 3D fetal US using a CNN that performs both plane classification and regression to estimate the transformation parameters. Di Vece et al. <cit.> improved the previous results obtained estimating the six-dimensional pose of arbitrary oriented US plane of the fetal brain with respect to a template normalized frame using a CNN regression network. Recently, reinforcement learning (RL) has shown great potential in addressing SP localization as a regression task <cit.>. Although RL approaches have achieved high performance, several issues remain to be addressed, such as the reliance of current studies on initial registration to ensure data orientation consistency, which can easily fail if the pre-registration process is unsuccessful. Moreover, unlike parameters regression models using a CNN, RL models simplify the problem of regressing the transformation parameters by considering a discrete action space and, in consequence, limiting the transformation that can be applied. 
To avoid dependence on pre-registration and ensure no limitations on the transformations that can be performed, we choose a classical yet effective parameters regression approach using a CNN. To help the network learn the structures that should be present in the SPs, we add a differentiable block that incorporates the image loss between estimated and GT planes as part of the minimization strategy. As a result, the number of parameters of the network is significantly reduced because no classification blocks are used.







§ METHOD





  §.§.§ Data & Pre-processing:

The dataset used consists of 1180 fetal facial 3D US images acquired between week 20 and 35 of gestation (26.56 ± 2.72) using a Voluson E8 RSA (BT-20) with a convex probe (4D-RAB6-D, 2–8 MHz) at two hospitals in Barcelona (Hospital de Mar and Hospital Universitari Dexeus) according to their Ethical Research Committee and the current legislation (Organic Law 15/1999). The study population comprises subjects from low-risk pregnancies, meaning without any pathology, or known family cases of craniofacial or syndromic pathologies, which were all carried to term. The data is divided into training, validation, and test sets, with 72%, 12%, and 16% of the data allocated to each set, respectively. To facilitate the use of deep learning, it is essential to standardize the input image size across the training, validation, and test sets. We implemented this through two steps: 1) down-sampling the 3D US image by a factor of two, to reduce the computational cost of the network; 2) symmetric zero-padding to the center of the 3D US to achieve a size of U ∈ℝ^C× H × W × D where H,W,D = 256 are the height, width, and depth dimension and C= 1. The latter is performed to ensure that no information is cropped-out during rotation and translation. The 2D initial planes are defined by I_0 = [I_s, I_c, I_a] with I_s = U(1,H/2, : , :), I_c = U(1,;, W/2 , :) and I_a = U(1,;, ; , D/2), corresponding to the sagittal, coronal, and axial planes, respectively.





  §.§.§ Ground Truth Standard Planes:

The facial GT SPs we are interested in locating are the axial, coronal and sagittal planes that define the canonical axes of the fetal face. They are obtained by minimizing the 3 orthogonal planes defined by 23 anatomical landmarks located by expert clinicians in the 3D US (see Appendix Fig. 1) and following the recommendations from the international 3D focus group <cit.>. We constrain the planes' normal vectors n_a ,n_c, and n_s to be orthonormal. The center of the planes c = (c_x, c_y, c_z)^T is defined as the intersection point of the 3 planes. The GT is obtained using custom code in 3D Slicer. The extracted normal vectors of the 3 orthogonal planes are used to compute the rotation matrix R_gt∈ℝ^3×3 needed to standardize the image axes to the estimated facial SPs. The GT rotation can be written as the change-of-basis matrix S_𝔹_1 →𝔹_2 from 𝔹_1 →𝔹_2. In our case, 𝔹_1 coordinates correspond to the canonical basis and 𝔹_2 coordinates are the estimated normal vectors. Thus,

    R_gt= S_𝔹_1 →𝔹_2  =  [ n_s^T; n_c^T; n_a^T ] = [ r_11 r_12 r_13; r_21 r_22 r_23; r_31 r_32 r_33;      ]

The rotation regression is performed in terms of quaternion representation for compactness, i.e., q_gt = (q_0, q_1, q_2, q_3)^T. The conversion from rotation matrix to quaternion follows



    q_gt = [ q_0; q_1; q_2; q_3;     ]
     =[ 1/2√(1 + r_11 - r_22 - r_33);            r_12 + r_21/4 q_0;            r_13 + r_31/4 q_0;            r_23 - r_32/4 q_0;                              ]

The intersection/center of the planes is the GT translation t_gt = c∈ℝ^3. Then, the GT transformation matrix is 



    θ_gt = [ R_gt ,   t_gt;        ] =[ r_11 r_12 r_13  t_x; r_21 r_22 r_23  t_y; r_31 r_32 r_33  t_z;      ]∈ℝ^3×4
 
To ensure compatibility with the spatial transformer block, the translation is expressed in image relative size. Thus, each component of t_gt is in the range [-1,1]. The 2D GT sagittal, coronal and axial SPs are obtained as I_gt = [I_s, I_c, I_a] where I_s = V_gt(1,H/2, : , :), I_c =  V_gt(1,;, W/2 , :) and I_a =  V_gt(1,;, ; , D/2), and V_gt is the transformed US using θ_gt.




  §.§.§ Feature Extractor Block:

The inputs of the feature extractor block (I_0 ∈ℝ^H× W × 3) are the 3 orthogonal planes located at the center of the 3D US image (sagittal, coronal and axial plane). Each branch uses the AG-SonoNet <cit.> as the backbone feature extractor, with the weights shared among the three branches. Then, each view has a specialized aggregation block that adds the attention information from multiple layers of the network to extract specialized features from each view. This information is concatenated and fed to the translation and rotation regression block. 




  §.§.§ Translation & Rotation Regression Block:
 It consists of two fully connected layers that convert the extracted features from the three orthogonal planes into the translation and rotation necessary to achieve the standardized axes/planes. The output is the regression vector z∈ℝ^7. The first three positions correspond to the translation vector t_es∈ℝ^3 where t_es= (t_x,t_y,t_z)^T with each component being in the range [-h_max,h_max]. The remaining 4 positions correspond to the quaternion representation of the rotation matrix q_es= (q_0,q_1,q_2,q_3)^T with each component being in the range [-1,1]. To represent a valid rotation, ||q_es|| needs to be 1. To ensure that this condition is satisfied, a normalization layer was added after the last fully connected layer. Given q_es, the estimated rotation can be expressed as:

    R_es = [ 1 - 2(q_2^2 + q_3^2)   2(q_1q_2 - q_0q_3)   2(q_1q_3 + q_0q_2);   2(q_1q_2 + q_0q_3) 1 - 2(q_1^2 + q_3^2)   2(q_2q_3 - q_0q_1);   2(q_1q_3 - q_0q_2)   2(q_2q_3 + q_0q_1) 1 - 2(q_1^2 + q_2^2) ] = [ r_11 r_12 r_13; r_21 r_22 r_23; r_31 r_32 r_33;      ]





  §.§.§ Spatial Transformer Block:
 It is a differentiable module capable of applying spatial transformations to the original 3D US image U ∈ R^H× W × D × C, resulting in a new standardized 3D US image V ∈ R^H × W × D × C and facial 2D SPs. The spatial transformer uses a differentiable 3D bi-linear sampling as defined in <cit.>. Each output value for pixel i can be written as

    V_i = ∑^H _n ∑^W_m ∑^D_l U^c_n,m,lmax(0,1-|x^s_i -m|) max(0,1-|y^s_i -n|) max(0,1-|z^s_i -l|)

where (x_i^inp, y_i^inp, z_i^inp) are the input coordinates that define the sampling points in the original 3D US U and U_nml is the value of U at location (n,m,l). The output coordinates (x_i^out, y_i^out, z_i^out) are defined to lie on a regular grid G = G_i of pixels G_i = (x^out_i , y^out_i , z^out_i) and can be obtained by a 3D affine transformation: 




    [ x^inp_i; y^inp_i; z^inp_i;         ] = 𝒯_θ(G) = θ(q,t) [ x^out_i; y^out_i; z^out_i;       1;         ] =[ r_11 r_12 r_13  t_x; r_21 r_22 r_23  t_y; r_31 r_32 r_33  t_z;      ][ x^out_i; y^out_i; z^out_i;       1;         ]

θ (q,t) is the 3D transformation matrix estimated. We use height, weight and depth normalized coordinates, such that x_i,y_i,z_i ∈ [-1,1].  The 2D estimated SPs are obtained  as I_es = [I_s, I_c, I_a] where I_s = V(1,H/2, : , :), I_c = V(1,;, W/2 , :) and I_a = V(1,;, ; , D/2).






  §.§.§ Cumulative Transformations & Initialization:

At initialization time (it=0), R_0 is a random rotation defined by the Euclidean angles α_x, α_y, α_z ∈ [-20,20] degrees and t_0  =(t_x,t_y,t_z)^T is random translation with each component being in the range [-0.05, 0.05]. To preserve image quality, we accumulate the transformations and perform a unique transformation to the input 3D image. If multiple steps of the network are applied, we define the transformations at step it as R_es^it = R_es^it-1 R_es and t_es^it= t_es^it-1 + t_es, whereas R_gt^it = R_es^-1 R_gt^it-1 and t_gt^it = - (R_es^it-1)^-1t_es + (R_es^it-1)^-1t_gt. Here, R_es, t_es correspond to the rotation and translation estimated by the CNN at the current iteration, whereas R_es^it, t_es^it denote the accumulated rotation and translation at iteration it. We found that 3 iterations are enough to improve performance by refining the SP estimates. However, there was no need to train the network in an iterative way.




  §.§.§ Network Loss:
 It
is defined as a combination of the mean absolute error (MAE) between the GT and estimated translation, the relative angle (SO3) between the GT and estimated rotation and the image loss between computed as the Frobenius norm of the difference between the GT and estimated SPs:





    ℒ =   β ||t_es - t_gt||_1 + γ   acos(0.5 * (Tr(R_es (R_gt)^T)-1))+ ||I_gt - I_es||_1

where Tr correspond to the trace, β and γ are the translation and rotation weights in the loss, I_gt and I_es are the 2D GT and estimated slices corresponding to the fetal facial sagittal, coronal and axial SPs.





§ EXPERIMENTS


The proposed method is compared to the inter-observer variability obtained from 3 different observers placing landmarks on the 3D US images as described in González-Aranceta et al. <cit.>, which are then used for estimating the planes as described in Section <ref>. Next, we compare to the state of the art method proposed by Li et al. <cit.>. For a fair comparison, the network used is their M1 baseline model with the addition of the differential spatial transformer model to be able to train with the image loss. The only task learned is the regression of t,qusing as input 3 orthogonal planes. The predicted SPs/axes are evaluated in the test set against the GT using the distance between the GT and the estimated translation, and the rotation angles between the GT and the estimated rotation. The image similarity of the planes is also measured using the peak to noise-ratio (PSNR) and structural similarity index (SSIM).



§ RESULTS & DISCUSSION






Table <ref> summarizes the results obtained. The proposed method outperforms the state of the art method from Li et al. <cit.> and even challenges inter-observer variability, producing smaller angular errors although larger translation errors.






Table <ref> shows the translation and rotation Euclidean angles error obtained per SP/axes. It highlights that the coronal and sagittal planes are more challenging to locate in terms of rotation than the axial plane. Our approach reduces the inter-observer rotation error of the sagittal, axial and coronal planes. The translation error obtained is around 6mm per plane/axes, higher than the inter-observer error. Fig. <ref> shows some qualitative examples and comparisons between the GT 2D planes and the estimated by the proposed method. Despite a translation error per plane of approximately 6 mm in patient 2 and 3, the estimated plane closely approximates the GT plane.



Thus, the proposed network correctly learns the plane angles but its translation error is higher than the inter-observer error. This could be due to the high variability defining the planes localization, as multiple slices could closely resemble each other <cit.>. However, it could also be that the most informative planes for the network are not the ones defined as GT planes and it sacrifices translation accuracy for rotation accuracy. Although the translation errors obtained with the proposed network are larger than the inter-observer, the mean translation error per plane is around 6 mm. Moreover, rotation can be regarded more important than translation, as the aim is to standardize the facial axes. This allows the sonographer to examine the standardized fetal 3D facial US, where the fetal pose/US probe rotations are removed, facilitating the evaluation of the fetal face. Furthermore, the proposed network is able to reduce the rotation inter-observer variability. This can be highly beneficial to homogenize the facial analysis evaluation criteria and to reduce the reliance on clinician expertise while reducing the time burden of manually locating the planes.  




§ CONCLUSIONS

We propose a network that estimates the transformation necessary to obtain the standard sagittal, coronal and axial facial US axes taking as input three 2D orthogonal planes. Evaluation on 184 US volumes shows that the network correctly standardizes the US 3D axes while reducing the rotation variability across observers. The method has the potential to be applied easily in a clinical setting due to its low computational cost. The standardization of the facial 3D US aims to facilitate the analysis of the facial biometric measurements to assess the presence of craniofacial abnormalities.
 

 

 



  §.§.§ 
 This work was partly supported by grants PID2020-114083GB-I00 and PRE2021-097544 funded by MICIU/AEI/10.13039/501100011033/ and under the ICREA Academia programme.


  §.§.§ 

None of the authors have any competing interests.







splncs03




