{"cells":[{"cell_type":"markdown","metadata":{},"source":"# RSNA Knee DICOM Protocol Audit\n\n**目的**: 提取所有 DICOM metadata，生成设备和协议统计表\n\n**输出**:\n- `rsna_dicom_metadata_full.csv`: 每个 series 的完整 metadata\n- `rsna_protocol_summary.json`: 统计摘要\n- `protocol_clusters.csv`: Protocol cluster 分配\n\n**用途**: 判断 fastMRI/MRNet 等外部数据与 RSNA 的相似度"},{"cell_type":"code","execution_count":null,"metadata":{},"outputs":[],"source":"import pandas as pd\nimport pydicom\nfrom pathlib import Path\nfrom collections import defaultdict, Counter\nimport json\nimport numpy as np\nfrom tqdm import tqdm\nimport matplotlib.pyplot as plt\nimport seaborn as sns\n\nprint(\"Libraries imported successfully\")"},{"cell_type":"code","execution_count":null,"metadata":{},"outputs":[],"source":"def extract_dicom_metadata(dcm_path):\n    \"\"\"提取单个 DICOM 文件的关键 metadata\"\"\"\n    try:\n        dcm = pydicom.dcmread(dcm_path, stop_before_pixels=True)\n\n        metadata = {\n            # 设备信息\n            'Manufacturer': getattr(dcm, 'Manufacturer', 'Unknown'),\n            'ManufacturerModelName': getattr(dcm, 'ManufacturerModelName', 'Unknown'),\n            'MagneticFieldStrength': getattr(dcm, 'MagneticFieldStrength', None),\n\n            # 序列参数\n            'RepetitionTime': getattr(dcm, 'RepetitionTime', None),\n            'EchoTime': getattr(dcm, 'EchoTime', None),\n            'FlipAngle': getattr(dcm, 'FlipAngle', None),\n            'ImagingFrequency': getattr(dcm, 'ImagingFrequency', None),\n\n            # 空间参数\n            'PixelSpacing': getattr(dcm, 'PixelSpacing', None),\n            'SliceThickness': getattr(dcm, 'SliceThickness', None),\n            'SpacingBetweenSlices': getattr(dcm, 'SpacingBetweenSlices', None),\n            'Rows': getattr(dcm, 'Rows', None),\n            'Columns': getattr(dcm, 'Columns', None),\n\n            # 序列描述\n            'SeriesDescription': getattr(dcm, 'SeriesDescription', 'Unknown'),\n            'SequenceName': getattr(dcm, 'SequenceName', 'Unknown'),\n            'ScanningSequence': getattr(dcm, 'ScanningSequence', 'Unknown'),\n            'SequenceVariant': getattr(dcm, 'SequenceVariant', 'Unknown'),\n            'ScanOptions': getattr(dcm, 'ScanOptions', 'Unknown'),\n\n            # Study/Series ID\n            'StudyInstanceUID': getattr(dcm, 'StudyInstanceUID', None),\n            'SeriesInstanceUID': getattr(dcm, 'SeriesInstanceUID', None),\n        }\n\n        # 处理 PixelSpacing（可能是 list）\n        if metadata['PixelSpacing'] is not None:\n            try:\n                ps = metadata['PixelSpacing']\n                metadata['PixelSpacing_0'] = float(ps[0])\n                metadata['PixelSpacing_1'] = float(ps[1])\n            except:\n                metadata['PixelSpacing_0'] = None\n                metadata['PixelSpacing_1'] = None\n\n        return metadata\n    except Exception as e:\n        print(f\"Error reading {dcm_path}: {e}\")\n        return None\n\nprint(\"extract_dicom_metadata() defined\")"},{"cell_type":"code","execution_count":null,"metadata":{},"outputs":[],"source":"# 读取 series metadata\nDATA_ROOT = Path('/kaggle/input/competitions/rsna-knee-abnormality-detection')\ntrain_series_csv = DATA_ROOT / 'train_series.csv'\ndicom_root = DATA_ROOT / 'train_series'\n\nseries_df = pd.read_csv(train_series_csv)\nprint(f\"Total series: {len(series_df)}\")\nprint(f\"\\nFirst few rows:\")\nprint(series_df.head())"},{"cell_type":"code","execution_count":null,"metadata":{},"outputs":[],"source":"# 提取每个 series 的一个代表性 DICOM\n# 为了速度，先采样 500 个 series 测试\nSAMPLE_SIZE = None  # 全量扫描所有 series\n\nif SAMPLE_SIZE:\n    series_sample = series_df.sample(n=min(SAMPLE_SIZE, len(series_df)), random_state=42)\n    print(f\"Sampling {len(series_sample)} series for quick audit\")\nelse:\n    series_sample = series_df\n    print(f\"Scanning all {len(series_sample)} series\")\n\nmetadata_list = []\n\nfor idx, row in tqdm(series_sample.iterrows(), total=len(series_sample), desc=\"Scanning DICOMs\"):\n    study_uid = row['StudyInstanceUID']\n    series_uid = row['SeriesInstanceUID']\n\n    # 构造 DICOM 路径\n    series_dir = dicom_root / study_uid / series_uid\n\n    if not series_dir.exists():\n        continue\n\n    # 读取该 series 的第一个 DICOM\n    dcm_files = list(series_dir.glob('*.dcm'))\n    if not dcm_files:\n        continue\n\n    metadata = extract_dicom_metadata(dcm_files[0])\n    if metadata:\n        # 添加 CSV 中的信息\n        metadata['Fluid_Sensitive'] = row['Fluid_Sensitive']\n        metadata['Fat_Suppression'] = row['Fat_Suppression']\n        metadata['Anatomical_Plane'] = row['Anatomical_Plane']\n        metadata_list.append(metadata)\n\nprint(f\"\\nSuccessfully extracted metadata from {len(metadata_list)} series\")"},{"cell_type":"code","execution_count":null,"metadata":{},"outputs":[],"source":"# 转换为 DataFrame\nmetadata_df = pd.DataFrame(metadata_list)\nprint(f\"Metadata DataFrame shape: {metadata_df.shape}\")\nprint(f\"\\nColumns: {list(metadata_df.columns)}\")\nprint(f\"\\nFirst few rows:\")\nmetadata_df.head()"},{"cell_type":"code","execution_count":null,"metadata":{},"outputs":[],"source":"# 生成统计报告\nreport = {}\n\n# 1. 设备厂商分布\nprint(\"=\"*60)\nprint(\"RSNA DICOM Protocol Summary\")\nprint(\"=\"*60)\n\nmanufacturer_counts = metadata_df['Manufacturer'].value_counts()\nprint(\"\\n📌 Manufacturer Distribution:\")\nfor mfr, count in manufacturer_counts.items():\n    pct = count / len(metadata_df) * 100\n    print(f\"  {mfr:30s}: {count:5d} ({pct:5.1f}%)\")\n    \nreport['Manufacturer'] = manufacturer_counts.to_dict()\n\n# 2. 磁场强度分布\nfield_strength_counts = metadata_df['MagneticFieldStrength'].value_counts()\nprint(\"\\n📌 Field Strength Distribution:\")\nfor fs, count in field_strength_counts.items():\n    if fs is not None:\n        pct = count / len(metadata_df) * 100\n        print(f\"  {fs}T{' ':28s}: {count:5d} ({pct:5.1f}%)\")\n        \nreport['FieldStrength'] = field_strength_counts.to_dict()\n\n# 3. 序列类型分布\nplane_counts = metadata_df['Anatomical_Plane'].value_counts()\nprint(\"\\n📌 Anatomical Plane Distribution:\")\nfor plane, count in plane_counts.items():\n    pct = count / len(metadata_df) * 100\n    print(f\"  {plane:30s}: {count:5d} ({pct:5.1f}%)\")\n\n# 4. Fluid Sensitive + Fat Suppression 组合\nprint(\"\\n📌 Sequence Type (Fluid + Fat):\")\nfs_fat_counts = metadata_df.groupby(['Fluid_Sensitive', 'Fat_Suppression']).size()\nfor (fluid, fat), count in fs_fat_counts.items():\n    label = f\"Fluid={fluid}, Fat={fat}\"\n    pct = count / len(metadata_df) * 100\n    print(f\"  {label:30s}: {count:5d} ({pct:5.1f}%)\")\n\n# 5. PixelSpacing 分布\nif 'PixelSpacing_0' in metadata_df.columns:\n    ps_values = metadata_df['PixelSpacing_0'].dropna()\n    print(\"\\n📌 PixelSpacing (mm):\")\n    print(f\"  Mean ± Std: {ps_values.mean():.3f} ± {ps_values.std():.3f}\")\n    print(f\"  Range: [{ps_values.min():.3f}, {ps_values.max():.3f}]\")\n    print(f\"  Median: {ps_values.median():.3f}\")\n    print(f\"  Q25, Q75: {ps_values.quantile(0.25):.3f}, {ps_values.quantile(0.75):.3f}\")\n    \n    report['PixelSpacing'] = {\n        'mean': float(ps_values.mean()),\n        'std': float(ps_values.std()),\n        'min': float(ps_values.min()),\n        'max': float(ps_values.max()),\n        'median': float(ps_values.median()),\n    }\n\n# 6. SliceThickness 分布\nif 'SliceThickness' in metadata_df.columns:\n    st_values = metadata_df['SliceThickness'].dropna()\n    print(\"\\n📌 SliceThickness (mm):\")\n    print(f\"  Mean ± Std: {st_values.mean():.2f} ± {st_values.std():.2f}\")\n    print(f\"  Range: [{st_values.min():.2f}, {st_values.max():.2f}]\")\n    print(f\"  Median: {st_values.median():.2f}\")\n    print(f\"  Q25, Q75: {st_values.quantile(0.25):.2f}, {st_values.quantile(0.75):.2f}\")\n    \n    report['SliceThickness'] = {\n        'mean': float(st_values.mean()),\n        'std': float(st_values.std()),\n        'min': float(st_values.min()),\n        'max': float(st_values.max()),\n        'median': float(st_values.median()),\n    }\n\n# 7. TR/TE 分布\nfor param in ['RepetitionTime', 'EchoTime']:\n    if param in metadata_df.columns:\n        values = metadata_df[param].dropna()\n        if len(values) > 0:\n            print(f\"\\n📌 {param} (ms):\")\n            print(f\"  Mean ± Std: {values.mean():.1f} ± {values.std():.1f}\")\n            print(f\"  Range: [{values.min():.1f}, {values.max():.1f}]\")\n            print(f\"  Median: {values.median():.1f}\")\n\nprint(\"\\n\" + \"=\"*60)"},{"cell_type":"code","execution_count":null,"metadata":{},"outputs":[],"source":"# 创建 Protocol Clusters\nmetadata_df['Protocol_Cluster'] = (\n    metadata_df['Manufacturer'].astype(str) + '_' +\n    metadata_df['MagneticFieldStrength'].astype(str) + 'T'\n)\n\ncluster_counts = metadata_df['Protocol_Cluster'].value_counts()\n\nprint(\"\\n📌 Protocol Clusters (Manufacturer_FieldStrength):\")\nfor cluster, count in cluster_counts.items():\n    pct = count / len(metadata_df) * 100\n    print(f\"  {cluster:40s}: {count:5d} ({pct:5.1f}%)\")"},{"cell_type":"code","execution_count":null,"metadata":{},"outputs":[],"source":"# 可视化\nfig, axes = plt.subplots(2, 3, figsize=(18, 10))\n\n# Manufacturer\nmanufacturer_counts.plot(kind='bar', ax=axes[0, 0], color='steelblue')\naxes[0, 0].set_title('Manufacturer Distribution', fontsize=14, fontweight='bold')\naxes[0, 0].set_ylabel('Count')\naxes[0, 0].tick_params(axis='x', rotation=45)\n\n# Field Strength\nfield_strength_counts.plot(kind='bar', ax=axes[0, 1], color='coral')\naxes[0, 1].set_title('Field Strength Distribution', fontsize=14, fontweight='bold')\naxes[0, 1].set_ylabel('Count')\n\n# Anatomical Plane\nplane_counts.plot(kind='bar', ax=axes[0, 2], color='mediumseagreen')\naxes[0, 2].set_title('Anatomical Plane Distribution', fontsize=14, fontweight='bold')\naxes[0, 2].set_ylabel('Count')\n\n# PixelSpacing histogram\nif 'PixelSpacing_0' in metadata_df.columns:\n    metadata_df['PixelSpacing_0'].dropna().hist(bins=50, ax=axes[1, 0], color='plum', edgecolor='black')\n    axes[1, 0].set_title('PixelSpacing Distribution', fontsize=14, fontweight='bold')\n    axes[1, 0].set_xlabel('PixelSpacing (mm)')\n    axes[1, 0].set_ylabel('Count')\n\n# SliceThickness histogram\nif 'SliceThickness' in metadata_df.columns:\n    metadata_df['SliceThickness'].dropna().hist(bins=50, ax=axes[1, 1], color='gold', edgecolor='black')\n    axes[1, 1].set_title('SliceThickness Distribution', fontsize=14, fontweight='bold')\n    axes[1, 1].set_xlabel('SliceThickness (mm)')\n    axes[1, 1].set_ylabel('Count')\n\n# Protocol Cluster\ncluster_counts.head(10).plot(kind='barh', ax=axes[1, 2], color='lightcoral')\naxes[1, 2].set_title('Top 10 Protocol Clusters', fontsize=14, fontweight='bold')\naxes[1, 2].set_xlabel('Count')\n\nplt.tight_layout()\nplt.savefig('rsna_dicom_protocol_summary.png', dpi=150, bbox_inches='tight')\nplt.show()\n\nprint(\"\\n✅ Visualization saved: rsna_dicom_protocol_summary.png\")"},{"cell_type":"code","execution_count":null,"metadata":{},"outputs":[],"source":"# 保存输出文件\noutput_dir = Path('./dicom_audit_output')\noutput_dir.mkdir(exist_ok=True)\n\n# 完整 metadata\nmetadata_df.to_csv(output_dir / 'rsna_dicom_metadata_full.csv', index=False)\nprint(f\"✅ Saved: {output_dir / 'rsna_dicom_metadata_full.csv'}\")\n\n# JSON 报告\nwith open(output_dir / 'rsna_protocol_summary.json', 'w') as f:\n    json.dump(report, f, indent=2, default=str)\nprint(f\"✅ Saved: {output_dir / 'rsna_protocol_summary.json'}\")\n\n# Protocol clusters\ncluster_df = metadata_df[['StudyInstanceUID', 'SeriesInstanceUID', 'Protocol_Cluster']]\ncluster_df.to_csv(output_dir / 'protocol_clusters.csv', index=False)\nprint(f\"✅ Saved: {output_dir / 'protocol_clusters.csv'}\")\n\nprint(\"\\n\" + \"=\"*60)\nprint(\"✅ DICOM Protocol Audit Complete!\")\nprint(\"=\"*60)\nprint(\"\\nNext Steps:\")\nprint(\"1. Download output files from Kaggle\")\nprint(\"2. Compare with fastMRI/MRNet protocol specifications\")\nprint(\"3. Decide which external data is most similar to RSNA\")\nprint(\"4. Design Leave-One-Protocol-Out validation strategy\")"}],"metadata":{"kernelspec":{"display_name":"Python 3","language":"python","name":"python3"},"language_info":{"name":"python","version":"3.10.0"}},"nbformat":4,"nbformat_minor":4}