Identification of eye diseases with small fundus image datasets using a new unsupervised deep learning approach

1 Introduction

Diabetic retinopathy (DR), age-related macular degeneration (AMD), glaucoma (GL), and hypertensive retinopathy are major retinal diseases that contribute significantly to visual impairment and blindness worldwide (Bilal et al., 2021). DR is also defined by chronic hyperglycemia progression that results in microvascular damage, which occurs in the form of microaneurysms, hemorrhages, hard exudates, and neovascularization in proliferative stages (Tsiknakis et al., 2021). AMD is a process that is characterized by the gradual degeneration of the macula in either geographic atrophy (dry AMD) or choroidal neovascularization (wet AMD) and causes permanent loss of central vision (Treder et al., 2018). GL refers to a set of optic neuropathies that have traits such as progressive retinal ganglion cell loss and optic nerve head cupping (Arrieta et al., 2025). These are mainly related to high intraocular pressure. These ailments are often comorbid in the geriatric population, making their diagnosis difficult and challenging to adopt multi-disease screening mechanisms over a complex diagnostic framework (Yang et al., 2025). The WHO has reported that retinal diseases account for a large proportion of the global population living with visual impairments (World Health Organization, 2019). According to recent epidemiological research, DR has an approximate prevalence of 103.12 million adults worldwide (Dinsa et al., 2025). AMD and GL have prevalences of 196 million and 76 million patients, respectively (Fleckenstein et al., 2024; Allison et al., 2020). These diseases are expected to grow substantially to 160.5 million (DR), 288 million (AMD), and 111.8 (GL) million by 2040 because of ageing populations and the reported alarming rise in diabetes. Early screening with retinal images can prevent vision loss. However, time-consuming inter-observer variability applies to manual diagnosis (Tsiknakis et al., 2021).

Supervised deep learning models have revolutionized the field of automated retinal analysis by achieving state-of-the-art performance in both classification and segmentation tasks (Khan et al., 2025; Kim and Tran, 2021; Rajarajeshwari and Selvi, 2024; Almustafa et al., 2023). Such models are based on different architectural paradigms, such as convolutional neural networks (CNNs), vision transformers (ViTs), and hybrid architectures to obtain hierarchical multi-scale features in fundus photographs and optical coherence tomography (OCT) images (Takahashi et al., 2024). Complex architectures like the U-Net variants of segmentation, ResNet-based models of classification, and attention-based refinements of features have exhibited impressive accuracy in identifying small quantities of pathological changes at multiple layers of the retinal component (Asadpour and Xie, 2025). Recent models such as MCG-Net (Lin et al., 2021) and MCGS-Net combine convolutional features with graphs or self-supervised modules to capture disease co-occurrence. However, these methods rely heavily on large-scale annotated datasets, which create challenges in medical domains owing to expert labeling costs and ethical constraints (Rani et al., 2024).

Unsupervised learning (Chen et al., 2022) is a viable alternative to supervised approaches. It allows models to learn meaningful visual representations from raw images without relying on labeled data. Unsupervised methods cover a wide variety of paradigms, such as clustering methods (k-means, hierarchical clustering, and DBSCAN), dimensionality reduction methods (PCA, t-SNE, and UMAP) (Fan et al., 2025; Rolf et al., 2025), as well as generative models (VAEs, GANs, and normalizing flow) and self-supervised learning models (Ren et al., 2025; Vafa et al., 2024) that have shown promise in uncovering latent disease patterns on unlabeled medical images. This is particularly valuable for sparsely annotated datasets in medical imaging, which are also subject to variability. Unsupervised methods are suitable for anomaly detection, clustering, and pretraining models without labels (Bengio et al., 2013; Chen et al., 2020). In the context of retinal image analysis, unsupervised models can automatically group images based on underlying pathological features (Jeyasri and Karthiyayini, 2025), image enhancement (Yu and Dai, 2024), and patient stratification (Li et al., 2024). They can discover types of diseases and their subgroups that are not recognized by usual supervised learning methods (Rani et al., 2024). The unsupervised approach of self-supervised learning is an important subgroup of unsupervised approaches. It has received significant academic interest largely due to the use of pretext tasks like rotation prediction (Yan et al., 2023), solving jigsaw puzzles (Konstantakos et al., 2025), masked image modeling (MAE) (Konstantakos et al., 2025), and contrastive learning models such as SimCLR and MoCo. These models enable the learning of robust and transferable representations from unlabeled data that can be fine-tuned for downstream diagnostic tasks. This is because the contrastive representation learning (Bengio et al., 2013) allows the learning of discriminative features without labels, and the models can perform even better than traditional pretraining.

Additionally, normalizing flow models, such as Glow (George et al., 2021), offer a promising direction for invertible latent modeling. They are useful in medical image enhancement and anomaly modeling because of their invertibility and suitable density estimation (Li et al., 2023). By applying flows with contrastive and clustering objectives, we can model retinal images that display different disease patterns. Therefore, we introduce an unsupervised framework that uses the advantages of structure-aware learning to enhance retinal images and identify common traits in retinal diseases. We evaluated our method on two fundus image datasets: STARE ¹ (multi-label) and a subset of IDRiD ² (binary class).

To the best of our knowledge, NF-HAE-UDL is the first framework to jointly combine a hybrid CNN-Transformer autoencoder, a normalizing flow module, and label-guided contrastive learning for retinal disease clustering on small-scale fundus datasets (N < 250 images). Unlike prior supervised methods that require thousands of annotated images and unlike simple unsupervised methods that operate on raw pixel features, NF-HAE-UDL learns structured and disease-discriminative latent representations only using image-level labels to guide contrastive pair construction during training. Three aspects define the core novelty of this work: (1) a dual-branch encoder that simultaneously captures complementary local texture and global structural features that single-branch models cannot provide; (2) the normalizing flow transforms the latent space into a structured Gaussian distribution by enabling invertible latent manipulation and explicit likelihood estimation that are unavailable in standard autoencoders; (3) the complete framework is validated on genuinely small medical imaging datasets where existing deep learning methods fail to generalize, demonstrating consistent performance across ten independent trials with controlled variance. The specific technical contributions of this study are as follows.

• Dual-branch hybrid encoder: a dual-branch encoder architecture is proposed that combines ResNet-50 for local convolutional texture features and Swin-Tiny for global window-attention structural context. The two branches are fused via a learned projection head into a unified latent vector z 0 ∈ R d, where d = 128 for STARE and d = 256 for IDRiD simultaneously capture complementary retinal cues such as microaneurysms, vessel morphology, and optic disc positioning. The detailed procedure is formalized in Section 3.1.
• Glow-based normalizing flow module: a K = 4-layer GlowFlow module is integrated to refine the fused latent vector through invertible ActNorm, linear, and learned-scale transformations (Section 3). This provides a structured Gaussian regularization of the latent space, enabling exact invertibility for latent manipulation and explicit density estimation for downstream anomaly detection, which are the benefits unavailable in standard variational or deterministic autoencoders.
• Label-guided contrastive training with multi-objective loss: a supervised contrastive loss constructs positive pairs using image-level disease labels, thus encouraging disease-specific cluster separation in the learned latent space. The total training objective (Section 3) combines reconstruction loss L rec, contrastive loss L con, and classification loss L cls in a unified multi-objective framework without requiring pixel-level annotations.
• Post-training clustering and evaluation: a post-training pipeline (Section 3) applies UMAP dimensionality reduction followed by DBSCAN density-based clustering on the learned latent codes, with ground-truth labels introduced only after training completes solely for computing evaluation metrics (ARI, NMI, and Silhouette). This guarantees clean separation between representation learning and evaluation.
• Post hoc Grad-CAM interpretability: a post hoc Grad-CAM protocol (Section 5) uses a lightweight linear probe trained on frozen latent codes to cluster pseudo-labels, enabling gradient-based saliency visualization without modifying the trained encoder. Expert-validated IoU scores ( 0.52 ± 0.11 ) confirm clinically meaningful attention to pathological regions.

