pip install ipympl numpy scipy scikit-image matplotlib pydicom pillow

!pip install ipympl numpy scipy scikit-image matplotlib pydicom pillow  #安裝套件

%matplotlib widget

import matplotlib as mpl
mpl.rcParams['font.sans-serif'] = [
    'Microsoft JhengHei',   # Windows（微軟正黑體）
    'PingFang TC',          # macOS
    'Heiti TC',             # macOS（舊版）
    'Noto Sans CJK TC',     # Linux（apt install fonts-noto-cjk）
    'SimHei',               # Windows（黑體）
    'Arial Unicode MS',
    'DejaVu Sans',          # fallback
]
mpl.rcParams['axes.unicode_minus'] = False

import numpy as np
import matplotlib.pyplot as plt
from matplotlib.widgets import Slider, Button, RadioButtons, CheckButtons, TextBox
from matplotlib.gridspec import GridSpec
from scipy import fft
from scipy.ndimage import map_coordinates, gaussian_filter1d
from skimage.transform import radon, iradon, rescale
from skimage.data import shepp_logan_phantom
from skimage import io, color

# （可選）檢查目前使用的字型：
# from matplotlib.font_manager import findfont, FontProperties
# print(findfont(FontProperties(family='Microsoft JhengHei')))

def normalize_image(image):
    image = image.astype(np.float64)
    lo, hi = image.min(), image.max()
    return (image - lo) / (hi - lo + 1e-8)

def load_grayscale_image(source, target_size=256):
    if isinstance(source, str):
        if source.lower().endswith(('.dcm', '.dicom')):
            import pydicom
            image = pydicom.dcmread(source).pixel_array
        else:
            image = io.imread(source)
    else:
        image = np.asarray(source)
    if image.ndim == 3:
        image = color.rgb2gray(image)
    image = normalize_image(image)
    return rescale(image, target_size / max(image.shape),
                   mode='reflect', anti_aliasing=True)

def ask_open_file_dialog():
    try:
        import tkinter as tk
        from tkinter import filedialog
        root = tk.Tk(); root.withdraw()
        path = filedialog.askopenfilename(
            filetypes=[("Image / DICOM", "*.png *.jpg *.jpeg *.tif *.tiff *.bmp *.dcm"),
                       ("All files", "*.*")])
        root.destroy()
        return path if path else None
    except Exception as e:
        print(f"無法開啟檔案對話框（{e}）。可改用 load_grayscale_image('path')")
        return None

