2
0

Delete article

Deleted articles cannot be recovered.

Draft of this article would be also deleted.

Are you sure you want to delete this article?

【コメントアウトで解説】PythonでCT被ばく線量を計算するコード

Last updated at Posted at 2025-10-17

1. はじめに

皆さんは、CT検査を受けたことがありますか? 大きなドーナツ状の機械に入っていく、アレです。CTは体の中を詳しく見ることができて非常に有用ですが、「放射線を使う」と聞いて、少し不安に思ったことがあるかもしれません。

「いったい、どれくらい被ばくするんだろう?」

そんな疑問に、Pythonコードを使って答えていくのがこの記事のテーマです。今回は、国際的な放射線防護の基準 (ICRP Publication 103) に基づき、「胸部CTを1回撮影した場合の被ばく線量」をシミュレーションするコードを、放射線の知識も交えながら、分かりやすく解説していきます。

  • 対象:Python初心者〜(診療放射線技師/医学物理士)
  • 環境:Windows + Anaconda + Jupyter(Mac/Linux可)
  • この記事で得られるもの
  1. CTの線量指標(CTDI、DLP、SSDE)の基本的な関係
  2. 吸収線量(Gy)から実効線量(Sv)への変換プロセス
  3. Pythonを使って、臓器ごとの線量を考慮したシミュレーションを行う方法
  4. 計算結果を可視化し、他の検査と比較して評価する方法

2. ライブラリのインポート

ライブラリを読み込まないと何も始まりません。
なので計算やデータ整理、グラフ化に必要な道具(ライブラリ)を読み込みます。

import math
import pandas as pd  
import matplotlib.pyplot as plt 
from IPython.display import display 
from io import BytesIO
from IPython.display import display, HTML

# グラフで日本語を使いたい場合に、どこかのセルで一度だけ実行してください。
!pip install japanize-matplotlib
import japanize_matplotlib
  • import math
    数学の複雑な計算(対数、三角関数、平方根など)を行うためのmathモジュールを読みこむ。

  • import pandas as pd
    データ分析でよく使われるpandasライブラリを「pd」という短い名前で使えるようにする。

  • import matplotlib.pyplot as plt
    グラフ作成ライブラリmatplotlibの中の描画機能(pyplot)を「plt」という名前で使えるようする。

  • from IPython.display import display
    Jupyter Notebookなどで、データや画像などを表示するためのdisplay関数を読みこむ。

  • pip install japanize-matplotlib
    matplotlibで日本語を表示できるようにする、japanize-matplotlibというライブラリをインストールするためのコマンド。

  • from io import BytesIO
    ioモジュールからBytesIOクラスをインポートする構文。これは、メモリ上でバイナリデータ(画像や音声など)をファイルのように扱えるようにするための機能。

  • from IPython.display import display, HTML
    HTMLコンテンツを直接表示するためにdisplay関数とHTML関数をインポートする。

3. 放射線被ばくの「モノサシ」を知ろう:Gy と Sv

3-1. 放射線の単位について

まず、放射線の量を表す重要な単位を押さえておきましょう。

  • 吸収線量 D [Gy]:物質が放射線から受け取ったエネルギーの量。物理的な量。
  • 等価線量 H [Sv]:放射線の種類による影響の違いを考慮した線量。
H_T=D_T×w_R(H_T:臓器Tの等価線量、D_T:臓器Tの吸収線量)
  • 実効線量 E [Sv]:等価線量に臓器別の w_Tを掛けて合計した全身のリスク指標。
E=\sum T(w_T×H_T)

3-2. ICRP103に基づく「定数」の定義

次に、線量計算の土台となる「係数」を定義します。これらは国際放射線防護委員会(ICRP)が定めた世界共通のルールブック(ICRP Publication 103)に基づいています。

① 放射線荷重係数 (w_R)
放射線にはX線、α線、中性子線など様々な種類があり、それぞれ"攻撃力"が違います。この攻撃力の違いを数値化したのが放射線荷重係数(w_R)です。

## 2-1. 放射線荷重係数 (w_R)# 放射線荷重係数(ICRP 103)
# CT(X線)は基準なので 1
W_R = {
    "xray": 1.0,      # X線/γ線
    "photon": 1.0,    # 同義
}

X線の場合、w_R=1なので、吸収線量と等価線量の数値は同じになりますが、単位と意味が違う(Gy → Sv)という点がポイントです。

def gy_to_sv(absorbed_dose_gy: float, radiation: str = "xray") -> float:
    """吸収線量(Gy)→等価線量(Sv)。X線では数値は同じ(単位が変わる)。"""
    wR = W_R.get(radiation, 1.0)
    return absorbed_dose_gy * wR

② 組織荷重係数(w_T) :臓器の "弱点度"
私たちの体にある臓器も、放射線に対する"弱さ(感受性)"が異なります。細胞分裂が活発な臓器ほど、影響を受けやすいとされています。この臓器ごとの弱点度を数値化したのが組織荷重係数(w_T)です。全身の合計が1になるように定められています。