The use of invertible networks is particularly valuable in medical imaging, where interpretability, reconstruction fidelity, and feature disentanglement are essential. Unlike standard autoencoders, flows enable exact likelihood estimation and invertible mapping (George et al., 2021). Their invertibility ensures traceability, reconstruction fidelity, feature disentanglement, and clinical relevance, particularly in ophthalmic applications (Li et al., 2023). Furthermore, integrating flow-based modeling with contrastive objectives promotes compact, structured embeddings that maintain inter-class separability and intra-class cohesion without labels (Chen et al., 2020; Bengio et al., 2013). This synergy is particularly beneficial for modeling complex retinal phenotypes and comorbidities, enabling unsupervised disease stratification and image enhancement within a unified framework.

This paper is organized as follows: Section 2 contains the technical background, an overview of contrastive representation learning, hybrid feature encoding, latent flow modeling, and the rationale of its integration. Section 3 presents the proposed approach, including hybrid feature extraction, normalizing flow architecture, image reconstruction, the unsupervised loss function, and clustering-based quality assessment. Section 4 discusses the experimental setup, including the description of the datasets, the parameter settings, the metrics, and the clustering results, such as the visualization of the latent features, the clustering performance, and the comparison results with and without normalizing flow as well as with the baseline methods. A discussion on model validation is presented in Section 5. Lastly, Section 6 concludes with a summary of contributions and future research directions.

2 Technical foundations

This section outlines the key components of our proposed method and our unsupervised retinal image enhancement elements that form the foundation of our framework.

2.1 Contrastive representation learning

To learn discriminative features in the absence of labels, we employed contrastive learning—specifically, InfoNCE loss (Oord et al., 2018). This loss encourages representations of augmented views of the same image (positives) to be similar while pushing apart representations of different images (negatives). This approach is particularly suitable for medical imaging where labeled data are scarce, and visual similarity serves as a critical cue for learning meaningful representations, such as retinal analysis. Given a batch of N images, we generate two augmentations per image, yielding 2 N views. Let z i , z j ∈ R d be the L 2-normalized latent representation of a positive pair ( i , j ) extracted from the encoder. The contrastive loss for a single positive pair ( i , j ) follows the NT-Xent formulation (Chen et al., 2020):In Equation 1, sim ( z a , z b ) = z a ⊤ z b / ( ‖ z a ‖ ‖ z b ‖ ) denotes cosine similarity, τ is a temperature hyperparameter controlling the concentration of the distribution, and ⊮ [ k ≠ i ] excludes only the anchor sample from the denominator. The denominator sums over all 2 N − 1 samples, excluding the anchor i, which includes the positive sample j and all 2 N − 2 negative samples. This formulation ensures proper normalization across all contrastive pairs in the batch.

L contrast z i , z j = − log exp sim z i , z j / τ ∑ k = 1 2 N ⊮ k ≠ i ⁡ exp sim z i , z k / τ . (1)

This loss was computed for all positive pairs in the batch and then averaged. Since contrastive learning uses similarity in visuals for grouping the data, it is very useful when annotations are scarce in medical imagery.

2.2 Hybrid feature encoding

To capture both local patterns and global context, we adopted a dual-branch encoder architecture that comprised a ResNet convolutional encoder and a Vision Transformer variant (Dosovit et al., 2020). Complementary local and global features were extracted using these convolutional and transformer encoders and fused into a unified latent vector z 0, where z 0 represents the fused latent vector combining local and global features. The detailed fusion mechanism, which leverages a learned projection combining ResNet and Swin Transformer features, is described in Section 3.1.

Let x denote an input image. The convolutional encoder f c ( ⋅ ) outputs a local feature vector f c ∈ R d c, while the transformer encoder f t ( ⋅ ) produces a global feature vector f t ∈ R d t. As shown in Equation 2, we extract complementary feature vectors from both encoders and project them into a unified latent space via a learned linear projection head:

f c = f c x ∈ R d c , f t = f t x ∈ R d t . (2)

The final hybrid latent representation is obtained by concatenating these vectors and applying a learned projection as in Equation 3:

z 0 = R e L U W e ⋅ f c ⊕ f t + b e ∈ R d . (3)

Full implementation details, including the specific backbone dimensions and projection specifications, are provided in Section 3.1. This fused embedding captures both detailed textures (e.g., microaneurysms, vessels) and broader structural cues (e.g., optic disc positioning), which are essential for fine-grained retinal analysis.

2.3 Latent flow modeling

To enable invertible image enhancement and facilitate the probabilistic interpretation of the latent space, we modeled z 0 using a normalizing flow.

Let f θ : R d → R d denote an invertible function parameterized by θ, implemented using a Glow-like architecture (Kingma et al., 2018). The flow transforms z 0 into a latent variable z = f θ ( z 0 ) such that the likelihood of z 0 under the model is

log ⁡ p z 0 = log ⁡ p z + log det ∂ f θ ∂ z 0 . (4)

In Equation 4, p ( z ) is a tractable prior (e.g., standard normal), and the Jacobian determinant accounts for the volume change under the transformation.

During training, we maximize this log-likelihood jointly with the contrastive loss. This dual-objective formulation enables the model to learn informative, discriminative, and generative latent representations that can be inverted for image enhancement or generation. In our framework, the contrastive loss is computed on the latent codes after refinement by normalizing flow to enhance discriminability in the transformed latent space.

2.4 Rationale for integrating normalizing flow

The integration of normalizing flow (NF) provides explicit likelihood-based regularization to the learning objective. Rather than solely maximizing contrastive similarity, NF encourages latent codes to follow a structured Gaussian distribution. An NF is a type of transformation that maps a complex distribution to a simpler probability distribution. Basically, NF is used to build invertible mapping f : X → Z. f maps the unknown distribution p ( x ) of the feature space to the latent space with a Gaussian distribution q ( z ). Then, for each sample point x in the distribution p ( x ), there will be a sample point z corresponding to it in q ( z ). Based on the change of variables formula, the probability distribution is defined as in Equation 5 (Papoulis and Pillai, 2001):where f denotes the mapping from x to z, J f refers the Jacobian matrix of the f operation, and d e t ( ⋅ ) denotes the determinant of the matrix.

p x = q z ⋅ | d e t ∂ z ∂ x | = q z ⋅ | d e t J f | , (5)

Multiple coupling layers are stacked in NF to improve the fitting ability of the network. The probability distribution of K coupling layers is as follows:

p 1 x = q z ⋅ | d e t J f 1 | p 2 x = q z ⋅ | d e t J f 1 | ⋅ | d e t J f 2 | ⋅ ⋅ ⋅ p K x = q z ⋅ | d e t J f 1 | ⋅ | d e t J f 2 | ⋯ | d e t J f K | , z = f 1 ⋯ f K x . (6)

Log-likelihood is used because the latent output z is modeled to fit a Gaussian (normal) distribution, so the log simplifies calculation. Based on Equation 6, we can estimate the log-likelihood of the image features by q ( z ). Thus, we denote q ( z ) as the probability density function, and the log-likelihood is defined then as in Equation 7 shown below:where D is the dimension of z.

log p K x = log ⁡ q z + ∑ k = 1 K log | d e t J f K | = log 1 2 π D / 2 exp − 1 2 ‖ f 1 ⋯ f K x ‖ 2 + ∑ k = 1 K log | d e t J f K | = − D 2 log 2 π − 1 2 ‖ f 1 ⋯ f K x ‖ 2 + ∑ k = 1 K log | d e t J f K | , (7)

During training, f maximizes the log-likelihood of the encoded features p ( z 0 ), providing an additional regularization signal that complements the contrastive loss. This incorporates two learning objective formulations, leading to modest but consistent improvements in latent space structure (see Section 4.2.5; Table 1 for empirical validation). Thus, we simplify the log-likelihood and formulate the learning objective as the minimization of the negative log-likelihood:

L f = − log p K x = 1 2 ‖ f 1 ⋯ f K x ‖ 2 − ∑ k = 1 K log | d e t J f K | . (8)

