This study was performed at Koneru Lakshmaiah University, Bachupally, Hyderabad 500043, Telangana, India, in accordance with applicable institutional ethical guidelines and the principles of the Declaration of Helsinki and its subsequent amendments. Ethical approval was obtained from the Institutional Ethics Committee (IEC), K L Deemed to be University, Bachupally Campus, Hyderabad, Telangana, India (Approval No. KLU/IEC/2026/014). The IEC waived the requirement for informed consent because the study involved retrospective analysis of anonymized MRI data with no personally identifiable information. The lumbar spine MRI data were de-identified before analysis, and all directly identifiable patient information was removed. Digital Imaging and Communications in Medicine (DICOM) headers were stripped of personal identifiers, and no patient names, identification numbers, contact information, facial reconstruction information, or institution-linking metadata were retained. The data were used solely for research purposes, with appropriate measures implemented to maintain patient confidentiality, privacy, and data security throughout model development, validation, testing, and reporting. Details of the software, computing infrastructure, and data resources used in this study are provided in the Table of Materials.
NOTE: Overview of the workflow. The complete pipeline proceeds in six stages, in the following order: (1) data acquisition and ingestion of de-identified lumbar spine MRI into HDFS; (2) preprocessing, including resizing, intensity normalization, and training-only augmentation; (3) convolutional feature extraction with a ResNet-50 backbone; (4) global context modeling and classification with a ViT encoder and a dense classification head; (5) distributed training with Apache Spark over the Hadoop cluster using a weighted cross-entropy objective; and (6) evaluation on the held-out test subset. The section numbering below follows this order, and Figure 1 shows the corresponding system architecture.
Dataset and partitioning
Expert-labeled sagittal T2-weighted lumbar spine MRI slices were assembled and annotated at each disc level according to the Pfirrmann grading system2. For this study, 5,000 slices from 625 patients were assembled, with eight slices obtained from each patient. Each slice carried a single slice-level label corresponding to the most degenerated disc visible in the slice. The ordinal Pfirrmann grades were mapped to a binary target once at the start of the study, and the same mapping was applied to every subset. Grades I and II were recorded as normal, whereas Grades III to V were recorded as degenerated. The mapping was stored in the manifest file, which contained the columns slice_id, patient_id, path, pfirrmann_grade, and label, to prevent alteration between runs.
The data were partitioned at the patient level rather than the slice level so that no patient contributed images to more than one subset. The stratified group-split function of the machine-learning library, StratifiedGroupKFold, was used with patient_id as the group key, label as the stratification key, and random_state set to 42 to produce a 70/15/15 training, validation, and test split while preserving the class proportion within each subset. The resulting dataset comprised 3,500 training slices from 438 patients, 750 validation slices from 94 patients, and 750 test slices from 93 patients. The partition was verified by confirming that the intersection of patient identifiers between any two subsets was empty before any training run was initiated. The pipeline was aborted if the intersection between any two subsets was non-empty. Nested evaluation subsets of 100, 200, and 300 slices were then drawn from the held-out test subset using a fixed random seed of 42 and stratification by class. These subsets were used to report performance at three test-subset sizes, and their class compositions were recorded in Table 1.
The manifest and evaluation-subset lists were stored in HDFS and reused without modification in all later analyses.
Preprocessing
Each slice was resized to 224 × 224 pixels, corresponding to the input resolution of the convolutional backbone, using bilinear interpolation. The single grayscale channel was replicated to three channels to enable the use of ImageNet-initialized weights. The intensity of each slice was normalized by subtracting the slice mean and dividing by the slice standard deviation, thereby placing the inputs on a common scale independent of scanner gain. This normalization stabilized optimization and improved convergence.
Data augmentation was applied only to the training subset. The training slices were subjected to random rotations within ±10°, horizontal flipping with a probability of 0.5, and additive Gaussian noise with a standard deviation of 0.02. No augmentation was applied to the validation or test subsets to avoid distortion of the reported evaluation metrics. The pre-processed slices were serialized in TFRecord format, with one record generated per slice containing the image tensor, slice_id, patient_id, and label, before the records were written to HDFS. This serialization reduced storage overhead and accelerated sequential reads during training.
Data ingestion into HDFS
The serialized dataset was written to HDFS using the command-line client with the command hdfs dfs -put train.tfrecord /ldd/train/. A replication factor of 3 and a block size of 128 MB were used. HDFS replicated each block across the data nodes, providing fault tolerance and enabling parallel read access32. The stored records were read into Spark DataFrames using PySpark with spark.read.format("tfrecord").load("hdfs:///ldd/train/"). The records were then converted to tensors on each worker so that input/output operations were distributed across the workers rather than being concentrated on a single node.
Convolutional feature extraction
A ResNet-50 backbone initialized with ImageNet weights was instantiated, with the classification top removed6. The residual connections enabled the network to be trained at greater depth while reducing the risk of vanishing gradients, and the pretrained weights provided edge and texture filters that could be transferred to MRI images.
Hierarchical feature maps were extracted from the convolutional layers. At layer l, the convolution operation was performed according to Equation 16:
F(l)i,j = σ (Σm=1M Σp=1P Σq=1Q W(l)m,p,q · X(l-1)i+p,j+q + b(l)) (1)
where F represented the feature map at layer l, W represented the filter weight at layer l, X represented the input from the preceding layer, b represented the bias term, and σ represented the rectified linear unit activation. Learnable filters were convolved with the input in a sliding-window manner to produce activation maps that emphasized disc contours, tissue density, and degenerative texture. The rectified linear unit then introduced the non-linearity required to represent these patterns. The output of the final convolutional block was a 7 × 7 × 2,048 feature map.
Vision Transformer encoding and classification
The final convolutional feature map was divided into 2 × 2 non-overlapping patches. Each patch was flattened into a one-dimensional vector, linearly projected into an embedding of 512 dimensions, and combined with a learned positional encoding according to Equation 218:
z(i)0 = E · x(i)p + E(i)pos (2)
where x represented the i-th patch, E represented the linear projection matrix, and Epos represented the positional encoding. This operation converted the feature map into a token sequence that could be processed by the transformer, while the positional encoding preserved the spatial arrangement of the discs that would otherwise have been lost during flattening.
Six transformer encoder layers were then applied, with each layer using eight-head self-attention over the token sequence according to Equation 3, followed by a two-layer feed-forward block with layer normalization and residual connections:
Attention (Q, K, V) = softmax ( QKT / √dk ) V (3)
where Q, K, and V represented the query, key, and value matrices calculated from the token embeddings, respectively, and dk represented the dimension of the key vectors. The attention operation computed the similarity between pairs of patches and re-weighted the value vectors accordingly, allowing the model to relate the appearance of one disc to the alignment and spacing of adjacent lumbar levels.
The transformer output was aggregated by global average pooling into a single feature vector and passed through a dense classification head with a softmax activation according to Equation 4:

