AI在医疗诊断中的应用与伦理挑战_03

import numpy as np
import pandas as pd
import matplotlib.pyplot as plt
import seaborn as sns
import os
import random
import cv2
import tensorflow as tf
from tensorflow.keras import layers, models, optimizers
from tensorflow.keras.preprocessing.image import ImageDataGenerator
from tensorflow.keras.callbacks import ModelCheckpoint, ReduceLROnPlateau, EarlyStopping
from tensorflow.keras.applications import ResNet50
from sklearn.model_selection import train_test_split
from sklearn.metrics import classification_report, confusion_matrix, roc_auc_score
from sklearn.utils.class_weight import compute_class_weight
import tensorflow_addons as tfa
import albumentations as A
from albumentations.tensorflow import ToTensorV2
import warnings
warnings.filterwarnings('ignore')

# 设置中文显示
plt.rcParams["font.family"] = ["SimHei", "WenQuanYi Micro Hei", "Heiti TC"]
plt.rcParams['axes.unicode_minus'] = False

# 随机种子设置，确保结果可复现
def set_seed(seed=42):
    random.seed(seed)
    np.random.seed(seed)
    tf.random.set_seed(seed)
    os.environ['PYTHONHASHSEED'] = str(seed)
    os.environ['TF_DETERMINISTIC_OPS'] = '1'

set_seed(42)

# 数据路径设置
DATA_DIR = "data/diabetic_retinopathy"
TRAIN_DIR = os.path.join(DATA_DIR, "train")
TEST_DIR = os.path.join(DATA_DIR, "test")
CSV_PATH = os.path.join(DATA_DIR, "train_labels.csv")

# 图像尺寸和训练参数
IMG_SIZE = 224
BATCH_SIZE = 32
EPOCHS = 30
LEARNING_RATE = 1e-4
NUM_CLASSES = 5  # DR分级：0-4级

# 加载并预处理数据
def load_data(csv_path, image_dir):
    df = pd.read_csv(csv_path)
    # 确保图像路径正确
    df['image_path'] = df['image_id'].apply(lambda x: os.path.join(image_dir, f"{x}.png"))
    # 检查图像文件是否存在
    df = df[df['image_path'].apply(os.path.exists)]
    return df[['image_path', 'diagnosis']]

# 数据增强策略 - 使用albumentations库
train_transform = A.Compose([
    A.RandomResizedCrop(height=IMG_SIZE, width=IMG_SIZE, scale=(0.8, 1.0)),
    A.HorizontalFlip(p=0.5),
    A.VerticalFlip(p=0.3),
    A.RandomRotate90(p=0.5),
    A.RandomBrightnessContrast(p=0.3, brightness_limit=0.2, contrast_limit=0.2),
    A.HueSaturationValue(p=0.2, hue_shift_limit=10, sat_shift_limit=10),
    A.GaussNoise(p=0.2),
    A.OneOf([
        A.MotionBlur(p=0.2),
        A.MedianBlur(p=0.1),
        A.GaussianBlur(p=0.1),
    ], p=0.2),
    A.OneOf([
        A.OpticalDistortion(p=0.3),
        A.GridDistortion(p=0.1),
    ], p=0.2),
    A.CLAHE(p=0.3),
    A.Normalize(),
    ToTensorV2()
])

val_transform = A.Compose([
    A.Resize(height=IMG_SIZE, width=IMG_SIZE),
    A.Normalize(),
    ToTensorV2()
])

# 自定义数据生成器
class DRDataset(tf.keras.utils.Sequence):
    def __init__(self, dataframe, transform=None, batch_size=32, shuffle=True):
        self.dataframe = dataframe
        self.transform = transform
        self.batch_size = batch_size
        self.shuffle = shuffle
        self.indices = self.dataframe.index.tolist()
        self.on_epoch_end()

    def __len__(self):
        return len(self.indices) // self.batch_size

    def __getitem__(self, index):
        batch_indices = self.indices[index * self.batch_size : (index + 1) * self.batch_size]
        batch_df = self.dataframe.loc[batch_indices]

        X = []
        y = []

        for _, row in batch_df.iterrows():
            image_path = row['image_path']
            diagnosis = row['diagnosis']

            # 读取图像
            image = cv2.imread(image_path)
            image = cv2.cvtColor(image, cv2.COLOR_BGR2RGB)

            # 应用变换
            if self.transform:
                image = self.transform(image=image)['image']

            X.append(image)
            y.append(diagnosis)

        return np.array(X), tf.keras.utils.to_categorical(y, num_classes=NUM_CLASSES)

    def on_epoch_end(self):
        if self.shuffle:
            random.shuffle(self.indices)

