rpotter6298
201 lines
7.6 KiB
Python
201 lines
7.6 KiB
Python
|
|
|
|
|
|
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")
|