rpotter6298
246 lines
8.5 KiB
Python
246 lines
8.5 KiB
Python
#!/usr/bin/env python3
|
||
"""
|
||
compare_clustering_methods.py
|
||
|
||
Compare two patient-clustering approaches against the v8 "reference":
|
||
1. Thumbnail-based (64×64 grayscale K-means — the manuscript's method)
|
||
2. Feature-based (PCA-50d of VGG16 features → K-means — our method)
|
||
|
||
Metrics (all label-invariant):
|
||
- Adjusted Rand Index (ARI)
|
||
- Normalized Mutual Information (NMI)
|
||
- V-measure (homogeneity + completeness)
|
||
|
||
Usage:
|
||
conda activate fundus_imaging
|
||
python scripts/compare_clustering_methods.py
|
||
"""
|
||
|
||
import os, sys, csv, re
|
||
import numpy as np
|
||
|
||
sys.path.insert(0, os.path.dirname(os.path.dirname(os.path.abspath(__file__))))
|
||
|
||
from sklearn.metrics import (
|
||
adjusted_rand_score,
|
||
normalized_mutual_info_score,
|
||
homogeneity_completeness_v_measure,
|
||
)
|
||
from sklearn.decomposition import PCA
|
||
from sklearn.cluster import KMeans
|
||
from PIL import Image
|
||
|
||
# ---------------------------------------------------------------------------
|
||
# Config
|
||
# ---------------------------------------------------------------------------
|
||
|
||
ROOT = os.path.dirname(os.path.dirname(os.path.abspath(__file__)))
|
||
DATASET_PATH = os.path.expanduser("~/Documents/data_leakage/The IQ-OTHNCCD lung cancer dataset")
|
||
MANIFEST_V8 = os.path.join(ROOT, ".archive", "results", "patient_manifest_v8.csv")
|
||
MANIFEST_SIMPLE = os.path.join(ROOT, "results", "simple_patient_manifest.csv")
|
||
FEATURES_PATH = os.path.join(ROOT, "features", "VGG16_features.npz")
|
||
|
||
PATIENT_COUNTS = {
|
||
"Bengin cases": 15,
|
||
"Malignant cases": 40,
|
||
"Normal cases": 55,
|
||
}
|
||
|
||
CLASS_ORDER = ["Bengin cases", "Malignant cases", "Normal cases"]
|
||
|
||
SEED = 42
|
||
THUMBNAIL_SIZE = (64, 64)
|
||
|
||
# ---------------------------------------------------------------------------
|
||
# Helpers
|
||
# ---------------------------------------------------------------------------
|
||
|
||
|
||
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
|
||
|
||
|
||
def load_v8_assignments(path):
|
||
"""Load v8 manifest, return {image_name: patient_id}."""
|
||
mapping = {}
|
||
with open(path, newline="") as f:
|
||
reader = csv.DictReader(f)
|
||
for row in reader:
|
||
for img in row["confirmed_images"].split(";"):
|
||
if img:
|
||
mapping[img] = row["patient_id"]
|
||
return mapping
|
||
|
||
|
||
def load_simple_assignments(path):
|
||
"""Load simple feature-based manifest, return {image_name: patient_id}."""
|
||
mapping = {}
|
||
with open(path, newline="") as f:
|
||
reader = csv.DictReader(f)
|
||
for row in reader:
|
||
for img in row["images"].split(";"):
|
||
if img:
|
||
mapping[img] = row["patient_id"]
|
||
return mapping
|
||
|
||
|
||
def run_thumbnail_clustering(dataset_path, filenames, labels):
|
||
"""Run 64×64 grayscale thumbnail K-means (manuscript method)."""
|
||
groups = np.array([None] * len(labels), dtype=object)
|
||
|
||
for class_name, n_clusters in PATIENT_COUNTS.items():
|
||
idx_class = np.where(labels == class_name)[0]
|
||
|
||
# Load thumbnails
|
||
X = []
|
||
for i in idx_class:
|
||
fname = filenames[i]
|
||
img_path = os.path.join(dataset_path, class_name, fname)
|
||
try:
|
||
img = Image.open(img_path).convert("L")
|
||
img = img.resize(THUMBNAIL_SIZE)
|
||
img_arr = np.array(img, dtype=np.float32) / 255.0
|
||
X.append(img_arr.flatten())
|
||
except Exception:
|
||
X.append(np.zeros(64 * 64, dtype=np.float32))
|
||
X = np.array(X, dtype=np.float32)
|
||
|
||
# PCA → 50d
|
||
n_pca = min(50, X.shape[0] - 1, X.shape[1])
|
||
pca = PCA(n_components=n_pca, random_state=SEED)
|
||
X_pca = pca.fit_transform(X)
|
||
|
||
# K-means
|
||
kmeans = KMeans(n_clusters=n_clusters, random_state=SEED, n_init=20)
|
||
clusters = kmeans.fit_predict(X_pca)
|
||
|
||
for i, cluster_id in zip(idx_class, clusters):
|
||
groups[i] = f"{class_name}_cluster_{cluster_id}"
|
||
|
||
sizes = np.bincount(clusters)
|
||
print(f" {class_name}: k={n_clusters}, sizes min={sizes.min()} max={sizes.max()} mean={sizes.mean():.1f}")
|
||
|
||
return groups
|
||
|
||
|
||
# ---------------------------------------------------------------------------
|
||
# Main
|
||
# ---------------------------------------------------------------------------
|
||
|
||
print("Loading VGG16 features for filename/label reference...")
|
||
data = np.load(FEATURES_PATH, allow_pickle=True)
|
||
all_filenames = data["filenames"]
|
||
all_labels = data["Y"]
|
||
X_vgg16 = data["X"]
|
||
|
||
print(f" {len(all_filenames)} images across {len(np.unique(all_labels))} classes")
|
||
|
||
# Load reference (v8) and our feature-based assignments
|
||
v8_map = load_v8_assignments(MANIFEST_V8)
|
||
simple_map = load_simple_assignments(MANIFEST_SIMPLE)
|
||
|
||
# Build per-image label arrays for all three methods, aligned by filename
|
||
# We need all images that exist in ALL three
|
||
v8_labels_list = []
|
||
simple_labels_list = []
|
||
thumb_labels_list = [] # filled after clustering
|
||
common_filenames = []
|
||
common_labels = []
|
||
|
||
# First, run thumbnail clustering
|
||
print("\nRunning thumbnail-based K-means (manuscript method)...")
|
||
thumb_groups = run_thumbnail_clustering(DATASET_PATH, all_filenames, all_labels)
|
||
|
||
# Build thumbnail mapping (using short names like B_001)
|
||
thumb_map = {}
|
||
for fname, group in zip(all_filenames, thumb_groups):
|
||
short = f2n(fname)
|
||
for cls in CLASS_ORDER:
|
||
if group.startswith(cls):
|
||
cluster_id = int(group.split("_cluster_")[-1])
|
||
prefix = {"Bengin cases": "Benign", "Malignant cases": "Malignant", "Normal cases": "Normal"}[cls]
|
||
thumb_map[short] = f"thumb_{prefix}_{cluster_id:02d}"
|
||
break
|
||
|
||
# Filter to images present in all three
|
||
# Feature filenames are like "Bengin case (1).jpg", manifests use "B_001"
|
||
for fname, label in zip(all_filenames, all_labels):
|
||
short_name = f2n(fname)
|
||
v8_id = v8_map.get(short_name)
|
||
simple_id = simple_map.get(short_name)
|
||
thumb_id = thumb_map.get(short_name)
|
||
if v8_id and simple_id and thumb_id:
|
||
common_filenames.append(short_name)
|
||
common_labels.append(label)
|
||
v8_labels_list.append(v8_id)
|
||
simple_labels_list.append(simple_id)
|
||
thumb_labels_list.append(thumb_id)
|
||
|
||
print(f"\nImages common to all three methods: {len(common_filenames)}")
|
||
|
||
# Convert to numpy arrays
|
||
v8_labels_arr = np.array(v8_labels_list)
|
||
simple_labels_arr = np.array(simple_labels_list)
|
||
thumb_labels_arr = np.array(thumb_labels_list)
|
||
common_labels_arr = np.array(common_labels)
|
||
|
||
# ---------------------------------------------------------------------------
|
||
# Compute agreement metrics — per-class and overall
|
||
# ---------------------------------------------------------------------------
|
||
|
||
|
||
def compute_metrics(ref, pred, name):
|
||
"""Compute clustering agreement metrics against reference."""
|
||
ari = adjusted_rand_score(ref, pred)
|
||
nmi = normalized_mutual_info_score(ref, pred)
|
||
h, c, v = homogeneity_completeness_v_measure(ref, pred)
|
||
return {"name": name, "ARI": ari, "NMI": nmi, "Homogeneity": h, "Completeness": c, "V_measure": v}
|
||
|
||
|
||
print("\n" + "=" * 80)
|
||
print("OVERALL AGREEMENT WITH v8 REFERENCE")
|
||
print("=" * 80)
|
||
|
||
results = []
|
||
for name, pred in [("Thumbnail (manuscript)", thumb_labels_arr), ("Feature-based (ours)", simple_labels_arr)]:
|
||
r = compute_metrics(v8_labels_arr, pred, name)
|
||
results.append(r)
|
||
|
||
print(f"\n{'Method':<30s} {'ARI':>8s} {'NMI':>8s} {'Homog':>8s} {'Compl':>8s} {'V_meas':>8s}")
|
||
print("-" * 72)
|
||
for r in results:
|
||
print(f"{r['name']:<30s} {r['ARI']:>8.4f} {r['NMI']:>8.4f} {r['Homogeneity']:>8.4f} {r['Completeness']:>8.4f} {r['V_measure']:>8.4f}")
|
||
|
||
# Per-class breakdown
|
||
print("\n" + "=" * 80)
|
||
print("PER-CLASS ARI WITH v8 REFERENCE")
|
||
print("=" * 80)
|
||
print(f"\n{'Class':<20s} {'Thumbnail':>10s} {'Feature-based':>15s}")
|
||
print("-" * 47)
|
||
for cls in CLASS_ORDER:
|
||
mask = common_labels_arr == cls
|
||
if mask.sum() < 2:
|
||
continue
|
||
thumb_ari = adjusted_rand_score(v8_labels_arr[mask], thumb_labels_arr[mask])
|
||
feat_ari = adjusted_rand_score(v8_labels_arr[mask], simple_labels_arr[mask])
|
||
better = " ←" if feat_ari > thumb_ari else ""
|
||
print(f"{cls:<20s} {thumb_ari:>10.4f} {feat_ari:>15.4f}{better}")
|
||
|
||
# Also: direct agreement between thumbnail and feature-based
|
||
print("\n" + "=" * 80)
|
||
print("THUMBNAIL vs FEATURE-BASED (direct agreement)")
|
||
print("=" * 80)
|
||
direct = compute_metrics(thumb_labels_arr, simple_labels_arr, "Thumb vs Feature")
|
||
print(f" ARI={direct['ARI']:.4f} NMI={direct['NMI']:.4f} V_measure={direct['V_measure']:.4f}")
|
||
|
||
print("\nDONE")
|