## 2-2. 組織荷重係数 (w_T)
# 体の臓器・組織は、それぞれ放射線に対する「弱さ(感受性)」が異なります。
# 細胞分裂が活発な臓器ほど影響を受けやすいとされます。
# この「臓器ごとの放射線感受性」を数値化したのが組織荷重係数です。
# 値が大きいほど、リスクへの寄与が大きいことを意味します。
W_T = {
    # 感受性が最も高いグループ (0.12)
    "red_bone_marrow": 0.12,  # 赤色骨髄(血液を作る場所)
    "colon":           0.12,  # 結腸
    "lung":            0.12,  # 肺
    "stomach":         0.12,  # 胃
    "breast":          0.12,  # 乳房
    "remainder":       0.12,  # その他13臓器の平均に対して適用
    
    # 感受性が次に高いグループ (0.08)
    "gonads":          0.08,  # 生殖腺(精巣・卵巣)
    
    # 感受性が比較的低いグループ (0.04)
    "bladder":         0.04,  # 膀胱
    "esophagus":       0.04,  # 食道
    "liver":           0.04,  # 肝臓
    "thyroid":         0.04,  # 甲状腺
    
    # 感受性が最も低いグループ (0.01)
    "bone_surface":    0.01,  # 骨表面
    "brain":           0.01,  # 脳
    "salivary_glands": 0.01,  # 唾液腺
    "skin":            0.01,  # 皮膚
}

# 上記の"remainder"(残りの臓器)に分類される13臓器のリストです。
REMAINDER_TISSUES_13 = [
    "adrenals", "extrathoracic_region", "gall_bladder", "heart", "kidneys",
    "lymphatic_nodes", "muscle", "oral_mucosa", "pancreas",
    "prostate_uterus", "small_intestine", "spleen", "thymus"
]

※以下、環境省のホームページから

image.png

出典:https://www.env.go.jp/chemi/rhm/current/02-03-05.html

3-3. 実効線量

私たちのゴールは、実効線量(Sv)をPythonで計算することです。

実効線量[Sv] = 吸収線量D [Gy] × 放射線荷重係数 (w_R) × 組織荷重係数 (w_T)

で算出することができます。

  • 対象はCTに限定(X線, (w_R=1))

  • X線(CT)の w_R=1 を使った Gy↔Sv変換が分かる/動く

  • ICRP103 の組織荷重係数 (w_T) を使い、等価線量 → 実効線量に変換

  • 医療被ばくは線量限度の適用外だが、公衆 1 mSv/年などと“目安比較”できる表も出力

Note: 医療被ばくには法令上の線量限度は適用されません。記事中の比較表は “感覚を掴むための目安” として利用してください。

4. CTDI, DLP, SSDEの計算

いよいよCTの線量指標の計算です。CT装置のコンソール画面や撮影レポートでよく見かけるこれらの指標の意味と関係性をコードで見ていきましょう。

  • CTDI (Computed Tomography Dose Index): CT装置が1スライスあたりにどれくらいの線量を出すかの指標。CTDI_vol [mGy] は、撮影ピッチ(ヘリカルスキャンの螺旋の間隔)を考慮した値。
  • DLP (Dose Length Product): 照射した範囲の「総線量」を表す指標。計算式は、
DLP=CTDI_{vol}×Scan Length[mGy・cm] 
  • SSDE (Size-Specific Dose Estimate): 患者さんの体格を考慮した吸収線量の推定値 [mGy]。同じ撮影条件でも、小柄な患者さんの方が中心部の線量は高くなる傾向があります。この体格差を補正するのがSSDEです。

4-1. 体格補正係数と実効直径の推定

SSDEを計算するには、AAPM (米国医学物理学会) が公開しているレポート(No. 204, 220)の補正係数を使います。このコードでは、胸部で標準的に使われる32cmファントム用の係数を辞書として定義しています。
出典:https://www.aapm.org/pubs/reports/rpt_204.pdf

# 4-1. AAPM 220(=204): Dwベース f_size(32 cmファントム, 8〜45 cm 抜粋)
F_DW_32CM = {
     8: 2.76,  9: 2.66, 10: 2.57, 11: 2.47, 12: 2.38, 13: 2.30, 14: 2.22, 15: 2.14, 16: 2.06,
    17: 1.98, 18: 1.91, 19: 1.84, 20: 1.78, 21: 1.71, 22: 1.65, 23: 1.59, 24: 1.53, 25: 1.48,
    26: 1.43, 27: 1.37, 28: 1.32, 29: 1.28, 30: 1.23, 31: 1.19, 32: 1.14, 33: 1.10, 34: 1.06,
    35: 1.02, 36: 0.99, 37: 0.95, 38: 0.92, 39: 0.88, 40: 0.85, 41: 0.82, 42: 0.79, 43: 0.76,
    44: 0.74, 45: 0.71,)
}

表にない中間の体格(例: 20.5 cm)にも対応できるよう、線形補間を行う関数 _interp_f_dw_32cm も用意しています。

患者さんの体格は、撮影画像から前後径(AP)と左右径(LAT)を計測して実効直径

(D_ {eff}=\sqrt{AP×LAT})

を計算するのが理想ですが、それが難しい場合は身長・体重からBMIを算出し、近似式で推定することもできます。

D_ {eff}=0.6414(cm・m^2/kg)×BMI+12.976(cm)

出典:Matheoud, R., Al-Maymani, N., Oldani, A. et al. The role of activity, scan duration and patient’s body mass index in the optimization of FDG imaging protocols on a TOF-PET/CT scanner. EJNMMI Phys 8, 35 (2021).
https://ejnmmiphys.springeropen.com/articles/10.1186/s40658-021-00380-9#ref-CR24