class CTReconstructionDemo:
    def __init__(self, image, name="Phantom", n_angles=180):
        self.n_angles = n_angles
        self.theta = np.linspace(0, 180, n_angles, endpoint=False)
        self.current_index = 0
        self.is_playing = False
        self.recon_vmax = 1.0

        self.fig, (self.ax_image, self.ax_sino, self.ax_recon) = plt.subplots(
            1, 3, figsize=(6, 3))
        plt.subplots_adjust(bottom=0.28, wspace=0.3)

        self._build_widgets()
        self.set_image(image, name)

        self.timer = self.fig.canvas.new_timer(interval=50)
        self.timer.add_callback(self._advance)
        plt.show()

    def _build_widgets(self):
        ax_slider = plt.axes([0.2, 0.10, 0.6, 0.03])
        self.slider = Slider(ax_slider, '已用角度', 0, self.n_angles - 1,
                             valinit=0, valstep=1)
        self.slider.on_changed(self._on_slider_change)

        ax_play = plt.axes([0.20, 0.03, 0.10, 0.05])
        self.btn_play = Button(ax_play, '播放 / 暫停', color='lightgreen')
        self.btn_play.on_clicked(self._toggle_play)

        ax_load = plt.axes([0.35, 0.03, 0.12, 0.05])
        self.btn_load = Button(ax_load, '載入圖片', color='lightblue')
        self.btn_load.on_clicked(self._load_dialog)

    def set_image(self, image, name):
        self.image = load_grayscale_image(image, target_size=256)
        self.name = name
        self.sinogram = radon(self.image, theta=self.theta)
        self.sinogram_vmax = max(float(self.sinogram.max()), 1e-6)
        recon_full = iradon(self.sinogram, theta=self.theta)
        self.recon_vmax = max(float(np.percentile(recon_full, 99.5)), 1e-6)

        self.ax_image.clear()
        self.ax_image.imshow(self.image, cmap='gray', vmin=0, vmax=1)
        self.ax_image.set_title(f'原始圖：{name}')
        self.ax_image.axis('off')

        self.ax_sino.clear()
        self.im_sino = self.ax_sino.imshow(
            np.zeros_like(self.sinogram), cmap='gray', aspect='auto',
            vmin=0, vmax=self.sinogram_vmax,
            extent=[0, 180, 0, self.sinogram.shape[0]])
        self.ax_sino.set_title('Sinogram（投影數據）')
        self.ax_sino.set_xlabel('角度（度）')
        self.indicator = self.ax_sino.axvline(0, color='red', lw=2)

        self.ax_recon.clear()
        self.im_recon = self.ax_recon.imshow(
            np.zeros_like(self.image), cmap='gray',
            vmin=0, vmax=self.recon_vmax)
        self.ax_recon.set_title('重建結果（FBP）')
        self.ax_recon.axis('off')

        self.current_index = 0
        self.slider.set_val(0)
        self._update_display(0)

    def _update_display(self, index):
        index = int(index)
        visible = np.zeros_like(self.sinogram)
        visible[:, :index + 1] = self.sinogram[:, :index + 1]
        self.im_sino.set_data(visible)
        self.indicator.set_xdata([self.theta[index], self.theta[index]])

        used_theta = self.theta[:index + 1]
        recon = iradon(self.sinogram[:, :index + 1], theta=used_theta)
        self.im_recon.set_data(recon)
        self.im_recon.set_clim(0, self.recon_vmax)
        self.ax_recon.set_title(f'重建結果（已用 {index + 1} 個角度）')
        self.fig.canvas.draw_idle()

    def _on_slider_change(self, val):
        self.current_index = int(val)
        self._update_display(self.current_index)

    def _toggle_play(self, _event):
        if self.is_playing:
            self.timer.stop()
        else:
            self.timer.start()
        self.is_playing = not self.is_playing

    def _advance(self):
        self.current_index = (self.current_index + 1) % self.n_angles
        self.slider.set_val(self.current_index)

    def _load_dialog(self, _event):
        self.timer.stop(); self.is_playing = False
        path = ask_open_file_dialog()
        if path:
            self.set_image(path, '自訂影像')

ct_demo = CTReconstructionDemo(shepp_logan_phantom(), name='Shepp-Logan Phantom')

def sense_1d_reconstruction(image, sens_maps, acceleration):
    Ny, Nx = image.shape
    n_coils = len(sens_maps)
    R = acceleration
    if R == 1:
        recon = np.zeros((Ny, Nx), dtype=complex)
        for s in sens_maps:
            recon += fft.ifft2(fft.fft2(image * s)) * np.conj(s)
        weight = sum(np.abs(s) ** 2 for s in sens_maps) + 1e-10
        return np.abs(recon / weight)
    if Ny % R != 0:
        raise ValueError(f"Ny ({Ny}) 必須是 acceleration ({R}) 的整數倍")
    Ny_alias = Ny // R

    aliased = []
    for s in sens_maps:
        kspace = fft.fft2(image * s)
        kspace_sampled = np.zeros_like(kspace)
        kspace_sampled[::R, :] = kspace[::R, :]
        aliased.append(fft.ifft2(kspace_sampled))

    recon = np.zeros((Ny, Nx), dtype=complex)
    for x in range(Nx):
        for y0 in range(Ny_alias):
            S = np.zeros((n_coils, R), dtype=complex)
            I = np.zeros(n_coils, dtype=complex)
            for c in range(n_coils):
                for p in range(R):
                    y_orig = y0 + p * Ny_alias
                    if y_orig < Ny:
                        S[c, p] = sens_maps[c][y_orig, x]
                I[c] = R * aliased[c][y0, x]
            try:
                rho, *_ = np.linalg.lstsq(S, I, rcond=None)
            except np.linalg.LinAlgError:
                rho = np.zeros(R, dtype=complex)
            for p in range(R):
                y_orig = y0 + p * Ny_alias
                if y_orig < Ny:
                    recon[y_orig, x] = rho[p]
    return np.abs(recon)


