1β₁ = 1e-4,  β_T = 0.02,  T = 1000 (線性排程)
2αₜ = 1 - βₜ
3ᾱₜ = ∏ᵢ₌₁ᵗ αᵢ  (cumulative product)

 1# 1. 建立虛擬環境（使用 uv 或 conda）
 2uv venv .venv --python 3.10
 3source .venv/bin/activate
 4
 5# 2. 安裝 PyTorch（依 CUDA 版本調整）
 6uv pip install torch torchvision --index-url https://download.pytorch.org/whl/cu121
 7
 8# 3. 取得程式碼
 9git clone https://github.com/junbopeng/conditional_DDPM.git
10cd conditional_DDPM
11
12# 4. 安裝其他依賴（程式碼內含的標準庫）
13uv pip install numpy pillow

 1<dataset_name>/
 2├── train/
 3│   ├── a/          # 目標影像（如 pCT）
 4│   │   ├── 001.npy
 5│   │   ├── 002.npy
 6│   │   └── ...
 7│   └── b/          # 條件影像（如 CBCT）
 8│       ├── 001.npy
 9│       ├── 002.npy
10│       └── ...
11└── test/
12    ├── a/          # 測試目標影像
13    └── b/          # 測試條件影像

 1import numpy as np
 2
 3# 假設原始 DICOM 已轉為 numpy array，HU 值範圍 [-1000, 3000]
 4def preprocess_slice(hu_array, hu_min=-1000, hu_max=3000):
 5    """將 HU 值正規化到 [-1, 1] 範圍"""
 6    normalized = (hu_array - hu_min) / (hu_max - hu_min)  # [0, 1]
 7    normalized = normalized * 2 - 1                         # [-1, 1]
 8    return normalized.astype(np.float32)
 9
10# 儲存為 .npy
11slice_ct = preprocess_slice(ct_dicom_array)
12slice_cbct = preprocess_slice(cbct_dicom_array)
13np.save(f"brain/train/a/{idx:04d}.npy", slice_ct)    # pCT
14np.save(f"brain/train/b/{idx:04d}.npy", slice_cbct)  # CBCT

 1# Train_condition.py — 條件式 DDPM 訓練
 2import os
 3import sys
 4import time
 5import datetime
 6import torch
 7import torch.optim as optim
 8from torch.utils.data import DataLoader
 9from torch.autograd import Variable
10
11from Diffusion_condition import GaussianDiffusionTrainer_cond
12from Model_condition import UNet
13from datasets import ImageDataset  # 通用版 dataset
14
15# ===== 超參數設定 =====
16dataset_name = "brain"       # 資料集資料夾名稱
17out_name     = "trial_1"     # 實驗名稱
18batch_size   = 2             # 批次大小（依 GPU VRAM 調整）
19T            = 1000          # 擴散步數
20ch           = 128           # 基礎通道數
21ch_mult      = [1, 2, 3, 4]  # 通道倍率
22attn         = [2]           # 加入 Attention 的層級索引
23num_res_blocks = 2           # 每層 ResBlock 數量
24dropout      = 0.3           # Dropout 比率
25lr           = 1e-4          # 學習率
26n_epochs     = 1000          # 訓練 epoch 數
27beta_1       = 1e-4          # 噪聲排程起始值
28beta_T       = 0.02          # 噪聲排程結束值
29grad_clip    = 1.0           # 梯度裁剪閾值
30
31# ===== 初始化 =====
32device = torch.device("cuda")
33Tensor = torch.cuda.FloatTensor
34save_weight_dir = f"./Checkpoints/{out_name}"
35os.makedirs(save_weight_dir, exist_ok=True)
36
37# 資料載入器
38train_dataloader = DataLoader(
39    ImageDataset(f"./{dataset_name}", transforms_=False, unaligned=True),
40    batch_size=batch_size,
41    shuffle=True,
42    num_workers=0,
43)
44
45# 模型與優化器
46net_model = UNet(T, ch, ch_mult, attn, num_res_blocks, dropout).to(device)
47optimizer = optim.AdamW(
48    net_model.parameters(),
49    lr=lr,
50    betas=(0.9, 0.999),
51    eps=1e-8,
52    weight_decay=0,
53)
54trainer = GaussianDiffusionTrainer_cond(net_model, beta_1, beta_T, T).to(device)
55
56# ===== 訓練迴圈 =====
57prev_time = time.time()
58for epoch in range(1, n_epochs + 1):
59    loss_save = 0
60    for i, batch in enumerate(train_dataloader, 1):
61        optimizer.zero_grad()
62        ct   = Variable(batch["a"].type(Tensor))    # 目標 CT（key 依 dataset 版本可能是 "a" 或 "pCT"）
63        cbct = Variable(batch["b"].type(Tensor))     # 條件 CBCT（"b" 或 "CBCT"）
64        x_0  = torch.cat((ct, cbct), 1)              # 2-channel 輸入
65        loss = trainer(x_0)
66        loss_save += loss / 65536                     # 正規化 loss
67        loss.backward()
68        torch.nn.utils.clip_grad_norm_(net_model.parameters(), grad_clip)
69        optimizer.step()
70    loss_save /= i
71
72    # 計時與日誌
73    duration  = datetime.timedelta(seconds=(time.time() - prev_time))
74    eta       = datetime.timedelta(seconds=(n_epochs - epoch) * (time.time() - prev_time))
75    prev_time = time.time()
76    sys.stdout.write(
77        f"\r[Epoch {epoch}/{n_epochs}] [ETA: {eta}] [Duration: {duration}] "
78        f"[MSE Loss: {loss_save.item():.6f}]"
79    )
80
81    # 每 10 epoch（100 之後）儲存 checkpoint
82    if epoch > 100 and epoch % 10 == 0:
83        torch.save(
84            net_model.state_dict(),
85            os.path.join(save_weight_dir, f"ckpt_{epoch}_.pt"),
86        )

 1# Test_condition.py — 條件式 DDPM 推論
 2import os
 3import numpy as np
 4import torch
 5from torch.utils.data import DataLoader
 6from torch.autograd import Variable
 7
 8from Diffusion_condition import GaussianDiffusionSampler_cond
 9from Model_condition import UNet