def estimate_deff_cm(height_cm: float, weight_kg: float, ap_cm: float = None, lat_cm: float = None):
    """
    実効直径 Deff[cm] を返す。AP/LATがあれば sqrt(AP×LAT) を採用。
    無ければ BMI から Deff ≈ 0.6414*BMI + 12.976 を近似(胸部での代用)。
    """
    h_m = height_cm / 100.0 if height_cm else None
    bmi = (weight_kg / (h_m**2)) if (h_m and weight_kg) else None
    if ap_cm and lat_cm:
        deff = math.sqrt(ap_cm * lat_cm)
        return float(deff), float(bmi) if bmi else None
    if bmi is None:
        raise ValueError("AP/LAT または 身長・体重のいずれかを指定してください。")
    deff = 0.6414 * bmi + 12.976
    return float(deff), float(bmi)

4-2. CT線量指標の計算

これらを使って、入力されたCTDIwCTDIvol、DLP、そして体格を考慮したSSDEまでを一気に計算する関数がこちらです。

def compute_ctdi_dlp_ssde(
    *,
    pitch: float = 1.2,
    scan_length_cm: float = 35.0,
    ctdiw_mGy: float = None,
    ctdi_vol_mGy: float = None,
    height_cm: float = 170.0,
    weight_kg: float = 65.0,
    ap_cm: float = None,
    lat_cm: float = None,
    phantom: str = "32cm",
) -> dict:
    """CTDIw/CTDIvol→DLP→SSDE(サイズ補正)→まとめ"""
    if ctdi_vol_mGy is None and ctdiw_mGy is None:
        raise ValueError("CTDIw または CTDIvol のどちらか一方を指定してください。")
    if ctdi_vol_mGy is None:
        ctdi_vol_mGy = ctdiw_mGy / pitch
    if ctdiw_mGy is None:
        ctdiw_mGy = ctdi_vol_mGy * pitch

    dlp_mGy_cm = ctdi_vol_mGy * scan_length_cm

    deff_cm, bmi = estimate_deff_cm(height_cm, weight_kg, ap_cm, lat_cm)
    if phantom != "32cm":
        raise NotImplementedError("本スクリプトは胸部=32 cmファントムのみ実装")
    f_size = _interp_f_dw_32cm(deff_cm)
    ssde_mGy = ctdi_vol_mGy * f_size

    return {
        "Pitch": pitch,
        "ScanLength_cm": scan_length_cm,
        "CTDIw_mGy": ctdiw_mGy,
        "CTDIvol_mGy": ctdi_vol_mGy,
        "DLP_mGy_cm": dlp_mGy_cm,
        "Height_cm": height_cm,
        "Weight_kg": weight_kg,
        "BMI": bmi,
        "Deff_cm": deff_cm,
        "f_size(32cm)": f_size,
        "SSDE_mGy": ssde_mGy,
    }

5. 胸部CTの被ばく線量シミュレーション・モデル

5-1. 撮影範囲のモデル化

ここがこのシミュレーションのコアです!
実際の胸部CT検査の状況を、簡単なモデルとして再現します。胸部CTのX線ビームに、各臓器がどのくらいの割合で含まれるかを推定した係数です。
胸部CTといっても、すべての臓器が100%X線に当たるわけではありません。肺や心臓はほぼ完全に撮影範囲に入りますが、肝臓のてっぺんが少しだけ、甲状腺の下の方が少しだけ、というように臓器によって撮影範囲に含まれる割合が異なります。

今回はこの割合を、以下のように仮定しました。

1.0: 100%撮影範囲に入る(例: 肺)
0.1: 10%だけ撮影範囲の端にかかる(例: 肝臓)
0.0: 全く撮影範囲に入らない(例: 脳、大腸)

※注意:これはあくまで簡易的なモデルです。実際の撮影範囲によって値は変動します。

CHEST_COVERAGE = {
    # 主要ターゲット臓器(胸部にほぼ全量含まれると想定)
    "lung": 1.0,
    "breast": 1.0,
    "esophagus": 1.0,
    # 分布臓器(胸郭内の一部が照射されると想定)
    "red_bone_marrow": 0.25,  # 胸骨や肋骨内の骨髄を概算
    "bone_surface":    0.25,
    "thyroid":         0.30,  # 撮影範囲の上端に少しだけ入ることを想定
    "skin":            0.50,  # 照射野内の皮膚面積を概算
    "liver":           0.10,  # 撮影範囲の下端に少しだけかかることを想定
    # 胸部から遠い臓器(完全に範囲外と想定)
    "stomach": 0.0, "colon": 0.0, "gonads": 0.0,
    "brain": 0.0, "salivary_glands": 0.0, "bladder": 0.0,
}

5-2. 臓器線量の計算と実効線量の合成

以上の係数を使って、各臓器の吸収線量D_Tを計算し、そこから等価線量H_T、そして組織荷重係数w_Tを掛け合わせて実効線量への寄与を求め、最後に全てを合計します。


E=\sum w_T(w_R×D_T×coverage_T)

