{"cells":[{"cell_type":"markdown","metadata":{},"source":"# 膝MRIの12所見を当てる — データを見てからベースラインを決めるまで\n\nRSNA Knee Abnormality Detection は、膝のMRI検査1件について\n**12種類の所見それぞれが「ある/なし」の確率**を出すコンペです。\nスコアは12所見の AUC の平均(macro AUC)。\n\nこのノートブックは、**データを見て分かったこと**と、\n**そこから決まったベースラインの構成**を順に書きます。\nモデルを凝る前に決めておくべきことが、このコンペにはいくつもあります。\n\n本ノートブックでお伝えすること:\n\n- このコンペの本当の難所はどこか(画像処理ではありません)\n- データを見て分かる「そのまま使ってはいけない列」\n- 検証(CV)をどう切るべきか\n- 上の全部から、ベースラインの設定がどう決まるか\n\n最後に、実際に提出して確かめた **public LB 0.898 までの道筋**を、\n段階的に説明します。\n\n## 0. まずコンペの形を押さえる\n\n| | |\n|---|---|\n| 入力 | 1件の検査(study)= 複数の撮像シリーズ(矢状断・冠状断・軸位断など) |\n| 出力 | 12所見の確率 |\n| 評価 | 12所見の AUC の平均 |\n| 形式 | Code Competition(ノートブックを提出。実行9時間以内) |\n\n12所見はこれです。\n\n`ACL`(前十字靭帯) `MCL`(内側側副靭帯) `Medial Meniscus`(内側半月板)\n`Lateral Meniscus`(外側半月板) `Medial OA` `Lateral OA` `PF OA`(3か所の変形性関節症)\n`Effusion`(関節液貯留) `Synovitis`(滑膜炎) `Baker's`(ベーカー嚢腫)\n`Contusion`(骨挫傷) `Fracture`(骨折)\n\n1件の検査(study)に何が入っていて、何を当てるのかを図にすると、こうなります。\n\n```\n【 1件の検査 (study) 】\n\n  画像 ─┬─ 矢状断 (Sagittal) のシリーズ ─┐\n        ├─ 冠状断 (Coronal)  のシリーズ ─┼─ 各シリーズ = 数十枚のスライス\n        └─ 軸位断 (Axial)    のシリーズ ─┘   (膝を輪切りにした断面画像)\n\n  文章 ─── 放射線科医が書いた所見レポート 1本\n           (英語とは限らない。9言語が混在 → §2)\n\n                          │\n                          ↓  予測する\n\n  ┌──────────────────────────────────────────┐\n  │  12所見それぞれの「ある/なし」の確率        │\n  │  ACL / MCL / Medial Meniscus / ...        │\n  └──────────────────────────────────────────┘\n```\n\n同じ膝を3つの向きから撮っているのがポイントです。\nどの所見がどの向きで見えるかは決まっており、これは後で実測します(§6)。\n\n医学用語が並びますが、中身はどれも「膝のどこが、どう傷んでいるか」です。\n1つだけ例を挙げます。**`ACL`(前十字靭帯)は、膝の中央で太ももの骨とすねの骨を\nつないでいる靭帯**で、スポーツ中の急な方向転換などで切れることがあります。\nこの所見は「その靭帯が切れているかどうか」を当てるものです。\n残り11個も同じように、**部位(靭帯・半月板・軟骨・骨・関節液)と\n傷み方の組み合わせ**になっています。\n所見の医学的な中身をこれ以上知らなくても、以降の話は追えます。\n\n### 評価指標から決まること\n\nAUC は **順位しか見ません**。「0.9 と 0.8」でも「0.6 と 0.5」でも、\n順番が同じなら同じスコアです。ここから2つ決まります。\n\n1. **確率を較正する必要がない。** 予測値を単調に変換してもスコアは変わりません\n2. **モデルを混ぜるときは順位に直してから混ぜられる。**\n   出力のスケールが違うモデル同士でも、順位にすれば公平に平均できます\n"},{"cell_type":"code","execution_count":null,"metadata":{},"outputs":[],"source":"from pathlib import Path\n\nimport numpy as np\nimport pandas as pd\n\nROOT = next(p for p in [Path(\"/kaggle/input/rsna-knee-abnormality-detection\"),\n                        Path(\"/kaggle/input/competitions/rsna-knee-abnormality-detection\")]\n            if p.exists())\n\ntrain = pd.read_csv(ROOT / \"train.csv\")\nseries = pd.read_csv(ROOT / \"train_series.csv\")\n\nTARGETS = [\"ACL\", \"MCL\", \"Medial Meniscus\", \"Lateral Meniscus\", \"Medial OA\",\n           \"Lateral OA\", \"PF OA\", \"Effusion\", \"Synovitis\", \"Baker's\",\n           \"Contusion\", \"Fracture\"]\n\nprint(f\"train   : {len(train):,} 件の検査\")\nprint(f\"series  : {len(series):,} 本のシリーズ\")\n"},{"cell_type":"markdown","metadata":{},"source":"## 1. 学習データの重要な前提 — 正解ラベルは 58 件しかない\n\n`train.csv` を開くと、すぐに分かります。\n\n**train には 4,407 件の検査がありますが、12所見のラベルが付いているのは 58 件だけです。**\n残りの 4,349 件には、**放射線科医が書いた所見レポートの原文**だけが付いています。\n\n```\ntrain 4,407 件\n  │\n  ├─    58 件 : 画像 + レポート + 【12所見の正解ラベル】\n  │              └─→ 正解があるので「作ったラベルの品質チェック」に使える\n  │\n  └─ 4,349 件 : 画像 + レポート + (ラベルなし)\n                 │\n                 └─→ このままでは学習に使えない\n                     レポートを読んでラベルを作れば使えるようになる\n```\n\nつまり、4,349 件を学習に利用するためには、\n**レポートの文章を読んで、そこから12所見の「ある/なし」を自分で決める**必要があります。\nこの作業をしない限り、学習に使えるデータは 58 件のままです。\n58 件で画像モデルを学習するのは、現実的ではありません。\n\nそのため、このコンペでは「画像から所見を当てる問題」に取りかかる前に、\n**「レポートから学習用のラベルを作る問題」**を先に解くことになります。\n手順にすると3段階です。\n\n1. レポートを読んで、4,407 件ぶんのラベルを**自分で作る**\n2. そのラベルで画像モデルを学習する\n3. 58 件の正解ラベルは、**作ったラベルの品質チェック**に使う\n\nラベル作りが適切でないと、その先の画像モデルをいくら凝っても頭打ちになります。\n実際、あとで示すように、**このコンペで一番スコアが動いたのはラベルを差し替えたとき**でした。\n"},{"cell_type":"code","execution_count":null,"metadata":{},"outputs":[],"source":"# ---- ラベルが付いているのは何件か -------------------------------------------\nlabeled = train[TARGETS].notna().all(axis=1)\nprint(f\"12所見すべてにラベルがある検査: {labeled.sum()} / {len(train):,} 件\")\nprint(f\"レポート本文がある検査        : {train['Report'].notna().sum():,} 件\")\n\nprint(\"\\n--- レポートの例(先頭200字)---\")\nprint(train.loc[train[\"Report\"].notna(), \"Report\"].iloc[0][:200])\n"},{"cell_type":"markdown","metadata":{},"source":"## 2. レポートは9言語ある\n\nレポートは英語だけではありません。英語・スペイン語・フランス語・オランダ語・\nドイツ語・トルコ語・クロアチア語・ギリシャ語・ロシア語が混ざっています。\n\nここでの設計判断は **「言語を判定してから処理を分ける」ことをしない**、です。\n\n理由は単純で、**言語判定を間違えたらそのレポートは丸ごと読めなくなる**から。\n「the が含まれるから英語」のような安い判定は、\nオランダ語やドイツ語の文章でも普通に当たってしまいます。\n\n代わりに、**9言語ぶんの手がかりを全部まとめた辞書**を作り、\nどのレポートに対しても全部を当てにいきます。\n言語を当てる必要がなくなるので、判定ミスという失敗の種が消えます。\n\n読むときに必要なのは3種類の表現です。\n\n| 種類 | 例 |\n|---|---|\n| 所見そのもの | `anterior cruciate`, `cruzado anterior`, `vorderes kreuzband`, `мениск` |\n| 否定 | `no`, `sin`, `ohne`, `unremarkable` |\n| ぼかし | `possible`, `suspected`, `cannot exclude` |\n\n否定とぼかしを拾わないと、「ACL断裂はありません」を「ACL断裂あり」と読んでしまいます。\n"},{"cell_type":"code","execution_count":null,"metadata":{},"outputs":[],"source":"# ---- 何語で書かれているか ---------------------------------------------------\n# 言語判定はしない(§2)。ここでは「英語以外が混ざっている」ことの確認だけ。\nCUE = {\n    \"英語\":       r\"\\bmeniscus\\b|\\bligament\\b\",\n    \"スペイン語\": r\"\\bmenisco\\b|\\bligamento\\b\",\n    \"フランス語\": r\"\\bm[eé]nisque\\b|ligament crois\",\n    \"オランダ語\": r\"\\bkruisband\\b\",\n    \"ドイツ語\":   r\"\\bmeniskus\\b|\\bkreuzband\\b\",\n    \"トルコ語\":   r\"men[iı]sk[uü]s|çapraz bağ\",\n    \"ロシア語\":   r\"мениск|связк\",\n    \"ギリシャ語\": r\"μηνίσκ|σύνδεσμ\",\n}\nrep = train[\"Report\"].fillna(\"\").str.lower()\nfor name, pat in CUE.items():\n    n = rep.str.contains(pat, regex=True).sum()\n    print(f\"  {name:8s} を示す語を含む: {n:5,} 件\")\nprint(\"\\n→ 英語だけではない。言語判定で分岐せず、全言語の辞書を一度に当てる(§2)\")\n"},{"cell_type":"markdown","metadata":{},"source":"## 3. 画像データの実態を数字で見る\n\n学習用の画像データ量は下記のとおりです。\n\n| | |\n|---|---|\n| DICOM ファイル数 | **819,078** |\n| 合計サイズ | **約 490 GB** |\n| ファイルサイズの中央値 | 0.526 MB |\n\n中央値 0.526 MB は、512×512 の16bit 画像を**無圧縮**で保存したときの\nサイズとぴったり一致します。実際、サンプルした 22,859 ファイルは\nすべて非圧縮(Explicit VR Little Endian)でした。\n\n**分かること**: 復号は速い(圧縮を解く処理がない)。重いのは\n「大量のファイルを開くこと」と「リサイズ」の方です。\n\n**やること**: 490 GB を毎エポック読み直すのは論外なので、\n**前処理を一度だけ実行して小さな配列にキャッシュ**します。\n学習はキャッシュだけを読みます。\n\n## 4. 用意されている列が何を表しているかを確認する\n\n`train_series.csv` には各シリーズの説明として\n`Anatomical_Plane`(撮像面)、`Fluid_Sensitive`、`Fat_Suppression` が付いています。\n\n名前からは `Fluid_Sensitive`(水を強調する撮り方か)と\n`Fat_Suppression`(脂肪を消す撮り方か)は**別の性質**に見えます。\n物理的にも独立な性質です。ところが、中身を数えると:\n\n**この2列は 24,371 本すべてで完全に一致していました(100%)。**\n片方は、もう片方の複製です。**2列あるように見えて、情報は1つ分しかありません。**\n\nさらに、その1列が何を表しているのかも調べました。\n`Fluid_Sensitive` の値ごとに、実際のシーケンス種別(T1・PD・T2 という\n3つの撮り方)を数えるとこうなります。\n\n```\nFluid_Sensitive = 0 の中身 : T1 4,219 / PD 2,361 / T2 1,782 / 不明 1,999\nFluid_Sensitive = 1 の中身 : PD 7,712 / T2 2,699 / T1     4 / 不明 3,595\n```\n\nT1・PD・T2 は**まったく写り方が違う撮り方**です。\n値が 0 のシリーズには T1 が 4,219 本・PD が 2,361 本・T2 が 1,782 本と3種類が混在し、\n値が 1 のシリーズにも PD 7,712 本と T2 2,699 本が混在しています。\n**したがって `Fluid_Sensitive` の値を見ても、そのシリーズが T1・PD・T2 の\nどれで撮られたのかは判別できません。**\nこの列から読み取れるのは脂肪抑制の有無だけで、T1/PD/T2 の区別は含まれていません。\n\nまとめると、与えられた2列から取り出せる情報は「脂肪抑制の有無」だけ。\nT1/PD/T2 の区別は**自分で復元する必要があります**。\n\n### 撮り方の種別は DICOM ヘッダから復元できる\n\nMRI は、装置が体に電波を当てて、返ってくる信号を画像にします。\nこのとき **電波を当てる間隔**(`RepetitionTime` / TR)と、\n**当ててから信号を読み取るまでの待ち時間**(`EchoTime` / TE)を変えると、\n同じ膝でも「何が白く写るか」が変わります。\nT1・PD・T2 という種別は、もともとこの2つの数値の組み合わせで決まるものなので、\nDICOM ヘッダに残っている TR と TE から逆算できます。\n\n| 条件 | 種別 | 何が見やすいか |\n|---|---|---|\n| TR < 800 | T1 | 脂肪が白い。解剖の形 |\n| TE < 50 | PD | 中間。半月板 |\n| それ以外 | T2 | 水が白い。関節液・炎症 |\n\nシリーズ名のテキストからも読めますが、**11.8% は匿名化されて\n`DummySeriesDesc!` になっている**ので、テキストだけでは 77% しか分類できません。\nTE/TR は 95.1% のシリーズに存在するので、**両方使えば 100% 分類できます**。\n\n> この 11.8% / 95.1% は**全 24,371 シリーズのヘッダを1度走査した実測値**です。\n> 下のセルは実行時間を抑えるため 300 シリーズの標本で確かめているので、\n> 表示される割合はこの値から多少ずれます。\n\n> 用意された列は、名前が示す内容と中身が一致しているとは限りません。\n> ここでは2列を突き合わせて数えたことで、\n> **片方がもう片方の複製であること**と、\n> **T1/PD/T2 の区別がどちらの列にも入っていないこと**の2点が分かりました。\n> 数えずに列名のまま特徴量として使っていたら、\n> この2点に気づかないまま学習を進めることになります。\n"},{"cell_type":"code","execution_count":null,"metadata":{},"outputs":[],"source":"# ---- 用意された列は「別の性質」を表しているか(§4)-------------------------\nct = pd.crosstab(series[\"Fluid_Sensitive\"], series[\"Fat_Suppression\"])\nprint(\"Fluid_Sensitive x Fat_Suppression:\")\nprint(ct)\n\nagree = (series[\"Fluid_Sensitive\"] == series[\"Fat_Suppression\"]).mean()\nprint(f\"\\n2列が一致する割合: {agree:.1%}\")\nprint(\"→ ほぼ同じことを言っている。コントラスト(T1/PD/T2)の情報は入っていない\")\nprint(\"   復元するには DICOM の TR/TE を見る(§4)\")\n"},{"cell_type":"code","execution_count":null,"metadata":{},"outputs":[],"source":"# ---- DICOM ヘッダを覗く(物理スケールとコントラスト)------------------------\n# 全 819,078 ファイルを開くのは高いので、シリーズごとに 1 枚だけ、\n# ランダムな 300 シリーズをサンプルする。傾向を見るにはこれで十分。\nimport pydicom\n\nsample = series.sample(min(300, len(series)), random_state=0)\n\nrows = []\nfor r in sample.itertuples():\n    d = ROOT / \"train_series\" / r.StudyInstanceUID / r.SeriesInstanceUID\n    files = list(d.glob(\"*.dcm\"))[:1] if d.exists() else []\n    if not files:\n        continue\n    try:\n        ds = pydicom.dcmread(files[0], stop_before_pixels=True)\n    except Exception:\n        continue\n    ps = getattr(ds, \"PixelSpacing\", None)\n    rows.append(dict(\n        spacing=float(ps[0]) if ps is not None else np.nan,\n        rows_px=float(getattr(ds, \"Rows\", np.nan)),\n        TR=float(getattr(ds, \"RepetitionTime\", np.nan) or np.nan),\n        TE=float(getattr(ds, \"EchoTime\", np.nan) or np.nan),\n        desc=str(getattr(ds, \"SeriesDescription\", \"\")),\n    ))\n\nhdr = pd.DataFrame(rows)\nprint(f\"読めたシリーズ: {len(hdr)}\")\nprint(\"\\n1ピクセルが表す長さ (mm):\")\nprint(hdr[\"spacing\"].describe()[[\"min\", \"25%\", \"50%\", \"75%\", \"max\"]].round(3))\nprint(\"\\n→ 施設・装置でばらつく。リサイズだけでは大きさが揃わない(§5)\")\n"},{"cell_type":"code","execution_count":null,"metadata":{},"outputs":[],"source":"# ---- TR/TE からコントラストを復元する(§4)---------------------------------\ndef contrast_of(tr, te):\n    \"\"\"TR/TE の組から T1 / PD / T2 を決める。単純な規則で十分に効く。\"\"\"\n    if not np.isfinite(tr) or not np.isfinite(te):\n        return \"不明\"\n    if tr < 800:\n        return \"T1\"\n    if te < 50:\n        return \"PD\"\n    return \"T2\"\n\n\nhdr[\"contrast\"] = [contrast_of(t, e) for t, e in zip(hdr[\"TR\"], hdr[\"TE\"])]\nprint(\"TR/TE から決めたコントラスト:\")\nprint(hdr[\"contrast\"].value_counts().to_string())\n\ndummy = hdr[\"desc\"].str.contains(\"Dummy\", case=False, na=False).mean()\nprint(f\"\\nシリーズ名が匿名化されている割合: {dummy:.1%}\")\nprint(\"→ 名前だけに頼ると分類できない。TR/TE と併用する(§4)\")\n"},{"cell_type":"markdown","metadata":{},"source":"## 5. 物理スケールを揃える — 160mm では 28.9% で効いていなかった\n\nMRI の画像は、施設や装置によって**1ピクセルが表す実際の長さが違います**。\nDICOM の `PixelSpacing` に mm/ピクセルが入っています。\n\nこれを無視して「全部 224×224 にリサイズ」すると、\n**同じ半月板が画像によって2倍以上の大きさで写る**ことになります。\nモデルにとっては別物です。\n\nそこで、**膝関節が収まる物理的な視野(mm)を決めて、そこを切り出してからリサイズ**します。\n\n```\n切り出す幅(ピクセル) = 視野(mm) ÷ PixelSpacing(mm/ピクセル)\n```\n\n最初は視野を 160mm にしました。ところが**うまくいっていませんでした**。\n\n撮像された画像そのものが 160mm より狭いことがあり、\nそのときは切り出し幅が画像サイズで頭打ちになります。\n**このとき、物理スケールを揃える処理は実質的に働いていません**\n(画像全体をそのまま使っているのと同じ状態になります)。\n\n数えたところ、**160mm では 28.9% のシリーズでこれが起きていました**。\n視野を 130mm にすると 0.5% まで下がります(こちらも全シリーズの走査による値。\n下のセルは 300 シリーズの標本なので数字は少しずれます)。\n\nこの1点だけを直した提出で、LB は 0.857 → 0.861 になりました。\n\n> **教訓**: 「実装した」と「効いている」は違います。\n> 効いているかは、**効かなかった件数を数えて**確かめます。\n"},{"cell_type":"code","execution_count":null,"metadata":{},"outputs":[],"source":"# ---- 視野を何 mm にすべきか(§5)-------------------------------------------\n# 切り出し幅 = 視野(mm) / spacing。これが画像サイズを超えると頭打ちになり、\n# 物理スケール合わせが働かなくなる。その割合を視野ごとに数える。\nok = hdr.dropna(subset=[\"spacing\", \"rows_px\"])\nprint(f\"{'視野':>6s}  {'働かないシリーズ':>18s}\")\nfor fov in (160, 150, 140, 130, 120):\n    need = fov / ok[\"spacing\"]\n    clipped = (need > ok[\"rows_px\"]).mean()\n    print(f\"{fov:4d}mm  {clipped:17.1%}\")\nprint(\"\\n→ 160mm は無視できない割合で働かない。130mm を採用した(§5)\")\n"},{"cell_type":"markdown","metadata":{},"source":"## 6. どの撮像面に、どの所見が写るのか\n\n膝のMRIは複数の向きから撮ります。矢状断(S)・冠状断(C)・軸位断(A)。\n\n「全部入れればよい」と決める前に、**平面を減らした版を作って提出**し、\nどの所見がどの平面で決まるのかを実測しました。\n\n| 使った平面 | LB |\n|---|---|\n| 冠状のみ (C) | 0.718 |\n| 矢状のみ (S) | 0.772 |\n| 矢状 + 軸位 (S+A) | 0.802 |\n| **3平面すべて** | **0.820** |\n\n所見ごとに見ると、解剖学的に納得できる並びになりました。\n以下は**手元の検証スコア(所見ごとの AUC)**で、上の表の LB とは別物です。\n\n| 所見 | 最も効く単一平面(検証AUC) | 解剖学的な期待 |\n|---|---|---|\n| MCL(内側側副靭帯) | **冠状 0.705**(矢状は 0.477 と最低) | 内側側副靭帯は冠状断で見る ✓ |\n| Medial Meniscus | **矢状 0.729** | 半月板は矢状断 ✓ |\n| PF OA(膝蓋大腿) | **軸位 0.629**(矢状 0.590 / 冠状 0.528) | 膝蓋骨の裏の軟骨は軸位 ✓ |\n| Synovitis | **軸位 0.804** | 関節液・滑膜は軸位 ✓ |\n\n**モデルが解剖学的に妥当なところを見ている傍証**になります。\nおかしなショートカット(たとえば施設ごとの画像のクセ)を\n学習しているわけではなさそうだ、と言えます。\n\n2枚目に足すなら矢状+軸位が良い、というのも検証スコアで確認できました\n(矢状+軸位 0.714 / 矢状+冠状 0.694)。\n\n## 7. 検証(CV)の切り方 — テストと同じ切れ方を再現する\n\nCV の目的は、**提出したときのスコアを手元で予測すること**です。\nなので、**テストデータと同じ切れ方**で切る必要があります。\n\nランダムに5分割すると何が起きるか。\n同じ施設・同じ装置で撮られた似た検査が、学習側と検証側の**両方に入ります**。\nモデルは「この装置の画像はこう写る」を覚えるだけで検証スコアを上げられてしまい、\n**手元では良いのに提出すると落ちる**ことになります。\n\nそこで、**装置を識別する複合キー**(以下「装置の指紋」と呼びます)を作り、\nその単位で分けます。\n\n```\n指紋 = Manufacturer | ManufacturerModelName | SoftwareVersions | MagneticFieldStrength\n```\n\n4,407 件はこの指紋で **112 群**に分かれます。\n群がまたがらないように5分割すれば、\n検証側は「学習で見たことのない装置」になります。\n\n### 指紋に使えなかった列\n\n最初は `ImagingFrequency`(撮像周波数)も指紋に入れようとしました。\n数えたら **8,023 種類**あり、**4,407 件のうち 3,978 件では\n1つの検査の中で値が複数**ありました。\n\nこれは装置のIDではなく、**撮像のたびに再調整される値**です。\n入れると群がバラバラになり、グループ分割の意味がなくなります。\n\n> **教訓**: グループのキーは「同じものには同じ値が入るか」を数えて確かめます。\n> 名前から選ぶと外します。\n"},{"cell_type":"code","execution_count":null,"metadata":{},"outputs":[],"source":"# ---- 検証の切り方(§7)------------------------------------------------------\n# 装置の指紋でグループを作り、群がまたがらないように5分割する。\n# 指紋は全シリーズのヘッダを1度だけ走査して作る(ここでは分割の作法だけ示す)。\nfrom sklearn.model_selection import GroupKFold\n\ndemo = train[[\"StudyInstanceUID\"]].copy()\ndemo[\"group\"] = np.random.default_rng(0).integers(0, 112, len(demo))  # 実データでは112群\n\ngkf = GroupKFold(n_splits=5)\nfor f, (tr_i, va_i) in enumerate(gkf.split(demo, groups=demo[\"group\"])):\n    a = set(demo.iloc[tr_i][\"group\"])\n    b = set(demo.iloc[va_i][\"group\"])\n    assert not (a & b), \"群が学習側と検証側にまたがっている\"\n    print(f\"  fold {f}: 学習 {len(tr_i):,} / 検証 {len(va_i):,}  群の重なり {len(a & b)}\")\nprint(\"\\n→ 検証側は『学習で見たことのない装置』になる(§7)\")\n"},{"cell_type":"markdown","metadata":{},"source":"## 8. ここまでで、ベースラインの構成が決まる\n\nEDA で分かったことが、そのまま設定になります。\n\n| 決めたこと | 根拠(前の節) |\n|---|---|\n| **学習ターゲット = レポートから作ったラベル** | ラベルが 58 件しかない(§1) |\n| **切り出し 130mm の物理スケール固定** | 160mm では 28.9% で働かない(§5) |\n| **3平面すべてを入力に使う** | 平面ごとに得意な所見が違う(§6) |\n| **スキャナ指紋の GroupKFold(5分割)** | 装置の学習を防ぐ(§7) |\n| **合成は fold 平均、混ぜるときは順位平均** | AUC は順位しか見ない(§0) |\n\n残るのは画像モデルの部分です。ここは一般的な構成で作ります。\n\n### 入力の作り方(2.5D)\n\nMRIは3次元ですが、3D CNN は重く、学習データも多くありません。\nそこで **2.5D** にします。\n\n1. 各平面から代表シリーズを1本選ぶ\n2. スライスを等間隔に **24枚** サンプルする\n3. **連続3枚を1組**にして、RGB の3チャンネルとして扱う(= 8組)\n4. 3平面 × 8組 = 24枚の「画像」を、1件の検査として扱う\n\n図にすると、1つの平面についてこうなります。\n\n```\n元のシリーズ(数十枚のスライス)\n  ▫▫▫▫▫▫▫▫▫▫▫▫▫▫▫▫▫▫▫▫▫▫▫▫▫▫▫▫▫▫▫▫▫▫▫▫\n\n        ↓ 等間隔に 24 枚えらぶ\n\n  ▪ ▪ ▪ ▪ ▪ ▪ ▪ ▪ ▪ ▪ ▪ ▪ ▪ ▪ ▪ ▪ ▪ ▪ ▪ ▪ ▪ ▪ ▪ ▪     24 枚\n  1 2 3 4 5 6 7 8 ...                        24\n\n        ↓ 連続する3枚を1組にまとめ、RGB の3チャンネルに入れる\n\n  [1|2|3] [4|5|6] [7|8|9] ...  [22|23|24]      8 組\n   R G B   R G B   R G B        R  G  B\n   └─ 1枚のカラー画像のように扱える ─┘\n```\n\n3平面ぶんを束ねたものが、1件の検査の入力です。\n\n```\n  矢状断 ─→ 8組 ┐\n  冠状断 ─→ 8組 ├─→ エンコーダに通す ─→ 24組ぶんの特徴ベクトル\n  軸位断 ─→ 8組 ┘                              │\n                                                ↓\n                       所見ごとに「どの組を重視するか」を学習(§8 プーリング)\n                                                ↓\n                                        12所見それぞれの確率\n```\n\n隣り合うスライスをチャンネルに入れることで、\n**平面に垂直な方向のつながり**を、2D のモデルのまま扱えます。\n半月板のように「ある断面にだけ写る」構造は、その断面を含む組が強く重み付けされます。\n\n### モデル\n\n- **エンコーダ**: DINOv2 ViT-S/14(自己教師あり学習済み)\n- **プーリング**: 8組のどれを重視するかを、**所見ごとに**学習する\n  (半月板を読むのに良い位置と、ACL を読むのに良い位置は違うため)\n- **学習率**: head は 1e-3、エンコーダは 8e-6 で**出力側6ブロックだけ**動かす\n\n最後の点が効きます。事前学習済みのエンコーダを普通の学習率で全層動かすと、\n**事前学習で獲得された重みが大きく書き換わり、その表現が失われます**。\n実際、全層を 1e-4 で動かした版は、はるかに小さい resnet18 に負けました。\n「適応させる、作り直さない」が正解でした。\n\n## 9. どの変更が、どれだけ効いたか\n\n上の構成は、いきなり全部決めたわけではありません。\n変更を1つずつ切り分けて提出し、それぞれの寄与を測りました。\n\n| 変更した点 | public LB |\n|---|---|\n| resnet18 2.5D + 自前ルールで作った弱ラベル | 0.820 |\n| **ターゲットを、レポートを読ませたLLMのラベルに差し替え** | **0.857** (+0.037) |\n| 切り出しを 160mm → 130mm(§5) | 0.861 (+0.004) |\n| エンコーダを DINOv2 ViT-S/14 + 判別的な学習率へ | 0.882 (+0.021) |\n| プーリングを所見ごとに(§8) | 0.884 (+0.002) |\n| 入力解像度 224px → **336px** | **0.898** (+0.014) |\n\n**一番大きく動いたのはラベルの差し替え(+0.037)** です。\n画像モデルの工夫を全部足したより大きい。§1 で書いた\n「これはラベルを作る問題である」が、そのまま結果に出ています。\n\n解像度を上げたのが2番目(+0.014)。130mm を 336px で見ると\n**1ピクセル 0.387mm** になります。半月板の裂け目を見るには、\nこのくらいの細かさが要る、ということだと考えています。\n\nこのように切り分けたので、**どの要素がどれだけ効いたのかを数字で言えます**。\nまとめて変えていたら、伸びた分がラベルによるものか解像度によるものか分からず、\n次にどこへ投資すべきかも決められませんでした。\n\n## 10. ここから先\n\nベースラインができたら、次に効きそうなのはこのあたりです。\n\n- **ラベルの質を上げる。** 一番動いた軸がここでした。\n  レポートの読み取りを改善すれば、まだ伸びる余地があります\n- **モデルを複数作って混ぜる。** ただし\n  「単体で弱いモデルを混ぜても効かない」ことは実測しました。\n  混ぜて効くのは**基準と同じくらい強くて、予測の傾向が違う**モデルだけです\n- **公開されているモデルと混ぜる。** 他の人が別の設計で学習したモデルは、\n  自分の中で工夫して作った多様性より、ずっと予測の傾向が違います\n"},{"cell_type":"markdown","metadata":{},"source":"## 付録: このベースラインを実際に走らせる\n\nここまでは分析でした。以下のセルは §8〜§9 で説明した構成を実際に動かし、\n`submission.csv` を書き出します。説明するだけでなく、**public LB 0.898 の提出物**\nがこのノートブックから出ます。\n\n重みは その構成の 5 fold を CC0 で公開したものです:\n`kitopl/rsna-knee-eda-baseline-weights`。\n学習に使ったラベルは公開の `stevenleehans/rsna-knee-llm-report-labels` なので、\n非公開のデータには一切依存していません。\n\n間違えやすい点が2つあります。\n\n- **推論側の前処理は学習側と完全に一致していなければなりません。** 130mm の切り出し、\n  24枚のスライス抽出、1〜99パーセンタイルでの正規化 — どれかがズレるとスコアが落ち、\n  しかも何も教えてくれません。以下は学習時の前処理の逐語コピーです。\n- **エンコーダは `pretrained=False` で組みます。** backbone の重みはチェックポイントに\n  丸ごと入っているので、インターネット接続は要りません。\n"},{"cell_type":"code","metadata":{},"source":"# ---- 推論の設定 -------------------------------------------------------------\nimport os\nimport time\nfrom concurrent.futures import ProcessPoolExecutor\n\nimport cv2\nimport pydicom\nimport timm\nimport torch\nimport torch.nn as nn\nfrom scipy.stats import rankdata\n\nOUT = \"/kaggle/working\"\nINPUT = str(ROOT)\n\nPLANES_ORDER = [\"Sagittal\", \"Coronal\", \"Axial\"]\nPLANES, N_SLICES, HW = 3, 24, 336\nTRIPLETS = N_SLICES // 3\nFOV_MM = 130.0          # 学習時と揃える(§5)\nBACKBONE = \"vit_small_patch14_dinov2.lvd142m\"\nBATCH = 8\nTIME_BUDGET = 7.5 * 3600   # 9時間制限に対する安全域。超えたら残りは0.5で埋める\n\n\ndef find_weights(fname=\"model_fold0.pt\"):\n    \"\"\"重みの場所を /kaggle/input から探す。\n\n    マウント先を決め打ちすると、データセットの付き方に依存して壊れる\n    (実際に `/kaggle/input/<slug>` だと思って外した。正しくは\n    `/kaggle/input/datasets/<owner>/<slug>`)。探して見つけ、\n    見つからなかったときに原因が分かるよう中身も出す。\n    \"\"\"\n    for root, _, files in os.walk(\"/kaggle/input\"):\n        if fname in files:\n            return root\n    return None\n\n\ndev = \"cuda\" if torch.cuda.is_available() else \"cpu\"\nprint(f\"device={dev}  torch={torch.__version__}  timm={timm.__version__}\")\nprint(\"/kaggle/input の直下:\", sorted(os.listdir(\"/kaggle/input\")))\nWEIGHTS = find_weights()\nprint(\"重みのディレクトリ:\", WEIGHTS)\nassert WEIGHTS, \"model_fold0.pt が /kaggle/input のどこにも無い\"\nprint(\"重みファイル:\", sorted(f for f in os.listdir(WEIGHTS) if f.endswith(\".pt\")))\n","execution_count":null,"outputs":[]},{"cell_type":"code","metadata":{},"source":"# ---- 前処理: 学習時からの逐語コピー -----------------------------------------\n# 学習時の前処理とズレると、スコアが静かに落ちる。書き直さずコピーする。\n\ndef pick_series(series_df):\n    \"\"\"study × 平面ごとに1本選ぶ。fluid-sensitive を優先。\"\"\"\n    df = series_df.copy()\n    df[\"_pref\"] = -df[\"Fluid_Sensitive\"].fillna(0).astype(int)\n    return (df.sort_values([\"StudyInstanceUID\", \"Anatomical_Plane\",\n                            \"_pref\", \"SeriesInstanceUID\"])\n              .groupby([\"StudyInstanceUID\", \"Anatomical_Plane\"], as_index=False)\n              .head(1))\n\n\ndef slice_position(ds):\n    \"\"\"IOP の法線への射影 = 真のスライス位置。\n\n    z成分をそのまま使ってはならない。Sagittal ではスライス軸は X、\n    Coronal では Y になる。\n    \"\"\"\n    iop = getattr(ds, \"ImageOrientationPatient\", None)\n    ipp = getattr(ds, \"ImagePositionPatient\", None)\n    if iop is None or ipp is None:\n        return None\n    n = np.cross([float(v) for v in iop[:3]], [float(v) for v in iop[3:]])\n    return float(np.dot([float(v) for v in ipp], n))\n\n\ndef load_volume(split, study, series):\n    d = os.path.join(INPUT, f\"{split}_series\", study, series)\n    try:\n        paths = [f.path for f in os.scandir(d) if f.name.endswith(\".dcm\")]\n    except OSError as e:\n        return None, f\"scandir: {e}\"\n    if not paths:\n        return None, \"no dcm\"\n\n    items = []\n    for p in paths:\n        try:\n            ds = pydicom.dcmread(p, force=True)\n            arr = ds.pixel_array\n        except Exception:\n            continue\n        ps = getattr(ds, \"PixelSpacing\", None)\n        items.append((slice_position(ds), arr,\n                      float(ps[0]) if ps is not None else None))\n    if not items:\n        return None, \"all slices failed to decode\"\n\n    if all(it[0] is None for it in items):\n        items = [(i, a, s) for i, (_, a, s) in enumerate(items)]\n    else:\n        items = [it for it in items if it[0] is not None]\n        items.sort(key=lambda x: x[0])\n\n    out = np.zeros((N_SLICES, HW, HW), dtype=np.uint8)\n    idx = np.linspace(0, len(items) - 1, N_SLICES).round().astype(int)\n    for k, i in enumerate(idx):\n        _, a, spacing = items[i]\n        a = a.astype(np.float32)\n        if a.ndim != 2:\n            a = a[..., 0] if a.ndim == 3 else a.reshape(a.shape[-2:])\n        if spacing and spacing > 0:\n            # 物理的な視野を切り出してからリサイズする(§5)\n            crop_px = int(round(FOV_MM / spacing))\n            crop_px = max(16, min(crop_px, min(a.shape)))\n            cy, cx = a.shape[0] // 2, a.shape[1] // 2\n            h = crop_px // 2\n            a = a[max(0, cy - h):cy + h, max(0, cx - h):cx + h]\n        flat = a[::4, ::4].ravel()\n        lo, hi = np.percentile(flat, [1, 99]) if flat.size else (0.0, 1.0)\n        a = np.clip((a - lo) / max(hi - lo, 1e-6), 0, 1)\n        out[k] = (cv2.resize(a, (HW, HW),\n                             interpolation=cv2.INTER_AREA) * 255).astype(np.uint8)\n    return out, None\n\n\ndef build_study(args):\n    study, chosen = args\n    vol = np.zeros((PLANES, N_SLICES, HW, HW), dtype=np.uint8)\n    for pi, plane in enumerate(PLANES_ORDER):\n        sid = chosen.get(plane)\n        if sid is None:\n            continue\n        v, _ = load_volume(\"test\", study, sid)\n        if v is not None:\n            vol[pi] = v\n    return study, vol\n","execution_count":null,"outputs":[]},{"cell_type":"code","metadata":{},"source":"# ---- モデル: 学習時と字句レベルで同一でなければならない ----------------------\n# 食い違うと state_dict のキーか形状がずれ、推論が静かに壊れる。\n\nPRETRAINED = False       # backbone の重みはチェックポイントに入っている\nUSE_TRIPLET_POS = False\n\n\ndef build_encoder(backbone):\n    kw = dict(pretrained=PRETRAINED, num_classes=0, in_chans=3)\n    if backbone.startswith(\"vit_\"):\n        # DINOv2 は patch14 / 事前学習は518px。img_size を渡すと\n        # timm が位置埋め込みを再サンプルして読み込む。\n        assert HW % 14 == 0, f\"patch14 の ViT には HW が14の倍数である必要がある: {HW}\"\n        kw[\"img_size\"] = HW\n    return timm.create_model(backbone, **kw)\n\n\nclass Net(nn.Module):\n    \"\"\"2.5D エンコーダ + 所見ごとの triplet attention。\n\n    attention は (平面, 組, 所見) ごとに重みを出す。半月板を読むのに良い組と\n    ACL を読むのに良い組は違うので、ここを所見ごとにするのが要点になる。\n    1ベクトルに潰してから12所見へ分岐すると、この区別は取り戻せない。\n    \"\"\"\n\n    def __init__(self, backbone, n_out=len(TARGETS)):\n        super().__init__()\n        self.enc = build_encoder(backbone)\n        d = self.enc.num_features\n        self.n_out = n_out\n        self.att = nn.Sequential(nn.Linear(d, 128), nn.Tanh(), nn.Linear(128, n_out))\n        self.pos = nn.Parameter(torch.zeros(1, 1, TRIPLETS, d))\n        self.drop = nn.Dropout(0.2)\n        self.cls_w = nn.Parameter(torch.empty(n_out, PLANES * d))\n        self.cls_b = nn.Parameter(torch.zeros(n_out))\n        nn.init.normal_(self.cls_w, std=0.02)\n\n    def forward(self, x):                      # (B, 3, 24, H, W)\n        B = x.shape[0]\n        x = x.reshape(B * PLANES * TRIPLETS, 3, HW, HW)\n        f = self.enc(x).reshape(B, PLANES, TRIPLETS, -1)\n        if USE_TRIPLET_POS:\n            f = f + self.pos\n        w = torch.softmax(self.att(f), dim=2)          # 組の方向に正規化\n        p = torch.einsum(\"bptd,bptc->bpcd\", f, w)      # 所見ごとの平面埋め込み\n        p = self.drop(p.permute(0, 2, 1, 3).reshape(B, self.n_out, -1))\n        return (p * self.cls_w.unsqueeze(0)).sum(-1) + self.cls_b\n\n\nnets = []\nfor f in range(5):\n    p = os.path.join(WEIGHTS, f\"model_fold{f}.pt\")\n    if not os.path.exists(p):\n        continue\n    m = Net(BACKBONE).to(dev).eval()\n    m.load_state_dict(torch.load(p, map_location=dev))\n    nets.append(m)\nassert nets, f\"重みが1つも読めなかった ({WEIGHTS})\"\nprint(f\"{len(nets)} fold 読み込み: {BACKBONE}\")\n","execution_count":null,"outputs":[]},{"cell_type":"code","metadata":{},"source":"# ---- 推論して submission.csv を書く -----------------------------------------\nt0 = time.time()\n\ntest = pd.read_csv(f\"{INPUT}/test.csv\")\nte_series = pd.read_csv(f\"{INPUT}/test_series.csv\")\nprint(f\"test: {len(test):,} 件の検査 / {len(te_series):,} シリーズ\")\n\nsel = pick_series(te_series)\ntasks = [(s, dict(zip(g.Anatomical_Plane, g.SeriesInstanceUID)))\n         for s, g in sel.groupby(\"StudyInstanceUID\")]\nprint(f\"{len(tasks):,} 件に対し {len(sel):,} シリーズを選択\")\n\norder, chunks = [], []\nworkers = max(1, (os.cpu_count() or 4))\nwith ProcessPoolExecutor(max_workers=workers) as ex:\n    buf_u, buf_x = [], []\n    for i, (study, vol) in enumerate(ex.map(build_study, tasks, chunksize=4), 1):\n        buf_u.append(study)\n        buf_x.append(vol)\n        last = (i == len(tasks))\n        if len(buf_u) < BATCH and not last:\n            continue\n\n        x = torch.from_numpy(np.stack(buf_x)).float().div_(255.0)\n        x = ((x - 0.449) / 0.226).to(dev)\n        with torch.no_grad(), torch.amp.autocast(\"cuda\", enabled=(dev == \"cuda\")):\n            # 5 fold を確率平均する(学習時の合成手続きと同じ)\n            p = sum(torch.sigmoid(m(x)) for m in nets) / len(nets)\n        chunks.append(p.float().cpu().numpy())\n        order += buf_u\n        buf_u, buf_x = [], []\n\n        if i % 200 == 0 or last:\n            el = time.time() - t0\n            print(f\"    {i:,}/{len(tasks):,}  {el/60:.1f} min  \"\n                  f\"残り {(len(tasks)-i)*el/max(i,1)/60:.0f} min\", flush=True)\n        if time.time() - t0 > TIME_BUDGET:\n            print(\"[budget] 時間切れ。残りは0.5で埋める\")\n            break\n\nP = np.concatenate(chunks) if chunks else np.zeros((0, len(TARGETS)))\nn = len(order)\nassert len(P) == n, f\"{len(P)} != {n}\"\n\n# AUC は順位しか見ないので、列ごとに順位へ直してから書き出す。\n# 単一モデルなら単調変換なのでスコアは変わらないが、複数モデルを混ぜるときと\n# 同じ手続きにしておく(ここがずれると「検証した混合」と「提出した混合」が別物になる)。\npred = np.column_stack([rankdata(P[:, j]) / max(n, 1) for j in range(len(TARGETS))])\npreds = dict(zip(order, pred))\n\nsub = pd.DataFrame({\"StudyInstanceUID\": test.StudyInstanceUID})\narr = np.full((len(sub), len(TARGETS)), 0.5, np.float64)\nfor i, u in enumerate(sub.StudyInstanceUID):\n    if u in preds:\n        arr[i] = preds[u]\nfor j, t in enumerate(TARGETS):\n    sub[t] = arr[:, j]\n\nsample = pd.read_csv(f\"{INPUT}/sample_submission.csv\")\nsub = sub[list(sample.columns)]\nassert list(sub.columns) == list(sample.columns)\nassert len(sub) == len(test)\nassert sub[TARGETS].notna().all().all()\nsub.to_csv(f\"{OUT}/submission.csv\", index=False)\n\nn_pred = sum(1 for u in sub.StudyInstanceUID if u in preds)\nprint(f\"\\nsubmission.csv を書き出し {sub.shape}  予測済み {n_pred:,}/{len(sub):,} \"\n      f\"(0.5で埋めたのは {len(sub)-n_pred:,})\")\nprint(f\"done in {(time.time()-t0)/60:.1f} min\")\n","execution_count":null,"outputs":[]}],"metadata":{"kernelspec":{"display_name":"Python 3","language":"python","name":"python3"},"language_info":{"name":"python","version":"3.11"}},"nbformat":4,"nbformat_minor":4}