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

import os
import numpy as np
import pandas as pd
import matplotlib.pyplot as plt
import tensorflow as tf
from tensorflow.keras.models import Sequential, Model
from tensorflow.keras.layers import Conv2D, MaxPooling2D, Dropout, Flatten, Dense, BatchNormalization, Input, concatenate
from tensorflow.keras.preprocessing.image import ImageDataGenerator
from tensorflow.keras.optimizers import Adam, SGD
from tensorflow.keras.callbacks import EarlyStopping, ModelCheckpoint, ReduceLROnPlateau
from sklearn.model_selection import train_test_split
from sklearn.metrics import classification_report, confusion_matrix, roc_curve, auc
import cv2
import random
import shutil
from pathlib import Path
import time

# 设置中文字体
plt.rcParams["font.family"] = ["SimHei", "WenQuanYi Micro Hei", "Heiti TC"]
plt.rcParams["axes.unicode_minus"] = False  # 解决负号显示问题

def load_and_preprocess_image(image_path, img_size=(224, 224)):
    """加载并预处理医学图像"""
    # 读取图像
    img = cv2.imread(image_path)
    # 转换为RGB格式（OpenCV默认读取为BGR）
    img = cv2.cvtColor(img, cv2.COLOR_BGR2RGB)
    # 调整图像大小
    img = cv2.resize(img, img_size)
    # 归一化像素值
    img = img / 255.0
    # 数据增强（可选）
    # 这里可以添加旋转、翻转、缩放等数据增强操作
    return img

def create_simple_cnn_model(input_shape=(224, 224, 3), num_classes=2):
    """创建一个简单的CNN模型用于医学图像分类"""
    model = Sequential([
        # 第一个卷积块
        Conv2D(32, (3, 3), activation='relu', padding='same', input_shape=input_shape),
        BatchNormalization(),
        Conv2D(32, (3, 3), activation='relu', padding='same'),
        BatchNormalization(),
        MaxPooling2D((2, 2)),
        Dropout(0.25),
        
        # 第二个卷积块
        Conv2D(64, (3, 3), activation='relu', padding='same'),
        BatchNormalization(),
        Conv2D(64, (3, 3), activation='relu', padding='same'),
        BatchNormalization(),
        MaxPooling2D((2, 2)),
        Dropout(0.25),
        
        # 第三个卷积块
        Conv2D(128, (3, 3), activation='relu', padding='same'),
        BatchNormalization(),
        Conv2D(128, (3, 3), activation='relu', padding='same'),
        BatchNormalization(),
        MaxPooling2D((2, 2)),
        Dropout(0.25),
        
        # 全连接层
        Flatten(),
        Dense(256, activation='relu'),
        BatchNormalization(),
        Dropout(0.5),
        Dense(num_classes, activation='softmax' if num_classes > 1 else 'sigmoid')
    ])
    
    # 编译模型
    optimizer = Adam(learning_rate=0.001)
    if num_classes > 1:
        model.compile(optimizer=optimizer, loss='categorical_crossentropy', metrics=['accuracy'])
    else:
        model.compile(optimizer=optimizer, loss='binary_crossentropy', metrics=['accuracy'])
    
    return model

def create_unet_model(input_shape=(256, 256, 1)):
    """创建U-Net模型用于医学图像分割"""
    # 输入层
    inputs = Input(input_shape)
    
    # 编码器部分
    def conv_block(x, filters, kernel_size=3, activation='relu'):
        x = Conv2D(filters, kernel_size, activation=activation, padding='same')(x)
        x = BatchNormalization()(x)
        x = Conv2D(filters, kernel_size, activation=activation, padding='same')(x)
        x = BatchNormalization()(x)
        return x
    
    # 下采样
    c1 = conv_block(inputs, 64)
    p1 = MaxPooling2D((2, 2))(c1)
    
    c2 = conv_block(p1, 128)
    p2 = MaxPooling2D((2, 2))(c2)
    
    c3 = conv_block(p2, 256)
    p3 = MaxPooling2D((2, 2))(c3)
    
    c4 = conv_block(p3, 512)
    p4 = MaxPooling2D((2, 2))(c4)
    
    # 瓶颈层
    c5 = conv_block(p4, 1024)
    
    # 解码器部分
    def up_conv_block(x, skip, filters, kernel_size=3, activation='relu'):
        x = tf.keras.layers.UpSampling2D((2, 2))(x)
        x = concatenate([x, skip])
        x = Conv2D(filters, kernel_size, activation=activation, padding='same')(x)
        x = BatchNormalization()(x)
        x = Conv2D(filters, kernel_size, activation=activation, padding='same')(x)
        x = BatchNormalization()(x)
        return x
    
    # 上采样
    u6 = up_conv_block(c5, c4, 512)
    u7 = up_conv_block(u6, c3, 256)
    u8 = up_conv_block(u7, c2, 128)
    u9 = up_conv_block(u8, c1, 64)
    
    # 输出层
    outputs = Conv2D(1, (1, 1), activation='sigmoid')(u9)
    
    # 创建模型
    model = Model(inputs=inputs, outputs=outputs)
    
    # 编译模型
    model.compile(optimizer=Adam(learning_rate=0.001), loss='binary_crossentropy', metrics=['accuracy'])
    
    return model

