import os
import gc
import glob
import shutil
import warnings
warnings.filterwarnings("ignore")
import numpy as np
import pandas as pd
from tqdm import tqdm
# 图像处理与特征
import cv2
from skimage import feature
from skimage.feature import canny
# 自动兼容旧版或新版 scikit-image
if hasattr(feature, 'graycomatrix'):
graycomatrix = feature.graycomatrix
graycoprops = feature.graycoprops
else:
graycomatrix = feature.greycomatrix
graycoprops = feature.greycoprops
# 统计学与评估
from scipy.stats import skew, kurtosis
# 机器学习
from sklearn.model_selection import train_test_split, StratifiedKFold
from sklearn.preprocessing import StandardScaler
from sklearn.decomposition import PCA
from sklearn.discriminant_analysis import LinearDiscriminantAnalysis as LDA
from sklearn.linear_model import LogisticRegression
from sklearn.metrics import (confusion_matrix, roc_curve, roc_auc_score,
precision_recall_fscore_support, brier_score_loss)
# 可视化（全部保存到 result1）
import matplotlib
matplotlib.use('Agg') # 不显示，只保存
import matplotlib.pyplot as plt
import os
import glob
# ======================
# 一、基础设置与工具函数
# ======================
# 结果输出目录
RESULT_DIR = "result1"
os.makedirs(RESULT_DIR, exist_ok=True)
# 可能的图像与标签目录（两种常见组织方式）
CANDIDATE_DIRS = [
# YOLO标准组织
{"img_dir": "images/train", "lab_dir": "labels/train"},
{"img_dir": "images/val", "lab_dir": "labels/val"},
# JPEGImages + labels 组织
{"img_dir": "JPEGImages", "lab_dir": "labels"},
]
# 允许的图片后缀
IMG_EXT = ["*.jpg", "*.jpeg", "*.png", "*.bmp"]
def list_images(img_dir):
"""列出目录下的所有图片路径（兼容多后缀）"""
paths = []
for ext in IMG_EXT:
paths.extend(glob.glob(os.path.join(img_dir, ext)))
return sorted(paths)
def yolo_label_path_for_image(img_path, lab_root):
"""根据图像路径与标签根目录，推断同名YOLO标签文件路径（.txt）"""
base = os.path.splitext(os.path.basename(img_path))[0]
return os.path.join(lab_root, base + ".txt")
def read_classes_map(classes_file="classes.txt"):
"""读取类别映射（若存在），仅用于信息提示；问题1只需二分类标签"""
if os.path.exists(classes_file):
try:
with open(classes_file, "r", encoding="utf-8") as f:
lines = [x.strip() for x in f.readlines() if x.strip()]
return {i: name for i, name in enumerate(lines)}
except Exception:
return None
return None
def imread_rgb(path, resize_to=None):
"""读取并返回RGB图像，必要时缩放；对异常图像进行鲁棒处理"""
img = cv2.imread(path, cv2.IMREAD_COLOR)
if img is None:
raise IOError(f"无法读取图像：{path}")
img = cv2.cvtColor(img, cv2.COLOR_BGR2RGB)
if resize_to is not None:
img = cv2.resize(img, resize_to, interpolation=cv2.INTER_AREA)
return img
def safe_minmax_norm(gray):
"""对灰度图做Min-Max归一化到[0,1]"""
gmin, gmax = np.min(gray), np.max(gray)
if gmax > gmin:
return (gray - gmin) / (gmax - gmin)
return gray * 0.0
def edge_density(gray_norm):
"""基于Canny的边缘密度（归一化灰度作为输入）"""
# 根据图像对比度自适应阈值（简单启发式）
v = np.median(gray_norm)
lower = max(0, (1.0 - 0.33) * v)
upper = min(1.0, (1.0 + 0.33) * v)
# skimage.feature.canny 接受 [0,1] 浮点图像
edges = canny(gray_norm, sigma=1.5, low_threshold=lower, high_threshold=upper)
return edges.mean(), edges # 返回比例与边缘图
def gradient_stats(gray_norm):
"""Sobel梯度模的统计量"""
g = (gray_norm * 255.0).astype(np.uint8)
gx = cv2.Sobel(g, cv2.CV_32F, 1, 0, ksize=3)
gy = cv2.Sobel(g, cv2.CV_32F, 0, 1, ksize=3)
mag = np.sqrt(gx * gx + gy * gy)
mag = mag.flatten()
return np.mean(mag), np.std(mag), skew(mag), kurtosis(mag, fisher=True)
def gray_hist_stats(gray_norm, bins=32):
"""灰度直方图统计：均值、方差、偏度、峰度、熵、能量等"""
g = (gray_norm * 255.0).astype(np.uint8)
hist = cv2.calcHist([g], [0], None, [bins], [0, 256]).ravel().astype(np.float32)
hist /= (hist.sum() + 1e-8)
vals = np.arange(bins, dtype=np.float32)
mean_ = (hist * vals).sum()
var_ = (hist * (vals - mean_) ** 2).sum()
sk_ = ((hist * (vals - mean_) ** 3).sum()) / (np.sqrt(var_) ** 3 + 1e-8)
ku_ = ((hist * (vals - mean_) ** 4).sum()) / (var_ ** 2 + 1e-8) - 3.0
entropy_ = -(hist * np.log(hist + 1e-12)).sum()
energy_ = (hist ** 2).sum()
return mean_, var_, sk_, ku_, entropy_, energy_
def glcm_features(gray_norm, distances=(1, 2, 3), angles=(0, np.pi / 4, np.pi / 2, 3 * np.pi / 4), levels=32):
"""GLCM纹理特征（多距离、多角度求均值）"""
g = (gray_norm * (levels - 1)).astype(np.uint8)
glcm = graycomatrix(g, distances=distances, angles=angles, levels=levels, symmetric=True, normed=True)
feats = {}
for prop in ["contrast", "dissimilarity", "homogeneity", "energy", "correlation", "ASM"]:
vals = graycoprops(glcm, prop).ravel()
feats[prop + "_mean"] = float(np.mean(vals))
feats[prop + "_std"] = float(np.std(vals))
return feats
def color_steady_features(img_rgb):
"""颜色稳态特征（HSV的H/S/V统计 + 对比度受限自适应直方图前后的亮度稳态差异）"""
hsv = cv2.cvtColor(img_rgb, cv2.COLOR_RGB2HSV)
H, S, V = hsv[..., 0], hsv[..., 1], hsv[..., 2]
# 归一化
Hn = H / 180.0
Sn = S / 255.0
Vn = V / 255.0
feats = {
"H_mean": float(Hn.mean()),
"H_std": float(Hn.std()),
"S_mean": float(Sn.mean()),
"S_std": float(Sn.std()),
"V_mean": float(Vn.mean()),
"V_std": float(Vn.std()),
}
# CLAHE 增强前后对比（度量光照鲁棒性）
clahe = cv2.createCLAHE(clipLimit=2.0, tileGridSize=(8, 8))
V_eq = clahe.apply((Vn * 255).astype(np.uint8)).astype(np.float32) / 255.0
feats["V_clahe_delta_mean"] = float((V_eq - Vn).mean())
feats["V_clahe_delta_std"] = float((V_eq - Vn).std())
return feats
def extract_features_for_image(img_path, resize_to=(640, 640)):
"""对单张图像提取特征，返回特征字典与中间过程（可选）"""
img = imread_rgb(img_path, resize_to=resize_to)
# 灰度与归一化
gray = cv2.cvtColor(img, cv2.COLOR_RGB2GRAY).astype(np.float32)
gray = safe_minmax_norm(gray)
# 直方图与统计
gh_mean, gh_var, gh_sk, gh_ku, gh_entropy, gh_energy = gray_hist_stats(gray)
# 边缘密度与梯度
ed_ratio, ed_map = edge_density(gray)
g_mean, g_std, g_sk, g_ku = gradient_stats(gray)
# GLCM纹理
glcm_feats = glcm_features(gray)
# 颜色稳态
color_feats = color_steady_features(img)
# 汇总
feats = {
"gray_hist_mean": gh_mean,
"gray_hist_var": gh_var,
"gray_hist_skew": gh_sk,
"gray_hist_kurt": gh_ku,
"gray_hist_entropy": gh_entropy,
"gray_hist_energy": gh_energy,
"edge_density": ed_ratio,
"grad_mean": g_mean,
"grad_std": g_std,
"grad_skew": g_sk,
"grad_kurt": g_ku,
}
feats.update(glcm_feats)
feats.update(color_feats)
return feats
def determine_binary_label(label_txt_path):
"""根据YOLO标签txt是否存在且非空，确定是否有破损(1=有；0=无)"""
if (not os.path.exists(label_txt_path)) or (os.path.getsize(label_txt_path) == 0):
return 0
# 存在且非空 → 有破损
return 1
# ======================
# 二、数据发现：定位真实数据目录
# ======================
def discover_dataset():
"""
聚合训练样本（图像路径 + 二分类标签）
优先使用 images/train + labels/train
若val存在，追加；若只有 JPEGImages + labels，则以此为准
"""
records = []
found_any = False
# 优先 YOLO 组织：train + val
for d in [{"img_dir": "images/train", "lab_dir": "labels/train"},
{"img_dir": "images/val", "lab_dir": "labels/val"}]:
if os.path.isdir(d["img_dir"]) and os.path.isdir(d["lab_dir"]):
imgs = list_images(d["img_dir"])
if len(imgs) > 0:
for ip in imgs:
lp = yolo_label_path_for_image(ip, d["lab_dir"])
y = determine_binary_label(lp)
records.append({"img_path": ip, "label_path": lp, "y": y})
found_any = True
# 若以上没有任何样本，再尝试 JPEGImages + labels
if (not found_any) and os.path.isdir("JPEGImages") and os.path.isdir("labels"):
imgs = list_images("JPEGImages")
if len(imgs) > 0:
for ip in imgs:
lp = yolo_label_path_for_image(ip, "labels")
y = determine_binary_label(lp)
records.append({"img_path": ip, "label_path": lp, "y": y})
found_any = True
# 收拢为 DataFrame
df = pd.DataFrame(records)
return df, found_any
# ======================
# 三、主流程：特征构建 → LDA → 概率校准（可选）
# ======================
def main():
# 读取类别映射（非必须，问题1只需要二分类标签）
classes_map = read_classes_map("classes.txt")
if classes_map is not None:
with open(os.path.join(RESULT_DIR, "classes_map.txt"), "w", encoding="utf-8") as f:
for k, v in classes_map.items():
f.write(f"{k}\t{v}\n")
# 发现数据
df, ok = discover_dataset()
if (not ok) or (len(df) == 0):
# 如无可用数据，给出友好提示并退出
msg = (
"未发现可用的图像与标签，请检查根目录数据组织是否为以下任一形式：\n"
" A) images/train + labels/train [+ images/val + labels/val]\n"
" B) JPEGImages + labels\n"
"同时确保图像为 jpg/png 等格式，标签为 YOLO 的 .txt 文件。"
)
print(msg)
# 保存提示到文件，便于排障
with open(os.path.join(RESULT_DIR, "DATA_NOT_FOUND.txt"), "w", encoding="utf-8") as f:
f.write(msg)
return
# 打乱并划分训练/验证（保留标签分布）
df = df.sample(frac=1.0, random_state=42).reset_index(drop=True)
train_df, val_df = train_test_split(df, test_size=0.2, random_state=42, stratify=df["y"])
# 提取特征（为节约时间，默认将图像缩放到 640x640；如需更快，可降到 480 或 384）
def build_feature_table(sub_df, tag):
feat_list = []
pbar = tqdm(sub_df.itertuples(index=False), total=len(sub_df), desc=f"提取特征[{tag}]")
for row in pbar:
feats = extract_features_for_image(row.img_path, resize_to=(640, 640))
feats["y"] = row.y
feats["img_path"] = row.img_path
feat_list.append(feats)
feat_df = pd.DataFrame(feat_list)
feat_df.to_csv(os.path.join(RESULT_DIR, f"features_{tag}.csv"), index=False, encoding="utf-8-sig")
return feat_df
train_feat = build_feature_table(train_df, "train")
val_feat = build_feature_table(val_df, "val")
# 特征矩阵与标签
feature_cols = [c for c in train_feat.columns if c not in ("y", "img_path")]
X_train = train_feat[feature_cols].values
y_train = train_feat["y"].values.astype(int)
X_val = val_feat[feature_cols].values
y_val = val_feat["y"].values.astype(int)
# 标准化 + PCA（保留95%方差）
scaler = StandardScaler()
X_train_std = scaler.fit_transform(X_train)
X_val_std = scaler.transform(X_val)
pca = PCA(n_components=0.95, svd_solver="full", random_state=42)
X_train_pca = pca.fit_transform(X_train_std)
X_val_pca = pca.transform(X_val_std)
# 训练 LDA（主模型）
lda = LDA()
lda.fit(X_train_pca, y_train)
# 验证集预测
val_decision = lda.decision_function(X_val_pca) # 连续分数
# 将决策分数线性映射到(0,1)近似为概率（便于画 ROC/Brier），也可换用逻辑回归做校准
# 这里保持“不过度工程化”的原则，直接做min-max归一化
p_min, p_max = val_decision.min(), val_decision.max()
if p_max > p_min:
val_prob = (val_decision - p_min) / (p_max - p_min)
else:
val_prob = (val_decision - p_min) # 常数，退化情况
y_pred = (val_prob >= 0.5).astype(int)
# 评估（混淆矩阵、ROC、分布可视化、PCA方差解释率、LDA系数解析）
cm = confusion_matrix(y_val, y_pred)
prc, rec, f1, _ = precision_recall_fscore_support(y_val, y_pred, average="binary", zero_division=0)
auc = roc_auc_score(y_val, val_prob) if len(np.unique(y_val)) == 2 else np.nan
brier = brier_score_loss(y_val, val_prob)
# 保存评估结果
with open(os.path.join(RESULT_DIR, "eval_ql.txt"), "w", encoding="utf-8") as f:
f.write("=== 问题1：有无破损二分类（LDA）评估 ===\n")
f.write(f"样本数（训练/验证）：{len(train_feat)} / {len(val_feat)}\n")
f.write("\n[混淆矩阵]\n")
f.write(pd.DataFrame(cm, index=["真0(无破损)", "真1(有破损)"], columns=["预测0", "预测1"]).to_string())
f.write("\n\n[指标]\n")
f.write(f"Precision: {prc:.4f}\n")
f.write(f"Recall: {rec:.4f}\n")
f.write(f"F1: {f1:.4f}\n")
f.write(f"ROC-AUC: {auc:.4f}\n")
f.write(f"Brier: {brier:.6f}\n")
f.write("\n[特征维度→PCA后维度]\n")
f.write(f"{X_train.shape[1]} → {X_train_pca.shape[1]}\n")
# ======================
# 四、可视化（全部保存到 result1）
# ======================
# 1）混淆矩阵热力图
plt.figure(figsize=(5, 4))
plt.imshow(cm, cmap="Blues")
plt.title("Confusion Matrix (LDA)")
plt.colorbar()
tick_labels = ["真0(无破损)", "真1(有破损)"]
plt.xticks([0, 1], ["预测0", "预测1"])
plt.yticks([0, 1], tick_labels)
for i in range(cm.shape[0]):
for j in range(cm.shape[1]):
plt.text(j, i, cm[i, j], ha="center", va="center", color="black")
plt.tight_layout()
plt.savefig(os.path.join(RESULT_DIR, "cm_lda.png"), dpi=200)
plt.close()
# 2）ROC 曲线
fpr, tpr, _ = roc_curve(y_val, val_prob)
plt.figure(figsize=(5, 4))
plt.plot(fpr, tpr, lw=2)
plt.plot([0, 1], [0, 1], linestyle="--")
plt.xlabel("False Positive Rate")
plt.ylabel("True Positive Rate")
plt.title(f"ROC Curve (AUC={auc:.3f})")
plt.tight_layout()
plt.savefig(os.path.join(RESULT_DIR, "roc_lda.png"), dpi=200)
plt.close()
# 3）概率直方图（校准感）
plt.figure(figsize=(5, 4))
plt.hist(val_prob[y_val == 0], bins=20, alpha=0.6, label="无破损")
plt.hist(val_prob[y_val == 1], bins=20, alpha=0.6, label="有破损")
plt.title("Predicted Score Distribution")
plt.xlabel("score")
plt.ylabel("count")
plt.legend()
plt.tight_layout()
plt.savefig(os.path.join(RESULT_DIR, "score_hist.png"), dpi=200)
plt.close()
# 4）PCA 累计方差解释率
plt.figure(figsize=(6, 4))
cumvar = np.cumsum(pca.explained_variance_ratio_)
plt.plot(np.arange(1, len(cumvar) + 1), cumvar, marker="o")
plt.xlabel("Number of Components")
plt.ylabel("Cumulative Explained Variance")
plt.title("PCA Cumulative Explained Variance")
plt.grid(alpha=0.3)
plt.tight_layout()
plt.savefig(os.path.join(RESULT_DIR, "pca_cumvar.png"), dpi=200)
plt.close()
# 5）LDA 投影分布
plt.figure(figsize=(5, 4))
proj = lda.transform(X_val_pca).ravel()
plt.hist(proj[y_val == 0], bins=30, alpha=0.6, label="无破损")
plt.hist(proj[y_val == 1], bins=30, alpha=0.6, label="有破损")
plt.title("LDA Projection (Validation)")
plt.xlabel("projection value")
plt.ylabel("count")
plt.legend()
plt.tight_layout()
plt.savefig(os.path.join(RESULT_DIR, "lda_projection.png"), dpi=200)
plt.close()
# 6）特征重要性近似可视化（以 LDA 在 PCA 空间的系数范数替代说明）
# 注：LDA 在 PCA 空间的系数 → 回映射到原特征并不直接，一般作为示意
coefs = lda.coef_.ravel()
plt.figure(figsize=(6, 3))
plt.plot(np.arange(len(coefs)), coefs, marker="o", linewidth=1)
plt.title("LDA Coefficients (in PCA space)")
plt.xlabel("component index")
plt.ylabel("coefficient")
plt.grid(alpha=0.3)
plt.tight_layout()
plt.savefig(os.path.join(RESULT_DIR, "lda_coefs.png"), dpi=200)
plt.close()
# 7）样例图像 + 边缘图叠加（增强可读性，随机抽取少量保存）
sample_show = val_df.sample(n=min(6, len(val_df)), random_state=7)
fig, axes = plt.subplots(2, 3, figsize=(10, 6))
axes = axes.ravel()
for ax, (_, row) in zip(axes, sample_show.iterrows()):
img = imread_rgb(row.img_path, resize_to=(640, 640))
gray = cv2.cvtColor(img, cv2.COLOR_RGB2GRAY).astype(np.float32)
gray = safe_minmax_norm(gray)
_, edges = edge_density(gray)
# 叠加显示
overlay = img.copy()
overlay[edges] = [255, 0, 0] # 用红色标记边缘
ax.imshow(overlay)
ax.set_title(f"y={row.y}")
ax.axis("off")
plt.tight_layout()
plt.savefig(os.path.join(RESULT_DIR, "samples_overlay_edges.png"), dpi=200)
plt.close()
# 保存最终验证集预测结果（供溯源）
out_df = val_feat[["img_path"]].copy()
out_df["y_true"] = y_val
out_df["y_prob"] = val_prob
out_df["y_pred"] = y_pred
out_df.to_csv(os.path.join(RESULT_DIR, "val_predictions.csv"), index=False, encoding="utf-8-sig")
print("训练与评估完成。结果文件与图片已保存到 ./result1/ 目录。")
if __name__ == "__main__":
main()
扫一扫
1万+
被折叠的 条评论
为什么被折叠?
到【灌水乐园】发言