def compute_chest_ct_organs(dose_gy: float) -> pd.DataFrame:
    """
    胸部CTの平均吸収線量 dose_gy [Gy] を入力とし、
    各臓器 D_T → H_T → w_T 合成(ICRP103)を計算。
    """
    rows = []

    # 主要臓器
    for organ, frac in CHEST_COVERAGE.items():
        D_T = dose_gy * frac
        H_T = gy_to_sv(D_T, radiation="xray")
        wT = W_T.get(organ, 0.0)
        E_T = wT * H_T
        rows.append({
            "臓器": organ,
            "吸収線量[mGy]": D_T * 1000,
            "等価線量[mSv]": H_T * 1000,
            "w_T(ICRP103)": wT,
            "臓器寄与E[mSv]": E_T * 1000,
        })

    # remainder(13臓器の平均)
    D_list_rem = []
    H_list_rem = []
    for t in REMAINDER_TISSUES_13:
        frac = REMAINDER_COVERAGE.get(t, 0.0)
        D_T_rem = dose_gy * frac       # [Gy]
        H_T_rem = gy_to_sv(D_T_rem)    # [Sv]
        D_list_rem.append(D_T_rem)
        H_list_rem.append(H_T_rem)

    D_mean_rem = sum(D_list_rem) / len(REMAINDER_TISSUES_13)  # [Gy]
    H_mean_rem = sum(H_list_rem) / len(REMAINDER_TISSUES_13)  # [Sv]
    E_rem = W_T["remainder"] * H_mean_rem

    rows.append({
        "臓器": "remainder(13臓器平均)",
        "吸収線量[mGy]": D_mean_rem * 1000,
        "等価線量[mSv]": H_mean_rem * 1000,
        "w_T(ICRP103)": W_T["remainder"],
        "臓器寄与E[mSv]": E_rem * 1000,
    })

6. シミュレーションの実行と結果表示

いよいよ最終段階です。これまでに作成した関数を使って、具体的な数値を入力し、計算を実行してみましょう。

6-1. 入力値の設定

ここに、実際のCT装置から得られた情報や患者さんの体格情報を入力します。

# ===================================================================
# 6. 入力(ここを自由に変更)
# ===================================================================
# --- 撮影条件 / 装置指標 ---
PITCH = 1.2
SCAN_LEN_CM = 35.0
CTDIW_INPUT_MGY   = 9.0      # 例:CTDIw(mGy)
CTDIVOL_INPUT_MGY = None     # CTDIwがあればNoneでOK

# --- 体格 ---
HEIGHT_CM = 170.0            # 身長(cm)
WEIGHT_KG = 65.0             # 体重(kg)

# --- DLP→E 換算係数(成人胸部の代表値) ---
K_CHEST_MSV_PER_MGY_CM = 0.014

6-2. 計算実行と結果表示

計算パートでは、これまで定義した関数を順番に呼び出していきます。

  1. compute_ctdi_dlp_ssde: CTDIwなどからDLPとSSDEを計算。
  2. compute_chest_ct_organs:SSDEを基に臓器ごとの線量を計算し、実効線量を合成。
  3. estimate_effective_dose_from_dlp_chest_adult : 比較のため、DLPに換算係数kを掛けるだけの簡易的な実効線量も計算。

そして、これらの結果をpandasのデータフレームやmatplotlibのグラフで見やすく表示します。

【実行結果の例】
まず、CTDIからSSDEまでの計算結果です。患者さんの体格(Deff_cm)で補正され、SSDEが計算されていることがわかります。

image.png

次に、臓器ごとの線量と、実効線量への寄与を計算した結果です。

image.png

そして、最終的な実効線量のまとめです。

image.png

臓器モデルから合成した方法と、DLPから簡易換算した方法で、少し結果が異なることがわかりますね。これは、簡易換算係数kが「標準的な体格・撮影範囲」を仮定しているのに対し、臓器モデルは「体格と簡易的な撮影範囲」の条件を反映しているためです。

臓器ごとの実効線量への寄与をグラフにすると、どの臓器が被ばくの主要因となっているかが一目瞭然です。

image.png

最後に、この線量が他の検査や日常の被ばくと比べてどの程度なのかを見てみましょう。

image.png

image.png

今回のシミュレーション結果(3.73 mSv)は、胸部X線写真の約187回分、自然放射線の約568日分に相当することが分かりました。

表示ユーティリティ

計算した線量がどの程度の大きさなのかを直感的に理解するために、様々な比較対象と並べて表示するユーティリティも用意しました。

バナナに含まれるカリウム40による内部被ばくを基準にしたバナナ等価線量 (BED) のような指標から、胸部X線写真、自然放射線まで、様々なスケールで比較できます。

ITEMS = [
    Item("banana", "🍌 バナナ等価線量(BED)", 0.0001, None, "1本 ≈ 0.1 μSv(BED。ネタ)"),
    Item("cxr", "🫁 胸部X線(正面1枚)", 0.02, (0.01, 0.10), "施設・体位・装置で±幅あり"),
    Item("dental", "🦷 デンタルX線(小)", 0.005, (0.002, 0.01), "撮影法で差あり"),
    Item("mammo", "🩺 マンモグラフィ(両側)", 0.4, (0.2, 0.6), "装置/圧迫/体格で変動"),
    Item("ct_head", "🧠 CT 頭部", 2.0, (1.0, 3.0), "プロトコル依存"),
    Item("ct_chest", "🫁 CT 胸部", 7.0, (4.0, 10.0), "低線量化で下がる場合あり"),
    Item("ct_ap", "🫓 CT 腹骨盤", 9.0, (8.0, 12.0), "撮影範囲で上下"),
    Item("flight_tyo_nyc", "✈️ 長距離フライト(東京↔NY)", 0.08, (0.05, 0.10), "高度・太陽活動・ルートで増減"),
    Item("background_day", "🌏 自然放射線(1日分)", 2.4/365.0, None, "世界平均2.4 mSv/年の概算"),
]