def train_model(model, train_data, val_data, epochs=50, batch_size=32, model_save_path='model.h5'):
    """训练深度学习模型"""
    # 设置回调函数
    early_stopping = EarlyStopping(monitor='val_loss', patience=10, restore_best_weights=True)
    model_checkpoint = ModelCheckpoint(model_save_path, monitor='val_loss', save_best_only=True)
    reduce_lr = ReduceLROnPlateau(monitor='val_loss', factor=0.5, patience=5, min_lr=1e-6)
    
    # 训练模型
    history = model.fit(
        train_data,
        validation_data=val_data,
        epochs=epochs,
        batch_size=batch_size,
        callbacks=[early_stopping, model_checkpoint, reduce_lr],
        verbose=1
    )
    
    return history

def evaluate_model(model, test_data, class_names=None):
    """评估模型性能"""
    # 获取预测结果
    y_pred = model.predict(test_data)
    
    # 处理预测结果
    if y_pred.shape[-1] > 1:
        # 多分类问题
        y_pred_classes = np.argmax(y_pred, axis=1)
    else:
        # 二分类问题
        y_pred_classes = (y_pred > 0.5).astype(int).flatten()
    
    # 获取真实标签
    y_true = np.concatenate([y for x, y in test_data], axis=0)
    if y_true.shape[-1] > 1:
        y_true_classes = np.argmax(y_true, axis=1)
    else:
        y_true_classes = y_true.flatten()
    
    # 计算性能指标
    print("分类报告:")
    print(classification_report(y_true_classes, y_pred_classes, target_names=class_names))
    
    # 计算混淆矩阵
    cm = confusion_matrix(y_true_classes, y_pred_classes)
    plt.figure(figsize=(10, 8))
    plt.imshow(cm, interpolation='nearest', cmap=plt.cm.Blues)
    plt.title('混淆矩阵')
    plt.colorbar()
    tick_marks = np.arange(len(class_names))
    plt.xticks(tick_marks, class_names, rotation=45)
    plt.yticks(tick_marks, class_names)
    
    # 在混淆矩阵上标注数值
    thresh = cm.max() / 2.
    for i in range(cm.shape[0]):
        for j in range(cm.shape[1]):
            plt.text(j, i, format(cm[i, j], 'd'),
                     horizontalalignment="center",
                     color="white" if cm[i, j] > thresh else "black")
    
    plt.ylabel('真实标签')
    plt.xlabel('预测标签')
    plt.tight_layout()
    plt.savefig('confusion_matrix.png')
    plt.close()
    
    # 计算ROC曲线和AUC
    if len(class_names) == 2:
        fpr, tpr, _ = roc_curve(y_true_classes, y_pred.flatten())
        roc_auc = auc(fpr, tpr)
        
        plt.figure(figsize=(10, 8))
        plt.plot(fpr, tpr, color='darkorange', lw=2, label='ROC曲线 (面积 = %0.2f)' % roc_auc)
        plt.plot([0, 1], [0, 1], color='navy', lw=2, linestyle='--')
        plt.xlim([0.0, 1.0])
        plt.ylim([0.0, 1.05])
        plt.xlabel('假阳性率')
        plt.ylabel('真阳性率')
        plt.title('接收器操作特征 (ROC) 曲线')
        plt.legend(loc="lower right")
        plt.tight_layout()
        plt.savefig('roc_curve.png')
        plt.close()
    
    return y_pred_classes