TABLE 1. Test set negative log-likelihood for learned flow models.

As seen in Equation 8, the minimization of the log likelihood function with normal distribution thus enables learning a deep model with small-sized data because the latent data z transformed from the training data then tend to fall into a small region on the support of p ( z ).

3 Proposed method: NF-HAE-UDL (normalizing flow–guided hybrid autoencoder-based unsupervised deep learning)

Building on the foundations introduced in Section 2, we propose a normalizing flow-guided hybrid autoencoder (NF-HAE) for unsupervised retinal image enhancement and clustering. As shown in Figure 1, the framework combines hybrid feature encoding, normalizing flow-based latent transformation, and clustering-guided reconstruction. This modular design enables both perceptual enhancement and semantic grouping in the absence of labels.

Let X ⊂ R ( C × H × W ) denote the space of input images and x ∼ P X. Our goal is to learn an invertible mapping F : X → Z and its inverse F − 1:

x ̂ = F − 1 G F E x ⏟ Normalizing Flow - Guided Transformation , F = G ◦ F E , (9)

where in Equation 9 F E is the hybrid encoder extracting complementary features, G is the normalizing flow module refining latent codes, and x ̂ is the reconstructed, enhanced image. This composition mirrors the architecture described in Section 2, where F E performs hybrid feature encoding and G represents the Glow-based normalizing flow.

3.1 Hybrid feature extraction

To capture both local texture and global semantic features, we extract complementary representations using a CNN and a transformer and fuse them into a latent vector z 0, as shown in Equation 10:

f res = ϕ CNN x = AvgPool ResNet - 50 x ∈ R 2048 , f trans = ϕ Trans x = CLS Swin - Tiny x ∈ R 768 , z 0 = ReLU W e ⋅ f res ⊕ f trans + b e ∈ R d . (10)

We use ResNet-50 (He et al., 2016) and Swin (Liu et al., 2021) as dual-bran-Tiny ch encoders, both initialized with ImageNet-1K pretrained weights. The ResNet-50 extracts local texture features via convolutional layers, outputting a 2048-dimensional feature vector from the final average pooling layer. The Swin-Tiny captures global context through hierarchical window-based self-attention, producing a 768-dimensional CLS token embedding. These complementary features are concatenated ( [ f res ⊕ f trans ] ∈ R 2816 ) and projected to a unified latent space via a learned linear transformation W e ∈ R ( d × 2816 ), where d = 128 for STARE and d = 256 for IDRiD. ReLU activation is applied after projection to introduce non-linearity. The complete step-by-step procedure for hybrid feature extraction, including input specifications and complexity analysis, is formalized in Algorithm 1.

Algorithm 1 Hybrid feature extraction.

3.2 Normalizing flow (glow)

We refine the fused latent vector z 0 using a Glow model (Kingma et al., 2018) composed of invertible transformations arranged in a multi-scale architecture with three levels and 32 flow steps per level. To prevent catastrophic forgetting of pretrained features on small datasets, we employ a two-stage protocol. (1) Warmup phase (first 50 epochs): encoder backbones are frozen; only the projection head ( W e , b e ) is trained with a learning rate of 1 × 1 0 − 4. (2) Fine-tuning phase (remaining epochs): encoders are unfrozen with a reduced learning rate of 1 × 1 0 − 5, while the projection head is maintained with 1 × 1 0 − 4. This stabilizes the training and mitigates overfitting on small-scale fundus datasets. Each flow step f k applies an affine coupling transformation, as shown in Equation 11:where ⊙ denotes element-wise multiplication, split ( ⋅ ) stands for the split operator, and CouplingNet is a three-layer MLP with hidden dimensions [ d , 2 d , d ] and Leaky ReLU ( α = 0.1 ) activations. The log-determinant of the Jacobian is efficiently computed as log | det ( J f k ) | = ∑ log ⁡ s. As in Equation 12, the complete transformation iswhere K = 96 in the full Glow architecture (3 levels × 32 steps) and K = 4 in the implemented simplified variant (Algorithm 2).

z a , z b = split z k − 1 , z a ′ = z a log ⁡ s , t = CouplingNet z a , z b ′ = z b ⊙ exp log ⁡ s + t z k = z a ′ , z b ′ , (11)

z = G z 0 = f K ◦ f K − 1 ◦ … ◦ f 1 z 0 , (12)

This architecture guarantees exact invertibility ( z 0 = G − 1 ( z ) ) and enables tractable likelihood evaluation. Note: the above formulation with K = 96 represents the full theoretical Glow architecture. In experiments, a computationally efficient simplified variant with K = 4 layers is used, as formalized in Algorithm 2, which retains the key properties of invertibility and likelihood estimation while reducing parameter count for small datasets:where p ( z ) = N ( 0 , I ) is a standard Gaussian prior.

log ⁡ p z 0 = log ⁡ p z + ∑ k = 1 K log | det J f k | = − 1 2 ‖ z ‖ 2 2 − d 2 log 2 π + ∑ k = 1 K log | det J f k | , (13)

For computational efficiency on small-scale datasets, the implemented GlowFlow module uses a simplified K = 4 layer architecture with ActNorm, affine and learned-scale transformations, as formalized in Algorithm 2. The full multi-scale Glow formulation described above serves as the theoretical foundation; the simplified variant retains approximate likelihood estimation while reducing parameter count. Note that the GELU activation in Algorithm 2 breaks exact invertibility. This simplified variant is used solely for regularization and density estimation, and exact invertibility is available only in the full K = 96 Glow formulation.

Algorithm 2 Normalizing flow forward pass (GlowFlow) simplified k = 4, layer implementation used in experiment.

3.3 Image reconstruction

To reconstruct enhanced images, we map the refined latent z back to image space using a progressive upsampling decoder with four transposed convolutional blocks. The process begins with a linear projection and a reshape (see Equation 14):where W up ∈ R ( 512 × 49 × d ) projects the latent code to a spatial feature map. Each upsampling block T i doubles spatial resolution via transposed convolution, as shown below by Equation 15:

h ̂ 0 = reshape W up ⋅ z + b up ∈ R 512 × 7 × 7 , (14)

T i h = ReLU BN ConvTranspose 2 d h ; stride = 2 , kernel = 4 . (15)

The complete decoder pipeline is

h ̂ 1 = T 1 h ̂ 0 ∈ R 256×14×14 h ̂ 2 = T 2 h ̂ 1 ∈ R 128×28×28 h ̂ 3 = T 3 h ̂ 2 ∈ R 64 × 56 × 56 h ̂ 4 = T 4 h ̂ 3 ∈ R 32 × 112 × 112 x ̂ = σ Conv 2 d Upsample h ̂ 4 , scale = 2 ; out_channels = 3 , kernel = 3 ∈ R 3 × 224 × 224 . (16)

In Equation 16, σ is the sigmoid activation constraining output pixels to [ 0,1 ]. All input images are resized to ( 224 × 224 ) during preprocessing.

3.4 Unsupervised loss function

Our total loss balances the reconstruction fidelity with latent discriminability, as formalized in Equation 17. We employ contrastive loss as defined in Section 2, computed over the transformed latent space z.where L NLL = − log ⁡ p ( z 0 ) is the negative log-likelihood from Equation 13, encouraging latent codes to follow a structured Gaussian distribution. We set λ rec = 1.0, λ con = 1.0, λ NLL = 0.1, and λ cls = 1.0. The classification loss L cls is computed using BCEWithLogitsLoss between predicted and true disease labels; it is used to guide contrastive pair construction during training as in Equation 18:where P ( i ) denotes the set of positive pairs for sample i (typically | P ( i ) | = 1 for single augmented views), z ̃ i = z i / ‖ z i ‖ 2 represents L 2-normalized latent codes, and 1 [ k ≠ i ] is an indicator function that equals 1 when k ≠ i and 0 otherwise. The denominator thus sums over all 2 N − 1 samples in the augmented batch except the anchor i itself, explicitly including the positive sample j. The denominator ∑ k = 1 2 N ⊮ [ k ≠ i ] sums over all samples in the batch except the anchor i. This includes the positive sample j, following the standard NT-Xent formulation, where the denominator represents the partition function over all possible pairs. This differs from some alternative formulations where the denominator sums only over negatives. Our formulation matches SimCLR (Chen et al., 2020) and ensures proper probabilistic normalization.