8. まとめ

今回は、Pythonを使って胸部CTの被ばく線量を評価する一連のプロセスを解説しました。

  • ICRPの定義に基づいて、吸収線量から実効線量へとステップを踏んで計算した。

  • SSDEを導入することで、患者さん一人ひとりの体格を考慮した線量評価の重要性を示した。

  • 臓器ごとの係数を用いて、より現実に近い形での実効線量合成を試みた。

  • 計算結果を様々な指標と比較することで、その線量の定性的な大きさを把握することができた。

もちろん、今回のコードは簡易的なモデルであり、実際の線量評価はモンテカルロシミュレーションなど、より高度な手法で行われます。しかし、線量評価の基本的な考え方や、各種指標がどのように関連しているのかを学ぶための第一歩としては、非常に有用なツールになったのではないでしょうか。

ぜひ、皆さんもこのコードを触ってみて、入力値を変えながら「こんな場合はどうなるんだろう?」とシミュレーションを楽しんでみてください!

この記事が、日々の業務や学習の助けになれば幸いです。最後までお読みいただき、ありがとうございました!

★配布用コード

初心者に分かりやすくコメントアウトで完全解説したフルコードです。
そのままコピペで張り付けてお使いください。

# ===================================================================
# 1. ライブラリのインポート
# ===================================================================
import math
import pandas as pd
import matplotlib.pyplot as plt
from IPython.display import display, HTML

# (任意)日本語フォントを使う場合は有効化してください
!pip install japanize-matplotlib
import japanize_matplotlib

# ===================================================================
# 2. ICRP103に基づく「定数」の定義
# ===================================================================

# 2-1. 放射線荷重係数 (w_R) - X線は1
W_R = {
    "xray": 1.0,
    "photon": 1.0,
}

# 2-2. 組織荷重係数 (w_T) - ICRP103
W_T = {
    # 0.12
    "red_bone_marrow": 0.12,
    "colon":           0.12,
    "lung":            0.12,
    "stomach":         0.12,
    "breast":          0.12,
    "remainder":       0.12,
    # 0.08
    "gonads":          0.08,
    # 0.04
    "bladder":         0.04,
    "esophagus":       0.04,
    "liver":           0.04,
    "thyroid":         0.04,
    # 0.01
    "bone_surface":    0.01,
    "brain":           0.01,
    "salivary_glands": 0.01,
    "skin":            0.01,
}

REMAINDER_TISSUES_13 = [
    "adrenals", "extrathoracic_region", "gall_bladder", "heart", "kidneys",
    "lymphatic_nodes", "muscle", "oral_mucosa", "pancreas",
    "prostate_uterus", "small_intestine", "spleen", "thymus"
]

# ===================================================================
# 3. 単位変換と便利ツール
# ===================================================================

def gy_to_sv(absorbed_dose_gy: float, radiation: str = "xray") -> float:
    """吸収線量(Gy)→等価線量(Sv)。X線では数値は同じ(単位が変わる)。"""
    wR = W_R.get(radiation, 1.0)
    return absorbed_dose_gy * wR

CT_ITEMS = [
    {"id": "ct_head",  "label": "🧠 CT 頭部",   "dose_msv": 2.0},
    {"id": "ct_chest", "label": "🫁 CT 胸部",   "dose_msv": 7.0},
    {"id": "ct_ap",    "label": "🫓 CT 腹骨盤", "dose_msv": 9.0},
]

def make_ct_equiv_table(total_msv: float) -> pd.DataFrame:
    rows = []
    for it in CT_ITEMS:
        unit = it["dose_msv"]
        cnt = (total_msv / unit) if unit > 0 else float("nan")
        rows.append({
            "項目": it["label"],
            "1回あたりの実効線量[mSv]": it["dose_msv"],
            "換算すると(回分)": f"{cnt:.2f}",
        })
    return pd.DataFrame(rows)

# ===================================================================
# 4. CTDIvol・DLP・SSDE(サイズ補正)の計算ユーティリティ(胸部)
#    ※胸部は32 cm ファントムを前提
# ===================================================================

# 4-1. AAPM 220(=204): Dwベース f_size(32 cmファントム, 8〜45 cm 抜粋)
F_DW_32CM = {
     8: 2.76,  9: 2.66, 10: 2.57, 11: 2.47, 12: 2.38, 13: 2.30, 14: 2.22, 15: 2.14, 16: 2.06,
    17: 1.98, 18: 1.91, 19: 1.84, 20: 1.78, 21: 1.71, 22: 1.65, 23: 1.59, 24: 1.53, 25: 1.48,
    26: 1.43, 27: 1.37, 28: 1.32, 29: 1.28, 30: 1.23, 31: 1.19, 32: 1.14, 33: 1.10, 34: 1.06,
    35: 1.02, 36: 0.99, 37: 0.95, 38: 0.92, 39: 0.88, 40: 0.85, 41: 0.82, 42: 0.79, 43: 0.76,
    44: 0.74, 45: 0.71,
}