def visualize_model_predictions(model, test_images, test_labels, class_names, num_samples=5):
    """可视化模型预测结果"""
    # 获取随机样本
    indices = random.sample(range(len(test_images)), min(num_samples, len(test_images)))
    
    plt.figure(figsize=(15, 3 * num_samples))
    
    for i, idx in enumerate(indices):
        # 预处理图像
        img = load_and_preprocess_image(test_images[idx])
        # 模型预测
        pred = model.predict(np.expand_dims(img, axis=0))
        
        # 处理预测结果
        if pred.shape[-1] > 1:
            pred_class = class_names[np.argmax(pred)]
            pred_prob = np.max(pred)
        else:
            pred_class = class_names[int(pred > 0.5)]
            pred_prob = pred[0][0] if pred[0][0] > 0.5 else 1 - pred[0][0]
        
        # 真实标签
        true_class = class_names[test_labels[idx]] if len(test_labels[idx].shape) == 0 else class_names[np.argmax(test_labels[idx])]
        
        # 绘制图像和预测结果
        plt.subplot(num_samples, 3, i*3 + 1)
        plt.imshow(img)
        plt.title(f"样本 {idx+1}")
        plt.axis('off')
        
        plt.subplot(num_samples, 3, i*3 + 2)
        plt.bar(range(len(class_names)), pred[0] if pred.shape[-1] > 1 else [1-pred[0][0], pred[0][0]])
        plt.xticks(range(len(class_names)), class_names, rotation=45)
        plt.title(f"预测概率")
        plt.ylim([0, 1])
        
        plt.subplot(num_samples, 3, i*3 + 3)
        result_text = f"预测: {pred_class} ({pred_prob:.2f})\n真实: {true_class}"
        plt.text(0.5, 0.5, result_text, fontsize=12, ha='center', va='center')
        plt.title("预测结果")
        plt.axis('off')
    
    plt.tight_layout()
    plt.savefig('predictions_visualization.png')
    plt.close()

