如果你听说研究癌症最好的“数字试验场”可能不是肿瘤模型而是视网膜第一反应多半是意外。但把视网膜数字孪生作为癌症模型正是 Simon Walker-Samuel 教授团队近年提出的研究方向通过眼底成像重建视网膜血管网络在虚拟空间里模拟血管生成、血流与氧供再把这些计算模型反哺给肿瘤抗血管生成治疗研究。这个方向听起来很“跨界”实际上它切中的是肿瘤研究里一个长期存在的老大难问题——新生血管很难在活体里被持续、无创地观察。这件事的意义不在于“多了一个 3D 模型”而在于它把数字孪生的对象从厂房、城市、电网这类确定性的物理系统切换到了充满不确定性的生命系统。生命系统的数字孪生有一个工业数字孪生很少遇到的难题你无法为了建一棵血管树把组织切开来量。视网膜刚好提供了一个几乎不可替代的窗口——它是人体内少数能无创、高分辨率、反复观察微血管动态的位置之一。所以视网膜数字孪生的本质是把肿瘤血管生成的核心机制搬到一双可以持续观测的“虚拟眼睛”里去做可控实验。这篇文章会从四个角度展开为什么视网膜能成为癌症模型生物医学数字孪生和工业数字孪生有什么本质差异从眼底影像到血管网络建模的完整技术链路以及用 Python 开源工具链跑通一条最小示例需要哪些步骤、会踩哪些坑。如果你从事医学影像算法、计算建模或者正在了解数字孪生技术如何进入生命科学下面的内容会有一份比较直接的参考价值。1. 这篇文章真正要解决的问题数字孪生Digital Twin是近两年搜索热度非常高的技术词但大多数讨论仍然停留在智慧工厂、设备预测性维护、城市仿真等工业场景。切换到医学影像和癌症研究之后这个概念很容易被误读成两种东西一种认为是“把器官做成立体可视化”另一种认为是“高精度的器官仿真”。这两种理解都不完整尤其对视网膜数字孪生来说核心价值既不在可视化也不在于把仿真做到多精细而在于能否用可观测的患者数据反推出不可直接测量的生理参数并让模型随着随访动态更新。真正的问题意识在于肿瘤血管生成研究已经积累了大量的分子生物学知识但是把它们转化成一个可以预测个体治疗反应的定量模型非常困难。传统的研究路径通常是取组织样本在切片上数血管、测密度、看形态这是一种“终点式”的观察每一次得到的信息都来自一个时间点、一个部位无法回答血管网络如何动态演化。视网膜数字孪生要解决的正是这个痛点它允许研究者在活体上反复观测完整的微血管网络再在虚拟空间里重建、仿真和验证。按照下面的示例你能跑通一条从眼底图像到血管网络建模的最小链路并用公开数据集验证每一步的输出。这套流程虽然离真正的临床级数字孪生还有距离但它能帮助你理解整个方向的数据流、技术栈和核心难点。如果你正好在考虑转做 AI 医疗、计算病理或生物医学数字孪生这篇文章可以作为第一份技术地图。2. 数字孪生不是“3D 模型”生物医学数字孪生的特殊性数字孪生的经典定义是利用物理实体在真实环境中产生的数据在数字空间中建立一个多物理场、多尺度、动态更新的映射模型并通过实时数据同步和反馈控制来优化物理实体。这里的关键词不是“3D”而是“映射”和“双向数据流”。一个只有静态几何的模型顶多叫数字存档不叫数字孪生。生物医学数字孪生和工业数字孪生有一个非常明显的区别工业系统的运行数据很容易通过传感器持续获取设备状态、温度、振动、负载都是可观测的变量而人体内部的生理过程大部分观测手段都是有创的或者只能通过间接信号推断。以视网膜血管为例眼底彩照只能提供二维形态OCTA 可以提供血流信号但血流速度、血管壁切应力、血管内外氧浓度这些东西几乎都无法直接测量。于是生物医学数字孪生的核心任务变成了“参数反演”——用有限的可观测数据去反推那些不可直接测量的内部状态。工业数字孪生与生物医学数字孪生还有一个差异就是尺度复杂性。数字孪生需要考虑的尺度跨度很大宏观上是血管网络的几何形态介观上是单个血管段的血流量和血管壁力学微观上是内皮细胞和血管周围细胞的信号通路。没有一个单一模型能同时覆盖所有尺度所以实际工程中必须采用多尺度建模策略。你在文献里看到 Agent-based model、血流动力学模型、氧输运模型都是这个多尺度链条中的一环。可以这样理解工业数字孪生化的是一个“长寿命、可观测、强约束”的系统生物医学数字孪生化的是一个“个体差异大、观测稀疏、边界模糊”的系统。后者不是不能做而是不能照搬工业数字孪生的工程流程。它需要更强调个体化初始化、模型降阶、不确定性量化和临床验证闭环这也是为什么目前大多数生物数字孪生还停留在研究阶段而不是临床产品。3. 为什么偏偏是视网膜癌症血管生成的“活体窗口”要理解视网膜数字孪生为什么可以用来研究癌症得先理解实体肿瘤的一个基本生存法则肿瘤长到一定体积之后单纯依靠扩散已经无法获得足够的氧气和营养必须诱导新生血管进入肿瘤内部。这个过程就是血管生成angiogenesis。肿瘤细胞会分泌大量血管内皮生长因子VEGF激活周围原本静止的内皮细胞让它们出芽、增殖、迁移最终形成新的血管网络。问题在于肿瘤新生血管和正常血管不一样。它们往往是扭曲、囊状扩张、充满盲端和动静脉短路的网络血管壁不完整周细胞覆盖减少渗漏严重。这种异常结构导致血流混乱肿瘤内部出现大量缺氧区域而缺氧又会进一步促进 VEGF 分泌形成恶性循环。临床上很多抗肿瘤药物尤其是抗 VEGF 药物就是试图切断这个循环但疗效个体差异很大。如果我们能在治疗前就预测某个患者的肿瘤血管网络会如何响应这无疑会极大改善治疗决策。视网膜之所以成为研究这类问题的理想窗口原因有三个。第一视网膜的微血管网络在结构上非常规则正常情况下呈分层排列浅层和深层血管丛具有清晰的空间分布当发生新生血管时会出现异常出芽、渗漏、无灌注区等特征这些形态学异常与肿瘤血管异常高度相似。也就是说视网膜提供了一个“结构规则、异常可识别”的血管网络样本。第二视网膜血管和肿瘤血管共享大量分子通路。VEGF 信号通路、缺氧诱导因子通路、周细胞调控机制在视网膜新生血管疾病和肿瘤血管生成中都有重要角色。因此研究者可以在视网膜上验证关于血管生成的定量模型再把模型迁移到肿瘤场景中。第三视网膜是体内少数几个能无创、高分辨率、反复观察微血管的地方。眼底彩照、荧光素血管造影、OCT 血管成像都可以在活体上获取微血管形态和血流信息即便在模式动物上也可以用共聚焦显微镜等技术观察小鼠视网膜血管。这种“可重复观测”能力恰恰是建立和验证数字孪生模型最需要的条件。所以视网膜数字孪生的准确表述不是“把眼睛做成一模一样的模型”而是“把视网膜血管网络变成一个可观测、可建模、可验证的癌症血管生成研究平台”。这个平台的产物是虚拟血管树以及在这棵树上的血流和血管生成仿真。理解了这一点后面所有技术细节就都有了落点。4. 视网膜数字孪生的技术架构与建模流程视网膜数字孪生的完整流程可以拆成五层数据采集层、图像分割层、几何与拓扑建模层、物理仿真层、生物学仿真层和验证回调层。这里的分层不是严格的软件架构而更像一条研究流水线每一层都有自己的技术难点和坑。4.1 数据采集层数据来源可以是临床视网膜影像也可以是模式动物的视网膜显微图像。临床场景中最常用的是眼底彩照、OCTA、OCT 和荧光素血管造影。眼底彩照适合获取二维血管形态OCTA 不需要注射造影剂就能分层显示血流信号但对血管细节的分辨率受设备限制荧光素血管造影能观察渗漏和灌注缺损但属于有创检查。在建立数字孪生初期至少需要一套能稳定提供血管分割结果的影像序列。如果要做动态孪生还需要同一患者的多次随访影像时间间隔可能是数周或数月。4.2 图像分割层图像分割的目标是把血管从背景中分离出来。传统方法用 Frangi 滤波器这类管状结构增强算子再配合阈值和形态学后处理深度学习方法则普遍使用 U-Net 及其变体。分割质量直接决定后面所有步骤的可靠性。一个常见的工程问题是没有做好光源不均校正眼底图像四周过暗或中心过亮导致分割结果在边缘区域大量丢失细血管。解决方式是在分割前做绿色通道提取、CLAHE 对比度增强和光照归一化。4.3 几何与拓扑建模层分割得到的是二值掩膜不能直接用来做仿真。要把掩膜转成可供计算的血管树通常需要三个步骤骨架化、管径估计和图构建。骨架化把血管中心线提取为单像素宽的路径管径估计可以基于距离变换或局部血管宽度求解图构建则把血管网络抽象成由节点和边组成的拓扑图节点代表分支点或端点边代表血管段每条边附带长度和平均管径。血管树的拓扑结构对血流计算影响巨大分支角度、管径之比、连通性这些拓扑参数一旦出错后面仿真结果就没有意义。4.4 物理仿真层物理仿真层最常用的是血流动力学模型。对于视网膜这类微循环系统可以先用零维或一维模型做集中参数分析将每段血管视为圆管用泊肃叶定律计算阻力用基尔霍夫定律在网络中分配流量。更精细的方法是在三维血管几何上做计算流体力学但计算成本很高短期内不适合大规模个体化应用。在血流模型之外还需要考虑氧输运模型氧气在血管内随血流对流在血管壁附近扩散到组织中并被光感受器等耗氧结构消耗。这一步开始涉及大量生理参数很多参数无法直接测量只能依靠文献值或参数反演。4.5 生物学仿真层与验证回调层生物学仿真层的目标是模拟血管生成行为常用 Agent-based model。每个 agent 可以代表一个内皮细胞、周细胞或巨噬细胞它们根据周围环境的氧浓度、VEGF 浓度和机械信号改变行为增殖、迁移、凋亡或出芽。Agent-based model 的优势是能表现细胞级的随机性和涌现行为缺点是参数众多、计算量大而且很难直接与宏观影像数据对比。因此完整的流程还必须包含验证回调层把模型的输出与真实随访影像对比判断血管密度、分支模式、无灌注区位置是否一致如果不一致就需要调整模型参数回到前面各层重新计算。整个数字孪生本质上是一个“反演-预测-验证”的循环。5. 环境准备与公开数据集在进入代码之前先准备环境。下面的依赖版本不要求写死建议以你当前的 Python 版本和 pip 解析结果为准。本文演示所用到的核心库是 OpenCV、scikit-image、networkx、numpy 和 matplotlib。这些库在医学影像算法圈子里非常常用安装方式也最简单。python -m venv reti-twin-env source reti-twin-env/bin/activate pip install opencv-python scikit-image networkx numpy matplotlib如果网络环境允许建议同时安装 jupyter方便逐步查看中间结果。pip install jupyter公开数据集方面最经典的眼底血管分割数据集是 DRIVE 和 STARE。DRIVE 包含 40 张眼底彩照其中训练集和测试集各 20 张每张图都有人工标注的血管掩膜适合做分割模型的训练和验证。STARE 包含 20 张眼底图像和两位标注者的血管标注标注倾向于保留更多细血管。两个数据集的图像尺寸不同直接使用前需要统一预处理流程。这里需要注意下载这些数据集前应查看各自的许可协议研究用途通常免费但商业用途可能需要单独联系数据提供方。建议的项目目录结构如下retina_twin/ ├── reti_twin/ │ ├── preprocess.py │ ├── segment.py │ ├── graph_builder.py │ ├── hemodynamics.py │ └── abm_demo.py ├── data/ │ └── DRIVE/ │ ├── training/ │ │ ├── images/ │ │ └── manual_fov/ │ └── test/ └── run_pipeline.py这样一来每一步都有独立模块最后用run_pipeline.py将流程串联起来。这也是一个比较适合继续扩展的工程结构分割模型更换、图算法增加、血流模型升级都不需要重写整个项目。6. 最小可行链路从眼底图像到血管网络的代码实现这一部分的目标是把一张眼底彩照变成一株可以用图结构表示的血管树再对血管段做简化血流阻力计算。由于只是教学演示代码不会追求 SOTA 分割精度但会保证流程完整、可运行、可验证。6.1 读取与预处理眼底图像的绿色通道里血管和背景的对比度最高因此第一步通常是提取绿色通道再做 CLAHE 对比度增强。下面是retina_twin/preprocess.py的代码# 文件路径retina_twin/preprocess.py import cv2 import numpy as np def load_retina_image(path: str) - np.ndarray: 读取眼底图像。DRIVE 数据集图像为 .tif 格式。 img cv2.imread(path, cv2.IMREAD_COLOR) if img is None: raise FileNotFoundError(f无法读取图像: {path}) return img def extract_green_channel(img: np.ndarray) - np.ndarray: 眼底彩照的绿色通道对比度最高适合血管分割。 return img[:, :, 1] def enhance_vessels(green_channel: np.ndarray, clip_limit: float 2.0) - np.ndarray: 用 CLAHE 增强血管对比度参数需要根据图像质量调整。 clahe cv2.createCLAHE(clipLimitclip_limit, tileGridSize(8, 8)) return clahe.apply(green_channel)这段代码有三个关键点。第一OpenCV 读取图像时默认通道顺序是 BGR所以绿色通道的索引是 1这里是针对彩色眼底照片。第二CLAHE 的clip_limit如果设得过大会放大背景噪声让后续分割出现大量假阳性如果设得过小细血管可能不明显。第三tileGridSize(8, 8)对 565×584 的 DRIVE 图像只是经验值如果换成更大尺寸的图像应适当增加网格大小。6.2 血管分割与骨架化分割环节使用 Frangi 滤波器增强管状结构。Frangi 滤波器通过分析海森矩阵的特征值来突出“管状”区域非常适合血管和视网膜这类线性结构。然后用 Otsu 阈值产生二值掩膜再用remove_small_objects去掉小块噪声。最后做骨架化得到单像素宽度的血管中心线。# 文件路径retina_twin/segment.py from skimage.filters import frangi, threshold_otsu from skimage.morphology import remove_small_objects, skeletonize import numpy as np def segment_vessels(enhanced: np.ndarray) - np.ndarray: 用 Frangi 滤波器增强血管再用 Otsu 阈值生成血管掩膜。 # sigmasrange(1, 5) 表示检测宽度约为 1~5 像素的血管 vessel_response frangi(enhanced, sigmasrange(1, 5)) binary vessel_response threshold_otsu(vessel_response) # 去掉面积小于 64 像素的噪声区域 binary remove_small_objects(binary, min_size64) return binary def get_skeleton(binary: np.ndarray) - np.ndarray: 对二值血管掩膜骨架化得到单像素宽的血管中心线。 return skeletonize(binary)这里最容易踩的坑是sigmas的范围。如果只设range(1, 3)粗血管中心区域可能被误判为“非管状”导致血管内部出现空洞如果设得太大比如range(1, 10)一些弯曲的背景纹理也可能被增强成血管。比较稳妥的做法是先跑一组小范围参数用可视化确认效果后再批量处理。6.3 骨架像素转血管网络图骨架化之后每个像素并不是最终需要的血管段。我把骨架像素转换成networkx.Graph其中每个节点对应一个骨架像素再写一个压缩函数把连续路径节点压缩成带长度属性的边只保留端点和分支点。这样得到的图才是真正适合做血流网络分析的血管树拓扑。# 文件路径retina_twin/graph_builder.py import networkx as nx import numpy as np def skeleton_to_graph(skeleton: np.ndarray) - nx.Graph: 将骨架像素转换为全像素图。 说明全像素图节点数量很多后续血流分析通常不直接使用 但它是压缩为血管段图之前的必要中间表示。 g nx.Graph() ys, xs np.nonzero(skeleton) pixels set(zip(ys.tolist(), xs.tolist())) for y, x in pixels: neighbors [] for dy in (-1, 0, 1): for dx in (-1, 0, 1): if dy 0 and dx 0: continue if (y dy, x dx) in pixels: neighbors.append((y dy, x dx)) kind endpoint if len(neighbors) 2: kind path elif len(neighbors) 2: kind branch g.add_node((y, x), kindkind, neighborslen(neighbors)) for nb in neighbors: g.add_edge((y, x), nb) return g def compress_graph(g: nx.Graph) - nx.Graph: 将连续 path 节点压缩成带长度属性的边只保留端点和分支点。 这是教学演示版本真实项目中建议先去除短分支再做压缩。 h nx.Graph() branch_nodes [n for n, d in g.nodes(dataTrue) if d[kind] ! path] for n in branch_nodes: h.add_node(n, kindg.nodes[n][kind]) for b_node in branch_nodes: for neighbor in g.neighbors(b_node): prev b_node cur neighbor length 1 # 沿着 path 节点走到下一个分支点或端点 while g.nodes[cur][kind] path: candidates [n for n in g.neighbors(cur) if n ! prev] if not candidates: break prev, cur cur, candidates[0] length 1 if cur in h and cur ! b_node: h.add_edge(b_node, cur, lengthlength) return h这个压缩过程在血管树有环或者分割不完整时会出现奇怪结果。比如两个分支点之间如果存在两条平行路径compress_graph会为每条路径分别建立边其中第一条可以被networkx保留第二条可能覆盖属性。真实项目中更可靠的做法是先用连通域分析把每段血管分离再做段级别的拓扑连接而不是在全局像素图上做遍历。6.4 简化血流阻力计算有了血管段图和管径信息后就可以用泊肃叶定律做最基础的阻力计算。血流在圆管中的阻力公式是 R 8ηL / (πr⁴)其中 η 是血液黏度L 是血管长度r 是血管半径。注意管径的四次方在分母位置意味着半径只要变化 10%阻力就会变化约 46%。这正是血管收缩、扩张和新生血管对局部血流影响巨大的原因。# 文件路径retina_twin/hemodynamics.py import math def poiseuille_resistance( radius_m: float, length_m: float, viscosity_pa_s: float 0.0035, ) - float: 血管段阻力 R 8ηL / (πr^4)单位 Pa·s/m^3。 return (8.0 * viscosity_pa_s * length_m) / (math.pi * radius_m**4) def estimate_segment_flow(pressure_drop_pa: float, resistance: float) - float: 流量 Q ΔP / R单位 m^3/s。 return pressure_drop_p 视网膜数字孪生：癌症血管生成研究的新试验场
如果你听说研究癌症最好的“数字试验场”可能不是肿瘤模型而是视网膜第一反应多半是意外。但把视网膜数字孪生作为癌症模型正是 Simon Walker-Samuel 教授团队近年提出的研究方向通过眼底成像重建视网膜血管网络在虚拟空间里模拟血管生成、血流与氧供再把这些计算模型反哺给肿瘤抗血管生成治疗研究。这个方向听起来很“跨界”实际上它切中的是肿瘤研究里一个长期存在的老大难问题——新生血管很难在活体里被持续、无创地观察。这件事的意义不在于“多了一个 3D 模型”而在于它把数字孪生的对象从厂房、城市、电网这类确定性的物理系统切换到了充满不确定性的生命系统。生命系统的数字孪生有一个工业数字孪生很少遇到的难题你无法为了建一棵血管树把组织切开来量。视网膜刚好提供了一个几乎不可替代的窗口——它是人体内少数能无创、高分辨率、反复观察微血管动态的位置之一。所以视网膜数字孪生的本质是把肿瘤血管生成的核心机制搬到一双可以持续观测的“虚拟眼睛”里去做可控实验。这篇文章会从四个角度展开为什么视网膜能成为癌症模型生物医学数字孪生和工业数字孪生有什么本质差异从眼底影像到血管网络建模的完整技术链路以及用 Python 开源工具链跑通一条最小示例需要哪些步骤、会踩哪些坑。如果你从事医学影像算法、计算建模或者正在了解数字孪生技术如何进入生命科学下面的内容会有一份比较直接的参考价值。1. 这篇文章真正要解决的问题数字孪生Digital Twin是近两年搜索热度非常高的技术词但大多数讨论仍然停留在智慧工厂、设备预测性维护、城市仿真等工业场景。切换到医学影像和癌症研究之后这个概念很容易被误读成两种东西一种认为是“把器官做成立体可视化”另一种认为是“高精度的器官仿真”。这两种理解都不完整尤其对视网膜数字孪生来说核心价值既不在可视化也不在于把仿真做到多精细而在于能否用可观测的患者数据反推出不可直接测量的生理参数并让模型随着随访动态更新。真正的问题意识在于肿瘤血管生成研究已经积累了大量的分子生物学知识但是把它们转化成一个可以预测个体治疗反应的定量模型非常困难。传统的研究路径通常是取组织样本在切片上数血管、测密度、看形态这是一种“终点式”的观察每一次得到的信息都来自一个时间点、一个部位无法回答血管网络如何动态演化。视网膜数字孪生要解决的正是这个痛点它允许研究者在活体上反复观测完整的微血管网络再在虚拟空间里重建、仿真和验证。按照下面的示例你能跑通一条从眼底图像到血管网络建模的最小链路并用公开数据集验证每一步的输出。这套流程虽然离真正的临床级数字孪生还有距离但它能帮助你理解整个方向的数据流、技术栈和核心难点。如果你正好在考虑转做 AI 医疗、计算病理或生物医学数字孪生这篇文章可以作为第一份技术地图。2. 数字孪生不是“3D 模型”生物医学数字孪生的特殊性数字孪生的经典定义是利用物理实体在真实环境中产生的数据在数字空间中建立一个多物理场、多尺度、动态更新的映射模型并通过实时数据同步和反馈控制来优化物理实体。这里的关键词不是“3D”而是“映射”和“双向数据流”。一个只有静态几何的模型顶多叫数字存档不叫数字孪生。生物医学数字孪生和工业数字孪生有一个非常明显的区别工业系统的运行数据很容易通过传感器持续获取设备状态、温度、振动、负载都是可观测的变量而人体内部的生理过程大部分观测手段都是有创的或者只能通过间接信号推断。以视网膜血管为例眼底彩照只能提供二维形态OCTA 可以提供血流信号但血流速度、血管壁切应力、血管内外氧浓度这些东西几乎都无法直接测量。于是生物医学数字孪生的核心任务变成了“参数反演”——用有限的可观测数据去反推那些不可直接测量的内部状态。工业数字孪生与生物医学数字孪生还有一个差异就是尺度复杂性。数字孪生需要考虑的尺度跨度很大宏观上是血管网络的几何形态介观上是单个血管段的血流量和血管壁力学微观上是内皮细胞和血管周围细胞的信号通路。没有一个单一模型能同时覆盖所有尺度所以实际工程中必须采用多尺度建模策略。你在文献里看到 Agent-based model、血流动力学模型、氧输运模型都是这个多尺度链条中的一环。可以这样理解工业数字孪生化的是一个“长寿命、可观测、强约束”的系统生物医学数字孪生化的是一个“个体差异大、观测稀疏、边界模糊”的系统。后者不是不能做而是不能照搬工业数字孪生的工程流程。它需要更强调个体化初始化、模型降阶、不确定性量化和临床验证闭环这也是为什么目前大多数生物数字孪生还停留在研究阶段而不是临床产品。3. 为什么偏偏是视网膜癌症血管生成的“活体窗口”要理解视网膜数字孪生为什么可以用来研究癌症得先理解实体肿瘤的一个基本生存法则肿瘤长到一定体积之后单纯依靠扩散已经无法获得足够的氧气和营养必须诱导新生血管进入肿瘤内部。这个过程就是血管生成angiogenesis。肿瘤细胞会分泌大量血管内皮生长因子VEGF激活周围原本静止的内皮细胞让它们出芽、增殖、迁移最终形成新的血管网络。问题在于肿瘤新生血管和正常血管不一样。它们往往是扭曲、囊状扩张、充满盲端和动静脉短路的网络血管壁不完整周细胞覆盖减少渗漏严重。这种异常结构导致血流混乱肿瘤内部出现大量缺氧区域而缺氧又会进一步促进 VEGF 分泌形成恶性循环。临床上很多抗肿瘤药物尤其是抗 VEGF 药物就是试图切断这个循环但疗效个体差异很大。如果我们能在治疗前就预测某个患者的肿瘤血管网络会如何响应这无疑会极大改善治疗决策。视网膜之所以成为研究这类问题的理想窗口原因有三个。第一视网膜的微血管网络在结构上非常规则正常情况下呈分层排列浅层和深层血管丛具有清晰的空间分布当发生新生血管时会出现异常出芽、渗漏、无灌注区等特征这些形态学异常与肿瘤血管异常高度相似。也就是说视网膜提供了一个“结构规则、异常可识别”的血管网络样本。第二视网膜血管和肿瘤血管共享大量分子通路。VEGF 信号通路、缺氧诱导因子通路、周细胞调控机制在视网膜新生血管疾病和肿瘤血管生成中都有重要角色。因此研究者可以在视网膜上验证关于血管生成的定量模型再把模型迁移到肿瘤场景中。第三视网膜是体内少数几个能无创、高分辨率、反复观察微血管的地方。眼底彩照、荧光素血管造影、OCT 血管成像都可以在活体上获取微血管形态和血流信息即便在模式动物上也可以用共聚焦显微镜等技术观察小鼠视网膜血管。这种“可重复观测”能力恰恰是建立和验证数字孪生模型最需要的条件。所以视网膜数字孪生的准确表述不是“把眼睛做成一模一样的模型”而是“把视网膜血管网络变成一个可观测、可建模、可验证的癌症血管生成研究平台”。这个平台的产物是虚拟血管树以及在这棵树上的血流和血管生成仿真。理解了这一点后面所有技术细节就都有了落点。4. 视网膜数字孪生的技术架构与建模流程视网膜数字孪生的完整流程可以拆成五层数据采集层、图像分割层、几何与拓扑建模层、物理仿真层、生物学仿真层和验证回调层。这里的分层不是严格的软件架构而更像一条研究流水线每一层都有自己的技术难点和坑。4.1 数据采集层数据来源可以是临床视网膜影像也可以是模式动物的视网膜显微图像。临床场景中最常用的是眼底彩照、OCTA、OCT 和荧光素血管造影。眼底彩照适合获取二维血管形态OCTA 不需要注射造影剂就能分层显示血流信号但对血管细节的分辨率受设备限制荧光素血管造影能观察渗漏和灌注缺损但属于有创检查。在建立数字孪生初期至少需要一套能稳定提供血管分割结果的影像序列。如果要做动态孪生还需要同一患者的多次随访影像时间间隔可能是数周或数月。4.2 图像分割层图像分割的目标是把血管从背景中分离出来。传统方法用 Frangi 滤波器这类管状结构增强算子再配合阈值和形态学后处理深度学习方法则普遍使用 U-Net 及其变体。分割质量直接决定后面所有步骤的可靠性。一个常见的工程问题是没有做好光源不均校正眼底图像四周过暗或中心过亮导致分割结果在边缘区域大量丢失细血管。解决方式是在分割前做绿色通道提取、CLAHE 对比度增强和光照归一化。4.3 几何与拓扑建模层分割得到的是二值掩膜不能直接用来做仿真。要把掩膜转成可供计算的血管树通常需要三个步骤骨架化、管径估计和图构建。骨架化把血管中心线提取为单像素宽的路径管径估计可以基于距离变换或局部血管宽度求解图构建则把血管网络抽象成由节点和边组成的拓扑图节点代表分支点或端点边代表血管段每条边附带长度和平均管径。血管树的拓扑结构对血流计算影响巨大分支角度、管径之比、连通性这些拓扑参数一旦出错后面仿真结果就没有意义。4.4 物理仿真层物理仿真层最常用的是血流动力学模型。对于视网膜这类微循环系统可以先用零维或一维模型做集中参数分析将每段血管视为圆管用泊肃叶定律计算阻力用基尔霍夫定律在网络中分配流量。更精细的方法是在三维血管几何上做计算流体力学但计算成本很高短期内不适合大规模个体化应用。在血流模型之外还需要考虑氧输运模型氧气在血管内随血流对流在血管壁附近扩散到组织中并被光感受器等耗氧结构消耗。这一步开始涉及大量生理参数很多参数无法直接测量只能依靠文献值或参数反演。4.5 生物学仿真层与验证回调层生物学仿真层的目标是模拟血管生成行为常用 Agent-based model。每个 agent 可以代表一个内皮细胞、周细胞或巨噬细胞它们根据周围环境的氧浓度、VEGF 浓度和机械信号改变行为增殖、迁移、凋亡或出芽。Agent-based model 的优势是能表现细胞级的随机性和涌现行为缺点是参数众多、计算量大而且很难直接与宏观影像数据对比。因此完整的流程还必须包含验证回调层把模型的输出与真实随访影像对比判断血管密度、分支模式、无灌注区位置是否一致如果不一致就需要调整模型参数回到前面各层重新计算。整个数字孪生本质上是一个“反演-预测-验证”的循环。5. 环境准备与公开数据集在进入代码之前先准备环境。下面的依赖版本不要求写死建议以你当前的 Python 版本和 pip 解析结果为准。本文演示所用到的核心库是 OpenCV、scikit-image、networkx、numpy 和 matplotlib。这些库在医学影像算法圈子里非常常用安装方式也最简单。python -m venv reti-twin-env source reti-twin-env/bin/activate pip install opencv-python scikit-image networkx numpy matplotlib如果网络环境允许建议同时安装 jupyter方便逐步查看中间结果。pip install jupyter公开数据集方面最经典的眼底血管分割数据集是 DRIVE 和 STARE。DRIVE 包含 40 张眼底彩照其中训练集和测试集各 20 张每张图都有人工标注的血管掩膜适合做分割模型的训练和验证。STARE 包含 20 张眼底图像和两位标注者的血管标注标注倾向于保留更多细血管。两个数据集的图像尺寸不同直接使用前需要统一预处理流程。这里需要注意下载这些数据集前应查看各自的许可协议研究用途通常免费但商业用途可能需要单独联系数据提供方。建议的项目目录结构如下retina_twin/ ├── reti_twin/ │ ├── preprocess.py │ ├── segment.py │ ├── graph_builder.py │ ├── hemodynamics.py │ └── abm_demo.py ├── data/ │ └── DRIVE/ │ ├── training/ │ │ ├── images/ │ │ └── manual_fov/ │ └── test/ └── run_pipeline.py这样一来每一步都有独立模块最后用run_pipeline.py将流程串联起来。这也是一个比较适合继续扩展的工程结构分割模型更换、图算法增加、血流模型升级都不需要重写整个项目。6. 最小可行链路从眼底图像到血管网络的代码实现这一部分的目标是把一张眼底彩照变成一株可以用图结构表示的血管树再对血管段做简化血流阻力计算。由于只是教学演示代码不会追求 SOTA 分割精度但会保证流程完整、可运行、可验证。6.1 读取与预处理眼底图像的绿色通道里血管和背景的对比度最高因此第一步通常是提取绿色通道再做 CLAHE 对比度增强。下面是retina_twin/preprocess.py的代码# 文件路径retina_twin/preprocess.py import cv2 import numpy as np def load_retina_image(path: str) - np.ndarray: 读取眼底图像。DRIVE 数据集图像为 .tif 格式。 img cv2.imread(path, cv2.IMREAD_COLOR) if img is None: raise FileNotFoundError(f无法读取图像: {path}) return img def extract_green_channel(img: np.ndarray) - np.ndarray: 眼底彩照的绿色通道对比度最高适合血管分割。 return img[:, :, 1] def enhance_vessels(green_channel: np.ndarray, clip_limit: float 2.0) - np.ndarray: 用 CLAHE 增强血管对比度参数需要根据图像质量调整。 clahe cv2.createCLAHE(clipLimitclip_limit, tileGridSize(8, 8)) return clahe.apply(green_channel)这段代码有三个关键点。第一OpenCV 读取图像时默认通道顺序是 BGR所以绿色通道的索引是 1这里是针对彩色眼底照片。第二CLAHE 的clip_limit如果设得过大会放大背景噪声让后续分割出现大量假阳性如果设得过小细血管可能不明显。第三tileGridSize(8, 8)对 565×584 的 DRIVE 图像只是经验值如果换成更大尺寸的图像应适当增加网格大小。6.2 血管分割与骨架化分割环节使用 Frangi 滤波器增强管状结构。Frangi 滤波器通过分析海森矩阵的特征值来突出“管状”区域非常适合血管和视网膜这类线性结构。然后用 Otsu 阈值产生二值掩膜再用remove_small_objects去掉小块噪声。最后做骨架化得到单像素宽度的血管中心线。# 文件路径retina_twin/segment.py from skimage.filters import frangi, threshold_otsu from skimage.morphology import remove_small_objects, skeletonize import numpy as np def segment_vessels(enhanced: np.ndarray) - np.ndarray: 用 Frangi 滤波器增强血管再用 Otsu 阈值生成血管掩膜。 # sigmasrange(1, 5) 表示检测宽度约为 1~5 像素的血管 vessel_response frangi(enhanced, sigmasrange(1, 5)) binary vessel_response threshold_otsu(vessel_response) # 去掉面积小于 64 像素的噪声区域 binary remove_small_objects(binary, min_size64) return binary def get_skeleton(binary: np.ndarray) - np.ndarray: 对二值血管掩膜骨架化得到单像素宽的血管中心线。 return skeletonize(binary)这里最容易踩的坑是sigmas的范围。如果只设range(1, 3)粗血管中心区域可能被误判为“非管状”导致血管内部出现空洞如果设得太大比如range(1, 10)一些弯曲的背景纹理也可能被增强成血管。比较稳妥的做法是先跑一组小范围参数用可视化确认效果后再批量处理。6.3 骨架像素转血管网络图骨架化之后每个像素并不是最终需要的血管段。我把骨架像素转换成networkx.Graph其中每个节点对应一个骨架像素再写一个压缩函数把连续路径节点压缩成带长度属性的边只保留端点和分支点。这样得到的图才是真正适合做血流网络分析的血管树拓扑。# 文件路径retina_twin/graph_builder.py import networkx as nx import numpy as np def skeleton_to_graph(skeleton: np.ndarray) - nx.Graph: 将骨架像素转换为全像素图。 说明全像素图节点数量很多后续血流分析通常不直接使用 但它是压缩为血管段图之前的必要中间表示。 g nx.Graph() ys, xs np.nonzero(skeleton) pixels set(zip(ys.tolist(), xs.tolist())) for y, x in pixels: neighbors [] for dy in (-1, 0, 1): for dx in (-1, 0, 1): if dy 0 and dx 0: continue if (y dy, x dx) in pixels: neighbors.append((y dy, x dx)) kind endpoint if len(neighbors) 2: kind path elif len(neighbors) 2: kind branch g.add_node((y, x), kindkind, neighborslen(neighbors)) for nb in neighbors: g.add_edge((y, x), nb) return g def compress_graph(g: nx.Graph) - nx.Graph: 将连续 path 节点压缩成带长度属性的边只保留端点和分支点。 这是教学演示版本真实项目中建议先去除短分支再做压缩。 h nx.Graph() branch_nodes [n for n, d in g.nodes(dataTrue) if d[kind] ! path] for n in branch_nodes: h.add_node(n, kindg.nodes[n][kind]) for b_node in branch_nodes: for neighbor in g.neighbors(b_node): prev b_node cur neighbor length 1 # 沿着 path 节点走到下一个分支点或端点 while g.nodes[cur][kind] path: candidates [n for n in g.neighbors(cur) if n ! prev] if not candidates: break prev, cur cur, candidates[0] length 1 if cur in h and cur ! b_node: h.add_edge(b_node, cur, lengthlength) return h这个压缩过程在血管树有环或者分割不完整时会出现奇怪结果。比如两个分支点之间如果存在两条平行路径compress_graph会为每条路径分别建立边其中第一条可以被networkx保留第二条可能覆盖属性。真实项目中更可靠的做法是先用连通域分析把每段血管分离再做段级别的拓扑连接而不是在全局像素图上做遍历。6.4 简化血流阻力计算有了血管段图和管径信息后就可以用泊肃叶定律做最基础的阻力计算。血流在圆管中的阻力公式是 R 8ηL / (πr⁴)其中 η 是血液黏度L 是血管长度r 是血管半径。注意管径的四次方在分母位置意味着半径只要变化 10%阻力就会变化约 46%。这正是血管收缩、扩张和新生血管对局部血流影响巨大的原因。# 文件路径retina_twin/hemodynamics.py import math def poiseuille_resistance( radius_m: float, length_m: float, viscosity_pa_s: float 0.0035, ) - float: 血管段阻力 R 8ηL / (πr^4)单位 Pa·s/m^3。 return (8.0 * viscosity_pa_s * length_m) / (math.pi * radius_m**4) def estimate_segment_flow(pressure_drop_pa: float, resistance: float) - float: 流量 Q ΔP / R单位 m^3/s。 return pressure_drop_p