10from datasets import ImageDataset
11
12# ===== 參數（須與訓練時一致） =====
13dataset_name = "brain"
14out_name     = "trial_1"
15T            = 1000
16ch           = 128
17ch_mult      = [1, 2, 3, 4]
18attn         = [2]
19num_res_blocks = 2
20dropout      = 0.3
21beta_1       = 1e-4
22beta_T       = 0.02
23
24device = torch.device("cuda:0")
25Tensor = torch.cuda.FloatTensor
26
27# ===== 載入訓練好的模型 =====
28net_model = UNet(T, ch, ch_mult, attn, num_res_blocks, dropout).to(device)
29ckpt = torch.load(f"./Checkpoints/{out_name}/ckpt_1000_.pt", map_location=device)
30net_model.load_state_dict(ckpt)
31net_model.eval()
32print("Model weights loaded.")
33
34# ===== 建立 Sampler =====
35sampler = GaussianDiffusionSampler_cond(net_model, beta_1, beta_T, T).to(device)
36
37# ===== 測試資料 =====
38test_dataloader = DataLoader(
39    ImageDataset(f"./{dataset_name}", transforms_=False, unaligned=True, mode="test"),
40    batch_size=1,
41    shuffle=False,
42    num_workers=0,
43)
44
45# ===== 推論迴圈 =====
46os.makedirs(f"test/{out_name}", exist_ok=True)
47results = []
48
49with torch.no_grad():
50    for idx, batch in enumerate(test_dataloader, 1):
51        print(f"Processing slice {idx}...")
52        cbct = Variable(batch["b"].type(Tensor))     # 條件 CBCT
53
54        # 從純噪聲開始
55        noisy_image = torch.randn(1, 1, 256, 256, device=device)
56        x_in = torch.cat((noisy_image, cbct), dim=1)  # [noise, CBCT] → 2-ch
57
58        # 反向擴散：T=1000 步去噪
59        x_out = sampler(x_in)
60
61        # 提取生成的合成 CT（第 0 channel）
62        synthetic_ct = x_out[:, 0:1, :, :]  # shape: [1, 1, 256, 256]
63
64        # 儲存結果
65        sct_np = synthetic_ct.cpu().numpy().squeeze()
66        np.save(f"test/{out_name}/sCT_{idx:04d}.npy", sct_np)
67        results.append(sct_np)
68
69print(f"Done. Generated {len(results)} synthetic CT slices.")

 1import torch
 2import numpy as np
 3from Diffusion_condition import GaussianDiffusionSampler_cond
 4from Model_condition import UNet
 5
 6def estimate_uncertainty(
 7    model_path: str,
 8    cbct_slice: np.ndarray,
 9    n_samples: int = 10,
10    device: str = "cuda:0",
11) -> dict:
12    """
13    對同一張 CBCT 進行多次取樣，計算 pixel-wise 不確定性。
14
15    Args:
16        model_path: 訓練好的 checkpoint 路徑
17        cbct_slice: 單張 CBCT 影像 (H, W)，已正規化至 [-1, 1]
18        n_samples: 取樣次數（越多越穩定，建議 10-20）
19        device: GPU 設備
20
21    Returns:
22        dict 包含 mean_sCT, std_map, all_samples
23    """
24    # 載入模型
25    T, ch, ch_mult, attn = 1000, 128, [1, 2, 3, 4], [2]
26    net_model = UNet(T, ch, ch_mult, attn, 2, 0.3).to(device)
27    net_model.load_state_dict(torch.load(model_path, map_location=device))
28    net_model.eval()
29
30    sampler = GaussianDiffusionSampler_cond(
31        net_model, beta_1=1e-4, beta_T=0.02, T=T
32    ).to(device)
33
34    # 準備 CBCT 條件
35    cbct_tensor = torch.from_numpy(cbct_slice).float().unsqueeze(0).unsqueeze(0).to(device)
36
37    # 多次取樣
38    all_samples = []
39    with torch.no_grad():
40        for i in range(n_samples):
41            noise = torch.randn(1, 1, 256, 256, device=device)
42            x_in = torch.cat((noise, cbct_tensor), dim=1)
43            x_out = sampler(x_in)
44            sct = x_out[:, 0:1, :, :].cpu().numpy().squeeze()
45            all_samples.append(sct)
46            print(f"  Sample {i+1}/{n_samples} done.")
47
48    all_samples = np.stack(all_samples, axis=0)  # (n_samples, H, W)
49
50    return {
51        "mean_sCT":    np.mean(all_samples, axis=0),   # 平均合成 CT
52        "std_map":     np.std(all_samples, axis=0),     # 不確定性地圖
53        "all_samples": all_samples,                      # 所有取樣結果
54    }
55
56# 使用範例
57# result = estimate_uncertainty("Checkpoints/trial_1/ckpt_1000_.pt", cbct_slice, n_samples=10)
58# high_uncertainty_mask = result["std_map"] > threshold  # 標記高不確定性區域

