Automated Quantitative Analysis of Pulmonary Vasculature in Congenital Diaphragmatic Hernia using Deep Learning

Summary

This protocol describes a fully automated deep learning pipeline for segmenting and analyzing the pulmonary vasculature in neonatal CT scans. The method enables quantitative morphometric evaluation of vascular development in congenital diaphragmatic hernia (CDH) and control subjects, supporting non-invasive characterization of pulmonary hypoplasia and vascular underdevelopment.

Abstract

Congenital diaphragmatic hernia (CDH) is characterized by pulmonary hypoplasia and vascular underdevelopment, leading to impaired gas exchange and high neonatal mortality. Accurate and quantitative assessment of pulmonary vasculature is crucial for understanding disease severity, but manual segmentation of three-dimensional vascular networks in medical images is time-consuming and operator-dependent.

This protocol presents a fully automated deep learning–based method for pulmonary vessel segmentation and morphometric analysis using postnatal computed tomography (CT) scans. The pipeline includes standardized preprocessing steps - conversion to Hounsfield units, windowing, isotropic resampling, and contrast-limited adaptive histogram equalization (CLAHE) - to normalize imaging data and enhance vascular visibility. A U-Net convolutional neural network (CNN) architecture is then trained to segment the pulmonary vasculature, followed by a three-dimensional skeletonization algorithm to quantify morphometric parameters such as branch number, mean branch length, and generational depth.

Representative results demonstrate that the proposed model achieves high segmentation accuracy, with the transfer learning configuration yielding the best performance. Quantitative morphometric analysis reveals markedly reduced vascular complexity in CDH compared with control lungs, consistent with the known pathological features of pulmonary hypoplasia.

This automated approach enables reproducible, quantitative, and non-invasive evaluation of pulmonary vascular morphology in CDH. The method can be adapted to other imaging modalities and applied to studies of fetal and neonatal lung development, facilitating translational research and future clinical integration.

Introduction

Congenital Diaphragmatic Hernia (CDH) is a life-threatening congenital anomaly characterized by a diaphragmatic defect, leading to the herniation of abdominal viscera into thorax^1,2. This physical compression severely impairs lung development, resulting in pulmonary hypoplasia and persistent pulmonary hypertension (PPHN), which are the primary drivers of morbidity and mortality. In addition to pulmonary hypoplasia and vascular underdevelopment, impaired cardiac development and ventricular dysfunction have also been shown to significantly influence clinical outcomes in neonates with CDH. The underlying pathophysiology involves abnormal development of both airways and the pulmonary vascular bed, leading to a reduced number of vessels, increased muscularization of arterioles, and consequently, elevated vascular resistance^3,4. Objective and quantitative biomarkers are needed to accurately stratify risk, guide interventions, and monitor treatment response in CDH patients⁵. One critical aspect of this evaluation is the detailed analysis of the pulmonary vasculature, which can provide insights into the extent of pulmonary hypoplasia and the functional capacity of the lungs. Advances in imaging techniques, particularly computed tomography (CT), have enhanced our ability to visualize and quantify the pulmonary vasculature in great detail^6,7.

While postnatal Computed Tomography (CT) provides high-resolution anatomical detail of the lungs, analysis of the intricate pulmonary vascular tree remains challenging. Existing methods for vascular segmentation often rely on traditional image processing techniques that require significant manual intervention, are susceptible to image artifacts, and may not be robust to the severe anatomical distortions in CDH^7,8,9,10. Deep learning, particularly convolutional neural networks (CNNs) such as the U-Net architecture, has achieved remarkable success in automated medical image segmentation. However, many existing models are trained on healthy subjects or on other disease contexts, limiting their applicability to congenital anomalies such as CDH^10,11,12,13.

Despite these advancements, there remain significant gaps in the literature. Many studies have focused on healthy individuals or specific pulmonary conditions, with limited attention to congenital anomalies like CDH¹². Additionally, while deep learning models have shown improved performance, they often require large, annotated datasets for training, which are not always available for rare conditions such as CDH. Furthermore, existing models have not fully addressed the challenge of distinguishing between different types of pulmonary vessels (e.g., arteries and veins) in the presence of severe anatomical distortions caused by CDH. This limitation underscores the need for further research to develop more robust models that can accurately segment and analyze the pulmonary vasculature in CDH patients.