L total = λ rec ⋅ ‖ x − x ̂ ‖ 2 2 + λ con ⋅ L contrast z + λ NLL ⋅ L NLL z 0 + λ cls ⋅ L cls z , (17)

L contrast = − ∑ i = 1 N 1 | P i | ∑ j ∈ P i log exp z ̃ i ⊤ z ̃ j τ ∑ k = 1 2 N ⊮ k ≠ i ⁡ exp z ̃ i ⊤ z ̃ k τ , (18)

In our implementation, we do not introduce an additional projection head g ( ⋅ ). Instead, we directly use the flow-refined latent representation z = G ( z 0 ) ∈ R d as the embedding for contrastive learning, where z 0 is the raw encoder output and G is the normalizing flow module. Prior to computing similarities, embeddings are ℓ 2 normalized to unit length as seen in Equation 19:and pairwise similarities are computed via cosine similarity (i.e., the dot product of normalized vectors), scaled by a temperature parameter τ. This normalization is applied immediately before constructing the similarity matrix implemented within the contrastive loss function, and no additional linear or nonlinear projection is used.

z ̃ = z ‖ z ‖ 2 (19)

3.5 Clustering-based quality assessment

We assessed the unsupervised enhancement quality by clustering the latent codes z using k-means or DBSCAN (Equation 20). In addition, DBSCAN does not require selection of the number of clusters a priori. We keep ε minPts fixed across all trials and do not tune them using validation or test labels. Clusters are reported exactly as produced by DBSCAN (including any noise points when present). Evaluation metrics are computed afterward solely for analytical purposes:

C z → Cluster Labels . (20)

We quantified the latent structure using the Silhouette ( S ) and Calinski–Harabasz ( C ) scores. Effective enhancement improves these metrics, as expressed in Equation 21:

M x ̂ > M x , where M ∈ S , C . (21)

All aspects of the training loop, including hybrid encoding, flow refinement, decoding, and classification, are formalized in Algorithm 3. The post-training clustering pipeline, which applies UMAP followed by DBSCAN on the learned latent codes without any label access, is detailed in Algorithm 4. Although the core representation learning is label-guided, clustering and all subsequent evaluation metrics are computed entirely without label access, ensuring a clean separation between training and evaluation. Images that have been reconstructed are only used to compute reconstruction loss and are not recycled into further training samples or reintroduced into evaluation. Clustering and classification measures are based solely on latent representations and model outputs, not on the reconstructions themselves. Although the main idea is unsupervised, we outline how to use labels in semi-supervised situations by introducing a classifier loss into the training procedure. Images that have been reconstructed are only used to compute reconstruction loss or error and are not recycled into further training samples or reintroduced into evaluation. The measures of clustering and classification are based on latent representations and model outputs but not on the reconstructions themselves.

Algorithm 3 Full training: NF-HAE-UDL.

Algorithm 4 Post-training clustering via UMAP and DBSCAN.

4 Experiment

4.1 Experimental setup

4.1.1 Datasets

We evaluate our method on two fundus image datasets: STARE (STARE dataset, 2025) (multi-label) and a subset of IDRiD (Porwal et al., 2018) (binary). The datasets are annotated in very different ways, the class labels have different levels of detail, and the patients included have many different conditions, all of which enables us to evaluate the method in various ways.

The STARE dataset contains 81 color fundus images, each annotated with one or more retinal disease labels of 11 possible categories (including “Normal Retina”). This dataset consists of 32 images of normal retinas and 49 images showing some pathology. Of a total ∼400 images, we considered the difficult pathological images for our resource-limited settings based on their optic discs (ODs), which are quite tricky to locate (Sattar et al., 2014) due to highly diseased retinas. Each STARE image is associated with a multi-hot binary vector of length 11, where multiple entries can be 1 simultaneously for concurrent conditions. For example, an image labeled “Background DR + Arteriosclerotic Retinopathy” has a vector [0,1,1,0,0,0,0,0,0,0,0]. Of the total images, 32 are labeled as normal (only index 0 active), while the remaining 49 images span diverse pathologies, with 15% exhibiting multiple concurrent diagnoses. These 49 span a diverse range of retinal pathologies, such as DR, vein/artery occlusions, and rare conditions like Coat’s disease and emboli. Since several images show various conditions together, STARE is beneficial for training models that are capable of distinguishing different diseases in multi-label learning, class imbalance, and fine-resolution separation. Sample annotations from STARE are shown in Table 2, with label definitions provided in Table 3 (left panel).

TABLE 2. Samples from multiclass multi-label STARE dataset (total: 81 images; normal: 32).

TABLE 3. Label definitions for STARE dataset (left) and sample entries from the IDRiD binary subset (right).

In contrast, the IDRiD dataset is focused solely on DR diagnosis, offering a binary classification task: “0” for normal and “1” for DR. From the original IDRiD dataset with 516 images (Porwal et al., 2018), we select a representative subset (n = 211) that includes both normal and DR diseased images from both training and testing datasets, retaining the original class distribution between the normal and DR images where all DR images are selected, including all DR grades to avoid the selection bias. Another 305 images referring to another disease (macular edema) are excluded for the sake of binary classification. Although limited in size for this experiment, IDRiD provides high-resolution, clinically validated images centered on DR progression; it is valuable for benchmarking performance in standard binary classification tasks. The subset used in our evaluation is illustrated in Table 3, whereas the whole subset can be seen in the Supplementary Appendix. The left panel shows STARE disease label definitions; the right panel shows IDRiD subset sample annotations. The distributions of the corresponding data samples with respect to two classes (normal (“0” label) and DR (“1” label)) and different DR grades (“1” (mild), “2” (moderate), “3” (severe), and “4” (proliferative DR)) are further depicted in Figures 2a,b, respectively. Our criteria to select the 211 images can be the inclusion for DR images of all grades (1, 2, 3, and 4), while their number decreases as the grade increases. We also make the class distribution imbalance more challenging. It can be noted that the images are patient-independent across splits.

The STARE dataset enables evaluation under a multi-label setting with overlapping and rare retinal pathologies, reflecting complex diagnostic scenarios. STARE is treated as a multi-label task. Each image is associated with a multi-hot vector of concurrent diagnoses, and the model predicts an 11-dimensional sigmoid probability vector. In contrast, the IDRiD subset provides a binary classification task focused on DR to be a representative of typical clinical screening workflows. As a result, the model’s performance can be measured in relation to the number of labels and complexity of tasks.

4.1.2 Parameter settings

Training strategy (label-guided contrastive framework): training only uses image-level disease labels to guide contrastive pair construction and the classification loss L cls. Ground-truth labels are never used for clustering assignments or evaluation metrics during training. To prevent catastrophic forgetting of pretrained features on small datasets, we adopt a two-stage training protocol: (1) warm-up phase (first 50 epochs), during which the encoder backbones are frozen and only the projection head ( W e , b e ) is trained using a learning rate of 1 × 1 0 − 4; (2) fine-tuning phase (remaining 200 epochs), in which the encoders are unfrozen and optimized with a reduced learning rate of 1 × 1 0 − 5, while the projection head retains a learning rate of 1 × 1 0 − 4. Ground-truth labels are never used to compute any training loss. The model is trained using the reconstruction loss L rec, contrastive similarity loss L con, classification loss L cls, and flow-based negative log-likelihood L NLL, as defined in Equation 17.

