OpenMP 数据并行
TIFF 按 z-slice 或 strip 分批；Zarr 按 chunk 分批。每个线程拿到互不重叠的数据区域，并行读取与解压。
00 · EXECUTIVE SUMMARY
研究背景与核心问题
PetaKit5D 主要解决两个效率问题：一是超大图像的读取、写入和任务执行难以扩展；二是传统 Richardson-Lucy 反卷积需要较多迭代，计算量过大。
$$v_{\mathrm{processing}} \gtrsim v_{\mathrm{acquisition}}$$
要实现在线处理，端到端处理速度至少应接近数据采集速度。如果前者较低，未处理数据会持续积压，也无法在实验过程中及时发现对焦、视野或采集参数错误。
01 · PARALLEL I/O & PROCESSING
创新一：多层级并行 I/O 与分布式处理
PetaKit5D 的并行方案包含三个层级：文件内部的多线程读写、超大体数据的分块处理，以及集群层面的任务编排。这三个层级分别处理 I/O、内存容量和任务可靠性问题。
StorageTIFF / Zarr
→
Read多线程解压
→
ProcessCPU / GPU
→
Write压缩 + 落盘
Split–Process–Merge
单文件超出内存时，按空间子体积分割。卷积类操作增加 overlap/buffer，计算后只保留核心区再合并。
Conductor 编排
将函数、输入、输出和资源需求变成任务；持续监控 Slurm worker，跳过已有结果，失败后自动重交。
01.1 · FROM FILE TO CLUSTER
并行 I/O 的具体实现
对单个处理批次而言，端到端延迟可以近似分解为读取、计算、写入和调度四部分。因此，只提高 GPU 计算速度并不能消除 I/O 瓶颈。
$$T_{\mathrm{total}}=T_{\mathrm{read}}+T_{\mathrm{compute}}+T_{\mathrm{write}}+T_{\mathrm{schedule}}$$
CPP-TIFF READER
把 z-slice / strip 平均分给线程
代码先读取 OpenMP 最大线程数，再计算 batchSize；每个 worker 独立打开 TIFF handle，读取自己的 directory 与 encoded strip。大二维图则进一步按 strip 并行。
Thread 1 · slice 1–250
Thread 2 · slice 251–500
Thread 3 · slice 501–750
Thread 4 · slice 751–1000
CPP-TIFF WRITER
压缩并行，容器写入串行
TIFF 是单容器，多个线程同时改目录和 byte offset 会冲突。因此实现让 OpenMP 线程并行压缩 strips，同时由一个异步 writer 按顺序写入。它承认格式约束，并把最重的压缩部分并行化。
4 threads · parallel compression
1 writer · serialized container
Conductor
→
Worker 1
Task 1
Task 1
Worker 2
Task 2
Task 2
Worker 3
Task 3
Task 3
Worker n
Task n
Task n
状态回传
←
完成则记录输出；失败则在重试上限内重交，并为 CPU 任务逐次增加核数。
Zarr 与超大规模分块
TIFF
一个大容器
擅长兼容传统显微镜数据，但并行写入受容器目录结构限制。局部访问也容易牵涉更多连续数据。
ZARR
很多独立 chunk
一个空间块就是一个独立文件/对象。不同任务能只读取感兴趣区域，并把不同结果块并行写回，天然对应 split–process–merge。
| 问题 | 大 chunk | 小 chunk | 设计结论 |
|---|---|---|---|
| 固定 I/O / metadata 开销 | 摊薄，顺序吞吐高 | 重复支付较多 | 不要无限切小 |
| 并行任务数 | 少，粗粒度 | 多，负载均衡好 | 需足够多的 chunk 喂满 worker |
| 局部区域访问 | 可能读入无用数据 | 只读邻近块 | 根据算法访问模式 rechunk |
| 边界正确性 | 卷积/反卷积不能硬切：每个子体积需带 overlap | decon buffer ≥ PSF 半尺寸 + 额外余量（默认 10） | |
吞吐模型
worker 数量与 I/O 上限
8 workers
并行计算能力4.0 TB/h
存储 I/O 上限8.0 TB/h
实际端到端吞吐4.0 TB/h
此时计算端仍是瓶颈；增加 worker 能提升总吞吐。
parallelreadtiff.cpp:146–207OpenMP
int32_t numWorkers = omp_get_max_threads(); int32_t batchSize = (z - 1) / numWorkers + 1; #pragma omp parallel for for (w = 0; w < numWorkers; w++) { // 每个线程只读取自己的 z-slice 批次 TIFFReadEncodedStrip(tif, strip, destination, bytes); }
02 · BACKWARD PROJECTION
创新二：OMW 反向投影器
显微镜成像可以表示为下式，其中 \(x\) 为真实样本，\(f\) 为点扩散函数（PSF），\(g\) 为观测图像，\(n\) 为噪声。RL 方法通过正向成像、误差计算和反向投影逐步估计 \(x\)。
$$g=f\ast x+n$$
① 当前估计 \(x^{(k)}\)初始化真实样本估计
→
② 正向投影\(\hat g=f\ast x^{(k)}\)
→
③ 比值误差\(r=g/\max(\hat g,\varepsilon)\)
④ 对误差比值进行反向投影，并将校正结果乘到当前估计上
$$x^{(k+1)}=x^{(k)}\odot\left[b\ast\frac{g}{\max\!\left(f\ast x^{(k)},\varepsilon\right)}\right]$$
RATIO = 1
预测正确
观测与预测一致，该位置的校正因子接近 1，当前估计基本不变。
RATIO > 1
预测偏暗
观测包含更多信号，误差反投影后会增加可能贡献到这里的样本位置。
RATIO < 1
预测偏亮
当前估计贡献过强，反投影后的乘法修正会抑制对应位置。
方法的关键：OMW 保留 RL 的正向模型 \(f\)，主要改变反向投影器 \(b\)。传统 matched projector 使用与 PSF 对应的反向核；OMW 则把正则化逆滤波与真实 OTF 支持域结合，使单次迭代能提供更强且受控的校正。
02.1 · OTF MASKED WIENER
OMW 反向投影器的构造
它不是直接使用 1/H。因为当某个频率的 OTF 幅值 |H| 很小时，噪声会被无限放大。OMW 用 Wiener 项限制放大，用符合真实 OTF 形状的窗函数 W 决定哪些频率值得相信。
$$H=\mathcal{F}\{f\},\qquad F_{\mathrm{Wiener}}=\frac{H^{*}}{|H|^{2}+\alpha},\qquad B=W\odot F_{\mathrm{Wiener}},\qquad b=\mathcal{F}^{-1}\{B\}$$
PSF 变换到频域
计算 \(H=\mathcal F\{f\}\)，用 \(|H|\) 描述显微镜对各空间频率的传递能力。
累计幅值阈值分割
将 \(|H|\) 从大到小排序，默认取累计和达到 90% 的位置作为阈值；代码另设 \(10^{-3}\) 的最低阈值以应对噪声。
清理并保留主支持域
形态学开闭运算去噪；常规空间保留中央主对象，skewed space 则处理三个主要分量。
Convex hull 补齐支持域
逐切面做凸包填充，避免噪声孔洞和凹缺让高频信息被不规则切断。
构造相对距离 \(D\)
从频谱中心沿射线归一化：中心为 0，真实支持域边界为 1，不再假设所有 PSF 都是椭球。
Hann 渐消窗 \(W\)
默认 \(D<0.8\) 全通过，\(0.8\leq D<1\) 平滑衰减，\(D\geq1\) 为零，以减少硬截断引起的振铃。
与 Wiener 滤波器相乘
计算 \(B=W\odot H^*/(|H|^2+\alpha)\)，再经逆 Fourier 变换得到迭代使用的 \(b\)。
$$m=\min\left\{j:\frac{\sum_{i=1}^{j}a_{(i)}}{\sum_{i=1}^{N}a_{(i)}}\geq0.9\right\},\qquad \tau=a_{(m)}$$
$$W(D)=\begin{cases}1,&D<l,\\[2pt]\cos^{2}\!\left(\dfrac{\pi(D-l)}{2(u-l)}\right),&l\leq D<u,\\[6pt]0,&D\geq u,\end{cases}\qquad l=0.8,\;u=1$$
WB：固定椭球假设边缘真实频率可能被截掉
→
OMW：跟随真实 support内部保留，边缘平滑衰减
参数示意
Wiener 参数 α 的作用
α = 0.005
|H|直接逆 1/|H|（噪声放大）Wiener 增益 |H|/(|H|²+α)
omw_backprojector_generation.m:158–163构造 b
OTF_bp_w = conj(OTF) ./ (abs_OTF.^2 + alpha); OTF_bp_omw = OTF_bp_om .* OTF_bp_w; b_omw = fftshift(real(ifftn(OTF_bp_omw)));
decon_lucy_omw_function.m:107–111一次 RL 更新
CX = real(ifftn(fftn(J) .* OTF_f)); CX = max(CX, eps); J = real(ifftn(fftn(I ./ CX) .* OTF_b)) .* J; J = max(J, 0);
| 方法 | 反向投影器 | 优点 | 主要问题 / 改进 |
|---|---|---|---|
| 传统 RL | matched projector | 准确、稳健、通用 | 通常要 10–100 次迭代，PB 数据计算代价高 |
| WB | Wiener–Butterworth | 少量迭代即可收敛 | 椭球支持域可能截断 lattice light-sheet 的近矩形高频区域 |
| OMW | OTF mask × Wiener | 1–2 次常可达到全分辨率；适配任意 PSF 支持域 | α 与迭代数仍需数据驱动选择；论文用 FSC 辅助确定 |
03 · OUR METHOD
从原论文 OMW，到我们的 Poisson 反传优化
下面才进入我们的工作。顺序是：先识别 RL 的精确反向传播，再分析 OMW 改写 backward projector 的收益与代价；随后用 Newton/Fisher 看清病态曲率，最后构造保留 \(H^{\mathsf T}\) 的频谱预条件与正则化更新。
1. RL：精确 Poisson 反传
输出层的 Poisson error signal 与传回图像域的梯度为：
$$\delta_k=\mathbf 1-\frac{y}{Hx_k},\qquad g_k=H^{\mathsf T}\delta_k$$
归一化 PSF 下，RL 可写成 \(x_{k+1}=x_k-\operatorname{Diag}(x_k)g_k\)。因此它是全批量 scaled gradient，而不是随机采样意义上的 SGD。
2. OMW：改写 backward
OMW 保持 forward 为 \(H\)，但用 inverse-like 的 \(B^{\mathsf T}\) 代替链式法则决定的 \(H^{\mathsf T}\)：
$$B(\omega)=W(\omega)\frac{H^*(\omega)}{|H(\omega)|^2+\alpha}$$
弱传递频率得到更强校正，因此前几轮很快；代价是 \(B^{\mathsf T}\delta\) 一般不再等于原 Poisson 目标的梯度，继续迭代也更容易放大噪声。
3. Newton/Fisher：解释加速，也暴露失稳
$$\nabla_x^2\mathcal L=H^{\mathsf T}\operatorname{Diag}\!\left(\frac{y}{(Hx)^2}\right)H,\qquad F=H^{\mathsf T}\operatorname{Diag}\!\left(\frac1{Hx}\right)H$$
完整 Newton 要反复求解空间变权的病态线性系统；更重要的是，小 Fisher 特征值正对应弱 OTF 与噪声主导方向，直接求逆会给予这些方向最大增益。
4. 我们的设计：反传后处理优化几何
不替换 \(H^{\mathsf T}\)。在正值 log-intensity 层继续反传，并在 OTF 支持域中做分数 Fisher 白化：
$$g_z=\frac{x}{\bar y}\odot H^{\mathsf T}\!\left(\mathbf 1-\frac{y}{Hx}\right),\qquad P_\beta=W\bigl(|H|^2+\lambda\bigr)^{-\beta}$$
\(\beta=0,\tfrac14,\tfrac12,1\) 分别连接 log-SGD、半白化、全白化和近似逆 Fisher；这使候选方法来自一条可解释的连续谱，而不是任意替换反向算子。
快速版本：连续逆曲率 + 显式曲率先验
SCR-RL 在早期使用较强的 OTF 对角逆曲率，随后几何退火回 RL 尺度；每轮再施加二阶 Sobolev 近端滤波，抑制不可辨识高频的 Poisson 噪声。
$$P_k=(1-\beta_k)I+\beta_kP_{\rm inv},\quad \beta_k=\beta_1\rho^{k-1},\qquad x_{k+1}=\operatorname{prox}_{\lambda\|\Delta x\|_2^2}\!\left[x_k-x_k\odot P_kg_k\right]$$
保留真实 error path所有新方向都先计算 Poisson score，再通过真实 \(H^{\mathsf T}\) 返回图像域。
用曲率决定频率步长OTF/Fisher 近似只改变梯度之后的优化度量，不改变 forward model。
用退火控制半收敛早期快速恢复可靠频率，后期减弱逆滤波增益，并显式压制噪声主导高频。
EQUAL-ITERATION VALIDATION · OMW BEST ITERATION
统一固定在第 4 次迭代比较恢复效果
验证采用 12 个独立 Poisson seeds。所有方法在每个 seed 内使用相同观测、相同初始化并统一比较第 4 轮；数值与置信区间由全部 seeds 计算，恢复图展示多种子结果中具有代表性的 seed 31，并对所有方法使用相同显示范围。
k = 4 for every method · n = 12 held-out seedsSCR-RL · 第 4 轮 NRMSE0.4802 ± 0.0071
PetaKit5D OMW · 第 4 轮0.6195 ± 0.0086
Matched RL · 第 4 轮0.6392 ± 0.0029
Half whitening · 第 4 轮0.6264 ± 0.0020
| 方法 | 统一轮次 | Held-out NRMSE ↓ | PSNR ↑ | 第 4 轮观察 |
|---|---|---|---|---|
| Matched RL | 4 | 0.6392 ± 0.0029 | 33.345 dB | 稳定，但细丝仍明显模糊 |
| PetaKit5D OMW | 4 | 0.6195 ± 0.0086 | 33.618 dB | 细丝分离更快，同时出现颗粒噪声 |
| Damped Newton | 4 | 0.6920 ± 0.0025 | 32.655 dB | NLL 下降未转化为恢复误差改善 |
| Log-SGD | 4 | 0.6312 ± 0.0156 | 33.460 dB | 正值参数化本身不足以解决病态频率 |
| Full Fisher whitening | 4 | 0.6577 ± 0.0108 | 33.099 dB | 保守线搜索使早期校正仍慢 |
| Geometric half whitening | 4 | 0.6264 ± 0.0020 | 33.519 dB | 早期优于 RL，但未超过 OMW |
| SCR-RL | 4 | 0.4802 ± 0.0071 | 35.830 dB | 逆曲率负责快速分离，先验抑制颗粒噪声 |
这张图回答的是“相同外层迭代预算下谁恢复得更好”。第 4 轮时，SCR-RL 相对 RL 和 OMW 的配对 NRMSE 分别改善 24.86% 与 22.48%，两项 Holm-adjusted \(p=0.00293\)。方法各自最佳轮次的冻结终点属于另一个问题，不再用于恢复图的横向视觉比较。