12017 Pix2Pix (paired GAN)
2 └─ 2018 CycleGAN (unpaired GAN)
3     └─ 2020 medSynthesis (medical GAN)
4         └─ 2020 DDPM (Ho et al.)
5             └─ 2023 conditional_DDPM (本專案：cDDPM for CBCT→CT) ◄── 你在這裡
6                 └─ 2024+ DenseDiffusion, SD-Messenger (進階擴散)
7                     └─ 2025+ 3D Volumetric Diffusion (下一個 blue ocean)

 1# 概念示範：將 cDDPM 邏輯遷移到 MONAI 框架
 2from monai.networks.nets import DiffusionModelUNet
 3from monai.networks.schedulers import DDPMScheduler
 4
 5# MONAI 風格的條件式 U-Net
 6model = DiffusionModelUNet(
 7    spatial_dims=2,
 8    in_channels=2,         # 與本專案相同：[noisy_CT, CBCT]
 9    out_channels=1,        # 預測噪聲
10    num_channels=(128, 256, 384, 512),  # 對應 ch_mult=[1,2,3,4]
11    attention_levels=(False, False, True, False),
12    num_res_blocks=2,
13)
14
15scheduler = DDPMScheduler(
16    num_train_timesteps=1000,
17    beta_start=1e-4,
18    beta_end=0.02,
19    schedule="linear",
20)

 1# 概念示範：2D → 3D 的核心修改點
 2# 1. Conv2d → Conv3d
 3self.head = nn.Conv3d(2, ch, kernel_size=3, stride=1, padding=1)
 4
 5# 2. 輸入維度：(B, 2, D, H, W) 取代 (B, 2, H, W)
 6x_0 = torch.cat((ct_volume, cbct_volume), dim=1)  # 5D tensor
 7
 8# 3. GroupNorm 不變（channel-wise 操作）
 9# 4. Attention 計算量急遽增加：D*H*W tokens
10#    建議改用 axial attention 或 windowed attention