{"metadata":{"kernelspec":{"language":"python","display_name":"Python 3","name":"python3"},"language_info":{"name":"python","version":"3.10.12","mimetype":"text/x-python","codemirror_mode":{"name":"ipython","version":3},"pygments_lexer":"ipython3","nbconvert_exporter":"python","file_extension":".py"},"kaggle":{"accelerator":"gpu","dataSources":[{"sourceId":36363,"databundleVersionId":4050810,"sourceType":"competition"},{"sourceId":10745238,"sourceType":"datasetVersion","datasetId":6663719}],"dockerImageVersionId":30887,"isInternetEnabled":true,"language":"python","sourceType":"notebook","isGpuEnabled":true}},"nbformat_minor":4,"nbformat":4,"cells":[{"cell_type":"code","source":"import os\nimport numpy as np\nimport pandas as pd\nimport tensorflow as tf\nimport SimpleITK as sitk\nimport cv2\nfrom tensorflow.keras import layers, models, callbacks\n\n# ====================================================\n# (A) 讀取原始 CSV 並自動衍生 is_cervical 標籤\n# ====================================================\n# 路徑請根據你的環境調整\ntrain_images_dir = \"/kaggle/input/rsna-2022-cervical-spine-fracture-detection/train_images\"\nlabels_csv_path = \"/kaggle/input/vertebrae-labels/vertebrae_labels.csv\"  # 原始 CSV，包含 Label1~Label19 及 Fracture 標籤\n\ndf = pd.read_csv(labels_csv_path)\nprint(\"原始資料筆數:\", len(df))\n\n# 利用 Label1～Label7（假設這幾個欄位代表頸椎各部位）衍生 is_cervical\nlabel_cols = [\"Label1\", \"Label2\", \"Label3\", \"Label4\", \"Label5\", \"Label6\", \"Label7\"]\ndf[\"is_cervical\"] = df[label_cols].max(axis=1)\nprint(\"is_cervical 分布:\")\nprint(df[\"is_cervical\"].value_counts())\n\n# ====================================================\n# (B) 過濾出實際存在的 DICOM 檔案\n# ====================================================\nvalid_rows = []\nfor i, row in df.iterrows():\n    spine_id = row[\"SpineID\"]\n    slice_num = int(row[\"SliceNumber\"])\n    dicom_path = os.path.join(train_images_dir, spine_id, f\"{slice_num}.dcm\")\n    if os.path.exists(dicom_path):\n        valid_rows.append(i)\ndf_filtered = df.loc[valid_rows].reset_index(drop=True)\nfiltered_csv_path = \"/kaggle/working/filtered_vertebrae_labels.csv\"\ndf_filtered.to_csv(filtered_csv_path, index=False)\nprint(\"過濾後筆數:\", len(df_filtered))\n\n# ====================================================\n# (C) Stage 1：建立頸椎分類器資料集\n# ====================================================\n# 取出所有影像與 is_cervical 標籤\nspine_ids_stage1 = df_filtered[\"SpineID\"].astype(str).values\nslice_nums_stage1 = df_filtered[\"SliceNumber\"].astype(np.int32).values\nis_cervical_labels = df_filtered[\"is_cervical\"].values.astype(np.float32)\n\n# 定義讀取影像函式（共用）\ndef load_dicom_image(spine_id, slice_number):\n    spine_id = spine_id.numpy().decode(\"utf-8\")\n    slice_number = int(slice_number.numpy())\n    dicom_path = os.path.join(train_images_dir, spine_id, f\"{slice_number}.dcm\")\n    try:\n        dicom_data = sitk.ReadImage(dicom_path)\n        image = sitk.GetArrayFromImage(dicom_data)[0]\n        image = cv2.normalize(image, None, 0, 255, cv2.NORM_MINMAX).astype(np.uint8)\n        image = cv2.resize(image, (256,256))\n        if len(image.shape) == 2:\n            image = cv2.cvtColor(image, cv2.COLOR_GRAY2RGB)\n    except:\n        image = np.zeros((256,256,3), dtype=np.uint8)\n    return image\n\ndef load_sample_stage1(spine_id, slice_number, is_cervical):\n    image = tf.py_function(func=load_dicom_image, inp=[spine_id, slice_number], Tout=tf.uint8)\n    image = tf.cast(image, tf.float32) / 255.0\n    image.set_shape([256,256,3])\n    is_cervical = tf.cast(is_cervical, tf.float32)\n    return image, is_cervical\n\n# 在 NumPy 層面切分 80% 訓練 / 20% 驗證\ndataset_size_stage1 = len(spine_ids_stage1)\ntrain_size_stage1 = int(0.8 * dataset_size_stage1)\n\ntrain_spine_ids_stage1 = spine_ids_stage1[:train_size_stage1]\ntrain_slice_nums_stage1 = slice_nums_stage1[:train_size_stage1]\ntrain_is_cervical = is_cervical_labels[:train_size_stage1]\n\nval_spine_ids_stage1 = spine_ids_stage1[train_size_stage1:]\nval_slice_nums_stage1 = slice_nums_stage1[train_size_stage1:]\nval_is_cervical = is_cervical_labels[train_size_stage1:]\n\ntrain_dataset_stage1 = tf.data.Dataset.from_tensor_slices((train_spine_ids_stage1, train_slice_nums_stage1, train_is_cervical))\ntrain_dataset_stage1 = train_dataset_stage1.map(load_sample_stage1, num_parallel_calls=tf.data.AUTOTUNE)\ntrain_dataset_stage1 = train_dataset_stage1.shuffle(1000).batch(32).prefetch(tf.data.AUTOTUNE)\n\nval_dataset_stage1 = tf.data.Dataset.from_tensor_slices((val_spine_ids_stage1, val_slice_nums_stage1, val_is_cervical))\nval_dataset_stage1 = val_dataset_stage1.map(load_sample_stage1, num_parallel_calls=tf.data.AUTOTUNE)\nval_dataset_stage1 = val_dataset_stage1.batch(32).prefetch(tf.data.AUTOTUNE)\n\nprint(\"Stage 1 - 訓練集 batch 數:\", tf.data.experimental.cardinality(train_dataset_stage1).numpy())\nprint(\"Stage 1 - 驗證集 batch 數:\", tf.data.experimental.cardinality(val_dataset_stage1).numpy())\n\n# --------------------------------------------------\n# Stage 1: 建立頸椎分類器模型 (二分類)\n# --------------------------------------------------\ndef build_cervical_classifier(input_shape=(256,256,3)):\n    inputs = layers.Input(shape=input_shape)\n    x = layers.Conv2D(16, (3,3), activation='relu', padding='same')(inputs)\n    x = layers.MaxPooling2D()(x)\n    x = layers.Conv2D(32, (3,3), activation='relu', padding='same')(x)\n    x = layers.MaxPooling2D()(x)\n    x = layers.Conv2D(64, (3,3), activation='relu', padding='same')(x)\n    x = layers.GlobalAveragePooling2D()(x)\n    outputs = layers.Dense(1, activation='sigmoid')(x)\n    model = models.Model(inputs, outputs)\n    return model\n\ncervical_model = build_cervical_classifier()\ncervical_model.compile(optimizer='adam', loss='binary_crossentropy', metrics=['accuracy'])\ncervical_model.summary()\n\n# 訓練 Stage 1 模型\ncervical_history = cervical_model.fit(train_dataset_stage1, epochs=10, validation_data=val_dataset_stage1)\n\n# --------------------------------------------------\n# (D) Stage 2：建立頸椎骨折預測模型資料集\n# --------------------------------------------------\n# 從 df_filtered 中只取出 is_cervical==1 的資料\ndf_fracture = df_filtered[df_filtered[\"is_cervical\"] == 1].reset_index(drop=True)\nprint(\"Stage 2: 頸椎骨折資料筆數:\", len(df_fracture))\n\ntarget_fracture_cols = [\"Fracture_C1\",\"Fracture_C2\",\"Fracture_C3\",\"Fracture_C4\",\"Fracture_C5\",\"Fracture_C6\",\"Fracture_C7\"]\nspine_ids_fracture = df_fracture[\"SpineID\"].astype(str).values\nslice_nums_fracture = df_fracture[\"SliceNumber\"].astype(np.int32).values\nfracture_labels = df_fracture[target_fracture_cols].values.astype(np.float32)\n\ndataset_size_fracture = len(spine_ids_fracture)\ntrain_size_fracture = int(0.8 * dataset_size_fracture)\n\ntrain_spine_ids_fracture = spine_ids_fracture[:train_size_fracture]\ntrain_slice_nums_fracture = slice_nums_fracture[:train_size_fracture]\ntrain_fracture_labels = fracture_labels[:train_size_fracture]\n\nval_spine_ids_fracture = spine_ids_fracture[train_size_fracture:]\nval_slice_nums_fracture = slice_nums_fracture[train_size_fracture:]\nval_fracture_labels = fracture_labels[train_size_fracture:]\n\ndef load_sample_fracture(spine_id, slice_number, label):\n    image = tf.py_function(func=load_dicom_image, inp=[spine_id, slice_number], Tout=tf.uint8)\n    image = tf.cast(image, tf.float32) / 255.0\n    image.set_shape([256,256,3])\n    label = tf.cast(label, tf.float32)\n    return image, label\n\ntrain_dataset_fracture = tf.data.Dataset.from_tensor_slices((train_spine_ids_fracture, train_slice_nums_fracture, train_fracture_labels))\ntrain_dataset_fracture = train_dataset_fracture.map(load_sample_fracture, num_parallel_calls=tf.data.AUTOTUNE)\ntrain_dataset_fracture = train_dataset_fracture.shuffle(1000).batch(32).prefetch(tf.data.AUTOTUNE)\n\nval_dataset_fracture = tf.data.Dataset.from_tensor_slices((val_spine_ids_fracture, val_slice_nums_fracture, val_fracture_labels))\nval_dataset_fracture = val_dataset_fracture.map(load_sample_fracture, num_parallel_calls=tf.data.AUTOTUNE)\nval_dataset_fracture = val_dataset_fracture.batch(32).prefetch(tf.data.AUTOTUNE)\n\nprint(\"Stage 2 - 訓練集 batch 數:\", tf.data.experimental.cardinality(train_dataset_fracture).numpy())\nprint(\"Stage 2 - 驗證集 batch 數:\", tf.data.experimental.cardinality(val_dataset_fracture).numpy())\n\n# --------------------------------------------------\n# Stage 2: 建立改進版 UNet 骨折預測模型 (7 標籤)\n# --------------------------------------------------\ndef build_improved_unet_classifier(input_shape=(256,256,3), num_classes=7):\n    inputs = layers.Input(shape=input_shape)\n    # Encoder Block 1\n    c1 = layers.Conv2D(32, (3,3), activation='relu', padding='same')(inputs)\n    c1 = layers.BatchNormalization()(c1)\n    c1 = layers.Conv2D(32, (3,3), activation='relu', padding='same')(c1)\n    c1 = layers.BatchNormalization()(c1)\n    p1 = layers.MaxPooling2D((2,2))(c1)\n    \n    # Encoder Block 2\n    c2 = layers.Conv2D(64, (3,3), activation='relu', padding='same')(p1)\n    c2 = layers.BatchNormalization()(c2)\n    c2 = layers.Conv2D(64, (3,3), activation='relu', padding='same')(c2)\n    c2 = layers.BatchNormalization()(c2)\n    p2 = layers.MaxPooling2D((2,2))(c2)\n    \n    # Encoder Block 3\n    c3 = layers.Conv2D(128, (3,3), activation='relu', padding='same')(p2)\n    c3 = layers.BatchNormalization()(c3)\n    c3 = layers.Conv2D(128, (3,3), activation='relu', padding='same')(c3)\n    c3 = layers.BatchNormalization()(c3)\n    p3 = layers.MaxPooling2D((2,2))(c3)\n    \n    # Encoder Block 4\n    c4 = layers.Conv2D(256, (3,3), activation='relu', padding='same')(p3)\n    c4 = layers.BatchNormalization()(c4)\n    c4 = layers.Conv2D(256, (3,3), activation='relu', padding='same')(c4)\n    c4 = layers.BatchNormalization()(c4)\n    p4 = layers.MaxPooling2D((2,2))(c4)\n    \n    # Bottleneck\n    c5 = layers.Conv2D(512, (3,3), activation='relu', padding='same')(p4)\n    c5 = layers.BatchNormalization()(c5)\n    c5 = layers.Conv2D(512, (3,3), activation='relu', padding='same')(c5)\n    c5 = layers.BatchNormalization()(c5)\n    \n    # Decoder Block 1\n    u6 = layers.UpSampling2D((2,2))(c5)\n    u6 = layers.Concatenate()([u6, c4])\n    c6 = layers.Conv2D(256, (3,3), activation='relu', padding='same')(u6)\n    c6 = layers.BatchNormalization()(c6)\n    c6 = layers.Conv2D(256, (3,3), activation='relu', padding='same')(c6)\n    c6 = layers.BatchNormalization()(c6)\n    \n    # Decoder Block 2\n    u7 = layers.UpSampling2D((2,2))(c6)\n    u7 = layers.Concatenate()([u7, c3])\n    c7 = layers.Conv2D(128, (3,3), activation='relu', padding='same')(u7)\n    c7 = layers.BatchNormalization()(c7)\n    c7 = layers.Conv2D(128, (3,3), activation='relu', padding='same')(c7)\n    c7 = layers.BatchNormalization()(c7)\n    \n    # Decoder Block 3\n    u8 = layers.UpSampling2D((2,2))(c7)\n    u8 = layers.Concatenate()([u8, c2])\n    c8 = layers.Conv2D(64, (3,3), activation='relu', padding='same')(u8)\n    c8 = layers.BatchNormalization()(c8)\n    c8 = layers.Conv2D(64, (3,3), activation='relu', padding='same')(c8)\n    c8 = layers.BatchNormalization()(c8)\n    \n    # Decoder Block 4\n    u9 = layers.UpSampling2D((2,2))(c8)\n    u9 = layers.Concatenate()([u9, c1])\n    c9 = layers.Conv2D(32, (3,3), activation='relu', padding='same')(u9)\n    c9 = layers.BatchNormalization()(c9)\n    c9 = layers.Conv2D(32, (3,3), activation='relu', padding='same')(c9)\n    c9 = layers.BatchNormalization()(c9)\n    \n    gap = layers.GlobalAveragePooling2D()(c9)\n    d1 = layers.Dense(64, activation='relu')(gap)\n    d1 = layers.Dropout(0.5)(d1)\n    outputs = layers.Dense(num_classes, activation='sigmoid')(d1)\n    \n    model = models.Model(inputs=inputs, outputs=outputs)\n    return model\n\nfracture_model = build_improved_unet_classifier()\nfracture_model.compile(optimizer='adam', loss='binary_crossentropy', metrics=[tf.keras.metrics.BinaryAccuracy(name='binary_accuracy')])\nfracture_model.summary()\n\n# 訓練 Stage 2 模型\nfracture_history = fracture_model.fit(train_dataset_fracture, epochs=50, validation_data=val_dataset_fracture)\n\n# --------------------------------------------------\n# (E) 最終應用：兩階段推論函式，並顯示準確率（機率）\n# --------------------------------------------------\ndef predict_image(spine_id, slice_number):\n    # 讀取影像\n    image = load_dicom_image(tf.constant(spine_id), tf.constant(slice_number))\n    image = image.astype(np.float32) / 255.0\n    image_exp = np.expand_dims(image, axis=0)\n    \n    # 第一階段：使用頸椎分類器判斷是否為頸椎\n    cervical_pred = cervical_model.predict(image_exp)\n    cervical_prob = cervical_pred[0,0]\n    is_cervical = cervical_prob > 0.5\n    print(f\"頸椎分類器信心：{cervical_prob:.2f}\")\n    \n    if is_cervical:\n        # 如果是頸椎，使用骨折預測模型\n        fracture_pred = fracture_model.predict(image_exp)\n        # fracture_pred 為 7 個標籤的機率\n        fracture_status = (fracture_pred[0] > 0.5).astype(int)\n        print(\"該影像為頸椎。\")\n        print(\"各椎體骨折預測 (機率)：\", np.round(fracture_pred[0], 2))\n        print(\"各椎體骨折狀態 (0/1)：\", fracture_status)\n        return fracture_status, fracture_pred[0]\n    else:\n        print(\"該影像並非頸椎影像。\")\n        return None, None\n\n# 測試推論 (請根據你的資料替換 SpineID 與 SliceNumber)\npredict_image(\"1.2.826.0.1.3680043.780\", 5)\n","metadata":{"_uuid":"8f2839f25d086af736a60e9eeb907d3b93b6e0e5","_cell_guid":"b1076dfc-b9ad-4769-8c92-a6c4dae69d19","trusted":true,"execution":{"iopub.status.busy":"2025-02-17T16:21:24.216511Z","iopub.execute_input":"2025-02-17T16:21:24.216797Z","iopub.status.idle":"2025-02-17T16:24:36.809873Z","shell.execute_reply.started":"2025-02-17T16:21:24.216774Z","shell.execute_reply":"2025-02-17T16:24:36.808481Z"}},"outputs":[],"execution_count":null}]}