def _interp_f_dw_32cm(dw_cm: float) -> float:
    """AAPM 220(32 cm, Dw)の線形補間。範囲外は端でクリップ。"""
    xs = sorted(F_DW_32CM.keys())
    if dw_cm <= xs[0]:
        return F_DW_32CM[xs[0]]
    if dw_cm >= xs[-1]:
        return F_DW_32CM[xs[-1]]
    for i in range(len(xs)-1):
        x0, x1 = xs[i], xs[i+1]
        if x0 <= dw_cm <= x1:
            y0, y1 = F_DW_32CM[x0], F_DW_32CM[x1]
            t = (dw_cm - x0) / (x1 - x0)
            return y0 + (y1 - y0) * t
    return F_DW_32CM[xs[-1]]

def estimate_deff_cm(height_cm: float, weight_kg: float, ap_cm: float = None, lat_cm: float = None):
    """
    実効直径 Deff[cm] を返す。AP/LATがあれば sqrt(AP×LAT) を採用。
    無ければ BMI から Deff ≈ 0.6414*BMI + 12.976 を近似(胸部での代用)。
    """
    h_m = height_cm / 100.0 if height_cm else None
    bmi = (weight_kg / (h_m**2)) if (h_m and weight_kg) else None
    if ap_cm and lat_cm:
        deff = math.sqrt(ap_cm * lat_cm)
        return float(deff), float(bmi) if bmi else None
    if bmi is None:
        raise ValueError("AP/LAT または 身長・体重のいずれかを指定してください。")
    deff = 0.6414 * bmi + 12.976
    return float(deff), float(bmi)

def compute_ctdi_dlp_ssde(
    *,
    pitch: float = 1.2,
    scan_length_cm: float = 35.0,
    ctdiw_mGy: float = None,
    ctdi_vol_mGy: float = None,
    height_cm: float = 170.0,
    weight_kg: float = 65.0,
    ap_cm: float = None,
    lat_cm: float = None,
    phantom: str = "32cm",
) -> dict:
    """CTDIw/CTDIvol→DLP→SSDE(サイズ補正)→まとめ"""
    if ctdi_vol_mGy is None and ctdiw_mGy is None:
        raise ValueError("CTDIw または CTDIvol のどちらか一方を指定してください。")
    if ctdi_vol_mGy is None:
        ctdi_vol_mGy = ctdiw_mGy / pitch
    if ctdiw_mGy is None:
        ctdiw_mGy = ctdi_vol_mGy * pitch

    dlp_mGy_cm = ctdi_vol_mGy * scan_length_cm

    deff_cm, bmi = estimate_deff_cm(height_cm, weight_kg, ap_cm, lat_cm)
    if phantom != "32cm":
        raise NotImplementedError("本スクリプトは胸部=32 cmファントムのみ実装")
    f_size = _interp_f_dw_32cm(deff_cm)
    ssde_mGy = ctdi_vol_mGy * f_size

    return {
        "Pitch": pitch,
        "ScanLength_cm": scan_length_cm,
        "CTDIw_mGy": ctdiw_mGy,
        "CTDIvol_mGy": ctdi_vol_mGy,
        "DLP_mGy_cm": dlp_mGy_cm,
        "Height_cm": height_cm,
        "Weight_kg": weight_kg,
        "BMI": bmi,
        "Deff_cm": deff_cm,
        "f_size(32cm)": f_size,
        "SSDE_mGy": ssde_mGy,
    }

def estimate_effective_dose_from_dlp_chest_adult(dlp_mGy_cm: float, k_mSv_per_mGy_cm: float = 0.014) -> float:
    """成人胸部の代表換算:E ≈ k × DLP"""
    return float(dlp_mGy_cm * k_mSv_per_mGy_cm)

# ===================================================================
# 5. 胸部CTの被ばく線量シミュレーション・モデル(臓器合成)
# ===================================================================

# 5-1. 撮影範囲のモデル化(カバレッジ係数)
CHEST_COVERAGE = {
    "lung": 1.0,
    "breast": 1.0,
    "esophagus": 1.0,
    "red_bone_marrow": 0.25,
    "bone_surface":    0.25,
    "thyroid":         0.30,
    "skin":            0.50,
    "liver":           0.10,
    "stomach": 0.0, "colon": 0.0, "gonads": 0.0,
    "brain": 0.0, "salivary_glands": 0.0, "bladder": 0.0,
}

REMAINDER_COVERAGE = {
    "heart": 1.0,
    "thymus": 1.0,
    "lymphatic_nodes": 1.0,
    "muscle": 0.5,
    "adrenals": 0.0, "gall_bladder": 0.0, "kidneys": 0.0, "pancreas": 0.0,
    "prostate_uterus": 0.0, "small_intestine": 0.0, "spleen": 0.0,
    "extrathoracic_region": 0.2, "oral_mucosa": 0.1,
}

