import os, sys, re, csv, json from collections import defaultdict import numpy as np import matplotlib import matplotlib.pyplot as plt from sklearn.decomposition import PCA from sklearn.manifold import TSNE from sklearn.cluster import KMeans import argparse #!/usr/bin/env python3 """ simple_patient_tsne.py The short path: 1. Load VGG16 features 2. K-means per class (15/40/55 patients) in PCA-50d space 3. Optional: iterative centroid refinement 4. Plot t-SNE colored by cluster, with centroid labels Usage: conda activate fundus_imaging python scripts/simple_patient_tsne.py """ matplotlib.use("Agg") sys.path.insert(0, os.path.dirname(os.path.dirname(os.path.abspath(__file__)))) PROJECT_ROOT = os.path.dirname(os.path.dirname(os.path.abspath(__file__))) FEATURES_DIR = os.path.join(PROJECT_ROOT, "features") PLOTS_DIR = os.path.join(PROJECT_ROOT, "plots") RESULTS_DIR = os.path.join(PROJECT_ROOT, "results") os.makedirs(PLOTS_DIR, exist_ok=True) os.makedirs(RESULTS_DIR, exist_ok=True) PATIENT_COUNTS = {"Bengin cases": 15, "Malignant cases": 40, "Normal cases": 55} CLASS_NAMES = {"Bengin cases": "Benign", "Malignant cases": "Malignant", "Normal cases": "Normal"} RANDOM_STATE = 42 # --------------------------------------------------------------------------- # 1. Load VGG16 features # --------------------------------------------------------------------------- print("Loading VGG16 features ...") data = np.load(os.path.join(FEATURES_DIR, "VGG16_features.npz"), allow_pickle=True) X, Y, filenames = data["X"], data["Y"], data["filenames"] def f2n(fname): m = re.search(r'\((\d+)\)', fname) num = int(m.group(1)) if m else None for cls_key, prefix in [("Bengin cases","B"),("Malignant cases","M"),("Normal cases","N")]: if fname.startswith(cls_key.rstrip("s")): return f"{prefix}_{num:03d}" if num else fname return fname img_nums = np.array([f2n(f) for f in filenames]) # --------------------------------------------------------------------------- # 2. K-means per class in PCA-50d space # --------------------------------------------------------------------------- print("Clustering in PCA-50d space ...") n_pca = min(50, X.shape[0] - 1, X.shape[1]) X_pca = PCA(n_components=n_pca, random_state=RANDOM_STATE).fit_transform(X) image_to_patient = {} patient_to_images = defaultdict(list) for class_name, k in PATIENT_COUNTS.items(): mask = Y == class_name X_class = X_pca[mask] idx_class = np.where(mask)[0] kmeans = KMeans(n_clusters=k, random_state=RANDOM_STATE, n_init=20) labels = kmeans.fit_predict(X_class) prefix = {"Bengin cases": "Benign", "Malignant cases": "Malignant", "Normal cases": "Normal"}[class_name] for i, cluster_id in enumerate(labels): pid = f"{prefix}_{cluster_id:02d}" img = img_nums[idx_class[i]] image_to_patient[img] = pid patient_to_images[pid].append(img) print(f" {len(patient_to_images)} patients, {len(image_to_patient)} images") # --------------------------------------------------------------------------- # 3. Iterative centroid refinement (optional, 3 passes) # --------------------------------------------------------------------------- print("Refining assignments (nearest-centroid, 5 passes) ...") for iteration in range(5): # Compute centroids centroids = {} for pid, imgs in patient_to_images.items(): idxs = [np.where(img_nums == img)[0][0] for img in imgs] centroids[pid] = X_pca[idxs].mean(axis=0) # Reassign moves = 0 for class_name in PATIENT_COUNTS: mask = Y == class_name for i in np.where(mask)[0]: img = img_nums[i] old_pid = image_to_patient[img] # Find nearest centroid in same class best_pid = old_pid best_dist = float('inf') for pid, c in centroids.items(): if pid.startswith(CLASS_NAMES[class_name]): d = float(np.linalg.norm(X_pca[i] - c)) if d < best_dist: best_dist = d best_pid = pid if best_pid != old_pid: patient_to_images[old_pid].remove(img) patient_to_images[best_pid].append(img) image_to_patient[img] = best_pid moves += 1 print(f" Pass {iteration+1}: {moves} moves") if moves == 0: break # --------------------------------------------------------------------------- # 4. t-SNE # --------------------------------------------------------------------------- print("Computing t-SNE ...") X_tsne = TSNE(n_components=2, perplexity=35, learning_rate="auto", init="pca", random_state=RANDOM_STATE).fit_transform(X_pca) # --------------------------------------------------------------------------- # 5. Plot — one figure per class # --------------------------------------------------------------------------- for class_name, display_name in CLASS_NAMES.items(): fig, ax = plt.subplots(1, 1, figsize=(14, 10)) ax.scatter(X_tsne[:, 0], X_tsne[:, 1], c="lightgray", s=3, alpha=0.15) mask = Y == class_name class_pids = sorted([p for p in patient_to_images if p.startswith(display_name)]) n_patients = len(class_pids) cmap = plt.cm.tab20 if n_patients <= 20 else plt.cm.gist_ncar for pi, pid in enumerate(class_pids): color = cmap(pi % 20) if n_patients <= 20 else cmap(pi / max(n_patients-1, 1)) pts_x, pts_y = [], [] for img in patient_to_images[pid]: i = np.where(img_nums == img)[0][0] pts_x.append(X_tsne[i, 0]) pts_y.append(X_tsne[i, 1]) ax.scatter(pts_x, pts_y, c=[color], s=18, alpha=0.8, label='_nolegend_') # Centroid diamond (no label) cx, cy = np.mean(pts_x), np.mean(pts_y) ax.scatter(cx, cy, c=[color], s=60, marker='D', edgecolors='black', linewidths=0.6, zorder=5, label='_nolegend_') # Legend elements from matplotlib.lines import Line2D legend_elements = [ Line2D([0], [0], marker='o', color='w', markerfacecolor='gray', markersize=8, label='Patient images (dots)'), Line2D([0], [0], marker='D', color='w', markerfacecolor='gray', markersize=8, label='Patient centroids (diamonds)'), ] ax.legend(handles=legend_elements, loc='lower right') ax.set_title(f"{display_name} — VGG16 t-SNE ({n_patients} patients)") ax.set_xlabel("t-SNE dim 1") ax.set_ylabel("t-SNE dim 2") plt.tight_layout() ap = argparse.ArgumentParser() ap.add_argument("--tag", default="", help="Append tag to filename") args = ap.parse_args() tag = f"_{args.tag}" if args.tag else "" out = os.path.join(PLOTS_DIR, "tsne", f"tsne_{display_name}{tag}.png") os.makedirs(os.path.dirname(out), exist_ok=True) plt.savefig(out, dpi=150) plt.close() print(f" Saved {out}") # --------------------------------------------------------------------------- # 6. Save assignments # --------------------------------------------------------------------------- manifest = [] for pid in sorted(patient_to_images.keys()): imgs = sorted(patient_to_images[pid]) manifest.append({"patient_id": pid, "class": pid.split("_")[0], "n_images": len(imgs), "images": ";".join(imgs)}) with open(os.path.join(RESULTS_DIR, "simple_patient_manifest.csv"), "w", newline="") as f: w = csv.DictWriter(f, fieldnames=["patient_id", "class", "n_images", "images"]) w.writeheader() w.writerows(manifest) print(f"\nSaved simple_patient_manifest.csv ({len(manifest)} patients, " f"{sum(m['n_images'] for m in manifest)} images)") print("DONE")