Data augmentation: following standard unsupervised contrastive learning protocols (He et al., 2020), positive pairs are generated by applying data augmentation to the same image. For each input image, two augmented views are created using the following transformations: (1) random horizontal flip with probability p = 0.5, (2) random rotation within ± 15 °, (3) color jitter with brightness ( ± 0.2 ), contrast ( ± 0.2 ), saturation ( ± 0.2 ), and hue ( ± 0.1 ), and (4) Gaussian blur with a ( 3 × 3 ) kernel and σ ∈ [ 0.1 , 2.0 ]. These augmentations are carefully tuned to preserve clinically relevant retinal features while providing sufficient variation for contrastive learning. Extreme transformations, such as large rotations ( > 20 ° ), elastic distortions, or aggressive color shifts, are deliberately excluded, as they may alter or obscure diagnostically relevant structures, including microaneurysms, hemorrhages, exudates, and vascular morphology. All images are resized to ( 224 × 224 ) during preprocessing.

Optimization: for all experiments, we set the batch size to 8 and used the Adam optimizer. The temperature parameter τ for the contrastive loss was set to 0.07. For the normalizing flow, we used a Glow architecture with three levels and 32 steps per level. The loss weights were empirically chosen as λ rec = 1.0, λ con = 1.0, λ NLL = 0.1, and λ cls = 1.0. All models were trained for 250 epochs (50 warmup + 200 fine-tuning) based on reconstruction error. The number of flow stacks was selected as four for STARE and one for IDRiD, giving the best clustering performance. Detailed architectural specifications and hyperparameters are provided in Table 4.

TABLE 4. Detailed architectural specifications and hyperparameters.

ReLU, rectified linear unit; BN, batch normalization.

Evaluation protocol: all results are reported over ten independent trials, each corresponding to a complete model retraining, with a different random seed affecting parameter initialization and minibatch ordering. Unless otherwise explicitly stated, all model hyperparameters and dataset splits remain fixed across trials. Ground-truth labels are never used during training or representation learning and are introduced only after training to compute evaluation metrics. Reported results correspond to the mean ± standard deviation across trials.

Data leakage prevention: we do not use validation or test labels as criteria for tuning clustering or model parameters. The cluster labels are not used in any way, since the unsupervised loss functions do not require ground-truth class information to derive all the training objectives. The clustering process, consisting of UMAP followed by DBSCAN, is employed only for monitoring and analytical purposes. Hyperparameters of DBSCAN (e.g., ε and minPts) and UMAP settings are fixed in advance and do not depend on the experimental run; they are not optimized on any held-out data. Ground-truth labels are introduced only after training is complete, solely to compute evaluation metrics such as the adjusted Rand index (ARI), normalized mutual information (NMI), and multi-label classification scores. These labels are never used to modify the model weights.

4.1.3 Computational cost and hardware requirements

All experiments were conducted on a single NVIDIA RTX 3090 GPU (24 GB VRAM) with an Intel Core i9-10900K CPU and 64GB RAM running PyTorch 2.0 on Ubuntu 22.04. Table 5 summarizes the training time and memory requirements for each dataset and model configuration.

TABLE 5. Computational cost and GPU memory requirements per dataset. Training time is reported as wall-clock time per epoch (mean over 250 epochs) and total training time for a single trial. Peak GPU memory is measured during the forward pass with batch size 8. Inference time is per image (single forward pass-through encoder and flow module only excluding the decoder).

The total parameter count of ∼58M is dominated by the two frozen pretrained backbones (∼52M), whereas only ∼6.7–6.9M parameters are actively trained on the retinal datasets, significantly reducing the risk of overfitting relative to the nominal model size. The GlowFlow module (K = 4 for STARE; K = 1 for IDRiD) introduces a training time overhead of less than 10% relative to the no-NF ablation baseline, consistent with our earlier claim in Section 4.2.5. Inference operates in real time (<25 ms per image), which is suitable for clinical screening workflows. All ten independent trials have been run sequentially on the same GPU, while the total compute time across all trials, are approximately 21 h (STARE) and 58 h (IDRiD).

4.1.4 Performance metrics

4.1.4.1 Adjusted Rand index (ARI)

The ARI measures the similarity between two partitions by considering all pairs of samples and counting pairs assigned consistently in both true and predicted labels. It is adjusted for chance and defined as in Equation 22.where n i j is the number of samples in the intersection of the ith true class and jth predicted cluster, a i = ∑ j n i j and b j = ∑ i n i j.

A R I = ∑ i j n i j 2 − ∑ i a i 2 ∑ j b j 2 n 2 1 2 ∑ i a i 2 + ∑ j b j 2 − ∑ i a i 2 ∑ j b j 2 n 2 , (22)

4.1.4.2 Normalized mutual information (NMI)

NMI quantifies the mutual dependence between predicted clusters and true labels, normalized to ensure values between 0 and 1. It is given by Equation 23.where I ( Y ; Y ̂ ) is the mutual information between ground truth labels Y and predicted clusters Y ̂, and H ( ⋅ ) denotes entropy.

N M I = 2 I Y ; Y ̂ H Y + H Y ̂ , (23)

4.1.4.3 Silhouette score

The Silhouette score evaluates cluster cohesion and separation for each sample. For a sample x, let a ( x ) denote the mean intra-cluster distance and b ( x ) the smallest mean distance to any other cluster. The Silhouette score for a dataset is:In Equation 24, a higher silhouette value indicates more compact and well-separated clusters.

S i l h o u e t t e = 1 n ∑ x b x − a x max a x , b x . (24)

4.2 Results

The evaluation and comparison results of the method are presented below.

4.2.1 Evaluation results through latent feature representation

A 2D pattern is generated for the latent features in terms of the reconstruction error and distortion, which is calculated according to Jin et al. (2014) using sampling entropy histograms. In Figures 3a,b, we show the discriminant ability of our latent features through nonoverlapping groups referring to various prominent classes for the STARE and IDRiD datasets.

4.2.2 Evaluation results through clustering metrics

The clustering scores (ARI (Rodriguez et al., 2019), NMI (Saitta et al., 2007), and Silhouette (Silhouettes, 1987)) over number of epochs are shown in Figure 4.

As in Figure 4, the Silhouette score plot has a few dips caused by the detection of data points as noise due to a mismatch between ground truth labels and predicted values.

The performance of the proposed method in terms of unsupervised clustering quality across ten independent trials is summarized in Table 6. These metrics assess how effectively similar disease patterns are grouped by the learned latent space without using ground-truth labels during training. We report the mean ± standard deviation of the adjusted Rand index (ARI, measuring post hoc agreement with true labels), normalized mutual information (NMI, measuring shared information between cluster assignments and labels), and the Silhouette score (measuring cluster compactness) computed independently for each trial.

TABLE 6. DBSCAN clustering scores (ARI, NMI, and Silhouette) over ten trials on the STARE and IDRiD datasets. Clustering is performed on the learned latent codes without using ground-truth labels during training. ARI and NMI metrics compare clusters to the true labels post hoc to assess clustering quality, while the Silhouette score evaluates internal cluster structure. Mean ± standard deviation values are reported across ten independent random seed trials. See Section 4.2.4 for 95% confidence intervals.

For the STARE dataset, the method achieves ARI, NMI, and Silhouette scores of 0.58 ± 0.10, 0.70 ± 0.06, and 0.79 ± 0.07, respectively. For the IDRiD dataset, the corresponding scores are 0.58 ± 0.06, 0.63 ± 0.05, and 0.66 ± 0.01. Results are reported for each trial.

The consistency of the scores of clustering on ten trials with random initializations (standard deviations of 0.10, 0.06, and 0.07 of ARI, NMI, and Silhouette on STARE) implies that the model does not learn localized latent representations but instead learns generalized latent representations. In the case of very high overfitting, high variance would be observed in cross-trial and very poor separation of unobserved samples.

4.2.3 Likelihood evaluation

To verify that the normalizing flow successfully transforms latent codes to structured distributions, we report the negative log-likelihood (in nats per dimension) in Table 1.