def compute_chest_ct_organs(dose_gy: float) -> pd.DataFrame:
    """
    胸部CTの平均吸収線量 dose_gy [Gy] を入力とし、
    各臓器 D_T → H_T → w_T 合成(ICRP103)を計算。
    """
    rows = []

    # 主要臓器
    for organ, frac in CHEST_COVERAGE.items():
        D_T = dose_gy * frac
        H_T = gy_to_sv(D_T, radiation="xray")
        wT = W_T.get(organ, 0.0)
        E_T = wT * H_T
        rows.append({
            "臓器": organ,
            "吸収線量[mGy]": D_T * 1000,
            "等価線量[mSv]": H_T * 1000,
            "w_T(ICRP103)": wT,
            "臓器寄与E[mSv]": E_T * 1000,
        })

    # remainder(13臓器の平均)
    D_list_rem = []
    H_list_rem = []
    for t in REMAINDER_TISSUES_13:
        frac = REMAINDER_COVERAGE.get(t, 0.0)
        D_T_rem = dose_gy * frac       # [Gy]
        H_T_rem = gy_to_sv(D_T_rem)    # [Sv]
        D_list_rem.append(D_T_rem)
        H_list_rem.append(H_T_rem)

    D_mean_rem = sum(D_list_rem) / len(REMAINDER_TISSUES_13)  # [Gy]
    H_mean_rem = sum(H_list_rem) / len(REMAINDER_TISSUES_13)  # [Sv]
    E_rem = W_T["remainder"] * H_mean_rem

    rows.append({
        "臓器": "remainder(13臓器平均)",
        "吸収線量[mGy]": D_mean_rem * 1000,
        "等価線量[mSv]": H_mean_rem * 1000,
        "w_T(ICRP103)": W_T["remainder"],
        "臓器寄与E[mSv]": E_rem * 1000,
    })

    df = pd.DataFrame(rows)
    df = df.sort_values("臓器寄与E[mSv]", ascending=False).reset_index(drop=True)
    return df

def compare_to_limits(total_effective_msv: float) -> pd.DataFrame:
    """(参考)線量限度との比較。※医療被ばくは適用外"""
    PUBLIC_LIMIT = 1.0
    WORKER_YEAR  = 50.0
    WORKER_5Y    = 100.0
    rows = [
        {"区分": "公衆(年間)", "基準[mSv]": PUBLIC_LIMIT, "今回の線量との比": total_effective_msv / PUBLIC_LIMIT},
        {"区分": "職業被ばく(年間)", "基準[mSv]": WORKER_YEAR, "今回の線量との比": total_effective_msv / WORKER_YEAR},
        {"区分": "職業被ばく(5年平均)", "基準[mSv]": WORKER_5Y / 5, "今回の線量との比": total_effective_msv / (WORKER_5Y / 5)},
    ]
    return pd.DataFrame(rows)

# ===================================================================
# 6. 表示ユーティリティ(比較・ALARAなど)
# ===================================================================
USV_PER_MSV = 1000.0

class Item:
    def __init__(self, id, label, dose_msv, range_msv=None, note=""):
        self.id = id
        self.label = label
        self.dose_msv = float(dose_msv)
        self.range_msv = range_msv
        self.note = note

ITEMS = [
    Item("banana", "🍌 バナナ等価線量(BED)", 0.0001, None, "1本 ≈ 0.1 μSv(BED。ネタ)"),
    Item("cxr", "🫁 胸部X線(正面1枚)", 0.02, (0.01, 0.10), "施設・体位・装置で±幅あり"),
    Item("dental", "🦷 デンタルX線(小)", 0.005, (0.002, 0.01), "撮影法で差あり"),
    Item("mammo", "🩺 マンモグラフィ(両側)", 0.4, (0.2, 0.6), "装置/圧迫/体格で変動"),
    Item("ct_head", "🧠 CT 頭部", 2.0, (1.0, 3.0), "プロトコル依存"),
    Item("ct_chest", "🫁 CT 胸部", 7.0, (4.0, 10.0), "低線量化で下がる場合あり"),
    Item("ct_ap", "🫓 CT 腹骨盤", 9.0, (8.0, 12.0), "撮影範囲で上下"),
    Item("flight_tyo_nyc", "✈️ 長距離フライト(東京↔NY)", 0.08, (0.05, 0.10), "高度・太陽活動・ルートで増減"),
    Item("background_day", "🌏 自然放射線(1日分)", 2.4/365.0, None, "世界平均2.4 mSv/年の概算"),
]

def as_count(total_msv: float, unit_msv: float) -> float:
    return float("nan") if unit_msv <= 0 else total_msv / unit_msv

def fmt_count(n: float) -> str:
    if n >= 1000:
        k = n/1000
        return f"約{round(k,1)}千回" if k < 10 else f"約{round(n/10000,1)}万回"
    if n >= 10:
        return f"{int(round(n))}回"
    if n >= 1:
        return f"{round(n,1)}回"
    if n >= 0.01:
        return f"{round(n,2)}回"
    return f"{n:.2e}回"

def fmt_range_text(total_msv: float, r):
    if not r:
        return ""
    lo = as_count(total_msv, r[0])
    hi = as_count(total_msv, r[1])
    return f"{fmt_count(lo)}{fmt_count(hi)}"

def make_dataframe(total_msv: float, show_ranges: bool=True) -> pd.DataFrame:
    rows = []
    for it in ITEMS:
        cnt = as_count(total_msv, it.dose_msv)
        rng = fmt_range_text(total_msv, it.range_msv) if show_ranges else ""
        rows.append({
            "項目": it.label,
            "1回 ≈ [mSv]": it.dose_msv,
            "換算(◯回分)": fmt_count(cnt),
            "換算レンジ": rng,
            "メモ": it.note,
        })
    return pd.DataFrame(rows)

def alara_badge(msv: float) -> str:
    if msv < 0.01:  return "🟢 Very Low"
    if msv < 1.0:   return "🟡 Low"
    if msv < 10.0:  return "🟠 Moderate"
    return "🔴 High"