where z represented the aggregated output vector of the ViT, Wc and bc represented the weights and bias of the classification layer, respectively, and
represented the predicted class probability. The resulting output represented the probability that a slice belonged to the degenerated class rather than the normal class.
Distributed training and optimization
The training data were sharded across the cluster, and partial gradients were computed independently on each node according to Equation 533:
Map(k, v) → (k', ∇Lk) (5)
where k indexed the data shard, v represented the corresponding input batch, and ∇Lk represented the gradient computed on that shard. The run was launched using spark-submit --num-executors 4 --executor-memory 48g train.py, with one executor assigned to each node and one GPU assigned to each executor. This data-parallel phase allowed each worker to evaluate the identical model on a distinct subset of the training data.
The partial gradients were then aggregated and synchronized according to Equation 633:
Reduce(k', [∇L1 ... ∇LK]) → ∇L = Σk=1K ∇Lk (6)
The worker gradients were summed, and the resulting weight update was broadcast to the workers so that all nodes remained synchronized. This distributed training procedure provided the same gradient update as single-node training with the corresponding effective batch size.
A weighted cross-entropy objective was optimized to compensate for class imbalance according to Equation 7:
L = - Σc=1C wc · yc · log yc , wc = N / (C · Nc) (7)
where C represented the number of classes, yc represented the reference label, yc represented the predicted probability for class c, N represented the total number of training samples, and Nc represented the number of samples in class c. The class weights increased the penalty associated with under-represented classes and reduced the tendency of the classifier to default to the majority class. For the 1,855 degenerated and 1,645 normal training slices, the class weights were 0.943 and 1.064, respectively.
The model was trained using the AdamW optimizer34 with an initial learning rate of 1 × 10⁻4, a weight decay of 1 × 10⁻2, and a batch size of 32 per node, resulting in an effective batch size of 128 across four nodes. Training was performed for a maximum of 50 epochs. The learning rate was reduced by a factor of 0.5 when the validation loss did not improve for four consecutive epochs, and training was stopped when the validation loss did not improve for eight consecutive epochs. The checkpoint with the lowest validation loss was retained. A random seed of 42 was fixed for weight initialization, data shuffling, and augmentation. The same configuration was applied to every model and every training run.
NOTE: The protocol can be paused after Train with the AdamW optimizer step. Checkpoints are written to HDFS at the end of every epoch and training resumes from the latest checkpoint.
Evaluation
The operating threshold was selected using the validation subset only, based on Youden's index according to Equation 8:
J = sensitivity + specificity - 1 (8)
The threshold that jointly maximized the true-positive and true-negative rates was identified using the validation data. The selected threshold for the proposed model was 0.47. The threshold was then fixed and was not adjusted after evaluation of the test subset, thereby preventing optimistic bias in the test estimate.
The trained model was evaluated on each held-out test subset generated in Draw nested evaluation subsets. Confusion-matrix counts, including true positives, false positives, false negatives, and true negatives, were recorded with the degenerated class designated as the positive class. Precision, recall, specificity, F1-score, and accuracy were calculated directly from these counts and were reported in Supplemental File 1. Precision represented the proportion of predicted degenerated discs that were truly degenerated, whereas recall represented the proportion of truly degenerated discs that were detected. The F1-score was calculated as the harmonic mean of precision and recall and was used as a summary measure when both types of classification errors carried clinical cost.
A 95% confidence interval (CI) was calculated for each proportion, including precision, recall, and accuracy, using the Wilson score method. These confidence intervals were used to characterize the sampling uncertainty associated with the test subsets and to facilitate comparison of model performance.
The comparison baseline models, AlexNet and VGG-16, were trained under the same experimental protocol. The same patient-level data split, augmentation strategy, input resolution, optimizer, learning-rate schedule, early-stopping rule, random seed, training budget, and evaluation pipeline were used for all models so that differences in performance could be attributed to the model architectures rather than differences in the experimental setup.
Model export and inference
The checkpoint selected in Train with the AdamW optimizer was exported to HDFS in the SavedModel format and served through the model-serving component behind a REST interface. Single-case inference was performed using the POST /v1/models/ldd:predict endpoint, with a 224 × 224 × 3 tensor provided in the request body.
Experimental environment and scalability measurement
Model development and testing were performed on a cluster of four identical nodes. Each node was equipped with a multi-core server processor, 64 GB of system memory, one data-center GPU with 32 GB of memory, and 2 TB of solid-state storage, and the nodes were connected through a 10 GbE network. The cluster was configured with one node serving as both the HDFS Name Node and Data Node and three nodes serving as Data Nodes. The same operating system, Python version, deep-learning framework, GPU toolkit, Hadoop version, and Spark version were installed and configured on every node.
To assess scalability, the identical model, preprocessing pipeline, hyperparameters, and evaluation procedure were used for training on both a single node and the four-node cluster at dataset sizes of 100, 200, 300, 1,000, and 5,000 images. For each dataset size and configuration, the end-to-end wall-clock training time was recorded from immediately before data loading and initialization until the final training epoch and validation had been completed. The recorded time included initialization, data communication, synchronization, and checkpointing overhead.
Throughput was calculated as the number of images processed per second, and speed-up was calculated as the ratio of the single-node training time to the distributed training time. Parallel efficiency was calculated as the speed-up divided by the number of participating nodes and was expressed as a percentage. The resulting scalability measurements were recorded in Table 2.