# 示例使用
if __name__ == "__main__":
    print("AI辅助医疗诊断系统演示开始...")
    
    # 注意：由于实际医学数据的敏感性和复杂性，这里使用模拟数据进行演示
    # 在实际应用中，应使用真实的医学数据进行训练和测试
    
    # 1. 准备模拟数据
    print("\n1. 准备模拟数据...")
    
    # 创建临时目录
    temp_dir = "temp_medical_data"
    os.makedirs(temp_dir, exist_ok=True)
    
    # 创建训练、验证和测试目录
    train_dir = os.path.join(temp_dir, "train")
    val_dir = os.path.join(temp_dir, "val")
    test_dir = os.path.join(temp_dir, "test")
    
    # 创建正常和异常类别的子目录
    for class_name in ["normal", "abnormal"]:
        os.makedirs(os.path.join(train_dir, class_name), exist_ok=True)
        os.makedirs(os.path.join(val_dir, class_name), exist_ok=True)
        os.makedirs(os.path.join(test_dir, class_name), exist_ok=True)
    
    # 生成模拟医学图像数据
    num_samples_per_class = 100
    img_size = (224, 224)
    
    # 生成正常样本
    for i in range(num_samples_per_class):
        # 创建模拟的正常肺部CT图像（简单的圆形代表肺部）
        img = np.zeros((img_size[0], img_size[1], 3), dtype=np.uint8)
        cv2.circle(img, (img_size[0]//3, img_size[1]//2), img_size[0]//4, (200, 200, 200), -1)
        cv2.circle(img, (2*img_size[0]//3, img_size[1]//2), img_size[0]//4, (200, 200, 200), -1)
        
        # 添加一些随机噪声使其更真实
        noise = np.random.normal(0, 10, img.shape).astype(np.uint8)
        img = cv2.add(img, noise)
        
        # 保存图像
        if i < 70:  # 70%用于训练
            cv2.imwrite(os.path.join(train_dir, "normal", f"normal_{i}.png"), cv2.cvtColor(img, cv2.COLOR_RGB2BGR))
        elif i < 85:  # 15%用于验证
            cv2.imwrite(os.path.join(val_dir, "normal", f"normal_{i}.png"), cv2.cvtColor(img, cv2.COLOR_RGB2BGR))
        else:  # 15%用于测试
            cv2.imwrite(os.path.join(test_dir, "normal", f"normal_{i}.png"), cv2.cvtColor(img, cv2.COLOR_RGB2BGR))
    
    # 生成异常样本（在正常样本基础上添加结节）
    for i in range(num_samples_per_class):
        # 创建模拟的异常肺部CT图像
        img = np.zeros((img_size[0], img_size[1], 3), dtype=np.uint8)
        cv2.circle(img, (img_size[0]//3, img_size[1]//2), img_size[0]//4, (200, 200, 200), -1)
        cv2.circle(img, (2*img_size[0]//3, img_size[1]//2), img_size[0]//4, (200, 200, 200), -1)
        
        # 添加结节（随机位置和大小）
        nodule_x = random.randint(img_size[0]//4, 3*img_size[0]//4)
        nodule_y = random.randint(img_size[1]//4, 3*img_size[1]//4)
        nodule_radius = random.randint(5, 15)
        cv2.circle(img, (nodule_x, nodule_y), nodule_radius, (100, 100, 100), -1)
        
        # 添加一些随机噪声使其更真实
        noise = np.random.normal(0, 10, img.shape).astype(np.uint8)
        img = cv2.add(img, noise)
        
        # 保存图像
        if i < 70:  # 70%用于训练
            cv2.imwrite(os.path.join(train_dir, "abnormal", f"abnormal_{i}.png"), cv2.cvtColor(img, cv2.COLOR_RGB2BGR))
        elif i < 85:  # 15%用于验证
            cv2.imwrite(os.path.join(val_dir, "abnormal", f"abnormal_{i}.png"), cv2.cvtColor(img, cv2.COLOR_RGB2BGR))
        else:  # 15%用于测试
            cv2.imwrite(os.path.join(test_dir, "abnormal", f"abnormal_{i}.png"), cv2.cvtColor(img, cv2.COLOR_RGB2BGR))
    
    print(f"模拟数据准备完成: {num_samples_per_class} 个正常样本和 {num_samples_per_class} 个异常样本")
    
    # 2. 创建数据生成器
    print("\n2. 创建数据生成器...")
    
    train_datagen = ImageDataGenerator(
        rescale=1./255,
        rotation_range=10,
        width_shift_range=0.1,
        height_shift_range=0.1,
        shear_range=0.1,
        zoom_range=0.1,
        horizontal_flip=True,
        fill_mode='nearest'
    )
    
    val_datagen = ImageDataGenerator(rescale=1./255)
    test_datagen = ImageDataGenerator(rescale=1./255)
    
    batch_size = 16
    
    train_generator = train_datagen.flow_from_directory(
        train_dir,
        target_size=img_size,
        batch_size=batch_size,
        class_mode='binary'
    )
    
    val_generator = val_datagen.flow_from_directory(
        val_dir,
        target_size=img_size,
        batch_size=batch_size,
        class_mode='binary'
    )
    
    test_generator = test_datagen.flow_from_directory(
        test_dir,
        target_size=img_size,
        batch_size=batch_size,
        class_mode='binary',
        shuffle=False
    )
    
    # 3. 创建和训练模型
    print("\n3. 创建和训练模型...")
    
    model = create_simple_cnn_model(input_shape=(img_size[0], img_size[1], 3), num_classes=1)  # 二分类问题
    
    # 打印模型结构
    model.summary()
    
    # 训练模型
    history = train_model(
        model,
        train_generator,
        val_generator,
        epochs=20,
        batch_size=batch_size,
        model_save_path='medical_diagnosis_model.h5'
    )
    
    # 绘制训练历史
    plt.figure(figsize=(12, 5))
    
    # 绘制准确率曲线
    plt.subplot(1, 2, 1)
    plt.plot(history.history['accuracy'], label='训练准确率')
    plt.plot(history.history['val_accuracy'], label='验证准确率')
    plt.title('模型准确率')
    plt.xlabel(' epoch')
    plt.ylabel('准确率')
    plt.legend(loc='lower right')
    
    # 绘制损失曲线
    plt.subplot(1, 2, 2)
    plt.plot(history.history['loss'], label='训练损失')
    plt.plot(history.history['val_loss'], label='验证损失')
    plt.title('模型损失')
    plt.xlabel('epoch')
    plt.ylabel('损失')
    plt.legend(loc='upper right')
    
    plt.tight_layout()
    plt.savefig('training_history.png')
    plt.close()
    
    # 4. 评估模型
    print("\n4. 评估模型性能...")
    
    # 加载最佳模型
    best_model = tf.keras.models.load_model('medical_diagnosis_model.h5')
    
    # 评估模型
    evaluate_model(best_model, test_generator, class_names=['正常', '异常'])
    
    # 5. 可视化预测结果
    print("\n5. 可视化预测结果...")
    
    # 获取测试图像路径和标签
    test_image_paths = []
    test_labels = []
    
    for class_name, class_idx in test_generator.class_indices.items():
        class_dir = os.path.join(test_dir, class_name)
        for img_file in os.listdir(class_dir):
            test_image_paths.append(os.path.join(class_dir, img_file))
            test_labels.append(class_idx)
    
    # 可视化预测结果
    visualize_model_predictions(best_model, test_image_paths, test_labels, class_names=['正常', '异常'], num_samples=5)
    
    print("\nAI辅助医疗诊断系统演示完成！")
    print("注意：此演示使用模拟数据，实际应用中应使用真实的医学数据进行训练和验证。")
    print("在实际应用中，请确保遵守相关的数据隐私和医疗伦理法规。")