def show_table_and_plot(total_msv: float):
    print(f"入力: {total_msv:.4f} mSv")
    print("ALARA:", alara_badge(total_msv))
    print("自然放射線(日)換算:", fmt_count(as_count(total_msv, 2.4/365.0)))
    df = make_dataframe(total_msv, show_ranges=True)
    display(df)

# ===================================================================
# 7. 入力(ここを自由に変更)
# ===================================================================

# --- 撮影条件 / 装置指標 ---
PITCH = 1.2
SCAN_LEN_CM = 35.0
CTDIW_INPUT_MGY   = 9.0     # 例:CTDIw(mGy)。CTDIvolしかない場合は None に
CTDIVOL_INPUT_MGY = None    # 例:CTDIvol(mGy)。CTDIwしかない場合は None に

# --- 体格(どちらか:AP/LAT か 身長・体重) ---
HEIGHT_CM = 170.0           # 身長(cm)
WEIGHT_KG = 65.0            # 体重(kg)
AP_CM = None                # scout等から測れれば入力
LAT_CM = None               # scout等から測れれば入力

# --- 臓器モデルへの入力モード ---
USE_SSDE_AS_INPUT = True    # True: SSDEを平均吸収線量に使用 / False: 手動で与える
MANUAL_MEAN_ABSORBED_DOSE_GY = 0.02  # Falseのときに使用(例:20 mGy)

# --- DLP→E 換算係数(成人胸部の代表値) ---
K_CHEST_MSV_PER_MGY_CM = 0.014

# ===================================================================
# 8. 計算本体
# ===================================================================

# (1) CTDI/DLP/SSDE を計算
dose_idx = compute_ctdi_dlp_ssde(
    pitch=PITCH, scan_length_cm=SCAN_LEN_CM,
    ctdiw_mGy=CTDIW_INPUT_MGY, ctdi_vol_mGy=CTDIVOL_INPUT_MGY,
    height_cm=HEIGHT_CM, weight_kg=WEIGHT_KG, ap_cm=AP_CM, lat_cm=LAT_CM,
    phantom="32cm",
)

# (2) 臓器モデルへ渡す平均吸収線量(Gy)
if USE_SSDE_AS_INPUT:
    EXAMPLE_CHEST_CT_DOSE_GY = dose_idx["SSDE_mGy"] / 1000.0
else:
    EXAMPLE_CHEST_CT_DOSE_GY = MANUAL_MEAN_ABSORBED_DOSE_GY

# (3) 臓器別→実効線量の合成
df_organs = compute_chest_ct_organs(EXAMPLE_CHEST_CT_DOSE_GY)
E_total_mSv = df_organs["臓器寄与E[mSv]"].sum()

# (4) 参考:DLP→E(成人胸部の代表換算)
E_from_DLP_mSv = estimate_effective_dose_from_dlp_chest_adult(
    dose_idx["DLP_mGy_cm"], k_mSv_per_mGy_cm=K_CHEST_MSV_PER_MGY_CM
)

# ===================================================================
# 9. 結果表示
# ===================================================================

print("=== 胸部CT: CTDI/DLP/SSDE の概算 ===")
display(pd.DataFrame([dose_idx]))

print("\n=== 臓器別の計算結果(ICRP103, 簡易カバレッジ) ===")
print(f"胸部CT 平均吸収線量(モデル入力): {EXAMPLE_CHEST_CT_DOSE_GY*1000:.1f} mGy")
display(df_organs[["臓器", "吸収線量[mGy]", "w_T(ICRP103)", "臓器寄与E[mSv]"]])

print("\n=== 実効線量(2系統) ===")
print(f"臓器合成 E ≈ {E_total_mSv:.2f} mSv")
print(f"DLP換算   E ≈ {E_from_DLP_mSv:.2f} mSv  (k={K_CHEST_MSV_PER_MGY_CM})")
if E_from_DLP_mSv > 0:
    diff = E_total_mSv - E_from_DLP_mSv
    pct  = diff / E_from_DLP_mSv * 100.0
    print(f"差分: {diff:+.2f} mSv  ({pct:+.1f} %)")

print("\n--- 年間線量限度との比較(参考:医療被ばくは適用外)---")
display(compare_to_limits(E_total_mSv))

print("\n--- 他のCT検査との比較(参考)---")
display(make_ct_equiv_table(E_total_mSv))

# グラフ(臓器寄与)
plot_df = df_organs[df_organs["臓器寄与E[mSv]"] > 0]
plt.figure(figsize=(10, 6))
plt.bar(plot_df["臓器"], plot_df["臓器寄与E[mSv]"])
plt.ylabel("実効線量寄与 [mSv]")
plt.title("胸部CT: 臓器別 実効線量寄与(ICRP103, SSDE入力, 簡易モデル)")
plt.xticks(rotation=60, ha="right")
plt.tight_layout()
plt.show()

# ===================================================================
# 10. 追加の参考表示(ALARA/換算表)
# ===================================================================
total_msv = E_total_mSv
show_table_and_plot(total_msv)
2
0
0

Register as a new user and use Qiita more conveniently

  1. You get articles that match your needs
  2. You can efficiently read back useful information
  3. You can use dark theme
What you can do with signing up
2
0

Delete article

Deleted articles cannot be recovered.

Draft of this article would be also deleted.

Are you sure you want to delete this article?