class MRIKSpaceDemo:
    N_COILS = 4
    def __init__(self, image, name="Brain Phantom"):
        self.size = 256
        self.use_parallel = False
        self.image = None
        self.kspace = None
        self.fig = plt.figure(figsize=(6, 4))
        self._build_widgets()
        self.set_image(image, name)
        plt.show()

    def _build_widgets(self):
        ax_slider = plt.axes([0.25, 0.15, 0.45, 0.03])
        self.slider = Slider(ax_slider, '採樣比例', 0.1, 1.0, valinit=0.4)
        self.slider.on_changed(lambda _v: self._refresh())

        ax_radio = plt.axes([0.05, 0.02, 0.12, 0.12], facecolor='#f8f8f8')
        self.radio = RadioButtons(ax_radio, ('中心密集', '均勻採樣', '隨機採樣'))
        self.radio.on_clicked(lambda _v: self._refresh())

        ax_check = plt.axes([0.20, 0.02, 0.22, 0.08], facecolor='#f8f8f8')
        self.check = CheckButtons(ax_check, ('開啟 SENSE 並行模式',), (False,))
        self.check.on_clicked(self._toggle_mode)

        ax_btn = plt.axes([0.85, 0.03, 0.10, 0.05])
        self.btn_load = Button(ax_btn, '載入圖片')
        self.btn_load.on_clicked(self._load_dialog)

    def set_image(self, image, name):
        self.image = load_grayscale_image(image, target_size=self.size)
        self.kspace = fft.fftshift(fft.fft2(self.image))
        self.name = name
        self._refresh()

    def _toggle_mode(self, _label):
        self.use_parallel = bool(self.check.get_status()[0])
        self._refresh()

    def _build_mask(self, ratio, mode):
        res = self.size
        mask = np.zeros((res, res), dtype=bool)
        if mode == '中心密集':
            r = max(2, int(res * ratio / 2))
            mask[res // 2 - r: res // 2 + r, res // 2 - r: res // 2 + r] = True
        elif mode == '均勻採樣':
            step = max(1, int(round(1 / ratio)))
            mask[::step, :] = True
        else:
            mask = np.random.rand(res, res) < ratio
        return mask

    def _refresh(self):
        widget_axes = [self.slider.ax, self.radio.ax,
                       self.check.ax, self.btn_load.ax]
        for ax in list(self.fig.axes):
            if ax not in widget_axes:
                ax.remove()
        ratio = float(self.slider.val)
        if self.use_parallel:
            self._draw_sense(ratio)
        else:
            mask = self._build_mask(ratio, self.radio.value_selected)
            aliased = np.abs(fft.ifft2(fft.ifftshift(self.kspace * mask)))
            self._draw_basic(mask, aliased)
        self.fig.canvas.draw_idle()

    def _draw_basic(self, mask, aliased):
        res = self.size
        low_mask = np.zeros((res, res), dtype=bool)
        low_mask[res // 2 - 12: res // 2 + 12,
                 res // 2 - 12: res // 2 + 12] = True
        low = np.abs(fft.ifft2(fft.ifftshift(self.kspace * low_mask)))
        high = np.abs(fft.ifft2(fft.ifftshift(self.kspace * (~low_mask))))
        panels = [
            (self.image, '原始解剖圖', 'gray'),
            (np.log(1 + np.abs(self.kspace)), 'K-space (log)', 'magma'),
            (mask, '採樣遮罩', 'gray'),
            (aliased, '零填充重建', 'gray'),
            (low, '低頻（對比）', 'gray'),
            (high, '高頻（邊緣）', 'gray'),
            (np.angle(self.kspace), '相位圖', 'hsv'),
            (np.abs(self.image - aliased), '重建誤差', 'gray'),
        ]
        gs = GridSpec(2, 4, figure=self.fig)
        plt.subplots_adjust(bottom=0.25, hspace=0.3)
        for i, (data, title, cmap) in enumerate(panels):
            ax = self.fig.add_subplot(gs[i // 4, i % 4])
            if i == 6:
                ax.imshow(data, cmap=cmap, vmin=-np.pi, vmax=np.pi)
            else:
                ax.imshow(data, cmap=cmap)
            ax.set_title(title, fontsize=10)
            ax.axis('off')

    def _draw_sense(self, ratio):
        res = self.size
        x, y = np.meshgrid(np.linspace(-1, 1, res), np.linspace(-1, 1, res))
        angles = np.linspace(0, 2 * np.pi, self.N_COILS, endpoint=False)
        sens_maps = [
            np.exp(-((x - 0.7 * np.cos(a)) ** 2 + (y - 0.7 * np.sin(a)) ** 2) / 0.3)
            for a in angles
        ]
        R = int(round(1 / ratio))
        R = max(1, min(self.N_COILS, R))
        if R == 1:
            recon = self.image.copy()
        else:
            recon = sense_1d_reconstruction(self.image, sens_maps, R)

        raw_mask = np.zeros((res, res), dtype=bool)
        raw_mask[::R, :] = True
        display_mask = fft.fftshift(raw_mask)

        coil_images = []
        for s in sens_maps:
            kspace = fft.fft2(self.image * s)
            ks = np.zeros_like(kspace)
            ks[::R, :] = kspace[::R, :]
            coil_images.append(np.abs(fft.ifft2(ks)))

        error = np.abs(self.image - recon)
        gs = GridSpec(3, 4, figure=self.fig)
        plt.subplots_adjust(bottom=0.25, hspace=0.45, wspace=0.3)

        row0 = [
            (self.image, '原始影像', 'gray'),
            (display_mask, f'欠採樣遮罩（R={R}）', 'bone'),
            (error, '重建誤差', 'hot'),
            (recon, f'SENSE 重建（R={R}）', 'gray'),
        ]
        for i, (data, title, cmap) in enumerate(row0):
            ax = self.fig.add_subplot(gs[0, i])
            ax.imshow(data, cmap=cmap)
            ax.set_title(title, fontsize=9)
            ax.axis('off')

        for i in range(self.N_COILS):
            ax_s = self.fig.add_subplot(gs[1, i])
            ax_s.imshow(sens_maps[i], cmap='viridis')
            ax_s.set_title(f'線圈 {i + 1} 敏感度', fontsize=9)
            ax_s.axis('off')
            ax_i = self.fig.add_subplot(gs[2, i])
            ax_i.imshow(coil_images[i], cmap='gray')
            ax_i.set_title(f'線圈 {i + 1} 折疊圖', fontsize=9)
            ax_i.axis('off')

    def _load_dialog(self, _event):
        path = ask_open_file_dialog()
        if path:
            self.set_image(path, '自訂影像')


def make_brain_phantom(size=256):
    x, y = np.meshgrid(np.linspace(-1, 1, size), np.linspace(-1, 1, size))
    img = np.zeros((size, size))
    img[(x ** 2 / 0.75 ** 2 + y ** 2 / 0.85 ** 2) < 1] = 0.2
    img[(x ** 2 / 0.60 ** 2 + y ** 2 / 0.75 ** 2) < 1] = 0.5
    img[((x - 0.2) ** 2 + (y - 0.3) ** 2) < 0.05] = 0.8
    img += np.random.normal(0, 0.01, (size, size))
    return np.clip(img, 0, 1)

mri_demo = MRIKSpaceDemo(make_brain_phantom(), name='Brain Phantom')

def tof_resolution_to_sigma(fwhm_ps, c_mm_per_ps=0.3):
    """時間解析度 FWHM (ps) → 空間 σ (mm)"""
    return c_mm_per_ps * fwhm_ps / (2 * 2.355)


def _build_tof_kernel(t_fine_px, t_centers_px, sigma_px):
    """
    (n_tof, n_fine) TOF 權重矩陣，列和 = 1。
    嚴格伴隨對的核心：forward 和 backward 用同一個矩陣。
    """
    diff = t_fine_px[None, :] - t_centers_px[:, None]  # (n_tof, n_fine)
    G = np.exp(-diff ** 2 / (2 * sigma_px ** 2))
    G[np.abs(diff) > 4 * sigma_px] = 0.0                # 截斷，加速
    col_sum = G.sum(axis=0, keepdims=True) + 1e-10
    return G / col_sum


def _sample_along_lors(image, theta, meta):
    """
    沿所有 LOR 做 fine 採樣（bilinear 插值）。
    回傳 (n_det, n_fine, n_angles)。
    """
    Ny, Nx = image.shape
    R_px, T_px = np.meshgrid(meta['r_indices_px'], meta['t_centers_fine_px'],
                              indexing='ij')
    center_y_px, center_x_px = (Ny - 1) / 2, (Nx - 1) / 2
    th_rad = np.deg2rad(theta)
    n_det, n_fine = R_px.shape
    n_angles = len(theta)

    profiles = np.zeros((n_det, n_fine, n_angles))
    for i, th in enumerate(th_rad):
        c, s = np.cos(th), np.sin(th)
        x_px = center_x_px + R_px * c - T_px * s
        y_px = center_y_px + R_px * s + T_px * c
        profiles[:, :, i] = map_coordinates(
            image, [y_px, x_px], order=1, mode='constant', cval=0)
    return profiles


def make_tof_sinogram(image, theta, sigma_tof_mm, fov_mm, step_px=4):
    """嚴格生成 TOF sinogram（與 tof_forward 一致）。"""
    Ny, Nx = image.shape
    px_mm = fov_mm / Ny
    max_r_px = int(np.ceil(fov_mm * np.sqrt(2) / 2 / px_mm))
    r_indices = np.arange(-max_r_px, max_r_px + 1)

    max_t_px = max_r_px
    t_fine_px = np.arange(-max_t_px, max_t_px + 1)
    t_centers_px = np.arange(-max_t_px + step_px // 2,
                              max_t_px + 1, step_px)

    sigma_px = sigma_tof_mm / px_mm
    G_ker = _build_tof_kernel(t_fine_px, t_centers_px, sigma_px)

    meta = {
        'r_indices_px':      r_indices,
        't_centers_fine_px': t_fine_px,       # ← 所有 fine 位置
        't_centers_px':      t_centers_px,    # ← 粗 bin 中心
        'px_mm':             px_mm,
        'fov_mm':            fov_mm,
        'sigma_px':          sigma_px,
        'step_px':           step_px,
        'max_t_px':          max_t_px,
        'n_fine':            len(t_fine_px),
        'n_tof':             len(t_centers_px),
        'bin_mm':            step_px * px_mm,
        'G_ker':             G_ker,           # ← TOF 權重矩陣
    }

    profiles = _sample_along_lors(image, theta, meta)     # (n_det, n_fine, n_angles)
    # 注意：einsum 索引 t 是 n_fine、k 是 n_tof
    tof_sino = np.einsum('kt,rti->kri', G_ker, profiles) * px_mm
    return tof_sino, meta


def tof_forward(image, theta, meta):
    """TOF forward：與 make_tof_sinogram 完全一致。"""
    profiles = _sample_along_lors(image, theta, meta)
    return np.einsum('kt,rti->kri', meta['G_ker'], profiles) * meta['px_mm']


def tof_backward(tof_sino, theta, image_shape, meta):
    """
    TOF backward：嚴格伴隨。
    先反投影回 fine 位置，再 scatter 到像素。
    """
    Ny, Nx = image_shape
    # 反投影到 fine：fine[t, r, θ] = Σ_k G_ker[k, t] · tof_sino[k, r, θ] · px_mm
    fine = np.einsum('kt,kri->tri', meta['G_ker'], tof_sino) * meta['px_mm']

    R_px, T_px = np.meshgrid(meta['r_indices_px'], meta['t_centers_fine_px'],
                              indexing='ij')
    center_y_px, center_x_px = (Ny - 1) / 2, (Nx - 1) / 2
    th_rad = np.deg2rad(theta)

    recon_flat = np.zeros(Ny * Nx)
    for i, th in enumerate(th_rad):
        c, s = np.cos(th), np.sin(th)
        x_px = center_x_px + R_px * c - T_px * s
        y_px = center_y_px + R_px * s + T_px * c
        xi = np.round(x_px).astype(int).ravel()
        yi = np.round(y_px).astype(int).ravel()
        # fine[:, :, i] 形狀 (n_fine, n_det)，.T 後 (n_det, n_fine) 與 R_px 對齊
        vals = fine[:, :, i].T.ravel()
        valid = (xi >= 0) & (xi < Nx) & (yi >= 0) & (yi < Ny)
        flat_idx = yi[valid] * Nx + xi[valid]
        recon_flat += np.bincount(flat_idx, weights=vals[valid],
                                   minlength=Ny * Nx)
    return recon_flat.reshape((Ny, Nx))


def radon_forward(image, theta, meta):
    """非 TOF forward：TOF forward 對 t 加總即得（因列和為 1）。"""
    profiles = _sample_along_lors(image, theta, meta)
    return profiles.sum(axis=1) * meta['px_mm']


def radon_backward(sino, theta, image_shape, meta):
    """非 TOF backward。"""
    Ny, Nx = image_shape
    R_px, T_px = np.meshgrid(meta['r_indices_px'], meta['t_centers_fine_px'],
                              indexing='ij')
    center_y_px, center_x_px = (Ny - 1) / 2, (Nx - 1) / 2
    th_rad = np.deg2rad(theta)
    n_fine = len(meta['t_centers_fine_px'])

    recon_flat = np.zeros(Ny * Nx)
    for i, th in enumerate(th_rad):
        c, s = np.cos(th), np.sin(th)
        x_px = center_x_px + R_px * c - T_px * s
        y_px = center_y_px + R_px * s + T_px * c
        xi = np.round(x_px).astype(int).ravel()
        yi = np.round(y_px).astype(int).ravel()
        vals = np.repeat(sino[:, i] * meta['px_mm'], n_fine)
        valid = (xi >= 0) & (xi < Nx) & (yi >= 0) & (yi < Ny)
        flat_idx = yi[valid] * Nx + xi[valid]
        recon_flat += np.bincount(flat_idx, weights=vals[valid],
                                   minlength=Ny * Nx)
    return recon_flat.reshape((Ny, Nx))

class PETIterativeDemo:
    """PET 迭代重建互動展示：MLEM / OSEM、衰減、TOF。"""

    def __init__(self, image, name='Phantom', n_angles=60,
                 image_size=128, fov_mm=256.0):
        self.size = image_size
        self.fov_mm = fov_mm
        self.n_angles = n_angles
        self.theta = np.linspace(0, 180, n_angles, endpoint=False)
        self.algorithm = 'MLEM'
        self.atten_coeff = 0.005
        self.tof_fwhm_ps = 400.0
        self.history = []

        self.fig, (self.ax_image, self.ax_sino, self.ax_recon) = plt.subplots(
            1, 3, figsize=(6, 4))
        plt.subplots_adjust(bottom=0.35, wspace=0.3)

        self._build_widgets()
        self.set_image(image, name)
        plt.show()

    # ---------------------------------------------------------------
    # UI
    # ---------------------------------------------------------------
    def _build_widgets(self):
        ax_check = plt.axes([0.05, 0.05, 0.16, 0.12], facecolor='#f9f9f9')
        self.check = CheckButtons(
            ax_check,
            ('模擬衰減', '衰減補償', '啟用 TOF'),
            (False, False, False))
        self.check.on_clicked(self._toggle_physics)

        ax_box = plt.axes([0.10, 0.20, 0.06, 0.03])
        self.text_box = TextBox(ax_box, 'μ (1/mm) ',
                                 initial=str(self.atten_coeff))
        self.text_box.on_submit(self._submit_atten)

        ax_tof_box = plt.axes([0.24, 0.20, 0.06, 0.03])
        self.tof_box = TextBox(ax_tof_box, 'TOF FWHM (ps) ',
                                initial=str(self.tof_fwhm_ps))
        self.tof_box.on_submit(self._submit_tof_fwhm)

        ax_radio = plt.axes([0.24, 0.05, 0.08, 0.10], facecolor='#f9f9f9')
        self.radio = RadioButtons(ax_radio, ('MLEM', 'OSEM'))
        self.radio.on_clicked(self._change_algo)

        ax_slider = plt.axes([0.36, 0.16, 0.50, 0.03])
        self.slider = Slider(ax_slider, '迭代步', 0, 1, valinit=0, valstep=1)
        self.slider.on_changed(self._on_slider)

        ax_run = plt.axes([0.36, 0.05, 0.14, 0.05])
        self.btn_run = Button(ax_run, '跑一次迭代', color='#FFD700')
        self.btn_run.on_clicked(self._run_iteration)

        ax_load = plt.axes([0.52, 0.05, 0.12, 0.05])
        self.btn_load = Button(ax_load, '載入圖片', color='#ADD8E6')
        self.btn_load.on_clicked(self._load_dialog)

    # ---------------------------------------------------------------
    # 影像 / sinogram
    # ---------------------------------------------------------------
    def set_image(self, image, name):
        self.image = load_grayscale_image(image, target_size=self.size)
        self.name = name
        self._compute_sinogram()
        self._reset_iterations()

    def _compute_sinogram(self):
        """生成 TOF sinogram（含衰減與 Poisson 噪聲）。"""
        sigma_tof_mm = tof_resolution_to_sigma(self.tof_fwhm_ps)
        tof_true, self.tof_meta = make_tof_sinogram(
            self.image, self.theta, sigma_tof_mm, self.fov_mm)

        # 衰減
        simulate_atten = self.check.get_status()[0]
        if simulate_atten:
            body = (self.image > 0.02).astype(float)
            path = radon_forward(body, self.theta, self.tof_meta) * self.atten_coeff
            self.atten_map = np.exp(-path)
            tof_obs = tof_true * self.atten_map[None, :, :]
        else:
            self.atten_map = np.ones((len(self.tof_meta['r_indices_px']),
                                       self.n_angles))
            tof_obs = tof_true

        # Poisson 噪聲
        scale = 100
        self.tof_sinogram = np.random.poisson(tof_obs * scale) / scale
        self.sinogram_2d = self.tof_sinogram.sum(axis=0)

        # 靜態圖
        self.ax_image.clear()
        self.ax_image.imshow(self.image, cmap='hot')
        self.ax_image.set_title(f'原始圖：{self.name}')
        self.ax_image.axis('off')

        self.ax_sino.clear()
        self.ax_sino.imshow(self.sinogram_2d, cmap='magma', aspect='auto')
        self.ax_sino.set_title('觀測 Sinogram（非 TOF，即 TOF 對時間軸加總）')
        self.ax_sino.set_xlabel('角度（度）')

    def _reset_iterations(self):
        self.history = [np.ones_like(self.image) * 0.5]
        self.slider.valmax = 1
        self.slider.ax.set_xlim(0, 1)
        self.slider.set_val(0)
        self._show_slice(0)

    # ---------------------------------------------------------------
    # 迭代
    # ---------------------------------------------------------------
    def _run_iteration(self, _event):
        status = self.check.get_status()
        simulate_atten, apply_ac, use_tof = status
        apply_correction = apply_ac and simulate_atten

        if use_tof:
            new_recon = self._run_tof_iteration(apply_correction)
        else:
            new_recon = self._run_2d_iteration(apply_correction)

        if new_recon is None:
            return
        if not np.isfinite(new_recon).all() or new_recon.max() > 1e10:
            print('警告：重建發散，請降低衰減係數或改用 MLEM。')
            return

        self.history.append(new_recon)
        new_max = len(self.history) - 1
        self.slider.valmax = new_max
        self.slider.ax.set_xlim(0, new_max)
        self.slider.set_val(new_max)

    def _run_2d_iteration(self, apply_correction):
        """非 TOF 的 MLEM / OSEM。"""
        recon = self.history[-1].copy()
        n_subsets = 4 if self.algorithm == 'OSEM' else 1

        for subset in range(n_subsets):
            idx = np.arange(subset, self.n_angles, n_subsets)
            theta_sub = self.theta[idx]
            obs_sub = self.sinogram_2d[:, idx]
            if apply_correction:
                ac_sub = self.atten_map[:, idx]
            else:
                ac_sub = np.ones((len(self.tof_meta['r_indices_px']), len(idx)))

            forward = radon_forward(recon, theta_sub, self.tof_meta) * ac_sub + 1e-7
            ratio = np.nan_to_num(obs_sub / forward)

            correction = radon_backward(ratio * ac_sub, theta_sub,
                                        recon.shape, self.tof_meta)
            sens = radon_backward(ac_sub, theta_sub,
                                  recon.shape, self.tof_meta) + 1e-7
            recon = np.clip(recon * correction / sens, 0, None)
        return recon

    def _run_tof_iteration(self, apply_correction):
        """
        TOF-MLEM / TOF-OSEM 通用更新。
        MLEM 是 N=1 特例；OSEM 用 N=4 子集，每個子集各自用子集 sens 歸一化
        （Hudson & Larkin 1994 原版）。
        """
        recon = self.history[-1].copy()
        ac = self.atten_map if apply_correction else np.ones_like(self.atten_map)
        ac_3d = ac[None, :, :]

        n_subsets = 4 if self.algorithm == 'OSEM' else 1

        for subset in range(n_subsets):
            idx = np.arange(subset, self.n_angles, n_subsets)
            theta_sub = self.theta[idx]
            obs_sub = self.tof_sinogram[:, :, idx]
            ac_sub = ac_3d[:, :, idx]

            # 前向投影（子集版 A_TOF）
            forward = tof_forward(recon, theta_sub, self.tof_meta) * ac_sub + 1e-7
            ratio = np.nan_to_num(obs_sub / forward)

            # 反投影
            correction = tof_backward(ratio * ac_sub, theta_sub, recon.shape,
                                      self.tof_meta)

            # 子集敏感度 s_n = A_TOFᵀ · AC
            sens_sub = tof_backward(ac_sub * np.ones_like(obs_sub),
                                    theta_sub, recon.shape, self.tof_meta) + 1e-7

            # MLEM 更新
            recon = np.clip(recon * correction / sens_sub, 0, None)

        return recon

    # ---------------------------------------------------------------
    # 顯示
    # ---------------------------------------------------------------
    def _on_slider(self, val):
        self._show_slice(int(val))

    def _show_slice(self, index):
        if index >= len(self.history):
            return
        data = self.history[index]
        vmax = np.percentile(data, 99.5) or 1
        use_tof = self.check.get_status()[2]

        algo_name = f'{self.algorithm} (TOF)' if use_tof else self.algorithm

        self.ax_recon.clear()
        self.ax_recon.imshow(data, cmap='hot', vmin=0, vmax=vmax)
        if use_tof:
            self.ax_recon.set_title(
                f'重建圖（{algo_name}，第 {index} 次迭代）\n'
                f'TOF FWHM = {self.tof_fwhm_ps:.0f} ps')
        else:
            self.ax_recon.set_title(
                f'重建圖（{algo_name}，第 {index} 次迭代）')
        self.ax_recon.axis('off')
        self.fig.canvas.draw_idle()

    # ---------------------------------------------------------------
    # 輸入 / 開關
    # ---------------------------------------------------------------
    def _submit_atten(self, text):
        try:
            self.atten_coeff = float(text)
            self._compute_sinogram()
            self._reset_iterations()
        except ValueError:
            print('請輸入數字')

    def _submit_tof_fwhm(self, text):
        try:
            self.tof_fwhm_ps = float(text)
            self._compute_sinogram()
            self._reset_iterations()
        except ValueError:
            print('請輸入數字')

    def _toggle_physics(self, _label):
        self._compute_sinogram()
        self._reset_iterations()

    def _change_algo(self, label):
        self.algorithm = label
        self._reset_iterations()

    def _load_dialog(self, _event):
        path = ask_open_file_dialog()
        if path:
            self.set_image(path, '自訂影像')

pet_demo = PETIterativeDemo(shepp_logan_phantom(), name='Shepp-Logan Phantom')

項目	量測到什麼	重建方法
CT 電腦斷層	各角度的 X 光投影	濾波反投影（FBP）
MRI 磁振造影	頻率域訊號（k-space）	逆傅立葉變換 + SENSE
PET 正子造影	各角度的光子計數 + 時間差	迭代重建（MLEM / OSEM），可選 TOF

醫學影像重建 — 互動式教學¶

Part 0｜共用工具¶

Part 1｜CT：投影重建¶

Part 2｜MRI：K-space 重建與 SENSE 並行成像¶

Part 3｜PET：迭代重建（MLEM / OSEM）+ TOF¶

理論背景¶

標準 MLEM 更新式¶

衰減（Attenuation）與 TOF-PET¶