Negative values indicate that the learned distribution assigns a higher probability density to the trained samples than in a standard Gaussian (log ⁡ p ( z ) = − d 2 log ( 2 π ) − 1 2 ‖ z ‖ 2 ≈ − 1.42 nats/dim for normalized data). All NLL values are computed per dimension by dividing the total negative log-likelihood by the latent dimension d, ensuring comparability across datasets with different d values (d = 128 for STARE, d = 256 for IDRiD). These results confirm successful latent structure learning and demonstrate that the normalizing flow effectively regularizes the latent space beyond simple contrastive alignment. The impact of this structured regularization on clustering performance is further analyzed in the ablation study (Section 4.2.5), where we show that the NF module provides consistent improvements in both clustering metrics and density estimation.

4.2.4 Significance testing

The 95% CIs (confidence intervals) of the clustering results for significance testing (Huang et al., 2015), are shown here.

STARE dataset: ARI→ CI: [0.37, 0.67]; NMI→ CI: [0.55, 0.71]; Silhouette→ CI: [0.73, 0.79]

IDRiD dataset: ARI→ CI: [0.46, 0.60]; NMI→ CI: [0.58, 0.71]; Silhouette→ CI: [0.61, 0.68]

Figure 5 illustrates the 95% CIs varying with epochs on various clustering scores for the STARE dataset.

For the plots in Figure 5, we calculate the mean-difference CI (MDCI) by Equation 25 as:where U i and L i are the upper and lower limits of CI for the ith trial and N is the total number of trials. The MDCI values are found as 0.3, 0.15, and 0.05 with ARI, NMI, and Silhouette scores, respectively, for the STARE dataset. Similarly, with the IDRiD dataset, the MDCI values are obtained as 0.14, 0.13, and 0.06 for the ARI, NMI, and Silhouette scores, respectively. All MDCI values are ≪ 1, indicating the robustness of the model.

MDCI = ∑ i = 1 N U i − L i N , (25)

4.2.5 Ablation study

In this ablation study, we examine the impact of the normalizing flow (NF) module on clustering performance. The results with and without NF (Table 7), show consistent but modest improvements (+ 0.03 ARI, + 0.06 NMI, + 0.07 Silhouette on STARE; + 0.05 ARI, + 0.02 NMI, and + 0.10 Silhouette on IDRiD). Although these gains are limited in magnitude, they are statistically significant given the confidence intervals, and they remain stable across multiple independent trials. Beyond clustering metrics, the NF module offers additional benefits, including explicit density estimation that enables anomaly detection via likelihood thresholding (test NLL: − 2.43 ± 0.15 nats/dim on STARE and − 2.18 ± 0.12 on IDRiD), invertibility that facilitates latent space manipulation (e.g., disease severity interpolation) without decoder-induced artifacts, and structured regularization through a Gaussian prior that promotes smooth and interpretable latent geometry. While a simpler autoencoder may suffice for clustering-only applications, incorporating the flow module enhances interpretability and controllability with minimal computational overhead, as a three-level Glow architecture increases training time by less than 10 %.

TABLE 7. DBSCAN clustering performance (mean ± std) with and without the normalizing flow (NF) module.

Furthermore, the NF is replaced by a VAE (variational autoencoder) to find the clustering scores for the complex STARE dataset, obtained as: ARI = 0.41, NMI = 0.55, and Silhouette = 0.87. Then, the overall clustering performance with a VAE is found to be lower than an NF, supporting the choice of the normalizing flow over the simpler latent regularization techniques.

4.2.6 Comparison results

The performances of the traditional baseline methods and our proposed method—for example, using the STARE and IDRiD datasets—are listed in Table 8, showing high performance compared to other baseline methods. The classical algorithms (k-means, GMM, and hierarchical clustering) have moderate performance. Usually, they do not reflect the pathological subtleties of small-scale fundus images. DEC enhances the performance both in joint feature learning and in clustering, yet it is sensitive to the limited data. DBSCAN has the best ARI, NMI, and Silhouette scores because it has the advantage of density-based clustering. Nevertheless, even DBSCAN cannot make use of complex latent structures to the fullest of their ability, which in turn leads to the given hybrid AE with normalizing flow, which produces much more expressive representations and yields improved clustering results on either of the datasets.

TABLE 8. Comparative clustering results on the STARE and IDRiD datasets evaluated using ARI, NMI, and Silhouette scores. Traditional clustering baselines (spectral clustering through DEC) operate on the same hybrid encoder latent embeddings ( z 0 ) as the proposed method, isolating the effect of normalizing flow and contrastive learning. Self-supervised baselines (DINO, BYOL, SimSiam, and MAE) use frozen features from their respective ImageNet-pretrained backbones without any domain-specific fine-tuning; the same UMAP + DBSCAN pipeline is applied to their representations for fair evaluation. All baseline hyperparameters are set to standard defaults without label-based tuning.

Traditional clustering baselines use the same hybrid encoder latent embeddings ( z 0 ), ensuring that observed differences reflect the clustering strategy rather than feature extraction. Self-supervised baselines (DINO, BYOL, SimSiam, and MAE) use their own pretrained backbone features; the same post hoc UMAP + DBSCAN pipeline is applied uniformly across all methods. NF-HAE-UDL applies DBSCAN on its flow-refined latent codes z; the ARI equality with standalone DBSCAN on z 0 is coincidental, while NMI and Silhouette improvements reflect genuine gains from normalizing flow regularization and contrastive training.

Each of the baseline approaches was tested based on the same experimental protocol as that utilized in the proposed approach. Precisely, the identical image preprocessing (resizing to ( 224 × 224 ) and then converting into a tensor) and datasets were retained, and where applicable, the identical trial definitions and splits upon which the proposed approach is based were kept. In the case of clustering baselines, feature extraction occurred in the same representation space we compare, whether (i) learned latent embeddings our encoder means that the comparison is at the representation level, or (ii) baseline-specific features when a method intrinsically forms its representation. Each clustering measure (ARI, NMI, and Silhouette score) was calculated with a common evaluation code and presented in a homogeneous format (mean ± standard deviation of all the trials and/or bootstrap confidence intervals). In order to protect against unfairness, clustering baselines (e.g., k-means, Gaussian mixture models, DBSCAN) were used on the same latent embeddings consumed by the encoder, therefore guaranteeing that the different observed performance results were due to the clustering phenomenon and not just differences in feature extraction (Table 8).

Baseline hyperparameters had not been optimized in terms of ground-truth labels—we did not look to select parameter values that would maximize cluster quality measures like ARI or NMI or even improve the confusion matrix. In the case of the DBSCAN method, constant configurations were used over the entire study: the dimensionality reduction step was uniform manifold approximation and projection (UMAP) with n neighbors = 15 and min_dist = 0.1, and the random state was fixed to ensure reproducibility. Further clustering was performed using DBSCAN with ε = 1.5 and min_samples = 3. Although ARI scores are identical between standalone DBSCAN and NF-HAE-UDL on both datasets, the improvements in NMI (+0.06 STARE) and Silhouette (+0.07 STARE, +0.10 IDRiD) confirm that the normalizing flow produces more cohesive and better-separated clusters even when the partition agreement metric remains unchanged. Any further baseline algorithms (e.g., HDBSCAN or DEC) were initialized with standard reference defaults or with non-labeled data criteria (e.g., objective convergence). The exact values of the hyperparameters of all baselines are clearly recorded in Table 8 to foster transparency.

The performance of the proposed hybrid autoencoder with normalizing flow was evaluated against a range of unsupervised and modified approaches. Traditional unsupervised methods like Vanilla AE + KMeans, Variational AE (VAE) + GMM, deep embedded clustering (DEC), and deep clustering network (DCN) primarily rely on latent-space clustering. These approaches achieve moderate performance on the datasets but fail to capture subtle variations in retinal features, particularly under conditions of limited data. Similarly, density-based models, such as InfoGAN + KMeans, Self-Organizing Map (SOM), t-SNE + KMeans, and UMAP + HDBSCAN, show reasonable clustering behavior. On the other hand, with DEC + AE and Hybrid VAE + DCN, such modified approaches improve representation learning by integrating autoencoder-based feature extraction with clustering. However, these methods still underperform on the IDRiD dataset compared to the proposed model, reflecting their limited generalization capability for smaller datasets.