# 加载数据
df = load_data(CSV_PATH, TRAIN_DIR)

train_df, val_df = train_test_split(df, test_size=0.2, random_state=42, stratify=df['diagnosis'])

train_dataset = DRDataset(train_df, transform=train_transform, batch_size=BATCH_SIZE)
val_dataset = DRDataset(val_df, transform=val_transform, batch_size=BATCH_SIZE, shuffle=False)

# 计算类别权重以处理不平衡数据
class_weights = compute_class_weight('balanced', classes=np.unique(df['diagnosis']), y=df['diagnosis'])
class_weights = {i: class_weights[i] for i in range(NUM_CLASSES)}

# 构建模型 - 基于ResNet50迁移学习
def build_model(input_shape=(IMG_SIZE, IMG_SIZE, 3), num_classes=NUM_CLASSES):
    # 加载预训练模型
    base_model = ResNet50(weights='imagenet', include_top=False, input_shape=input_shape)

    # 冻结基础模型初始层
    for layer in base_model.layers[:-4]:
        layer.trainable = False

    # 添加分类头部
    model = models.Sequential([
        base_model,
        layers.GlobalAveragePooling2D(),
        layers.BatchNormalization(),
        layers.Dense(512, activation='relu'),
        layers.BatchNormalization(),
        layers.Dropout(0.5),
        layers.Dense(256, activation='relu'),
        layers.BatchNormalization(),
        layers.Dropout(0.3),
        layers.Dense(num_classes, activation='softmax')
    ])

    # 编译模型
    optimizer = optimizers.Adam(learning_rate=LEARNING_RATE)
    model.compile(
        optimizer=optimizer,
        loss='categorical_crossentropy',
        metrics=[
            'accuracy',
            tf.keras.metrics.AUC(multi_label=True, name='auc'),
            tfa.metrics.F1Score(num_classes=num_classes, average='macro', name='f1_score')
        ]
    )

    return model

# 创建模型
model = build_model()
model.summary()

# 定义回调函数
callbacks = [
    ModelCheckpoint(
        'best_dr_model.h5',
        monitor='val_f1_score',
        mode='max',
        save_best_only=True,
        save_weights_only=False,
        verbose=1
    ),
    ReduceLROnPlateau(
        monitor='val_loss',
        factor=0.5,
        patience=5,
        min_lr=1e-6,
        verbose=1
    ),
    EarlyStopping(
        monitor='val_f1_score',
        mode='max',
        patience=10,
        restore_best_weights=True,
        verbose=1
    )
]

# 训练模型
history = model.fit(
    train_dataset,
    validation_data=val_dataset,
    epochs=EPOCHS,
    class_weight=class_weights,
    callbacks=callbacks,
    verbose=1
)

# 绘制训练曲线
def plot_training_history(history):
    plt.figure(figsize=(15, 5))

    # 准确率
    plt.subplot(1, 3, 1)
    plt.plot(history.history['accuracy'], label='训练准确率')
    plt.plot(history.history['val_accuracy'], label='验证准确率')
    plt.title('模型准确率')
    plt.xlabel('Epoch')
    plt.ylabel('准确率')
    plt.legend()

    # 损失
    plt.subplot(1, 3, 2)
    plt.plot(history.history['loss'], label='训练损失')
    plt.plot(history.history['val_loss'], label='验证损失')
    plt.title('模型损失')
    plt.xlabel('Epoch')
    plt.ylabel('损失')
    plt.legend()

    # F1分数
    plt.subplot(1, 3, 3)
    plt.plot(history.history['f1_score'], label='训练F1分数')
    plt.plot(history.history['val_f1_score'], label='验证F1分数')
    plt.title('模型F1分数')
    plt.xlabel('Epoch')
    plt.ylabel('F1分数')
    plt.legend()

    plt.tight_layout()
    plt.savefig('training_history.png')
    plt.close()

plot_training_history(history)