This study aims to address these gaps by developing and validating a fully automated deep learning framework to segment the pulmonary vasculature and extract quantitative morphometric features from CT scans. A key innovation of our approach is training our model on a combined dataset of CDH and control patients, enabling it to learn a robust representation of both normal and pathological vascular patterns. While CT imaging involves ionizing radiation, making it unsuitable for routine longitudinal screening, this study serves as a crucial proof-of-concept. The primary goal of this study is to establish that automated radiological quantification of vascular structure is feasible and can reliably differentiate CDH patients from controls. Success in this domain provides the necessary validation to adapt this quantitative framework to radiation-free imaging modalities, such as Magnetic Resonance Imaging (MRI), for future clinical applications.

Protocol

All procedures involving human participants were conducted in accordance with the guidelines of the institutional human research ethics committee and were approved by the Institutional Review Board (IRB #2017-6361). The study was performed in compliance with applicable regulatory standards. Patient data were retrospectively collected and de-identified prior to analysis. An overview of the developed system is shown in Figure 1.

1. Image preprocessing

1. Load raw CT scans of neonatal patients in DICOM format into the working environment.
   1. Open the medical image analysis software and create a new project workspace.
   2. Click File > Import > DICOM series.
   3. Navigate to the folder containing the neonatal thoracic CT scan in DICOM format and select the complete image series.
   4. Verify that all slices in the series are correctly detected and ordered based on acquisition metadata (e.g., Instance Number).
   5. Confirm image dimensions, voxel spacing, and slice thickness in the metadata panel to ensure consistency across scans.
   6. Load the series into the workspace and visually inspect axial, coronal, and sagittal views to confirm correct orientation and absence of loading errors.
2. Convert DICOM files to NIfTI format using open-source software to enable easier manipulation and analysis.
   1. Open a terminal window within the scientific computing environment.
   2. Navigate to the directory containing the DICOM image series using the cd command.
   3. Execute the DICOM-to-NIfTI conversion tool using the following command structure: dcm2niix -z y -f output_filename -o /output_directory /input_directory
   4. Ensure that compression is enabled (-z y) to generate a compressed .nii.gz file.
   5. Verify that the output file is successfully generated in the specified output directory.
   6. Open the converted NIfTI file in a medical image viewer and visually inspect axial, coronal, and sagittal planes to confirm correct spatial orientation and image integrity.
   7. Confirm voxel spacing and image dimensions to ensure consistency across all subjects prior to preprocessing.
3. Convert image intensities to Hounsfield Units (HU) using scanner-specific metadata or standard formulas.
   1. Extract the DICOM metadata parameters Rescale Slope and Rescale Intercept from the image header for each scan.
   2. For each voxel intensity value (I_raw), compute the corresponding Hounsfield Unit (HU) using the following formula: HU = (I_raw × Rescale Slope) + Rescale Intercept.
   3. Apply the conversion to the entire 3D image volume using element-wise matrix operations within the scientific computing environment.
   4. Confirm correct conversion by verifying that air regions approximate −1000 HU and soft tissue regions fall within expected physiological ranges.
   5. Save the converted volume as a new NIfTI file to preserve the original data.
4. Apply windowing to highlight the lung and soft tissues. Set window level (WL) to -400 HU and window width (WW) to 1500 HU.
   1. Open the converted NIfTI file in the image visualization software.
   2. Navigate to the image display or intensity settings panel.
   3. Select the Window/Level adjustment option.
   4. Manually set the Window Level (WL) value to -400 Hounsfield Units (HU).
   5. Set the Window Width (WW) value to 1500 HU.
   6. Confirm and apply the settings to update the image visualization.
   7. Verify that lung parenchyma and pulmonary vessels are clearly distinguishable from bone and mediastinal structures before proceeding to further preprocessing steps.
5. Perform isotropic resampling to ensure voxel dimensions are uniform (e.g., 1 mm × 1 mm × 1 mm) using trilinear interpolation.
   1. Load the NIfTI image into the Python environment using a medical image processing library.
   2. Extract the original voxel spacing from the image header metadata.
   3. Define the target isotropic spacing as (1.0, 1.0, 1.0) mm.
   4. Compute the new image dimensions using the formula: new_size = original_size × (original_spacing / target_spacing)
   5. Initialize a resampling object.
   6. Set the interpolation method to trilinear interpolation.
   7. Assign the target voxel spacing (1.0 mm × 1.0 mm × 1.0 mm).
   8. Set the computed new image size.
   9. Preserve the original image direction and origin metadata.
   10. Execute the resampling operation.
   11. Save the resampled image in NIfTI format for subsequent preprocessing steps.
   12. Verify isotropic spacing by checking the updated voxel dimensions in the image header before proceeding.
6. Apply Contrast-Limited Adaptive Histogram Equalization (CLAHE) to enhance contrast and improve the visibility of vascular structures.
   1. Import the required image processing library into the Python environment.
   2. Convert the resampled 3D CT volume to 8-bit grayscale format if necessary, using linear intensity normalization to map the selected HU range to 0–255.
   3. Process the CT volume slice-by-slice in the axial plane to apply CLAHE in 2D.
   4. Initialize the CLAHE object using the following parameters: clipLimit = 2.0, tileGridSize = (8, 8).
   5. For each axial slice, apply the CLAHE function to enhance local contrast.
   6. Reconstruct the processed slices back into a 3D volume after CLAHE application.
   7. Save the contrast-enhanced volume in NIfTI format for subsequent segmentation.
   8. Visually verify that pulmonary vessels are more distinguishable from surrounding parenchyma without excessive noise amplification before proceeding.
7. Visually inspect a subset of preprocessed images to ensure quality and consistency across the dataset.
   1. Randomly select at least 10% of the total dataset for manual quality assessment.
   2. Open each selected preprocessed volume in the medical image viewer.
   3. Inspect axial, coronal, and sagittal planes to verify the following: Correct spatial orientation; Absence of truncation artifacts; Proper application of windowing parameters; Successful isotropic resampling (uniform voxel spacing); Adequate contrast enhancement after CLAHE.
   4. Confirm that pulmonary vessels are clearly distinguishable from adjacent parenchymal structures without excessive noise amplification.
   5. Compare preprocessed images with original HU-converted volumes to ensure that preprocessing steps did not introduce distortion or anatomical inconsistencies.
   6. Document any preprocessing errors and repeat preprocessing for affected cases if necessary.

2. Manual annotation

1. Select a representative subset of CT scans from both control and CDH patients to be used for manual annotation. Ensure a balanced distribution of anatomical variability.
   1. Identify all eligible preprocessed CT scans from both control and CDH cohorts.
   2. Exclude scans with severe motion artifacts or incomplete lung coverage.
   3. Randomly select a predefined number of cases from each group to ensure balanced representation.
   4. Ensure inclusion of cases demonstrating a range of anatomical variability, including differences in lung volume, vascular density, and mediastinal shift severity.
   5. Confirm that selected CDH cases represent varying degrees of pulmonary hypoplasia, where available.
   6. Document the selected case identifiers prior to proceeding with manual annotation.
2. Load the preprocessed NIfTI images into a 3D medical image annotation tool.
   1. Open the 3D medical image annotation software.
   2. Create a new project or segmentation session.
   3. Click File > Open Image (or equivalent import option).
   4. Navigate to the directory containing the preprocessed NIfTI file (.nii or .nii.gz) and select the image.
   5. Confirm successful loading by verifying correct image orientation in axial, coronal, and sagittal views.
   6. Adjust display settings if necessary to optimize visualization.
   7. Create a new segmentation label or mask layer to store manual vascular annotations.
   8. Save the project file before initiating manual annotation.
3. Using axial, coronal, and sagittal views, manually segment the pulmonary vasculature by outlining vascular structures in each relevant slice.
   1. Activate the segmentation label layer created in step 2.2.
   2. Select the manual drawing or brush tool within the annotation software.
   3. Adjust the brush size dynamically according to vessel diameter to ensure accurate boundary tracing.
   4. Using axial slices as the primary reference plane, manually outline visible pulmonary vascular structures, including both arterial and venous branches.
   5. Exclude non-vascular structures such as bronchi, airway walls, and mediastinal tissues.
   6. Scroll slice-by-slice through the entire lung volume to ensure continuous annotation of each vascular branch.
   7. Cross-validate each annotated region in coronal and sagittal views to confirm anatomical consistency and avoid discontinuities.
   8. Include vessels down to the smallest visually distinguishable branches while avoiding over-segmentation of noise artifacts.
   9. Periodically render a 3D preview of the segmentation to verify the spatial continuity of the vascular tree.
   10. Save the completed segmentation mask in NIfTI format before proceeding to model training.
4. Annotate only the pulmonary vessels, excluding the heart, bronchi, and major non-pulmonary structures.
   1. Identify pulmonary arteries and veins within the lung parenchyma using axial slices as the primary reference.
   2. Include intraparenchymal vascular branches originating from the main pulmonary arteries and extending distally within the lung fields.
   3. Exclude the heart chambers, atria, ventricles, and great vessels.
   4. Exclude airway structures, including bronchi and bronchial walls, by distinguishing them from vessels based on morphology and luminal characteristics.
   5. Avoid labeling mediastinal soft tissues, pleura, and chest wall structures.
   6. Use multiplanar views (axial, coronal, sagittal) to confirm that annotated structures follow expected vascular continuity and branching patterns.
   7. When vessel-airway differentiation is uncertain, verify continuity across adjacent slices to confirm vascular trajectory before labeling.
   8. Perform a final 3D rendering of the segmented volume to ensure that only the pulmonary vascular tree has been included.
5. If available, consult clinical imaging experts to validate ambiguous regions during annotation.
   1. Identify regions in which vessel boundaries are uncertain due to low contrast, anatomical distortion, or proximity to bronchi or mediastinal structures.
   2. Flag these regions within the annotation software using a temporary label or comment tool.
   3. Present the flagged regions to a clinical imaging expert.
   4. Review the axial, coronal, and sagittal views together to determine whether the structure represents pulmonary vasculature.
   5. Modify the segmentation mask based on expert consensus.
   6. Document any corrected regions before finalizing the ground-truth mask.
6. Perform inter-rater validation by having at least two independent annotators review and refine each segmentation. In case of disagreement, reach a consensus through discussion or arbitration.
   1. Assign each selected CT scan to two independent annotators with experience in thoracic imaging analysis.
   2. Ensure that annotators perform segmentation independently and are blinded to each other’s results.
   3. After completion of independent annotations, compare the segmentation masks using a quantitative overlap metric (e.g., Dice similarity coefficient).
   4. Identify regions of disagreement by computing voxel-wise differences between masks.
   5. Review discrepant regions jointly in axial, coronal, and sagittal planes.
   6. Reach consensus through structured discussion.
   7. If disagreement persists, involve a third senior reviewer to arbitrate the final decision.
   8. Save the consensus segmentation mask as the final ground-truth label for model training.
7. Save the annotated vessel masks in the same resolution and space as the original CT images. Store them in NIfTI format using consistent naming conventions.
   1. Ensure that the final consensus segmentation mask is stored as a binary label map, where vascular voxels are assigned a value of 1 and background voxels are assigned a value of 0.
   2. Confirm that the segmentation mask retains the same voxel spacing, image dimensions, origin, and orientation matrix as the corresponding preprocessed CT image.
   3. Export the segmentation mask in NIfTI format (.nii or .nii.gz) using the annotation software’s export function.
   4. Use a consistent naming convention structured as follows: SubjectID_Group_VesselMask.nii.gz; (e.g., CDH_012_VesselMask.nii.gz).
   5. Store the masks in a dedicated directory parallel to the image dataset to maintain pairing consistency.
   6. Perform a final verification by reloading both the CT image and its corresponding mask to confirm perfect spatial alignment before model training.
8. Use these manually labeled segmentations as the ground truth for model training and evaluation.

3. Model training and validation

1. Organize the dataset into three subsets: training, validation, and test sets. In this study, use 35 control and 20 CDH cases for training, 5 control and 10 CDH cases for validation, and hold out the remaining for independent testing.
   1. Compile all preprocessed CT images and their corresponding consensus vessel masks into a single dataset.
   2. Perform dataset splitting at the patient level to prevent data leakage between subsets.
   3. Randomly assign 35 control cases and 20 CDH cases to the training set.
   4. Randomly assign 5 control cases and 10 CDH cases to the validation set.
   5. Assign all remaining cases to an independent hold-out test set that is not accessed during model training or hyperparameter tuning.
   6. Ensure that each CT image and its corresponding segmentation mask remain paired throughout the splitting process.
   7. Verify class distribution in each subset to maintain representation of both control and CDH cases.
   8. Document the final allocation of subject identifiers for reproducibility.
2. Normalize the intensity values of all images between 0 and 1 to improve neural network convergence during training.
   1. For each CT volume, first restrict intensity values to a predefined Hounsfield Unit range to remove extreme outliers.
   2. Apply intensity clipping such that values below −1000 HU are set to −1000 HU and values above 500 HU are set to 500 HU.
   3. Perform min–max normalization independently for each volume using the following transformation: Normalized value = (I − I_min) / (I_max − I_min), where I_min and I_max correspond to the clipped minimum and maximum intensity values of the volume.
   4. Apply the normalization to the full 3D volume using element-wise operations within the computational environment.
   5. Verify that all voxel intensities are within the interval [0, 1] before inputting the data into the neural network.
   6. Ensure that normalization parameters are derived independently for each image to prevent information leakage between training, validation, and test sets.
3. Implement a U-Net convolutional neural network architecture using a deep learning framework such as PyTorch or TensorFlow.
   1. Create a new project in a deep learning development environment and set a fixed random seed for reproducibility.
   2. Define a 2D U-Net architecture (Figure 2) for binary segmentation with an encoder-decoder structure and skip connections.
   3. Set the model input as single-channel CT slices (grayscale) and set the model output as a single-channel probability map representing vessel likelihood.
   4. Configure the final layer to use a sigmoid activation function to produce values in the range [0, 1].
   5. Define the loss function for binary segmentation (e.g., binary cross-entropy) and initialize an optimizer.
   6. Specify training hyperparameters, including batch size, number of epochs, and learning rate, and record these values for reporting.
   7. Prepare data loaders to feed the training and validation datasets into the model with consistent shuffling and batching.
   8. Save the complete model definition and hyperparameter configuration for reproducibility.
4. Configure the model with an encoder-decoder structure, incorporating skip connections and batch normalization layers to improve segmentation accuracy.
   1. Define an encoder consisting of repeated convolutional blocks. Ensure each block includes:
      Two consecutive 2D convolution layers (kernel size 3 × 3, padding = 1)
      Batch normalization applied after each convolution
      Rectified Linear Unit (ReLU) activation
   2. Apply 2 × 2 max pooling with a stride of 2 after each encoder block to progressively reduce spatial resolution.
   3. Double the number of feature channels after each downsampling step.
   4. Define the decoder using transposed convolution (2 × 2 kernel, stride 2) for upsampling.
   5. Concatenate feature maps from the corresponding encoder layer to the decoder layer via skip connections to preserve spatial information.
   6. Apply two convolutional layers with batch normalization and ReLU activation after each concatenation step.
   7. Use a final 1 × 1 convolution layer to map features to a single-channel output.
   8. Apply a sigmoid activation function to produce a voxel-wise probability map for vessel segmentation.
5. Train three configurations of the model:
   1. Baseline model (training from scratch)
      1. Initialize the U-Net weights randomly (e.g., He initialization).
      2. Use only the manually annotated artery-vein dataset as input.
      3. Resize input images to 512 × 512 pixels.
      4. Normalize intensity values to the range [0,1].
      5. Set batch size to 8 (or maximum allowed by GPU memory).
      6. Use Adam optimizer with learning rate = 1 × 10⁻⁴.
      7. Use Binary Cross-Entropy (for binary segmentation) or Cross-Entropy loss (for artery-vein classification).
      8. Train for 30 epochs.
      9. Monitor validation loss after each epoch.
      10. Save the model weights corresponding to the lowest validation loss.
      11. In PyTorch, follow steps 3.5.1.13–3.5.1.14.
      12. Define optimizer: torch.optim.Adam(model.parameters(), lr=1e-4)
      13. Define loss: torch.nn.BCEWithLogitsLoss() or torch.nn.CrossEntropyLoss()
      14. Use model.train() during training and model.eval() during validation.
   2. Segmentation-input model (auxiliary vessel mask input)
      1. Modify the input layer to accept two channels:
         Channel 1: CT image
         Channel 2: Binary vessel segmentation mask
      2. Concatenate the CT image and vessel mask along the channel dimension before feeding them into the network.
      3. Keep the architecture identical to the baseline model.
      4. Use the same optimizer, learning rate, batch size, and epoch number as in step 3.5.1.
      5. Save the best-performing weights based on validation F1 score.
   3. Transfer learning model (pretrained initialization)
      1. Load pretrained U-Net weights trained on a generic vessel segmentation dataset (as described in the progress report).
      2. Freeze encoder layers for the first 5 epochs (optional stabilization step).
      3. Unfreeze all layers and continue fine-tuning for the remaining epochs.
      4. Use reduced learning rate = 5 × 10⁻⁵ during fine-tuning.
      5. Train for 30 epochs total.
      6. Apply horizontal flipping as data augmentation during training.
      7. Save the model with the highest validation F1 score.
6. Use binary cross-entropy loss and the Adam optimizer with an initial learning rate of 0.001. Reduce the learning rate adaptively if validation loss plateaus.
   1. Defining the loss function
      1. For binary vessel segmentation, use Binary Cross-Entropy loss with logits.
      2. In PyTorch, define: criterion = torch.nn.BCEWithLogitsLoss()
      3. If performing multi-class artery-vein classification, use: criterion = torch.nn.CrossEntropyLoss()
   2. Defining the optimizer
      1. Use the Adam optimizer with an initial learning rate of 0.001.
      2. In PyTorch: optimizer = torch.optim.Adam(model.parameters(), lr=0.001)
   3. Enabling adaptive learning rate reduction
      1. Implement a learning rate scheduler to reduce the learning rate when validation loss stops improving.
      2. Use ReduceLROnPlateau scheduler.
      3. In PyTorch: scheduler = torch.optim.lr_scheduler.ReduceLROnPlateau(
         optimizer,
         mode='min',
         factor=0.5,
         patience=5,
         verbose=True
         )
      4. After each validation phase, call: scheduler.step(validation_loss)
   4. Training loop configuration.
      1. For each epoch, follow steps 3.6.4.2–3.6.4.7:
      2. Set model to training mode: model.train()
      3. Perform a forward pass.
      4. Compute loss.
      5. Backpropagate: loss.backward()
      6. Update weights: optimizer.step()
      7. Zero gradients before next iteration: optimizer.zero_grad()
      8. After the training phase, switch to evaluation mode using model.eval() and compute validation loss.
   5. Early stopping (optional but recommended)
      1. If validation loss does not improve for 10 consecutive epochs, stop training to prevent overfitting.
7. Apply data augmentation techniques such as random rotations, flips, and elastic deformations to increase robustness and reduce overfitting.
   1. Defining an augmentation pipeline
      1. Apply augmentations only to the training dataset.
      2. Apply identical spatial transformations to both the CT image and its corresponding segmentation mask.
   2. Random rotations.
      1. Apply random in-plane rotations between -15° and +15°.
      2. In PyTorch (using torchvision or Albumentations): RandomRotation(degrees=15)
   3. Horizontal and vertical flips.
      1. Apply a horizontal flip with probability p = 0.5.
      2. Apply a vertical flip with probability p = 0.5 (if anatomically acceptable).
   4. Elastic deformation.
      1. Apply an elastic transformation to simulate anatomical variability.
      2. Use small deformation parameters to avoid unrealistic distortion.
   5. Implementation example (PyTorch + Albumentations)
      1. Define transformation pipeline before training:
         transform = A.Compose([
         A.Rotate(limit=15, p=0.5),
         A.HorizontalFlip(p=0.5),
         A.ElasticTransform(alpha=1, sigma=50, alpha_affine=10, p=0.3)
         ])
      2. Apply transformation inside the dataset __getitem__() method to ensure synchronized transformation of image and mask.
   6. Validation and test sets
      1. Do not apply augmentation to validation or test datasets.
      2. Use only normalization and resizing for these sets.
   7. Quality control.
      1. Visually inspect augmented samples before training to ensure anatomical plausibility.
      2. Verify that masks remain aligned with the transformed images.
8. Train the model for a fixed number of epochs or until convergence, monitoring performance on the validation set after each epoch.
   1. Defining training duration
      1. Set the maximum number of training epochs to 30.
      2. Alternatively, continue training until convergence criteria are met (see Early Stopping below).
   2. Training loop
      1. For each epoch, follow steps 3.8.2.2–3.8.2.10:
      2. Set the model to training mode: model.train()
      3. Iterate over all mini-batches in the training set.
      4. Load batch of CT images and corresponding masks.
      5. Perform a forward pass.
      6. Compute loss using the defined loss function.
      7. Backpropagate gradients: loss.backward()
      8. Update model weights: optimizer.step()
      9. Reset gradients before next batch: optimizer.zero_grad()
      10. Compute average training loss for the epoch.
   3. Validation step (after each epoch)
      1. Switch model to evaluation mode: model.eval()
      2. Disable gradient computation: with torch.no_grad():
      3. Iterate over the validation dataset.
      4. Perform a forward pass.
      5. Compute validation loss.
      6. Compute performance metrics (F1 score, precision, sensitivity, DICE score).
      7. Record validation loss and metrics.
   4. Learning rate adjustment
      1. Update scheduler after validation step: scheduler.step(validation_loss)
   5. Convergence criteria
      1. Stop training if validation loss does not improve for 10 consecutive epochs (early stopping), or if performance metrics plateau.
   6. Model checkpointing
      1. Save model weights whenever the validation F1 score improves.
      2. Retain the best-performing model for final evaluation on the test set.
   7. Logging
      1. Store training and validation loss values per epoch.
      2. Plot learning curves (loss vs. epoch) to verify convergence behavior.
9. Select the best-performing model based on the highest F1 score on the validation dataset.
   1. Validation metric tracking
      1. After each epoch, compute the F1 score on the validation dataset.
      2. Store the F1 score together with the corresponding epoch number.
   2. Model comparison
      1. Compare validation F1 scores across all epochs.
      2. Identify the epoch that achieved the highest validation F1 score.
   3. Model checkpointing
      1. During training, save model weights whenever the validation F1 score improves.
   4. Final model selection
      1. After completion of training, load the weights corresponding to the highest validation F1 score:
         model.load_state_dict(torch.load("best_model.pth"))
   5. Independent testing
      1. Evaluate the selected model only once on the independent test dataset.
      2. Do not use test set performance for model selection.
   6. Reproducibility
      1. Record the selected epoch number and corresponding validation metrics.
      2. Fix random seed to ensure reproducibility.
10. Save the trained model weights and configurations for downstream segmentation tasks.

4. Segmentation and skeletonization

1. Load the trained U-Net model and apply it to the preprocessed CT scans in the independent test set.
2. For each scan, generate a binary segmentation mask of the pulmonary vasculature by applying a threshold to the model’s probabilistic output.
3. Visually inspect the segmentation results to confirm anatomical plausibility, particularly in regions affected by CDH-induced distortion.
4. Convert the segmented 3D binary masks into skeletonized representations using a 3D thinning algorithm implemented in software.
5. Label the vascular skeleton using a breadth-first search algorithm, designating the pulmonary trunk as the root node and assigning generation levels to each branch based on connectivity.
6. Remove small, disconnected components or spurious branches that are likely due to noise or segmentation errors, using a minimum voxel size or branch length threshold.
7. Save the skeletonized structures in 3D mesh or graph-compatible formats (e.g., VTK or SWC) for further morphometric analysis.

5. Morphometric feature extraction

1. Load the skeletonized vascular graphs generated from the segmented CT images.
2. Identify all individual branches by traversing the graph structure between bifurcation points and terminal nodes.
3. Calculate the total number of branches by summing all identified segments within the vascular graph.
4. Compute the length of each branch by summing the Euclidean distances between connected voxels along the skeleton.
5. Determine the mean branch length by averaging the lengths of all branches in each subject’s skeleton.
6. Assign generation levels to each branch starting from the main pulmonary artery as generation 0 and incrementing by one at each bifurcation using breadth-first traversal.
7. Calculate the maximum number of distal generations by identifying the longest path from the root node to any terminal branch.
8. Calculate the maximum number of proximal generations by identifying the longest path from any peripheral branch back to the root.
9. Store all morphometric features in a structured spreadsheet or database format (e.g., CSV or SQL) with subject identifiers, diagnosis group, and extracted metrics.
10. Visually inspect a subset of skeleton graphs and corresponding morphometric features to confirm accuracy and biological plausibility.

6. Classification analysis

1. Import the structured dataset containing morphometric features (e.g., total branch count, mean branch length, distal and proximal generation depth) and corresponding diagnostic labels (CDH or control).
2. Split the dataset into training and testing sets using stratified sampling to maintain the class distribution.
3. Implement Random Forest and Decision Tree classifiers using a standard machine learning library.
4. Configure model parameters such as the number of estimators (e.g., 100 trees for Random Forest) and maximum depth based on cross-validation performance.
5. Train each classifier using the training subset of the data.
6. Evaluate classification performance on the hold-out test set using accuracy as the primary metric.
7. Generate confusion matrices to assess true positive, true negative, false positive, and false negative rates for CDH detection.
8. Compare the performance of Random Forest and Decision Tree classifiers, and select the model with the highest accuracy for reporting representative results.
9. Visualize classification results using bar charts, ROC curves, or decision trees if applicable, and save the outputs for inclusion in figures.
10. Document all model parameters, performance metrics, and any preprocessing applied to the data to ensure reproducibility.

Results

The deep learning model trained with transfer learning achieved the highest segmentation performance among all configurations, with a precision of 0.714, sensitivity of 0.706, and F1 score of 0.672 on the independent test set. The baseline model trained from scratch showed reduced performance (precision: 0.703, sensitivity: 0.589, F1 score: 0.551), while the segmentation-input model performed moderately (F1 score: 0.630).

Following segmentation, morphometric analysis revealed a significant red...

Discussion

Several steps in this protocol are critical to achieving accurate and reproducible results. The image preprocessing stage must be carefully executed, particularly the conversion to Hounsfield Units and the application of windowing (WL = -400 HU, WW = 1500 HU), as these parameters determine the visibility of vascular structures^14,15,16,17. Skeletonization also requires precise parameter tuning t...

Materials

List of materials used in this article

Name	Company	Catalog Number	Comments
Analyze 12.0	AnalyzeDirect	https://analyzedirect.com/	Used for CT image review, editing, and measurements
CT Scan Data (postnatal thoracic scans)	Institutional Archive	N/A	Retrospective dataset of CDH and control neonates
dcm2niix	Open-source (GitHub)	https://github.com/rordenlab/dcm2niix	Converts DICOM to NIfTI format
ITK-SNAP	Open-source (http://www.itksnap.org)	http://www.itksnap.org	3D medical image annotation tool
NetworkX	Open-source	https://networkx.org	Used for graph-based vascular tree analysis
OpenCV (CLAHE function)	Open-source	https://opencv.org	Used for image contrast enhancement
Pandas, Numpy	Open-source	https://pandas.pydata.org, https://numpy.org	Data management and numerical operations
Python 3.8+	Python Software Foundation	https://www.python.org	Programming language for analysis and model development
PyTorch 1.13+	Meta AI	https://pytorch.org	Deep learning framework for U-Net implementation
Scikit-image	Open-source	https://scikit-image.org	Image processing library used for skeletonization
Scikit-learn	Open-source	https://scikit-learn.org	Machine learning library for classification
Ubuntu 20.04 LTS	Canonical	https://ubuntu.com	Operating system used throughout all processing
Workstation with NVIDIA RTX 3090 GPU	NVIDIA	https://www.nvidia.com	Required for training deep learning models