To position NF-HAE-UDL against current state-of-the-art representation learning techniques, we additionally compared against four self-supervised baselines: DINO (Caron et al., 2021), BYOL (Grill et al., 2020), SimSiam (Chen and He, 2021), and MAE (He et al., 2022). For each baseline, we extract frozen features from the respective pretrained backbone (without any fine-tuning on the retinal datasets) and apply the same UMAP + DBSCAN clustering pipeline used in our method. These methods represent powerful general-purpose representation learners but were developed for natural image datasets with large-scale pretraining corpora. Despite their strong performance on standard benchmarks, they are underperformed by NF-HAE-UDL on both STARE and IDRiD (Table 8), which we attribute to three factors: (1) the domain gap between natural images and retinal fundus photography, where pathological features such as microaneurysms and hemorrhages are absent from ImageNet pretraining; (2) the absence of any disease-specific contrastive objective DINO, BYOL, SimSiam, and MAE optimize for general invariance rather than disease-discriminative structure; (3) the extremely small dataset size ( N < 250 ), which limits the effectiveness of direct domain fine-tuning for these large-parameter architectures. SimSiam achieves the strongest performance among the self-supervised baselines, likely because masked image modeling encourages learning of local structural patterns that partially transfer to retinal images. Nevertheless, NF-HAE-UDL outperforms all four baselines across both datasets and all three clustering metrics, confirming that our combined hybrid encoding, normalizing flow regularization, and label-guided contrastive training is specifically advantageous for small-scale medical imaging datasets.

5 Discussion

5.1 Model validation through explainable AI analysis

To interpret which regions of the input image contributed most to cluster assignments, we employed gradient-weighted class activation mapping (Grad-CAM) (Suara et al., 2023). Our framework is unsupervised. Grad-CAM is applied post hoc after all representation learning completes and does not influence encoder training in any way. Standard Grad-CAM cannot be applied directly without supervised class logits, so we adopt the following protocol.

Step 1 (post-training): after unsupervised training completes, each image is assigned to a cluster using DBSCAN on the learned latent codes z.

Step 2 (interpretability probe): a lightweight linear classification head is trained to map latent codes z to cluster pseudo-labels, with encoder weights frozen as given by Equation 26:where c i is the cluster label assigned by DBSCAN and ℓ CE is cross-entropy loss. This classifier is trained for 250 epochs solely to enable visualization; it is never used in the unsupervised pipeline.

W cls * , b cls * = arg min W cls , b cls ∑ i = 1 N ℓ CE softmax W cls z i + b cls , c i , (26)

Step 3 (Grad-CAM application): the scalar objective for Grad-CAM is the predicted cluster logit and calculated as given by Equation 27.where ( ⋅ ) c denotes the cth element. Gradients ∂ score c / ∂ A ( k ) are backpropagated through the frozen encoder to the final convolutional layer feature maps A ( k ), producing saliency heatmaps via Equation 29.

score c = W cls z + b cls c , (27)

This approach is standard in self-supervised learning (Chen et al., 2020), where linear probes assess representation quality without invalidating unsupervised learning.

Scalar objective for Grad-CAM: the differentiable score used for Grad-CAM is the predicted cluster logit from the post hoc linear classifier. For an image assigned to cluster c by DBSCAN, the scalar objective is calculated using Equation 28.where ( ⋅ ) c denotes the cth element of the logit vector. Grad-CAM then computeswhere A ( k ) ∈ R ( H × W ) represents the kth feature map from the final convolutional layer of the encoder, and α k c are the importance weights obtained by global average pooling of the gradients.

score c = W cls z + b cls c , (28)

α k c = 1 H W ∑ i , j ∂ score c ∂ A i j k , L Grad - CAM c = ReLU ∑ k α k c A k , (29)

As shown in Figure 6, the heat maps demonstrate that our model focuses on clinically relevant regions for different retinal conditions. For the normal class, attention is distributed relatively uniformly across the fundus, with slight emphasis on the optic macula reflecting the absence of focal pathology. In cases of background DR, the model concentrates on areas with microaneurysms and hemorrhages, typically visible as small red lesions. For non-isolated telangiectasia, activation intensifies around vascular abnormalities and areas of dilated blood vessels. Finally, in central retinal vein occlusion (C.R.V. Occlusion), the model highlights regions with extensive hemorrhages and cotton wool, spotting characteristic of this severe vascular disorder. The main reason for using Grad-CAM is that it gives a clear explanation of how the model can distinguish between clusters and ensure that the model’s attention fits the pathological characteristics of integration in the STARE and IDRiD datasets. The whole process on post hoc Grad-CAM visualization is presented in Algorithm 5.

We acknowledge that the relatively small dataset size (81 images for STARE) raises legitimate concerns about overfitting. However, several factors mitigate this risk: (1) contrastive learning operates on augmented pairs, effectively increasing sample diversity, (2) the normalizing flow provides explicit regularization through likelihood maximization on a Gaussian prior, and (3) the consistent clustering performance across multiple trials suggests learned features generalize beyond memorizing individual training samples.

Algorithm 5 Post hoc Grad-CAM visualization.

5.1.1 Clinical validation of Grad-CAM explanations

To ensure clinical relevance beyond visual inspection, we performed the following sanity checks.

• Perturbation analysis: we systematically masked critical anatomical regions (optic disc, major vessels, lesion areas) using ( 32 × 32 ) patches and measured the change in cluster prediction confidence. Masking regions with high Grad-CAM activation caused an average confidence drop of 0.34 ± 0.08, while masking low-activation regions caused only 0.09 ± 0.05 drop, confirming that highlighted regions are indeed influential.
• Expert comparison: an experienced ophthalmologist (15+ years) annotated pathological regions in 20 randomly selected STARE images. We computed intersection-over-union (IoU) between Grad-CAM heatmaps (thresholded at top 25% activation) and expert annotations, achieving IoU = 0.52 ± 0.11; this indicates moderate spatial agreement with clinical ground truth.
• Negative control: Grad-CAM applied to randomly initialized encoders (no training) produced diffuse, non-localized attention maps with high spatial entropy ( 0.87 ± 0.04 ) compared to trained models (0.43 ± 0.06, p < 0.001); this confirmed that learned features drive focused attention rather than architectural biases.
• Class discriminability: for multi-label images (e.g., “Background DR + Arteriosclerotic Retinopathy”), separate Grad-CAMs for each disease cluster highlighted anatomically distinct regions (e.g., DR → microaneurysms near macula; ASR → vessel tortuosity), thus demonstrating class-specific explanations.
• Interpretation: these validation experiments confirm that the clinical relevance of Grad-CAM attention patterns arises from learned disease-specific features rather than architectural biases or random chance. The moderate IoU score ( 0.52 ± 0.11 ) reflects both meaningful alignment with expert annotations and the inherent variability in pathology localization across observers.

5.2 Model validation through auxiliary classification

This subsection presents auxiliary classification validation to assess whether learned representations preserve disease-discriminative information. We trained a lightweight linear classifier on frozen encoder outputs. This classifier is not part of the unsupervised framework but serves solely as an interpretability probe. All results are in-sample (evaluated on the same data used for representation learning) and represent representation quality assessment, not diagnostic performance. This approach follows standard practice in self-supervised learning (Chen et al., 2020) for validating learned features.

It can be noted that the clusters are mapped to disease labels using Hungarian matching (Chen et al., 2025) to calculate the classification accuracies by learning the full dataset. Here, the Hungarian algorithm (Hashemi et al., 2021), a powerful tool for bipartite graph matching, is used to find the best match between clusters and labels. The weighted correlation distance is applied to subsets obtained in each matching to sort the clusters.