# 评估模型
def evaluate_model(model, val_dataset):
    # 获取预测结果
    y_pred = []
    y_true = []

    for x, y in val_dataset:
        pred = model.predict(x, verbose=0)
        y_pred.extend(np.argmax(pred, axis=1))
        y_true.extend(np.argmax(y, axis=1))

    # 计算评估指标
    print("\n分类报告:")
    print(classification_report(y_true, y_pred))

    # 计算AUC
    auc = roc_auc_score(y_true, model.predict(val_dataset), multi_class='ovr')
    print(f"AUC: {auc:.4f}")

    # 混淆矩阵
    cm = confusion_matrix(y_true, y_pred)
    plt.figure(figsize=(10, 8))
    sns.heatmap(cm, annot=True, fmt='d', cmap='Blues', cbar=False)
    plt.title('糖尿病视网膜病变分级混淆矩阵')
    plt.xlabel('预测类别')
    plt.ylabel('真实类别')
    plt.savefig('confusion_matrix.png')
    plt.close()

    return {
        'accuracy': np.mean(np.array(y_true) == np.array(y_pred)),
        'auc': auc,
        'classification_report': classification_report(y_true, y_pred, output_dict=True)
    }

# 评估最佳模型
best_model = tf.keras.models.load_model('best_dr_model.h5')
evaluation_results = evaluate_model(best_model, val_dataset)
print(f"最佳模型验证准确率: {evaluation_results['accuracy']:.4f}")
print(f"最佳模型AUC: {evaluation_results['auc']:.4f}")

# 模型解释性 - Grad-CAM可视化
def grad_cam(model, img_array, layer_name='conv5_block3_out', eps=1e-8):
    # 创建模型，输出最后一个卷积层和预测层
    grad_model = tf.keras.models.Model(
        [model.inputs], [model.get_layer(layer_name).output, model.output]
    )

    # 计算梯度
    with tf.GradientTape() as tape:
        conv_outputs, predictions = grad_model(img_array)
        class_idx = tf.argmax(predictions[0])
        loss = predictions[:, class_idx]

    # 获取梯度
    output = conv_outputs[0]
    grads = tape.gradient(loss, conv_outputs)[0]

    # 全局平均池化
    pooled_grads = tf.reduce_mean(grads, axis=(0, 1))

    # 加权组合
    heatmap = output @ pooled_grads[..., tf.newaxis]
    heatmap = tf.squeeze(heatmap)

    # ReLU激活
    heatmap = tf.maximum(heatmap, 0) / (tf.reduce_max(heatmap) + eps)

    return heatmap.numpy()

# 可视化样本预测和Grad-CAM
def visualize_predictions(model, val_df, num_samples=5):
    sample_df = val_df.sample(num_samples, random_state=42)

    plt.figure(figsize=(15, 3*num_samples))

    for i, (_, row) in enumerate(sample_df.iterrows()):
        img_path = row['image_path']
        true_label = row['diagnosis']

        # 读取和预处理图像
        img = cv2.imread(img_path)
        img = cv2.cvtColor(img, cv2.COLOR_BGR2RGB)
        img_resized = cv2.resize(img, (IMG_SIZE, IMG_SIZE))
        img_array = val_transform(image=img_resized)['image']
        img_array = np.expand_dims(img_array, axis=0)

        # 预测
        pred_probs = model.predict(img_array)[0]
        pred_label = np.argmax(pred_probs)

        # 生成Grad-CAM热图
        heatmap = grad_cam(model, img_array)
        heatmap = cv2.resize(heatmap, (img.shape[1], img.shape[0]))
        heatmap = np.uint8(255 * heatmap)
        heatmap = cv2.applyColorMap(heatmap, cv2.COLORMAP_JET)

        # 叠加热图到原图
        superimposed_img = heatmap * 0.4 + img
        superimposed_img = np.uint8(superimposed_img)

        # 绘制结果
        plt.subplot(num_samples, 3, i*3+1)
        plt.imshow(img)
        plt.title(f'原始图像\n真实标签: {true_label}')
        plt.axis('off')

        plt.subplot(num_samples, 3, i*3+2)
        plt.imshow(superimposed_img)
        plt.title(f'Grad-CAM可视化\n预测标签: {pred_label} (置信度: {pred_probs[pred_label]:.2f})')
        plt.axis('off')

        plt.subplot(num_samples, 3, i*3+3)
        plt.bar(range(NUM_CLASSES), pred_probs)
        plt.title('预测概率分布')
        plt.xlabel('DR等级')
        plt.ylabel('概率')
        plt.xticks(range(NUM_CLASSES))

    plt.tight_layout()
    plt.savefig('dr_predictions_visualization.png')
    plt.close()

# 可视化预测结果
visualize_predictions(best_model, val_df)

print("模型训练和评估完成！")
print(f"最佳模型保存在: best_dr_model.h5")
print(f"评估指标: {evaluation_results}")