Table 9 summarizes the trial-wise classification (Raini et al., 2024) of average accuracies (macro-level) over the classes. These results validate representation quality, not diagnostic capability. The classifier is trained on the training set latent codes and evaluated on both the training and test datasets as they randomly split with an 80/20 ratio to assess how well unsupervised representations preserve label-discriminative information by limiting the number of epochs and deploying pretrained backbones to prevent overfitting. The proposed model achieves high performance and stability across the classes since, on STARE, it has an average accuracy of 0.88 with a standard deviation of 0.03. Based on this, the probability vector of 11 dimensions in our model is obtained as the result of a sigmoid activation, which is trained using a multi-label loss (BCEWithLogitsLoss). Multi-label performance is assessed by applying subset (exact-match) and example-based evaluation, computed against the entire label set. Metrics are computed using the auxiliary classifier on frozen latent codes. At convergence (epoch 250 of classifier training, not encoder training), the classifier achieves subset accuracy of 0.7654 and a sample accuracy of 0.8817. It also attains a micro-F1 of 0.9000, a macro-F1 of 0.7083, a weighted-F1 of 0.8911, and a Hamming loss of 0.0224, indicating strong preservation of multi-label structure in the learned representations.

TABLE 9. Auxiliary classification validation across ten trials. A linear classifier trained on frozen latent codes assesses representation quality through in-sample accuracy. These metrics measure preservation of disease-discriminative information, not diagnostic performance.

To analyze the latent space, we also evaluate latent space clustering using UMAP embeddings and DBSCAN. The cluster measures at epoch 250 are as follows: ARI = 0.6585 (95% CI: [0.5028, 0.8345]), NMI = 0.7423 (95% CI: [0.6678, 0.8506]), and Silhouette score = 0.6314 (95% CI: [0.5766, 0.7073]).

Lastly, we incorporate a standard confusion matrix as a visual representation. To construct this, we reduce the multi-label ground truth and predictions to one representative class by performing an argmax operation. This provides an interpretable one-label view but does not substitute for the main multi-label assessment indicated above. In the case of IDRiD, the mean accuracy is 0.98 with a small variance, which indicates the high predictive power of the model with even small cases of disease and health. These findings affirm that the suggested technique can be used to capture meaningful attributes with small-scale fundus datasets better than the conventional unsupervised and modified baseline methods.

The confusion matrices for the STARE and IDRiD datasets are shown in Figures 7a,b, respectively, for illustration.

Based on the results in Figures 7a,b, the F1 measures are obtained as: micro-F1 = 0.95 and macro-F1 = 0.77 for the STARE dataset and micro-F1 = 0.99 and macro-F1 = 0.98 for the IDRiD dataset.

5.3 Limitations and future directions

Limitations. While our method demonstrates strong clustering performance and high auxiliary classification accuracy (88% for STARE and 98% for IDRiD), several limitations should be noted. First, the classification results represent in-sample validation rather than held-out diagnostic benchmarks, which may overestimate generalization. Second, the relatively small dataset sizes (81 images for STARE and 211 for IDRiD, with only 16 and 42 test images, respectively) limit the statistical power of our generalization assessment. Third, both datasets originate from single institutions with controlled imaging conditions, so cross-dataset evaluation on external databases (e.g., Messidor and APTOS) is needed to strengthen generalization claims. The high IDRiD test accuracy (94%) may indicate that the model exploits dataset-specific characteristics. Additionally, the normalizing flow module provides consistent but incremental improvements (+0.03 ARI and +0.06 NMI on STARE); while statistically significant, these gains are modest, with primary value in likelihood estimation (for anomaly detection) and invertibility (for latent manipulation) rather than dramatic clustering improvements. Grad-CAM visualizations require post hoc linear probes, and the moderate IoU ( 0.52 ± 0.11 ) with expert annotations suggests room for improvement in aligning learned attention with clinical reasoning. Auxiliary classification validation confirms representation quality but does not constitute clinical validation for diagnostic deployment. Lastly, other public datasets include a limited number of classes of retinal disorders, which restrain us to study subset selection or data bias regarding the reproducibility of the experimental results.

Future Work. To address these limitations, future studies should focus on: (1) external validation using larger, multi-center datasets with diverse imaging conditions to strengthen generalization claims; (2) semi-supervised or fine-tuning approaches to enable direct diagnostic performance while preserving data efficiency; (3) improving interpretability by developing end-to-end attention mechanisms or more robust alignment with clinical annotations; (4) prospective studies and regulatory approval necessary before clinical deployment, including the development of interactive visualization tools (e.g., Gradio-based interfaces) to facilitate translation to practical clinical workflows; (5) merging multiple datasets to mitigate the data bias problem.

6 Conclusion

Diagnosing retinal pathology using fundus images is a vital clinical problem, especially in resource-constrained locations where trained annotators are few and where it is not feasible to label large data sets. This study presents NF-HAE-UDL, a framework that automatically groups retinal images by disease type without requiring pixel-level annotations. The system combines three main components: a dual-branch hybrid encoder (ResNet-50 and Swin-Tiny) that extracts complementary local texture and global structural features from retinal images, a Glow-based normalizing flow module that organizes the learned features into a structured mathematical space to enable invertible analysis and anomaly detection, and a label-guided contrastive training objective that encourages the separation of disease-specific patterns in the feature space. Following training, disease groupings are produced by applying UMAP dimensionality reduction and DBSCAN clustering to the learned features, with ground-truth labels used only afterward to measure how well the groupings match clinical diagnoses. Tests on two publicly available datasets indicate that NF-HAE-UDL cluster the disease images better than any of the compared baseline methods. On the STARE dataset (81 multi-label fundus images spanning 11 disease categories), the model achieved consistent clustering across ten independent trials, outperforming classical methods such as k-means and hierarchical clustering by substantial margins. The same strong performance was observed on the IDRiD subset (211 diabetic retinopathy images), confirming that the framework generalizes across both multi-label and binary disease classification scenarios. The normalizing flow module was also found to have more practical value in clustering as it allows explicit density estimation that can be applied to anomaly detection and invertible manipulation of disease representations. To verify the clinical relevance of the learned features, Grad-CAM heatmaps were compared against annotations from an experienced ophthalmologist, confirming that the model correctly focuses on disease-relevant anatomical regions such as microaneurysms, hemorrhages, and vascular abnormalities. A simple linear classifier trained on the frozen learned features without any retraining of the main model achieved 88% accuracy on STARE and 98% on IDRiD, confirming that the features carry strong disease-discriminative information even when the model was never explicitly told what diseases to look for during training. The main limitations of this study are the small dataset sizes, the use of single-institution imaging data, and the fact that classification results reflect in-sample validation rather than independent diagnostic benchmarking. Future research will address these through external validation on larger multi-center datasets such as Messidor, APTOS, and DRHAGIS, semi-supervised fine-tuning to enable direct diagnostic performance benchmarking with regulatory validation, the development of end-to-end attention mechanisms that align more closely with clinical expert reasoning, and the creation of interactive visualization tools to make the framework accessible for clinical deployment. In general, NF-HAE-UDL is an exciting and feasible future for interpretable, data-efficient, and scalable automated retinal image analysis.

Generative AI statement

The author(s) declared that generative AI was not used in the creation of this manuscript.

Any alternative text (alt text) provided alongside figures in this article has been generated by Frontiers with the support of artificial intelligence and reasonable efforts have been made to ensure accuracy, including review by the authors wherever possible. If you identify any issues, please contact us.

Publisher’s note

All claims expressed in this article are solely those of the authors and do not necessarily represent those of their affiliated organizations, or those of the publisher, the editors and the reviewers. Any product that may be evaluated in this article, or claim that may be made by its manufacturer, is not guaranteed or endorsed by the publisher.

Supplementary material

The Supplementary Material for this article can be found online at: https://www.frontiersin.org/articles/10.3389/frsip.2026.1769553/full#supplementary-material

Supplementary material Supplementary material at the source

Footnotes

Edited by: Manish Sharma, Institute of Infrastructure, Technology, Research and Management, India

Reviewed by: Muhammad Hammad Malik, University of Science and Technology Beijing, China

Samantha Rajapaksha, Sri Lanka Institute of Information Technology, Sri Lanka

Sathishkumar Rajamani, Manakula Vinayagar Institute of Technology, India

Footnotes

STARE dataset: https://cecas.clemson.edu/∼ahoover/stare

IDRiD dataset: https://paperswithcode.com/dataset/idrid