CN115032196B - 一种全划片高通量彩色病理成像分析仪器及方法 - Google Patents
一种全划片高通量彩色病理成像分析仪器及方法 Download PDFInfo
- Publication number
- CN115032196B CN115032196B CN202210961097.7A CN202210961097A CN115032196B CN 115032196 B CN115032196 B CN 115032196B CN 202210961097 A CN202210961097 A CN 202210961097A CN 115032196 B CN115032196 B CN 115032196B
- Authority
- CN
- China
- Prior art keywords
- image
- color
- light source
- brightness
- full
- Prior art date
- Legal status (The legal status is an assumption and is not a legal conclusion. Google has not performed a legal analysis and makes no representation as to the accuracy of the status listed.)
- Active
Links
- 238000003384 imaging method Methods 0.000 title claims abstract description 199
- 238000004458 analytical method Methods 0.000 title claims abstract description 89
- 230000001575 pathological effect Effects 0.000 title claims abstract description 33
- 238000000034 method Methods 0.000 title claims description 145
- 230000007170 pathology Effects 0.000 claims abstract description 65
- 230000006870 function Effects 0.000 claims description 139
- 238000005286 illumination Methods 0.000 claims description 71
- 238000004422 calculation algorithm Methods 0.000 claims description 68
- 210000001747 pupil Anatomy 0.000 claims description 57
- 238000010586 diagram Methods 0.000 claims description 48
- 239000011159 matrix material Substances 0.000 claims description 47
- 230000003287 optical effect Effects 0.000 claims description 36
- 230000004075 alteration Effects 0.000 claims description 30
- 238000011156 evaluation Methods 0.000 claims description 29
- 238000013519 translation Methods 0.000 claims description 27
- 238000009826 distribution Methods 0.000 claims description 26
- 238000012546 transfer Methods 0.000 claims description 25
- 238000012937 correction Methods 0.000 claims description 24
- 238000010521 absorption reaction Methods 0.000 claims description 23
- 238000001228 spectrum Methods 0.000 claims description 22
- 230000003044 adaptive effect Effects 0.000 claims description 20
- 230000010287 polarization Effects 0.000 claims description 20
- 239000003086 colorant Substances 0.000 claims description 18
- 238000004364 calculation method Methods 0.000 claims description 17
- 238000006243 chemical reaction Methods 0.000 claims description 17
- 238000011084 recovery Methods 0.000 claims description 15
- 238000005070 sampling Methods 0.000 claims description 15
- 238000006073 displacement reaction Methods 0.000 claims description 12
- 238000005259 measurement Methods 0.000 claims description 12
- 238000005457 optimization Methods 0.000 claims description 12
- 238000012634 optical imaging Methods 0.000 claims description 11
- 230000009977 dual effect Effects 0.000 claims description 10
- 230000003595 spectral effect Effects 0.000 claims description 10
- 238000010186 staining Methods 0.000 claims description 10
- 238000013507 mapping Methods 0.000 claims description 9
- 230000002829 reductive effect Effects 0.000 claims description 9
- 230000004907 flux Effects 0.000 claims description 8
- 238000012545 processing Methods 0.000 claims description 7
- 238000011049 filling Methods 0.000 claims description 6
- 238000003860 storage Methods 0.000 claims description 5
- 239000000126 substance Substances 0.000 claims description 5
- 230000002194 synthesizing effect Effects 0.000 claims description 5
- 229910052736 halogen Inorganic materials 0.000 claims description 4
- 150000002367 halogens Chemical class 0.000 claims description 4
- 230000000630 rising effect Effects 0.000 claims description 4
- 238000002922 simulated annealing Methods 0.000 claims description 4
- 238000013473 artificial intelligence Methods 0.000 claims description 3
- 230000015572 biosynthetic process Effects 0.000 claims description 3
- 230000006835 compression Effects 0.000 claims description 3
- 238000013500 data storage Methods 0.000 claims description 3
- 230000003993 interaction Effects 0.000 claims description 3
- 238000010606 normalization Methods 0.000 claims description 3
- 238000000926 separation method Methods 0.000 claims description 3
- 238000010200 validation analysis Methods 0.000 claims description 3
- 238000013144 data compression Methods 0.000 claims description 2
- 238000013508 migration Methods 0.000 claims description 2
- 238000003786 synthesis reaction Methods 0.000 claims description 2
- 238000010295 mobile communication Methods 0.000 claims 1
- 230000007547 defect Effects 0.000 abstract description 9
- 238000010827 pathological analysis Methods 0.000 abstract description 5
- 238000012106 screening analysis Methods 0.000 abstract description 2
- 239000000523 sample Substances 0.000 description 94
- 230000008569 process Effects 0.000 description 28
- 230000000875 corresponding effect Effects 0.000 description 26
- 230000008859 change Effects 0.000 description 14
- 230000035945 sensitivity Effects 0.000 description 13
- 230000011218 segmentation Effects 0.000 description 12
- 238000010187 selection method Methods 0.000 description 11
- 239000000306 component Substances 0.000 description 10
- 238000005516 engineering process Methods 0.000 description 10
- 238000000799 fluorescence microscopy Methods 0.000 description 10
- 238000001514 detection method Methods 0.000 description 9
- 238000012216 screening Methods 0.000 description 9
- 241000282414 Homo sapiens Species 0.000 description 8
- 230000003902 lesion Effects 0.000 description 8
- 230000009466 transformation Effects 0.000 description 8
- 206010028980 Neoplasm Diseases 0.000 description 7
- 210000004027 cell Anatomy 0.000 description 7
- 230000009194 climbing Effects 0.000 description 7
- 238000011160 research Methods 0.000 description 7
- 201000011510 cancer Diseases 0.000 description 6
- 238000013527 convolutional neural network Methods 0.000 description 6
- 238000013461 design Methods 0.000 description 5
- 238000003745 diagnosis Methods 0.000 description 5
- 238000001914 filtration Methods 0.000 description 5
- 238000012804 iterative process Methods 0.000 description 5
- 231100000915 pathological change Toxicity 0.000 description 5
- 230000036285 pathological change Effects 0.000 description 5
- 238000012549 training Methods 0.000 description 5
- 230000005540 biological transmission Effects 0.000 description 4
- 239000003795 chemical substances by application Substances 0.000 description 4
- 238000013135 deep learning Methods 0.000 description 4
- 230000000694 effects Effects 0.000 description 4
- 239000010408 film Substances 0.000 description 4
- 238000011176 pooling Methods 0.000 description 4
- 230000001629 suppression Effects 0.000 description 4
- 230000009897 systematic effect Effects 0.000 description 4
- 230000000007 visual effect Effects 0.000 description 4
- OAICVXFJPJFONN-UHFFFAOYSA-N Phosphorus Chemical compound [P] OAICVXFJPJFONN-UHFFFAOYSA-N 0.000 description 3
- 238000013459 approach Methods 0.000 description 3
- 150000001875 compounds Chemical class 0.000 description 3
- 230000001276 controlling effect Effects 0.000 description 3
- 238000005520 cutting process Methods 0.000 description 3
- 238000011161 development Methods 0.000 description 3
- 230000018109 developmental process Effects 0.000 description 3
- 201000010099 disease Diseases 0.000 description 3
- 208000037265 diseases, disorders, signs and symptoms Diseases 0.000 description 3
- 230000005484 gravity Effects 0.000 description 3
- 230000006872 improvement Effects 0.000 description 3
- 238000000386 microscopy Methods 0.000 description 3
- 239000000203 mixture Substances 0.000 description 3
- 230000036961 partial effect Effects 0.000 description 3
- 238000010845 search algorithm Methods 0.000 description 3
- 238000012360 testing method Methods 0.000 description 3
- 238000012935 Averaging Methods 0.000 description 2
- 206010008342 Cervix carcinoma Diseases 0.000 description 2
- 208000006105 Uterine Cervical Neoplasms Diseases 0.000 description 2
- 238000013528 artificial neural network Methods 0.000 description 2
- 230000002146 bilateral effect Effects 0.000 description 2
- 239000012472 biological sample Substances 0.000 description 2
- 201000010881 cervical cancer Diseases 0.000 description 2
- 239000008358 core component Substances 0.000 description 2
- 230000002596 correlated effect Effects 0.000 description 2
- 230000003247 decreasing effect Effects 0.000 description 2
- 239000000428 dust Substances 0.000 description 2
- 230000009843 endothelial lesion Effects 0.000 description 2
- 238000002474 experimental method Methods 0.000 description 2
- 230000036541 health Effects 0.000 description 2
- 238000012417 linear regression Methods 0.000 description 2
- 238000012423 maintenance Methods 0.000 description 2
- 230000036210 malignancy Effects 0.000 description 2
- 238000004519 manufacturing process Methods 0.000 description 2
- 230000007935 neutral effect Effects 0.000 description 2
- 210000003463 organelle Anatomy 0.000 description 2
- 239000002245 particle Substances 0.000 description 2
- 230000002265 prevention Effects 0.000 description 2
- 238000011158 quantitative evaluation Methods 0.000 description 2
- 230000005855 radiation Effects 0.000 description 2
- 238000009877 rendering Methods 0.000 description 2
- 230000004044 response Effects 0.000 description 2
- 230000000717 retained effect Effects 0.000 description 2
- 239000000758 substrate Substances 0.000 description 2
- 210000001519 tissue Anatomy 0.000 description 2
- 241000234282 Allium Species 0.000 description 1
- 235000002732 Allium cepa var. cepa Nutrition 0.000 description 1
- 101100001671 Emericella variicolor andF gene Proteins 0.000 description 1
- 206010014561 Emphysema Diseases 0.000 description 1
- 241001289753 Graphium sarpedon Species 0.000 description 1
- 101100533306 Mus musculus Setx gene Proteins 0.000 description 1
- 241000276498 Pollachius virens Species 0.000 description 1
- 230000002411 adverse Effects 0.000 description 1
- 101150117391 artQ gene Proteins 0.000 description 1
- 230000001174 ascending effect Effects 0.000 description 1
- 230000002238 attenuated effect Effects 0.000 description 1
- 238000005311 autocorrelation function Methods 0.000 description 1
- 230000009286 beneficial effect Effects 0.000 description 1
- 230000008901 benefit Effects 0.000 description 1
- 238000001574 biopsy Methods 0.000 description 1
- 230000000903 blocking effect Effects 0.000 description 1
- 230000015556 catabolic process Effects 0.000 description 1
- 210000003855 cell nucleus Anatomy 0.000 description 1
- 210000002421 cell wall Anatomy 0.000 description 1
- 230000001427 coherent effect Effects 0.000 description 1
- 230000004456 color vision Effects 0.000 description 1
- 239000002131 composite material Substances 0.000 description 1
- 238000007906 compression Methods 0.000 description 1
- 230000008094 contradictory effect Effects 0.000 description 1
- 238000001816 cooling Methods 0.000 description 1
- 238000005314 correlation function Methods 0.000 description 1
- 238000007405 data analysis Methods 0.000 description 1
- 238000006731 degradation reaction Methods 0.000 description 1
- 230000001419 dependent effect Effects 0.000 description 1
- 238000009792 diffusion process Methods 0.000 description 1
- 239000006185 dispersion Substances 0.000 description 1
- 210000001339 epidermal cell Anatomy 0.000 description 1
- 238000000605 extraction Methods 0.000 description 1
- 230000004438 eyesight Effects 0.000 description 1
- 230000002349 favourable effect Effects 0.000 description 1
- 238000002073 fluorescence micrograph Methods 0.000 description 1
- 238000010191 image analysis Methods 0.000 description 1
- 238000003365 immunocytochemistry Methods 0.000 description 1
- 230000002055 immunohistochemical effect Effects 0.000 description 1
- 238000003364 immunohistochemistry Methods 0.000 description 1
- 238000007689 inspection Methods 0.000 description 1
- 230000001678 irradiating effect Effects 0.000 description 1
- 239000004973 liquid crystal related substance Substances 0.000 description 1
- 238000011068 loading method Methods 0.000 description 1
- 238000007620 mathematical function Methods 0.000 description 1
- 239000002184 metal Substances 0.000 description 1
- 230000000877 morphologic effect Effects 0.000 description 1
- 239000000843 powder Substances 0.000 description 1
- 238000007781 pre-processing Methods 0.000 description 1
- 238000002360 preparation method Methods 0.000 description 1
- 238000000513 principal component analysis Methods 0.000 description 1
- 238000012887 quadratic function Methods 0.000 description 1
- 238000011002 quantification Methods 0.000 description 1
- 238000013139 quantization Methods 0.000 description 1
- 238000005295 random walk Methods 0.000 description 1
- 230000009467 reduction Effects 0.000 description 1
- 238000012827 research and development Methods 0.000 description 1
- 230000000452 restraining effect Effects 0.000 description 1
- 238000012552 review Methods 0.000 description 1
- 239000004576 sand Substances 0.000 description 1
- 230000007480 spreading Effects 0.000 description 1
- 238000003892 spreading Methods 0.000 description 1
- 206010041823 squamous cell carcinoma Diseases 0.000 description 1
- 230000003068 static effect Effects 0.000 description 1
- 238000000547 structure data Methods 0.000 description 1
- 239000010409 thin film Substances 0.000 description 1
- 230000007704 transition Effects 0.000 description 1
- 238000012800 visualization Methods 0.000 description 1
Images
Classifications
-
- G—PHYSICS
- G01—MEASURING; TESTING
- G01N—INVESTIGATING OR ANALYSING MATERIALS BY DETERMINING THEIR CHEMICAL OR PHYSICAL PROPERTIES
- G01N21/00—Investigating or analysing materials by the use of optical means, i.e. using sub-millimetre waves, infrared, visible or ultraviolet light
- G01N21/84—Systems specially adapted for particular applications
-
- G—PHYSICS
- G02—OPTICS
- G02B—OPTICAL ELEMENTS, SYSTEMS OR APPARATUS
- G02B21/00—Microscopes
- G02B21/06—Means for illuminating specimens
-
- G—PHYSICS
- G02—OPTICS
- G02B—OPTICAL ELEMENTS, SYSTEMS OR APPARATUS
- G02B21/00—Microscopes
- G02B21/36—Microscopes arranged for photographic purposes or projection purposes or digital imaging or video purposes including associated control and data processing arrangements
- G02B21/365—Control or image processing arrangements for digital or video microscopes
- G02B21/367—Control or image processing arrangements for digital or video microscopes providing an output produced by processing a plurality of individual source images, e.g. image tiling, montage, composite images, depth sectioning, image comparison
-
- G—PHYSICS
- G06—COMPUTING; CALCULATING OR COUNTING
- G06F—ELECTRIC DIGITAL DATA PROCESSING
- G06F17/00—Digital computing or data processing equipment or methods, specially adapted for specific functions
- G06F17/10—Complex mathematical operations
-
- G—PHYSICS
- G06—COMPUTING; CALCULATING OR COUNTING
- G06T—IMAGE DATA PROCESSING OR GENERATION, IN GENERAL
- G06T3/00—Geometric image transformations in the plane of the image
- G06T3/40—Scaling of whole images or parts thereof, e.g. expanding or contracting
- G06T3/4046—Scaling of whole images or parts thereof, e.g. expanding or contracting using neural networks
-
- G—PHYSICS
- G06—COMPUTING; CALCULATING OR COUNTING
- G06T—IMAGE DATA PROCESSING OR GENERATION, IN GENERAL
- G06T3/00—Geometric image transformations in the plane of the image
- G06T3/40—Scaling of whole images or parts thereof, e.g. expanding or contracting
- G06T3/4053—Scaling of whole images or parts thereof, e.g. expanding or contracting based on super-resolution, i.e. the output image resolution being higher than the sensor resolution
- G06T3/4076—Scaling of whole images or parts thereof, e.g. expanding or contracting based on super-resolution, i.e. the output image resolution being higher than the sensor resolution using the original low-resolution images to iteratively correct the high-resolution images
-
- G—PHYSICS
- G06—COMPUTING; CALCULATING OR COUNTING
- G06T—IMAGE DATA PROCESSING OR GENERATION, IN GENERAL
- G06T7/00—Image analysis
- G06T7/0002—Inspection of images, e.g. flaw detection
- G06T7/0012—Biomedical image inspection
-
- G—PHYSICS
- G06—COMPUTING; CALCULATING OR COUNTING
- G06T—IMAGE DATA PROCESSING OR GENERATION, IN GENERAL
- G06T7/00—Image analysis
- G06T7/60—Analysis of geometric attributes
- G06T7/66—Analysis of geometric attributes of image moments or centre of gravity
-
- G—PHYSICS
- G06—COMPUTING; CALCULATING OR COUNTING
- G06T—IMAGE DATA PROCESSING OR GENERATION, IN GENERAL
- G06T7/00—Image analysis
- G06T7/90—Determination of colour characteristics
-
- G—PHYSICS
- G06—COMPUTING; CALCULATING OR COUNTING
- G06V—IMAGE OR VIDEO RECOGNITION OR UNDERSTANDING
- G06V10/00—Arrangements for image or video recognition or understanding
- G06V10/20—Image preprocessing
- G06V10/25—Determination of region of interest [ROI] or a volume of interest [VOI]
-
- G—PHYSICS
- G06—COMPUTING; CALCULATING OR COUNTING
- G06V—IMAGE OR VIDEO RECOGNITION OR UNDERSTANDING
- G06V10/00—Arrangements for image or video recognition or understanding
- G06V10/40—Extraction of image or video features
- G06V10/46—Descriptors for shape, contour or point-related descriptors, e.g. scale invariant feature transform [SIFT] or bags of words [BoW]; Salient regional features
- G06V10/467—Encoded features or binary features, e.g. local binary patterns [LBP]
-
- G—PHYSICS
- G06—COMPUTING; CALCULATING OR COUNTING
- G06V—IMAGE OR VIDEO RECOGNITION OR UNDERSTANDING
- G06V10/00—Arrangements for image or video recognition or understanding
- G06V10/70—Arrangements for image or video recognition or understanding using pattern recognition or machine learning
- G06V10/764—Arrangements for image or video recognition or understanding using pattern recognition or machine learning using classification, e.g. of video objects
-
- G—PHYSICS
- G06—COMPUTING; CALCULATING OR COUNTING
- G06V—IMAGE OR VIDEO RECOGNITION OR UNDERSTANDING
- G06V10/00—Arrangements for image or video recognition or understanding
- G06V10/70—Arrangements for image or video recognition or understanding using pattern recognition or machine learning
- G06V10/82—Arrangements for image or video recognition or understanding using pattern recognition or machine learning using neural networks
-
- G—PHYSICS
- G06—COMPUTING; CALCULATING OR COUNTING
- G06T—IMAGE DATA PROCESSING OR GENERATION, IN GENERAL
- G06T2207/00—Indexing scheme for image analysis or image enhancement
- G06T2207/30—Subject of image; Context of image processing
- G06T2207/30168—Image quality inspection
-
- G—PHYSICS
- G06—COMPUTING; CALCULATING OR COUNTING
- G06V—IMAGE OR VIDEO RECOGNITION OR UNDERSTANDING
- G06V2201/00—Indexing scheme relating to image or video recognition or understanding
- G06V2201/03—Recognition of patterns in medical or anatomical images
Landscapes
- Engineering & Computer Science (AREA)
- Physics & Mathematics (AREA)
- General Physics & Mathematics (AREA)
- Theoretical Computer Science (AREA)
- Computer Vision & Pattern Recognition (AREA)
- Multimedia (AREA)
- Evolutionary Computation (AREA)
- Health & Medical Sciences (AREA)
- General Health & Medical Sciences (AREA)
- Analytical Chemistry (AREA)
- Medical Informatics (AREA)
- Databases & Information Systems (AREA)
- Chemical & Material Sciences (AREA)
- Artificial Intelligence (AREA)
- Software Systems (AREA)
- Data Mining & Analysis (AREA)
- Computing Systems (AREA)
- Mathematical Physics (AREA)
- Optics & Photonics (AREA)
- Mathematical Optimization (AREA)
- Life Sciences & Earth Sciences (AREA)
- Mathematical Analysis (AREA)
- Algebra (AREA)
- Geometry (AREA)
- Pure & Applied Mathematics (AREA)
- General Engineering & Computer Science (AREA)
- Computational Mathematics (AREA)
- Biochemistry (AREA)
- Immunology (AREA)
- Pathology (AREA)
- Nuclear Medicine, Radiotherapy & Molecular Imaging (AREA)
- Radiology & Medical Imaging (AREA)
- Quality & Reliability (AREA)
- Microscoopes, Condenser (AREA)
Abstract
本发明公开了一种全划片高通量彩色病理成像分析仪器及方法,涉及一整套全划片高通量彩色病理成像分析仪器及方法,能够同时实现5‑10mm左右直径的视场大小、200‑1000nm空间成像分辨率,景深深度50‑300μm,高效数据采集4s,全彩色重构1min内,病理分析1min内,初筛分析置信度98%,自动批量扫描重构4‑1000片。有效解决了传统数字病理扫描仪诸多缺陷,如拼接伪影、拼接重影、成功率低、效率低、场景受限等,具有更广的应用场景、临床潜力与更低的仪器成本。
Description
技术领域
本发明属于光学仪器,具体属于光学成像仪器、装置与技术领域。进一步地涉及全划片高通量彩色病理成像分析***及方法,其可被应用于血液病理学、组织病理学、细胞病理学、免疫组织化学、免疫细胞化学等分析、筛查、教学和研究。
背景技术
病理学作为癌症及癌症早期筛查明确诊断的“金标准”,是连接基础研究和临床应用的“桥梁”。近年来,随着人民生活水平的提高,生命医疗领域的临床应用及科学研究,尤其是数字病理学领域对高端光学显微镜的市场需求仍以较大幅度增长。GrandViewResearch发布了全球显微镜市场研究报告,报告显示:2013年的全球显微镜市场容量为56.8亿美元,从2014年至2020年的年均复合增长率预计为7.7%,到2020年全球显微镜市场容量预计将达到95.4亿美元。新兴的发展中国家,特别是中国,随着产业转型升级的需求,高通量的数字病理扫描仪器的发展对高端显微镜的需求增速会更快。据西部证券预测,我国病理行业的潜在市场超 300亿元,其中组织病理市场规模20-30亿元,细胞病理中***筛查潜在市场超200亿,免疫组化病理潜在市场空间超40亿,分子病理潜在检验空间超50亿元。
面向国民经济主战场和人民生命健康需求,当前我国人民对幸福健康生活的需求不断提升,然而病理医生人员素质要求高、培养周期长,导致全球范围内病理医生相对短缺;相比于检验、影像科室,病理科的自动化数字化设备少且昂贵,从而综合导致诊断时间长,通常需要5-10个工作日。尤其是术中病理,对时间非常敏感。其整个诊断流程如图1所示,分为三或四步。
制备样本通常需要20-30分钟,取决于制片人的熟练度,之后传统病理学依赖人眼通过显微镜目镜观察读片,需要来回切换不同倍率物镜、平移视场,以获取病理特征。举例来说,对一块厚度为1 mm~1 cm的待测组织活检,通常需用切片机以每10 μm的厚度切100~1000片,使用人眼观察无疑会存在易漏判、误判、主观、效率低下等问题,通常需要耗费数天时间进行观察与分析。针对传统的基于人眼的读片与分析环节中存在的效率与准确率低下的问题,得益于数字成像器件的发展,现代数字病理学将原先的一个人工步骤分割为两个步骤,采用数字化成像配合机械扫描的方式,先获得高分辨率再获取大视场,实现对单一切片的全视场(视野)成像,这类技术称为全划片成像技术,与之对应的成像仪器称为数字病理扫描仪,后续的分析步骤则通过数字病理分析仪或软件来实现。
但现有的数字病理扫描仪仍存在四个由扫描方案带来的本质(内生)缺陷,如图2所示,具体技术缺陷如下:
其一,成像质量与成功率不高:电子器件漂移易导致拼接边界出现伪影,如图2(a)所示;通常的解决方案是对相邻图像采用特征匹配方法,但对于稀疏样本而言,由于存在大量空白区域且缺乏相邻特征,往往导致无法拼接,即拼接失败,如图2(b)所示为一典型的肺气肿切片拼接失败结果,拼接成功率与样本特性直接相关,且易发生特征匹配错误导致重影(Double Image),如图2(c)所示;
其二,应用场景有限:由于采用高数值孔径物镜,景深浅,样本必须按照规定要求如在10 μm厚度内进行切片,如图2(d)所示,因而仅适用于部分组织病理学和血液病理学的应用,仅横向扫描无法用于存在一定厚度的细胞病理学,而上述市场规模中提到细胞病理学在病理学中占据目前最大的市场;添加z轴扫描则会导致装置更为复杂,成本及维护难度都大大增加;
其三,效率相对还不是很高:平移、自动聚焦、成像、拼接四步工作流程中,大部分时间被用于机械运动,成像时间占比短,探测器大部分时间在空转,对于一般应用,效率还不是最突出的问题,但对于术中病理,时间非常敏感;
其四,仪器造价与维护昂贵:依赖于高精密电动平移台、高数值孔径物镜、脉冲光源等制式结构,成本高,且因制式设计相应的零部件不易更换,维护困难。
针对以上存在的四个大的方面问题,本发明基于傅里叶叠层显微成像技术,提出一种全划片高通量彩色病理成像分析仪器及方法,能够同时实现5-10 mm左右直径的视场大小、200-1000 nm空间成像分辨率,景深深度50-300 μm,高效数据采集4s,全彩色重构1min内,自动批量扫描重构4-1000片。此外,本仪器同时还针对了图1分析过程中主观、低效率等问题,提出基于AI辅助的快速初筛方案,初筛分析置信度98%,能够完成病理分析1 min内。另外,该仪器还可定量恢复出样本相位信息,可用于高通量定量生物学分析研究,如测量癌症的病理结构,癌症自动分阶段等。
傅里叶叠层显微成像术(Fourierptychographicmicroscopy, FPM) 是新一代的计算光学成像技术,由美国加州理工学院ChanghueiYang教授等人发明于2013年[ZhengG,HorstmeyerR, YangC. Wide-field, high-resolutionFourierptychographicmicroscopy[J]. NaturePhotonics, 2013, 7(9): 739-745.]。与传统显微成像中使用高NA物镜扫描拼接获得高分辨率大视场不同,FPM使用低NA物镜从而天生具有大视场,它融合了微波合成孔径和光学相位恢复的概念,在傅立叶频域中迭代拼接低分辨率图像,突破了现代显微成像***中视场范围与空间分辨率之间的权衡关系,同时实现了高分辨率、大视场、定量相位成像。相同参数下其效率远高于传统扫描拼接方案,属于十亿像素级高通量成像技术;更重要的是,其无需扫描获得大视场兼具长景深优势,从而从本质上避免了扫描拼接带来的成像质量问题,成功率和应用场景得到了极大的提升。它将原先的阿贝远场衍射极限公式λ/2NA改进为λ/(NA illu +NA obj ),其中NA illu 和NA obj 分别代表了照明和物镜数值孔径(Numericalaperture, NA),并被写入领域知名专家Goodman教授的《傅里叶光学导论》(第四版)中,相当于该技术原理已获得了同行广泛的认可。根据该公式,传统非相干照明分辨率只能提升2倍,改进照明后分辨率可以提升约7-8倍,最高可提升10倍水平。
专有名词解释
其中“全划片”准确英文名称为“WholeSlideImaging”,相关扫描拼接仪器研究发展可参见文献综述[FarahaniN, ParwaniA, PantanowitzL.Wholeslideimaginginpathology: advantages, limitations,andemergingperspectives[J]. PathologyandLaboratoryMedicineInternational 2015,7: 23-33.]一文。中文翻译存在多种名称如全玻片、全视野数字切片、全划片等,这里我们统一选用全划片这一专业词汇,指的是对一张切片的全视场(视野)成像。
高通量包含两层定义,其一,对光学***而言,时间分辨率、空间分辨率和视场大小三者是矛盾权衡关系,构成“闭环三角”,学术上常采用单位时间里可分辨的像素个数,即时空间带宽积或有效像素数来描述通量,其单位为像素数。高通量定义1即单位时间里可分辨更多有效像素数。上述特征又包含定义2,单批次可自动成像切片数量,即“批量(自动化)扫描”,英文准确术语为“BatchScanning”,高通量即单次可实现更多的切片数量,市场规格常有4片、8片、100片、200片、1000片等,其主要与仪器自动化输运存储样本有关,与光学技术与方法无关。
傅里叶叠层显微成像技术,英文准确名称为“Fourierptychographicmicroscopy”,简称FPM技术,由美国加州理工学院Yang等人发明于2013年[ZhengG, HorstmeyerR, YangC. Wide-field, high-resolutionFourierptychog raphicmicroscopy[J]. NaturePhotonics, 2013, 7(9): 739-745.],中文翻译存在多种名称如傅立叶重叠关联成像、傅里叶叠层显微成像技术、傅里叶叠层成像术、傅里叶叠层技术等,本文统一采用傅里叶叠层成像技术这一专业术语。
初筛置信度指的是对病理样本进行病理特征的判断过程,分为有病变与无病变两种结果,有病变的结果又会细分为多个阶段,部分疾病的病理特征已经有相应金标准。通过对一定数量的已知结果样本(同时包含有病变与无病变样本,且通常无病变占绝大多数)测试所训练的模型,其准确率即初筛置信度。置信度98%可以简单理解为,100片待测样本,98片的检测结果都是与事实相符的。而进一步地判断有病变的样本属于哪一病理阶段,其置信度在不同场景中是完全不同的,一方面与具体病种有关,很多疾病即使人类也未必认全其全部特征形态;另一方面与病理特征识别难易度、可训练样本数量与成像质量有关。因此属于行业难点,尚缺乏相关标准,故不作为具体仪器技术指标提供。
发明内容
与本仪器技术方案相关的物理底层技术思想是“先获得大视场后获得高分辨率”、“多角度照明提高分辨率”。然而如何将该相关技术思想,稳定可持续地将其从实验室迈入医院、乡镇卫生院、第三方实验室、科研院所、中小学等日常生活等应用场景,这对科学仪器的可重复、稳定等提出了更高的要求。随着观察样本增多,应用需求更具体化、多样化等导致突破该技术瓶颈,不单是是要解决原先未解决的部分技术难题,还需要通盘***化考虑原先未解决的***工程难题。在现实的研发过程中我们遇到的具体难题包括:
硬件层面:
1、技术问题方面:光源亮度不均匀对成像质量有明显影响。传统方案采用算法补偿[BianZ, DongS, ZhengG. AdaptivesystemcorrectionforrobustFourierptychograph icimaging[J]. OpticsExpress, 2013, 21(26): 32400-32410.]或者手动反馈调节[PanA, ZhangY, WenK, etal. SubwavelengthresolutionFourierptychographywithhem isphericaldigitalcondensers[J]. OpticsExpress, 2018, 26(18): 23119-23131.],在工程应用方面缺乏自动化与自适应能力;
2、技术问题方面:FPM技术为主动相干成像技术。受日光等杂散光影响导致实验室级的原理样机只能严格放置在暗室中工作[KondaPC, LoetgeringL, ZhouKC, etal.Fourierptychography: currentapplicationsandfuturepromises[J]. OpticsExpress,2020, 28(7): 9603-9630.];
3、技术问题方面:缺乏偏振成像、荧光成像等单台机多模态成像手段[PanA,ZuoC, YaoB. High-resolutionandlargefield-of-viewFourierptychographicmicrosco pyanditsapplicationsinbiomedicine[J]. ReportsonProgressinPhysics 2020, 83(9):096101.];
4、工程应用方面:作为产品而言,稳定、可靠、可重复是非常重要的,通常用户并不需要像实验室级样机那样太多的调节自由度,相关***最优参数的选取缺乏相关理论层面的分析支撑;
5、工程应用方面:传统实验室级的原理样机,最大只获取物镜视场数(Fieldnumber, FN)大小的视场,若样本尺寸超过本技术所提供的最大视场大小,缺乏相应解决方案;
软件层面:
1、工程应用方面:传统实验室级的原理样机,只针对一张切片样本,当面对批量样本时,缺乏自动化手段,缺乏***工程思考下的方案,导致其重构算法需要重复对每一片样本进行数字化像差补偿[OuX, ZhengG, YangC. EmbeddedpupilfunctionrecoveryforFou rierptychographicmicroscopy[J]. OpticsExpress, 2014, 22(5): 4960-4972.]与波矢校正[SunJ, ChenQ, ZhangY, etal. Efficientpositionalmisalignmentcorrectionmet hodforFourierptychographicmicroscopy[J]. BiomedicalOpticsExpress, 2016, 7(3):1336–1350.];
2、工程应用方面:缺乏针对大批量样本的数据采集与图像重构的协同方案,易浪费计算机资源,影响效率;
3、技术问题方面:其数字化补偿方案由于采用分块思想,缺乏对全局像差的约束[SongP, JiangS, ZhangH, etal. Full-fieldFourierptychography (FFP): spatially varyingpupilmodelinganditsapplicationforrapidfield-dependentaberrationmetrology[J]. APLPhotonics, 2019, 4(5): 050802.],同时其计算结果稳定性不足,易受栅格噪声影响[GuoK, DongS, NandaP, etal. Optimizationofs amplingpatternandthedesignofFourierptychographicilluminator[J].OpticsExpress, 2015, 23(5): 6171-6180.];
4、技术问题方面:传统需要依赖手动调焦,缺乏自动聚焦能力;
5、技术问题方面:传统全彩色化依赖RGB三色合成,精度高,但效率还不够快,且易受灰尘等造成的背景色差影响[GaoY, ChenJ, WangA, etal. High-throughputfastfull-colordigitalpathologybasedonFourierptychographicmicroscopyviacolortransfer[J]. ScienceChina-PhysicsMechanics & Astronomy 2021, 64(11),:114211.];
6、技术问题方面:传统重构算法易受相机噪声影响[YehL-H, DongJ, ZhongJ,etal. ExperimentalrobustnessofFourierptychographyphaseretrievalalgorithms[J].OpticsExpress, 2015, 23(26): 33214–33240.],对探测器灵敏度有严格要求,影响了仪器成本;
7、技术问题方面:缺乏获得全彩色高通量图像后的分析与辅助诊断能力。
为了解决现有技术相关问题或技术缺陷,本发明提供一种全划片高通量彩色病理成像分析仪器及方法。
本发明的技术方案,分别解决了如下技术问题:
硬件层面:
1、设计与制造了自适应强度可变角度照明器,从而确保不同角度照明下将各个入射角度光源到达样本平面亮度自适应达到一致;
2、设计加载在光源附近处的第一偏振片,用以消除如日光等杂散光或者散射光影响,结合高亮度的照明器件,从而使仪器无需在严格暗室下正常工作;
3、针对单台机多模态成像需求,提出可选的基于双偏振片组的偏振成像方案,提供用于激发荧光的紫外灯作为选配,从而实现可选的宽场荧光成像;
4、设计提出大量程x-y扫描位移台,对于超过视场大小的样本,通过结合传统的扫描拼接思想,实现每次大视场扫描拼接,由于扫描步长非常大,达到毫米量级时,反而做到仪器对位移台精度要求极低,工程结构极易实现;同样大小的拼接交叠率下,视场更大所交叠的内容会更多,拼接特征更多的前提下,不易存在拼接失败等问题。
软件层面:
1、提出快速分析仪器的初始化参数标定方法,通过本发明的***误差矫正算法矫正分析仪器***光波矢误差和像差影响,完成分析仪器初始化参数标定。***误差矫正算法是在经典高斯牛顿重构算法基础上,依次添加了基于模拟退火、非线性回归波矢矫正及全局像差约束的三个矫正模块予以实现。该方法只需要一次***校正,校正后的***参数将直接用于后续重构算法当中,相比于传统每一个样本都需要用算法同时校正参数与重构样本,提升了效率、重构精度和稳定性。同时只需要采集少量图像即可实现该操作,初始化速度快,且有效地抑制栅格噪声影响;
2、提出自动聚焦成像方案,通过图像反馈配合高精度电动z轴位移台控制样本或者物镜上下移动实现;
3、提出数字重聚焦方法,通过在重构的初始化光瞳函数中添加额外的离焦相位项实现景深延拓,能够有效解决样本超过一定厚度后因景深限制无法实现全视场自动聚焦技术问题。
4、提出了光学窗口的判别方法,用于选取最佳单波长照明,便于后续的颜色传递与重构。
5、提出改进的ADMM-FPM重构算法,用于实现高通量灰度图像重构,进一步消除对成像探测器灵敏度的要求。
6、提出高通量真彩色图像重构算法,能够实现快速真彩色图像重构,便于后续病理学切片的观测和识别。
7、提出了基于深度学习的方法进行图像病理分析,实现对高通量病理图像的初筛与辅助诊断,能够筛选出有病变的病理切片,并给出可选的病理区域与置信度。
附图说明
图1是基于人工母舰或者自动化仪器的病理诊断工作流程图;
图2是读片过程所使用的传统病理扫描仪缺陷示意图;图2(a)为电子器件漂移导致拼接边界出现伪影缺陷;图2(b)为稀疏样本由于存在大量空白区域缺乏相邻特征无法拼接缺陷;图2(c)为特征匹配错误导致重影缺陷;图2(d)为景深狭小导致无法聚焦的缺陷。
图3是非暗室环境下本仪器硬件连接图;
图4是仪器的原理示意图;图4(a)为光机电***及控制器内部光学部件的结构示意图;图4(b)为样品的高频信息可以移入光学***的低频通带内以低分辨率的信息表达出来;图4(c)为高分辨率、大视场、定量相位的图像信息。
图5是本发明的仪器硬件结合示意图;
图中各符号表示:1、自适应强度可变角度照明器,1-1、可变角度照明器控制器;2、第一偏振片;3、样本;4、样本存储器;5、x-y轴电动载物台;6、载物台控制器;7、显微物镜;7-1显微物镜控制器(可选);8、z轴电动平移台;9、平移台控制器;10、显微镜镜架;11、第二偏振片(可选);12、筒透镜;13、目镜(可选);14 、卤素灯、紫外灯(可选);15 、滤光片(可选);16、成像探测器;17 、处理器;18 、存储器;19、分析器;20、显示器。
图6是仪器装置各组件硬件(包括可选)功能结构连接图;
图7是自适应强度可变角度照明器的多种结构示意图;其中图7(a)为矩形阵列;图7(b)为方形阵列;图7(c)为圆环形阵列;图7(d)-7(f)为异形结构阵列。
图8是傅立叶频域交叠率定义示意图;
图9是自适应强度可变角度照明器亮度反馈原理图;其中图9(a)为平板LED照明模型图;图9(b)为平板LED归一化照明强度分布图;图9(c)为平板LED归一化照明强度分布的俯视图;图9(d)为数字半球形聚光镜照明模型图;图9(e)为数字半球形聚光镜的归一化照明强度分布图;图9(f)数字半球形聚光镜的归一化照明强度分布的俯视图,图9(g)为可变角度照明器控制器1-1的控制方法逻辑图。
图10是本发明的中视场大小和空间成像分辨率标定的示意图;
图11是本发明添加可选的偏振成像后的多模态成像原理图;其中图11(a)为偏振多模态成像的简易装置图;图11(b)为偏振多模态成像中使用的四种偏振配置;图11(c)为琼斯矩阵重构结果;图11(d)为物函数重构结果;图11(e)为CTF重构结果;
图12是不同规格的样本存储器实施例附图。图12(a)4片装置;图12(b)8片装置;图12(c)40片装置;图12(d)200片装置;
图13是显微镜内部光路结构设计图;
图14是本发明添加可选的宽场荧光成像后的多模态成像实施例附图。图14(a)宽场荧光成像结果;图14(b)宽场荧光成像去卷积后的结果;图14(c)宽场荧光成像去卷积后叠加样本相位信息的结果;
图15是本发明的仪器工作与方法流程图;
图16是初始化中***误差矫正的算法流程图,从而为后续重构提供固定的全视场像差和不同角度的准确光波矢;
图17是初始化中随视场变化的光学***数字像差恢复算法流程图;
图18是初始化中随视场变化的光学***全视场数字像差恢复与补偿实施例附图;
图19是初始化中光波矢误差矫正算法实施例附图;
图20是自动聚焦方法原理图与实施例附图;其中图20(a1)为清晰度评价值与焦距的函数关系曲线图;图20(a2)为清晰度比率的定义曲线图;图20(a3)为陡峭峰的定义曲线图;图20(a4)为灵敏度因子的定义曲线图;图20(a5)为爬山搜索算法示意图;图20(b1)为样本在正向离焦14μm处的图像;图20(b2)为样本在负向离焦14μm处的图像;图20(b3)为样本在焦距处的图像。
图21是本发明自动聚焦方法流程图;
图22是本发明中通过数字重聚焦方法进行景深延拓的实施例附图。
图23是本发明的快速全彩色方案中光学窗口的选择方法;(a)已知染色剂及其光谱吸收曲线示意图;(b)未知染色剂及其光谱吸收曲线的处理方案;
图24是本发明的图像重构算法流程图;
图25是图像重构算法的实施例附图;
图26是快速全彩色化的彩色传递原理图;
图27是本发明的彩色传递思想示意图。(a)颜色匹配: 传递色调信息; (b)彩色传递: 传递彩色纹理信息;
图28是本发明的彩色传递算法流程图;
图29是彩色空间建立与跨设备显示解决方案原理图;
图30是本发明的特征金字塔网络FPN结构图;
图31是本发明的区域建议网络RPN结构图;
图32是目标检测中确定每个目标所在具***置的示意图;
图33是非极大值抑制思路逻辑图;
图34是双线性插值的原理图;
图35是ROIPooling操作原理图;
图36是ROIAlign操作原理图;
图37是本发明针对某一病理的数据分析算法模型训练流程图;
图38是以***筛查为实施例经过本***分析后的结果实例图;
图39是批量样品检测协同并行处理时序逻辑示意图。
具体实施方式
下面结合附图和具体实施方式对本发明技术内容作进一步详细说明。
图3展示了本发明仪器装置必备的4个部件。图4展示了该发明的基本物理原理,其中图4(a)为本发明仪器的简化原理结构图,通用设备及一些不涉及物理原理的光学元件被隐去。使用可变角度照明器替代传统光学显微镜的聚光镜,依次点亮不同角度照明器改变照射样品的光线角度。通常,物镜具有一定的收集光的能力,由物镜的NA表征。当照明角度在物镜NA之内时,成像探测器记录明场图像;当照明角度大于物镜NA范围时,光与物质(细胞器)相互作用,会发生瑞利或者米氏散射(取决于细胞器尺寸),从而有部分散射光也会被收集,成像探测器记录暗场散射图像。根据傅里叶变换的平移定理,多角度照明的实质是将样品的频谱在傅立叶频域中进行大小不同的平移,因而样品的高频信息可以转移到光学***的通频带内以低分辨率(lowresolution, LR)图像表达出来,如图4(b)所示。与传统的使用高NA物镜扫描拼接获得高分辨率大视场不同,FPM使用低NA物镜从而天生具有大视场,在傅立叶频域中迭代拼接不同角度照明下获取的低分辨率图像,如图4(c)所示,进而获得高分辨率、大视场的定量相位图像。它同时兼具相位恢复和合成孔径的思想,无需机械扫描和重聚焦就能获得十亿像素级图像。与传统FPM***不同[ZhengG, HorstmeyerR, YangC.Wide-field, high-resolutionFourierptychographicmicroscopy[J].NaturePhotonics, 2013, 7(9): 739-745.],在硬件方面,本装置采用自适应亮度可变角度照明器解决了工程上自动化与自适应的需求;增加了偏振片1与偏振片2(可选),该偏振片均为线偏振,配合高亮度的可变角度照明器,使得曝光时间极短,从而让仪器装置能够摆脱暗室的苛刻要求。进一步地,当两偏振片旋转至多种偏振方向组合时,还可实现偏振成像的多模态成像功能。
具体地,参见图5,本发明为一种全划片高通量彩色病理成像分析仪器,其核心部件至少包括:自适应强度可变角度照明器1,用于实现在不同的时间以多个入射角度照明样本3,最终各入射角度的光源到达样本平面时的亮度保持一致;
光机电***及其控制器,至少包括一个偏振片2、一个显微物镜7、一个扫描位移台,利用偏振片结合高亮度光源实现快速曝光和杂散光抑制,再利用扫描位移台实现对多个样本扫描,进一步通过显微物镜7收集经过不同扫描样本发出的光;所述的多个入射角度和显微物镜7在傅立叶频域中的多个重叠区域相对应;
成像探测器16,用于获取光机电***及其控制器输出的多个强度图像,每个强度图像分别对应自适应强度可变角度照明器生成的各个入射角度;
处理器17,至少用于实现重建比探测器输出的强度图像分辨率更高、图像精度更高的可重复图像,定量恢复样本相位信息。
进一步地,本发明的处理器17还可以实现分析仪器的初始化参数标定、自动聚焦成像。当处理器在实现自动聚焦成像失败时,还能通过数字重聚焦方法进行景深延拓;处理器还能用于实现快速高精度全彩色成像。
一、全划片高通量彩色病理成像分析仪器硬件部分
以下从4个方面阐述本发明的硬件部分:
1、关于自适应亮度可变角度照明器及其控制器
自适应强度可变角度照明器是本发明的核心部件,其光源是由多个高亮度独立光源单元间隔排列而形成的阵列式高亮度光源,从而无需机械扫描实现各个高亮度独立光源单元的照明角度切换。参见图7,阵列式高亮度光源排布形式为平板矩阵列排布、平板圆环阵列排布、曲面半球体阵列排布或异形变体阵列排布。
高亮度独立光源单元为R/G/B三波长光源或者R/G/B/W四波长光源,其中R为红光,中心波长范围为600-740 nm;G为绿光,中心波长范围为492-590 nm;B为蓝光,中心波长范围为400-490 nm;W为白光,中心波长范围为540-580 nm;不同波长的峰值半高宽范围均优选为20-40 nm。本实例均控制在30 nm内。
自适应强度可变角度照明器的光源采用LED、激光二极管、OLED、LCD或miniLED;所述的自适应强度可变角度照明器的高亮度独立光源单元的每个通道刷新帧率优选范围为20-1000Hz,电功率优选范围为0.2-3W。
例如:将自适应强度可变角度照明器的光源选择为LED时,考虑到曝光时间与成像效率,推荐刷新帧率为500Hz及以上。本实例中采用500Hz,将最短曝光时间控制为2ms,在2ms及以上的曝光时间下不会发现LED闪烁现象,从而能够保持LED亮度的稳定。为了控制成像效率且与成像探测器相匹配,单张图像曝光时间应控制在1ms-100ms为佳,本实例中采用2 ms。
另外,LED的电功率与亮度、曝光时间存在一定关系,电功率越高,亮度越高,曝光时间越短。对于常规的商用探测器,单颗单通道电功率需满足0.2-3 W,本实施例采用的单颗单通道电功率为1W。
当阵列式高亮度光源排布形式为平板矩阵列排布或平板圆环阵列排布时,考虑到照明高度以及傅立叶频域交叠率,独立光源单元的间距d优选为3 mm~6 mm,本实施例中选择了4 mm间距。
阵列式光源距样本的高度h优选为45-80 mm。距离太低导致傅立叶频域交叠率过低(低于40%),成像质量受影响,并且每颗LED的照明光无法看作平面波,产生明显的波矢失配;距离太远导致许多能量被无意义耗散掉,仪器的体积重量和成本都会增加,且傅立叶频域交叠率过高(高于80%)也会引发反作用。因此光源工作高度为45-80 mm是较为合理的区间。光源器件可选配备(自动/手动)升降装置,将傅立叶频域交叠率控制在40%-80%以内,以获得较好的成像结果。本实施例为获取接近最大的成像通量,且保证成像结果的稳定性,因此不添加升降装置,将高度固定为70 mm,傅立叶频域交叠率控制在55%。
高亮度独立光源单元为平板矩阵列排布时,阵列排布数量范围为3×3~64×64;高亮度独立光源单元为平板圆环阵列排布时,阵列排布数量范围为4~1000。
阵列式高亮度光源的照明数值孔径(NA)计算方式为:
1)当阵列式光源排布形式为平板矩阵列排布或平板圆环阵列排布时,其照明NA为:
其中m、n分别为阵列式排布光源分别沿行、列排布的独立光源单元的个数,d为各个独立光源单元的间距,η为折射率,h为阵列式光源距样本的高度;
2)当阵列式光源排布形式为曲面半球体阵列排布或异形变体阵列排布时:
其中η为折射率,θ max 为最大入射角;
所述的阵列式高亮度光源的照明NA优选为0.3-1。
阵列式高亮度光源的傅立叶频域交叠率计算方式如下:
1)当阵列式排布光源的各个高亮度独立光源单元为均匀分布时,该光源的傅立叶频域交叠率计算如下:
其中S t 是NA的增加步长,NA obj 为物镜的数值孔径;
2)当阵列式排布光源的各个高亮度独立光源单元为非均匀分布时,该光源的傅立叶频域交叠率通过加权平均计算获得:
其中S t,i 是第i圈高亮度独立光源单元NA的增加步长,NA obj 为物镜的数值孔径,m为高亮度独立光源单元总数,n为第i圈高亮度独立光源单元的个数;
参见图8,由于物体频谱与CTF之间是相对运动的关系,FPM的多角度照明可以理解为:物体频谱不动,CTF在物体频谱上进行扫描。图8展示了物体频谱上两个相邻的扫描CTF的位置关系,傅立叶频域交叠率则定义为相邻两CTF重叠面积与单孔径面积的比值。阵列式高亮度光源的傅立叶频域交叠率优选40%-85%。
本实施例以19×19平板方阵列为例,该平板方阵列易制作,满足稳定、可量产等工程应用要求。对于任一阵列或者异形变体而言,LED个数可选,如3×3、5×5、32×32阵列等。LED的排布方式决定了最大照明角度和交叠率,从而决定成像分辨率与成像质量。经测试,19×19、均匀间隔4 mm平板方阵是能够获得最大分辨率的尺寸,此时最大有效照明NA约为0.6-0.7。
参见图9,LED照明强度不均匀会导致重建质量的下降。图9(a)为平板LED照明的模型图,通过将采集到的明场图像取平均强度作为该LED的照明强度,归一化后绘制了图9(b)所示的照明强度曲线图和图9(c)所示的俯视图。照明强度的能量衰减与cos 4 θ的黑色曲线基本吻合。在照明角度为34度左右时,能量衰减一半,此时强度波动为50%,而在照明角度达到50度时,有效能量强度已经低于20%,强度波动达到80%,所造成的成像伪影会非常严重。图9(d-f)同理为异形LED阵列结果,其能量衰减正比于cosθ,这两种照明方式有着cos 3 θ的能量差别,这意味着在40度入射角时平面LED相比半球形聚光镜能量低50%,77度入射角时能量低99%,几乎无光,在50度入射角时半球形聚光镜是平板LED的6.5倍,根据实验测算,平板LED真实可用的照明角度为39度(照明NA为0.63),约为能量衰减至36%时的水平,对应的,半球形数字聚光镜需要在照明角度为68度时,能量才会衰减至相同水平。又因为,实际中半球形数字聚光镜的正入射LED能量就比平板LED高,所以实际归一化后,半球形数字聚光镜的极限还会有小幅度提升,达到所设计的71.8度(0.95NA)照明角。
自适应强度可变角度照明器1,还包括与处理器17相连的可变角度照明器控制器1-1,用于根据成像探测器的反馈控制每个高亮度独立光源单元,阵列式高亮度光源到达样品平面亮度自适应地保持一致。
参见图9(g),可变角度照明器控制器1-1的控制方法如下:
1)选择大于阵列式高亮度光源的NA,确保成像探测器拍摄任意角度照明下均为明场图像;
2)固定曝光时间,取灰度平均值作为每一角度抵达样本平面的亮度值;
3)以正入射光源为基准,计算亮度测量系数矩阵;
4)根据亮度测量系数矩阵反馈控制每个角度的到达样本平面亮度,使其均匀一致;所述的反馈控制通过基于亮度测量系数矩阵调节单独光源的电功率;
所述的亮度测量系数矩阵定义为:
其中,φ e m,n 表示低分辨估计光场,I c m,n (x,y)表示第x行y列高亮度独立光源单元得到的图像。
自适应强度可变角度照明器1的强度校正方法描述如下:不放置任何样本,通过使用高NA物镜(如100×/0.9NA),完全覆盖(大于)照明NA,使成像探测器拍摄任意角度照明下均为明场图像。取灰度平均值作为每一角度抵达样本平面的亮度值,如图9(b)和图9(c)所示。根据该亮度测量系数矩阵反馈控制每个角度的达到亮度。以中央LED为基准,调节其余LED亮度至与中央LED亮度相一致。实施例中,每个角度可有65535种(216)亮度(电压值)可选,从而可以精细自适应地调节每一个角度的亮度情况,在确保亮度误差在±1%内后停止,其极限亮度误差可控制至±0.1%以内。实际中,考虑到算法具有一定鲁棒性,根据经验,亮度误差在±10%以内都是可以接受的。
本发明的自适应强度可变角度照明器1还有一个重要的技术特点在于,它还能够通过数字化调节具体投入使用的高亮度独立光源单元个数从而控制照明NA实现整个成像分析仪器分辨率可变,该成像分析仪器的横向空间分辨率设定方式如下:
其中NA syn 是成像分析仪器的合成NA,NA obj 为物镜的NA。
根据本发明的横向空间分辨率设定方式,参见图10,横向空间分辨率最大可接近200nm光学衍射极限,故横向空间分辨率优选范围200nm-1000nm,同时本发明的视场能够达到4倍物镜下的极限视场5mm左右。
可变角度照明器控制器1-1,还能够作为信号发生器输出上升沿或下降沿信号,用于硬触发成像探测器完成自动序列采集图像,所输出的上升沿或下降沿电压范围优选3V-12V。本实施例采用5 V的硬触发模式。
2、关于光机电***及其控制器
如图4(a)及图6所示,光机电***及其控制器至少包括第一偏振片2、显微物镜7、一个扫描位移台,利用偏振片结合高亮度光源实现快速曝光和杂散光抑制,再利用扫描位移台实现对多个样本扫描,进一步通过显微物镜7收集经过不同扫描样本发出的光;所述的多个入射角度和显微物镜7在傅立叶频域中的多个重叠区域相对应;
所述的扫描位移台为单轴或双轴x-y轴电动载物台5及载物台控制器6,用于实现大尺寸样本的扫描拼接重构和多个样本的扫描;所述的大尺寸样本为10mm×10mm视场范围以上,多个样本为2个及以上。x-y轴电动载物台5及载物台控制器6用于实现高通量装载样本,样本存储器4用于更替样本。
光机电***及其控制器还包括样本3及样本存储器4、z轴电动平移台8以及升降台控制器9、筒透镜12;其中样本存储器4用于存储大批量待测或已测的样本3,该样本个数已超过直接由载物台装载个数极限,样本3及样本存储器4通过z轴电动平移台8实现交互;z轴电动平移台8及平移台控制器9,通过自动聚焦成像方法获取具体离焦量后调整z轴位置至实焦面实现自动聚焦。
自适应强度可变角度照明器1照明样本,经过第一偏振片2和高亮度光源抑制杂散光后到达样本平面,显微物镜7收集经过扫描样本发出的光,由筒透镜12反馈至成像探测器16,然后通过z轴电动平移台8及平移台控制器9实现自动聚焦,借助单轴或双轴x-y轴电动平移台台及其控制器实现大尺寸样本的扫描拼接重构和通过与样本存储器交互实现样本批量扫描。
处理器17还可以分别与单轴或双轴x-y轴电动载物台5及载物台控制器6、显微物镜控制器7-1、z轴电动平移台控制器9相连,显微物镜控制器7-1与显微物镜7相连。
光机电***及其控制器内,在样品与探测器之间还设置有第二偏振片11,用于实现偏振成像等多模态成像。
两块可插拔式偏振片(2、11),第一偏振片2用于产生线偏振光,抑制环境光等杂散光影响,配合高亮的照明器,获得极短的曝光时间,突破暗室限制,使仪器可于自然环境下工作。第二偏振片11可以为同方向的线偏振片,从而形成常规的“起偏器-检偏器”模式,进一步抑制杂散光。实际中,一个偏振片已经能够起到足够效果。该技术方案原理是自然光等环境光为非偏振光,漫反射光线为部分偏振光,通过偏振片的偏振方向选取能够起到滤除光线的作用。
第一偏振片2位置可位于样品至光源间任意处,实施例中采用薄膜偏振片,紧贴光源阵列放置。
第二偏振片11位置位于样品后至探测器前任意位置,实例中放置于筒透镜至成像探测器间,如图11(a)所示。
此外,当旋转两偏振片偏振方向至不同组合时,第二偏振片11还可用于配合第一偏振片2实现偏振成像等多模态成像,其原理图与实施例如图11(b)所示。
偏振成像的具体步骤如下:
其中,为阵列式高亮度光源第个LED发射的倾斜平面波,为经
过第一个偏振片后的矢量场,为复值标量,其中,
是样本频谱的琼斯矩阵,是光瞳函数的琼斯矩阵,是第二个偏振片,其中,是二维傅里叶变换矩阵,是各个高亮度独立
光源单元在像面产生的矢量场。
参见图11,图11(a)所示的偏振FPM光学装置由标准偏振光显微镜和LED阵列光源
组成,使用起偏器偏振光对样本进行多角度照明,散射光通过检偏器后被图像传感器记录
下来;图11(b)展示了每个LED照明下采集的图像所对应的四种偏振配置,其中起偏器和检
偏器的角度分别设置为;图11(c)展示了蚕豆根截面的基于4
通道琼斯矩阵的偏振FPM重建结果;图11(d)是利用特征向量分析从琼斯矩阵中提取的各种
与偏振相关的样品属性;图11(e)展示了偏振FPM对CTF的恢复结果。
参见图12,本发明的样本3为标准生物切片样本,尺寸为25×75mm载玻片,170 μm厚盖玻片(可选)。为了实现高通量技术特点2,样本存储器4用于批量(自动)成像,可为4、8、20、40、100、200等规格。
筒透镜12还通过显微镜架10连接有目镜13,可用于人眼的观察。显微物镜7还通过滤光片15连接卤素灯或紫外灯14,用于实现非相干成像和荧光多模态成像。
显微物镜7、显微物镜控制器7-1,放大倍率通常选用2×-20×,数值孔径范围为0.08NA-0.75NA。商业物镜的视场数通常为22-27mm,因此在2×物镜下实现本仪器的标称的10 mm直径视场是绰绰有余的。另外,可通过显微物镜控制器7-1实现物镜塔台自动旋转,本实施例只采用一颗4×/0.1NA物镜,因此无需显微物镜控制器。
z轴电动平移台8及平移台控制器9,用于实现自动聚焦。自动聚焦可以采用调节物镜高度或者调节载物台高度实现。本实例采用调节物镜高度实现自动聚焦。100-500 μm的量程和1 μm及以下的精度已经足够,本实例采用300μm量程和1μm精度。
显微镜镜架10,属于机械支撑结构,包含倒置和正置两种,本实施例选用三目倒置显微镜结构,其内部结构设计如图13所示。其中可机械切换反射镜用于切换光路是目镜观察还是相机观测,仪器光束出口处放置光学窗口,用以密封整个仪器,防止进入灰尘等影响仪器性能,光学窗口外为常规的c型接口用于连接相机。相机前装有1×或者2×的适配器,本实施例选择1×,从而使得相机法兰距满足物镜的要求,且物镜是严格的标称倍率,如物镜标称4×,则实际测量结果也是4×,误差控制在±0.1倍以内。
根据用户需求可以选择不同焦距制式,由于不同厂商制式不同,显微物镜7和筒透镜12需要匹配设计。如图13所示,根据镜架的机械设计,本实例采用200 mm焦距一英寸消色差筒透镜,使得筒镜到相机的距离为200 mm,符合尼康等显微厂商制式。用于人眼直接观察的目镜13,通常有10×、20×等几种,本实施例选择20×。
参见图13,侧边光源为可选,可用于单台机实现非相干成像、荧光成像等多模态成像。如果要实现反射式非相干成像,则选配卤素灯。如果要实现宽场荧光成像,则选配紫外灯,分光镜需替换为二项色镜。图14展示了对HeLa细胞的荧光多模态成像结果;图14(a)的宽场荧光成像可以提供更好的观测对比度;图14(b)在前者的基础上加入去卷积操作,图像清晰度得到进一步提升;图14(c)将荧光图像与FPM重构的相位信息结合,便于观察细胞核的结构与位置。
3、关于成像探测器
其中,λ为光波波长,Δx为成像探测器的像元尺寸,M为放大倍率;所述的成像探测器的像元尺寸Δx优选3-7μm,成像探测器大小优选为对角线大于等于20mm。
像元尺寸通常有3.75μm、4.5μm、5μm、5.5μm、6μm、6.5μm等若干种;动态范围位数通常为8位、12位、14位、16位,相关参数之间是互相权衡的制约关系,针对本仪器应用需求,需要仔细合理选择。同时该***可扩展至X-ray或者电子领域,换为辐射计量器即可。本实施例采用6.5μm像元,2048×2048 像素,滨松FlashV3_4.0,16位sCMOS相机。
除USB3.0接口(最大带宽320MB/s)相机外,其它相机均需要采集卡(可选),如千兆以太网(GigabitEthernet, 最大带宽100MB/s)、CameraLink(单根最大带宽850MB/s)、CoaXPress(单根最大带宽6.25GB/s)等。本实施例选择Camera Link采集卡。
4、关于通用仪器
参见图6,处理器17上还可以连接存储器18、分析器19和显示器20,构成通用仪器;所述的存储器18用于实现数据存储及压缩处理,分析器19利用人工智能识别、分类图像中的病理信息,输出图像中的病理区域与置信度,显示器20实现图像显示。显示器20还可以连接通信器实现屏幕共享或远程医疗。
处理器17可以是笔记本、台式机、图形工作站、服务器、嵌入式集成***如FPGA等,本实例采用台式图形工作站。存储器用于数据存储,考虑数据量巨大,推荐配置大容量硬盘,并采用相应压缩处理。本实施例采用标准I5特定压缩格式。分析器19与处理器17并行计算,实现图像分类、识别分析等辅助诊断。
二、全划片高通量彩色病理成像分析仪器软件部分
参见图15,本发明还公开了一种全划片高通量彩色病理成像分析方法,该方法包括如下步骤:
步骤一:成像分析仪器初始化参数标定;
步骤二:装填待测样本;
步骤三:点亮自适应强度可变角度照明器正中央的高亮度光源单元,依次单通道正入射照明样本,成像探测器采集低分辨率图像并存储,合成低分辨率彩色图像;
步骤四:自适应强度可变角度照明器自动采用特定波长扫描样本,通过成像探测器采集序列低分辨率图像并存储,直至完成待测样本采集;
步骤五:处理器重构高通量灰度图像;
步骤六:处理器将低分辨率真彩色信息传递给高通量灰度图像,获得高通量真彩色图像;
步骤七:显示器输出高通量真彩色图像。
下面详细介绍本发明的方法相关步骤的涉及的详细技术内容:
1、关于成像分析仪器初始化参数标定
参见图16,步骤一:成像分析仪器初始化参数标定,是通过***误差矫正算法矫正分析仪器***光波矢误差和像差影响,完成分析仪器初始化参数标定;***误差矫正算法是在经典高斯牛顿重构算法[YehL-H, DongJ, ZhongJ, etal. Experimentalrobustness ofFourierptychographyphaseretrievalalgorithms[J]. OpticsExpress, 2015, 23(26): 33214–33240.]基础上,依次添加了基于模拟退火、非线性回归波矢矫正及全局像差约束的三个矫正模块予以实现。
***误差矫正算法包括如下步骤:
步骤1:放置用于标定的已知样本,计算每个高亮度独立光源单元对应的傅立叶频域位置信息,通过成像探测器捕获一组低分辨率强度图像,根据高斯牛顿重构算法的重构结果判断是否存在光波矢误差,若存在则启用***误差矫正算法;若不存在,则认为照明波矢位置准确,记录照明波矢位置,用于设定后续每一待测样本重构算法中波矢位置,同时将高斯牛顿重构算法输出的全视场数字化像差保存,用于后续每一待测样本重构算法中的数字化像差补偿;
步骤2:计算阵列式排布光源的各个高亮度独立光源单元对应的频谱信息;
步骤3:在阵列式高亮度光源照明条件下,在多个方向随机进行傅立叶频域步长搜索,更新待测目标频谱信息;
步骤4:通过频谱信息误差矩阵,更新每个高亮度独立光源单元对应的频谱信息;
步骤5:借助非线性规划计算阵列式高亮度光源的整体偏差参数;
步骤6:通过自适应步长多次迭代,实现全局像差约束,得到阵列式高亮度光源全视场像差值。
我们以15×15LED阵列式高亮度光源为例,对于每个LED m,n (m行n列)及其照明波矢量(u m,n ,v m,n ),成像探测器捕获LR(低分辨率,lowresolution)强度图像,其公式表达如下:
其中是傅里叶逆变换算子,O是样品透射函数o的傅里叶谱,P(u,v)是光瞳函
数,用作成像***的低通滤波器,(u,v)是傅里叶平面中相对于(x,y)的二维空间频率坐标。
入射波矢量(u m,n ,v m,n )可以表示为:
其中(x 0 ,y 0 )是每个小段的中心位置,x m,n 和y m,n 表示LED元件在第m行第n列的位置,λ是照明波长,h是LED阵列和样品之间的距离。
然后,对样本估计对应的频谱区域进行如下更新:
其中α和β是更新的步长,通常取值为α=β=1。δ 1、δ 2为正则化常数,防止分母为零,i为迭代次数,Δφ i,m,n 为更新过程的辅助函数。为了获得最佳的鲁棒性和收敛效率,我们在程序中设置δ 1=1,δ 2=1000。整个迭代过程重复i次,直到获得收敛解,算法收敛通过每次迭代的误差度量值来判断,误差度量值的定义如下:
在理想的FPM设置中,***参数通常是准确的,但在实际情况下,***参数会出现各种形式的偏差。在此模型中,设置三个全局变量来修正误差:LED阵列旋转角度θ、中心LED沿x轴和y轴的位移Δx、Δy。每个LED高亮度独立光源单元的位置可以表示为:
其中d LED 表示相邻LED独立光源单元的间距。
首先,给出样本频谱O o (u,v)的和光瞳函数P o (u,v)的初始猜测来启动算法。其次,定义每次迭代的LED更新范围S i 。通常情况下,所有低分辨率图像交替迭代,以更新样本频谱和光瞳函数。然而,在迭代重建过程中,低频分量占主导地位。因此,在迭代的初始阶段,模块1和模块2对具有低照明NA(数值孔径)的明场图像进行多次初始迭代,以纠正这些低频光圈的参数。
模块一在前十次迭代中,仅使用中央5×5的明场图像,每两次迭代为一组,奇数次迭代更新强度系数,偶数次迭代重置初始位置参量(Δθ,Δx,Δy)。在对明场图像进行5组10次迭代后,将采集的所有图像带入后续的迭代过程,以获得所有LED单元的精确位置坐标。其中,每次迭代的LED更新范围S i 定义为:
令R表示算法的搜索方向,(Δu,Δv)分别为模拟退火算法的降温速度,即搜索步长。通常R取8,即逐次向8个方向搜索,可以在精度和效率上权衡最优的结果。在LED m,n 照明下,物体频谱信息更新如下:
其中(Δu m,n,R ,Δv m,n,R )表示第R 次搜索过程中采用的傅立叶频域搜索步长。定义一个额外的误差矩阵为E(R)用于判别是否找到了更优的位置:
随着迭代次数的增加,像面出射光场的强度将越来越接近于相机实际拍得的图像强度,即E(R)会越来越小。令q表示一维误差矩阵E(R)取最小值时所对应的搜索方向,则LED m,n 在频谱面的坐标可更新为:
随后的更新步骤与原始FPM算法相同。这里配合前十次迭代将每个偶数次迭代过程中的搜索步长(Δu,Δv)压缩为原来的一半,减至两像素时不再减少,迭代第十二次开始减至1 像素随后保持不变,直至算法迭代次数的增多,重构结果趋于稳定,(Δu,Δv)取0,LED m,n 的傅立叶频域坐标逐渐接近于真实值。(Δu,Δv)的初始值需要根据具体误差严重程度设定,通常设定在4 像素或者8 像素。考虑到目标函数同时存在振幅和强度两种,式(9)可能也存在另一种形式,如下:
模块二由LED阵列位置误差参数(Δθ,Δx,Δy)组成的目标函数T:
其中u(Δθ,Δx,Δy)和v(Δθ,Δx,Δy)表示LED m,n 的傅立叶频域坐标(u,v)随(Δθ,Δx,Δy)的变化关系,即:
根据非线性规划原理,当目标函数T取最小值时对应的参量值
即为更新后的LED阵列整体偏差参数,即:
(Δθ,Δx,Δy) update = argmin[T(Δθ,Δx,Δy)] (25)
模块三在每次迭代中,令光场更新为:
其中下标m, j分别表示第m个全视场分块和第j个LED照明角度,样本频谱通过以下公式更新:
3)对所有像差系数进行更新,以系数b 1为例,更新方式为b 1= b 1-αΔb 1,
其中参数α为更新步长。
4)重复上述步骤进行多次迭代,最终全局像差的系数应为:
参见图17,使用***误差矫正算法恢复光学***数字化像差的可视化流程。对输入的光瞳函数初始猜测进行迭代更新,将最终输出的全视场数字化像差保存,用于后续每一待测样本重构算法中的数字化像差补偿。
参见图18,使用***误差矫正算法对光学***数字化像差进行优化补偿的前后对比结果。由于采用周期性LED阵列照明,CTF会出现明显的栅格噪声。使用全局像差优化后,噪声伪影消除,像差得到精确重构。
参见图19,使用***误差矫正算法校正光波矢误差,该算法的x-y精度达到0.1mm,角度精度约为0.003rad(0.17°)。
2、关于装填待测样本
初始化完成后,样本存储器自动将标准尺寸为25×75mm载玻片样本装填到x-y轴电动载物台上。实施例所用载物台一次性可装载四片载玻片,当LED阵列扫描完样本后,电动载物台会将下一个样本移到中心位置。
3、自动聚焦方法
步骤二:装填待测样本之后,还包括步骤A自动聚焦待测样本;如果自动聚焦成像失败,通过数字重聚焦方法进行景深延拓。
所述的步骤A自动聚焦待测样本,具体通过聚焦窗口的选择、图像清晰度评价、搜索策略确定实现自动聚焦成像。
现有技术就聚焦窗口的选择进行了大量研究,提出了多种区域划分方法,比较经典的主要有:中心区域聚焦窗口选择法、多区分布聚焦窗口选择法和一阶矩聚焦窗口选择法。
1)聚焦窗口的选择
1.1)中心区域聚焦窗口选择法
中心区域聚焦窗口选择法是基于图像中心区域的聚焦窗口选择法,就是一种在图像的中心部分按照一定比例进行区域划分,并以此作为聚焦区域的一种窗口选择方法,对分辨率为M×N的图像,中央区域调焦窗口选择法的公式为:
中心区域聚焦窗口选择法采样较为简单,但当图像主体并不处于中心位置时就会使得聚焦结果产生严重偏差,且缺乏自适应调节能力。
1.2)多区分布聚焦窗口选择法
中心区域聚焦窗口选择法的采样区域固定在中心位置,自适应性能较差,当图像主体偏离中心区域则会严重影响聚焦效果。因此,在中心区域法的基础上,拓展出一类多区域聚焦窗口采样方法,主要包括黄金比例窗口采样法和倒T型窗口采样法。
1.2.1)黄金比例窗口采样法
基于黄金比例,摄影者在拍摄大型场景时会使用“ 三分法” 进行构图,以避免对称性构图,使主体景物更加自然;在拍摄较小型物体时,常采用“井分法”进行构图,将主体景物置于横竖分割线的交叉点,在突出主体的同时能以此为锚点引导观众的目光向周围扩散。根据美学和摄影艺术理论,除了中心点,图像中的四个黄金分割点也是人眼的视觉注意点。将四个黄金分割点围成的区域作为聚焦窗口过于狭隘,为了将图像信息放于聚焦区域内,可以将经过中心点和四个黄金分割点的黄金分割十字窗口作为聚焦窗口,计算公式如下:
1.2.2)倒T型窗口采样法
根据摄影者的拍摄习惯,一般倾向于把图像主体置于整体的下半区,因此产生一种倒T形状的区域划分方式。根据合适的比例把图像划分为多块,并取中下部的四个区域组成聚焦窗口。对应表达式如下所示:
相比单窗口的中心区域聚焦窗口选择法,多区分布的聚焦窗口选择方法通过综合考虑图片的多个部分,一定程度上减小了主体偏离图像中心位置造成的影响。
1.3)一阶矩聚焦窗口选择法
静态窗口采样法通常会以图像的几何中心为原点向周围扩散来划分聚焦区域,当主体部分不在图像几何中心附近时就会使聚焦结果产生较大的偏差。即通过计算图像的灰度一阶矩,得到主体所处的实际“重心”位置,再以此为基础划分窗口区域。
对于普通灰度图像,图像的实际“重心”位置可由如下所示的公式求得:
式中M和N表示图像的尺寸,f (x,y)表示该坐标点处的像素值,⌊⌋表示对数值取整。考虑到灰度图像易受不同程度的亮度的影响,该方法提出改进,先进行边缘图像的梯度化,仅提取边缘图像的细节信息,以忽略亮度对灰度图像的影响。对梯度值图像,图像重心位置可由下面的表达式得到:
图像一阶矩的取窗方式既保证了成像主体处于调焦窗口中,又减小了背景对聚焦效果的负面影响。图像一阶矩的取窗方式对成像主体的位置具有一定的自适应能力,适用于成像主体较小且偏离中心的情况。
2)图像清晰度评价函数
图像清晰度和图像锐度正相关,图像锐度与对比度有关,因此我们可以使用对比度值作为图像质量的评价值。目前,许多对比度聚焦法己用于自动聚焦***。对比度评价函数主要有四类:空域梯度函数、基于信息统计的函数、基于相关性度量的函数和基于变换域的函数等。F是不同的评价函数,f (x,y)是图像(x,y)处的灰度值。
2.1)空域梯度函数
基于图像分散的特点,离焦图像通常在黑暗和明亮的区域之间的差异较轻微,因为点扩散函数将一个像素的强度分散到几个像素之间,模糊它们并平均它们的灰度值。清晰图像相邻像素间的灰度差相对大,模糊图像的相邻像素经过扩散后灰度差变小。空域梯度函数利用图像邻域像素间的灰度变化的大小表示图像的清晰度,空域梯度函数又叫作基于图像锐度的评价函数。
2.1.1)Brenner函数
2.1.2)Tenengrad函数
2.1.3)能量梯度函数(Energy)
能量梯度函数类似于Tenengrad函数,不同的是能量梯度函数不是使用Sobel算子来计算图像的近似梯度,而是使用沿着扫描线的相邻像素之间的强度差来计算近似梯度。能量梯度函数能较好地解决图像噪声造成的问题,但也带来了大量的计算,增加算法用时。
空域梯度函数是公认的复杂度最低、操作简单的函数
2.2)基于信息统计的函数
2.2.1)方差函数(Variance)
M和N分别代表图像的长度和宽度。方差函数使用图像数据的方差为图像评估值。基于方差的方法是快速和强大的,基本思想是计算图像强度的方差,当方差达到最大时,图像准确聚焦。方差函数值取决于整幅图像灰度的概率分布而不是其空间分布。
2.2.2)熵函数
直方图显示图像整体的灰度分布情况,图像熵是基于灰度直方图计算的,反映像素灰度值的差异程度。
对于光学像来说,像素点灰度的分配是离散随机的,焦距不同得到图像的清晰度不同,图像的熵值也就不同。离焦时,像素点模糊造成图像灰度分布较均匀,相应的图像熵值也较大;正焦时,图像清晰,像素灰度分布差异大,图像熵值最小,我们将图像熵的相反数作为图像的评价值。
2.3)基于相关性度量的函数
2.3.1)自相关函数(Vollaths)
2.3.2)邻域相关函数(Correlation)
离焦时,由于图像点扩散的原理,像素在图像接收器上形成圆形斑点。由每个像素形成的圆形斑点彼此叠加并相互影响,这就导致像素与其相邻像素之间具有相关性,离焦度越高,相关性越大。
2.4)傅立叶频域函数
傅立叶变换将图像的空域灰度分布映射为傅立叶频域频率分布。在傅立叶频域中,频率的高低对应着原图像中灰度变化情况,若原图像中灰度分布没有太大变化,则反映到傅立叶频域里是低频分量;相应的原图像中灰度分布变化较大,反映到傅立叶频域里就是高频分量。基于傅立叶频域的评价函数通过计算图像的频谱得到:
R(u,v)和I(u,v)分别是傅里叶变换后的实部和虚部。当图像的尺寸较大时,要处理的数据量增大,傅里叶变换的计算量会降低***实时性。
3)图像清晰度评价指标
评价函数包括定性和定量两种性能评价指标。定性指标比较直观,可以第一时间判断出这条曲线是否适合用作聚焦评价函数;定量指标比较客观,当从定性指标上分辨不出区别时定量指标可以帮助我们更好的甄选评价函数。
3.1)定性评价指标
标准的评价函数曲线参见图20(a),曲线必须满足以下几点:
(1)单峰性。评价函数曲线必须具有单峰性,只存在一个最值,此类函数曲线在最值搜索时不会陷入局部最值的陷阱。
(2)无偏性。只有像平面与焦平面重合才能使评价函数达到最值,最值的尖锐结果重复性好。
(3)抗干扰性。即使受到一定噪声的干扰,或环境的各个参数(照明条件)发生变化,评价函数的功能也可以保证***正确检测,完成聚焦。
(4)适用性。评价函数要对不同类型的图像均适用,不能限于一些特殊类型。
3.2)定量评价指标
为了客观的衡量各种不同算法的性能,除了以上定性评价指标外,我们在实验中还使用了定量评价指标。两类评价指标结合,方便我们依据不同情况择取相应适合的算法。
3.2.1)清晰度比率(DefinitionRatio)
F max 是评价函数最大值,F min 是评价函数最小值,我们定义F max 与F min 的比值为清晰度比率Rat。一个标准的评价函数经过归一化处理后,F max 的值为1,F min 的值应接近于0 。Rat越大,代表正焦图像与离焦图像间的差别越大,聚焦的准确性就会有所保障,参见图20(a2)。
3.2.2)陡峭峰(Steepness)
根据某个范围内函数值的变化程度,我们将函数曲线分为平缓区和陡峭区。在平缓区内,曲线的函数值不会发生明显变化;在陡峭区内,焦距较小的调整也会引起曲线函数值较大的变化。平缓区分左、右两部分,两平缓区中间即是陡峭区,参见图20(a3)。
左平缓区与陡峭区的分界点为左临界点(p max , F max )。左临界点应符合公式要求如下:
从数学角度来说,曲线上某点的陡峭度等于该点处的函数微分值。陡峭区内,函数值变化量与陡峭区宽度的比值就近似为该区间的陡峭度,我们将左、右陡峭区的平均值定义为评价函数的陡峭度。好的陡峭度对后续的极值搜索至关重要,陡峭度较大的函数可以在平缓区使用较大步长,节省爬山时间,更好的找到最大值点。曲线的左、右陡峭度分别为:
曲线的平均陡峭度为:
3.2.2)灵敏度因子(SensitivityFactor)
灵敏度因子表征曲线在最大值附近的变化程度,灵敏度高代表曲线峰值附近的差异大;当灵敏度较小时,***最后很有可能会成像于正确焦点附近的虚假焦平面上,***受到噪声的干扰时,此类状况会加剧。图20(a4)是不同曲线的灵敏度示意图。F 1曲线和F 2曲线相比,在最大值左右两侧,横坐标改变大小为ε,函数值发生相应的变化,F 2曲线的改变值明显大于F 1曲线的改变值。灵敏度因子对聚焦过程的最后一步影响非常大,灵敏度因子的大小决定了峰值搜索时的准确度。
灵敏度定义如下:
p max 是函数最大值的位置,对应着函数的最大值F max ;最大值点偏移ε大小后位置为(p max +ε),对应的函数值为F(p max +ε)。将曲线左右两边灵敏度的平均值定义为灵敏度因子:
(4)平缓区波动(FlatFluctuation)
4)搜索策略
聚焦***的评价函数曲线部分为增、部分为减,近似对称,增减区间的转折点就是函数的峰值点,也就是函数的最大值点。爬山搜索是一种经典的极值查找算法,除此之外还有黄金分割搜索和曲线拟合法等。
4.1)爬山搜索
在峰值搜索的区间内任取两个点x 1和x 2(x 1与x 2间的距离小于搜索步长),若F(x 1)<F(x 2),则x 1到x 2为正确的搜索方向,沿此方向继续以预设的步长运行镜头寻找峰值。在运行过程中,当曲线函数值开始减小时,说明己经跨过峰值,此时应当改变搜索方向,并适当调小步长。接下来就是不断的重复以上步骤,当镜头运行到正确焦点位置搜索过程结束。爬山搜索通常使用自适应的搜索步长,首先用大步长“粗扫”整个聚焦区间,不断的缩小搜索范围,然后调整为小步长,“细扫”极值区间搜寻正确的聚焦位置。参见图20(a5)。
爬山搜索算法只能根据目前的坡向来判断山峰的方向,搜索一旦到达极值附近,对搜索方向的判断容易产生失误,会产生随机走动。同时,爬山搜索在峰值两侧有较大可能出现不断折返的情况,这大大地降低了算法效率
4.2)黄金分割搜索
假设搜索范围为[a,b],搜索精度要求为t,t>0。在黄金分割搜索中,分割系数T=0.618。选定两个点x 1=a+(1-T)(b-a)和x 2=b-(1-T)(b-a),若F(x 1)≤F(x 2),则搜索区间变为[x 1,b],令a=x 1,得到新的搜索区间[a,b];若F(x 1)>F(x 2),则搜索区间变为[a, x 2],令b=x 2,得到新的搜索区间[a,b]。不断重复以上步骤,当区间长度小于给定的搜索精度t,搜索停止,此时将搜索区间的中点x*=(a+b)/2作为函数最大值的解。
黄金分割搜索和Fibonacci搜索的原理一样,都是以二分查找为基础理论,但前者相对简单。黄金分割搜索和爬山搜索一样,在峰值的两侧容易发生摆动,浪费搜索资源。
4.3)曲线拟合搜索
直接根据曲线拟合法做出函数值的近似曲线,然后得到曲线的函数式,通过计算函数最大值点就可以一步到达焦点位置。一般搜索过程分为粗精度搜索和细精度搜索两个过程,曲线拟合多用于后段的细精度搜索过程。
曲线拟合是构建数学函数曲线的过程,函数曲线中导数为零的点即是函数的极值点,二次函数拟合就是一种较简单的拟合方式。此类算法要求曲线具备单峰性和对称性,这样才能得到好的拟合效果。聚焦函数的曲线远焦区较平缓,近焦区较陡峭,但是拟合曲线的收敛性不是很好,这就造成聚焦失败的可能性较大。曲线拟合法对拟合值的依赖性强,对数据波动的变化敏感,受曲线平滑度的影响较大。
选用改进的三点比较爬山搜索算法作为极值搜索策略。在函数曲线上选择一个初始点a 0。得到函数值F 0,接着按照***设置好的方向和步长移动两次到点a 1和a 2,得到评价函数值F 1和F 2,将三个函数值进行比较。若F 0> F 1> F 2,说明此时处于下坡状态或者已经越过极值点,极值点位于区间[a 0, a 1];若F 0<F 1< F 2,说明此时处于上坡状态或者己经越过极值点,极值点位于区间[a 1, a 2];若F 0<F 1且F 1>F 2或F 0>F 1且F 1<F 2,说明搜索有可能受到局部极值的干扰,改变步长重新取值再次作出判断。
如果自动聚焦成像失败,通过数字重聚焦方法进行景深延拓。数字重聚焦方法包括如下步骤:
步骤1:通过光机电***及其控制器获取阵列式高亮度光源的分布情况和待测物体的低分辨强度图像集;
步骤2:针对离焦待测物体计算原始强度图像的偏移量;
步骤3:将低分辨强度图像集的图像进行最大强度归一化处理;
步骤4:根据偏移量的大小分别进行偏移,最终合成一张低分辨率图像。
原始强度图像的偏移量如下:
其中(x j ,y j )是LED的位置坐标,H和h分别表示LED阵列式高亮度光源与样品的距离以及样品的离焦距离。
参见图22,传统科勒照明中4×/0.1NA物镜在红光照明下的景深为63µm,当观测样本的厚度超过物镜的景深时,传统的自动聚焦已经无法让样品处于实焦状态。对于切片机而且最多只能切100um即0.1mm,如果超过了这个值,就意味着不管怎么切,怎么制作样本都不容易发生失焦问题。FPM可以重构出复振幅光波场,从而能够通过上述数字重聚焦的方法计算出样品在一定深度范围内的实焦图,实现景深延拓。如图22,FPM在使用4×/0.1NA物镜时的景深延拓能力为300µm。
4、判别光学窗口
步骤三:点亮自适应强度可变角度照明器正中央高亮度光源单元,依次单通道正入射照明样本,成像探测器采集低分辨率图像并存储,合成低分辨率彩色图像;成像探测器还开启HDR功能获取高动态范围图像,并且依次通过步骤B:判别光学窗口,选取最佳单波长照明,通过步骤C:读取一维或二维条码信息。
所述的步骤B:判别光学窗口,选取最佳单波长照明,具体包括如下步骤:
1)获取样本染色的染色剂吸收曲线;
2)分析获取得到的染色剂吸收曲线是否是已知的染色剂吸收曲线;
3)当为染色剂吸收曲线为已知的染色剂吸收曲线时,通过其在已知光谱吸收曲线上的位置与红绿蓝三色光的光谱曲线进行对照,确定使用哪种单色光重构高分辨率灰度图像;
4)若染色剂的吸收曲线为未知,拍摄一张低分辨率彩色图片,得到红绿蓝三通道的灰度图片,将低分辨率彩色图片灰度化,通过计算红绿蓝三通道与灰度图的RMSE值,取RMSE值最小时对应的通道作为高分辨率重构色彩通道。
光学窗口的判别通常有两种情况,一种是已知染色剂吸收曲线,另一种是未知染色剂吸收曲线。若样本染色的染色剂吸收曲线已知,可以简单通过其在光谱吸收曲线上的位置与红绿蓝三色光的光谱曲线来确定使用哪种单色光重构高分辨率灰度图像,以便后续颜色传递重构,如图23(a)。若染色剂的吸收曲线未知,可拍摄一张低分辨率彩色图片,得到红绿蓝三通道的灰度图片,如图23(b2)-(b4),将低分辨率彩色图片灰度化,如图23(b1)。通过计算红绿蓝三通道与灰度图的RMSE值,取RMSE最小的通道作为高分辨率重构通道。
5、重构高通量灰度图像
参见图24,步骤五:处理器重构高通量灰度图像,具体包括如下步骤:
步骤1:读入初始化参数标定的光波矢位置及全视场像差值,构建光瞳约束条件和相位恢复图像重建模型;
步骤2:通过采集待测物体的低分辨强度图像,然后初始化待测物体的物函数,同时使用已标定好成像***的光波矢位置及光瞳函数;
步骤3:将经过步骤2初始化的物函数和标定的光瞳函数输入相位恢复图像重建模型,进一步求解光瞳约束值、图像的复振幅分布、光学成像***光瞳函数、光瞳约束优化条件,实现FPM高分辨率图像重建;其中,光学成像***的光瞳函数需要在每次迭代重构完毕后,用已知标定好的光瞳函数覆盖重建得到的光瞳函数及更新光瞳约束值;
步骤4:计算原始误差和对偶误差,采用原始误差衡量FPM高分辨率图像重建结果与待恢复图像之间的偏差;采用对偶误差评价光瞳约束值、图像的复振幅分布、光学成像***的光瞳函数、光瞳约束优化条件的结合与相位恢复图像重建模型得到的数值之间的偏差;
步骤5:设置原始误差和对偶误差判断容限,当误差容限满足设定要求时,输出图像重建结果,当误差容限不满足设定要求时,返回步骤3,直到输出满足误差设定要求的图像重建为止。
需要说明的,上述步骤1构建光瞳函数约束条件,引入光瞳函数约束中间变量为:
其中diag(·)是矩阵相乘运算的数学运算符,p表示光学***的光瞳函数,Q j 表示对应于第j个LED照明位置的采样矩阵,s表示待测物体;
相位恢复图像重建模型为:
其中“:=”是对“定义为”的符号表示,表示二维傅里叶逆变换,||·||2表示矩
阵的2-范数运算,q j 为光瞳函数约束中间变量,为待测物体的低分辨强度图像集,s表示
待测物体,N表示照明阵列中LED的数量,δ是不含任何相位信息的平面光场在傅立叶频域内
的表达,γ用来度量施加正则化项的强弱程度。
所述的步骤3中光瞳约束值求解方式如下:
图像的复振幅分布求解方式如下:
光学成像***光瞳函数求解方式如下:
光瞳约束优化条件的求解方式如下:
上述公式中,q j ,s,p,λ j 分别表示光瞳约束值、图像的复振幅分布、光学成像***的光瞳函数以及光瞳约束优化条件,携带的上标表明当前的迭代次数;“←”是赋值的符号表示,angle(·)表示求解复数相位的运算。
步骤3中光瞳约束值求解迭代过程如下:
图像的复振幅分布求解迭代过程:
光学成像***光瞳函数求解迭代过程:
光瞳约束优化条件的求解迭代过程:
其中和分别表示二维傅里叶变换及其逆变换,表示矩阵的复共轭;“←”
是赋值的符号表示,angle(·)表示求解复数相位的运算。Q j 和Q j H 分别表示对应于第j个LED
照明位置的采样矩阵及其逆矩阵,为待测物体的低分辨强度图像集,q j ,s,p,ω j 分别表
示光瞳约束值、图像的复振幅分布、光学成像***的光瞳函数以及光瞳约束优化条件,携带
的上标表明当前的迭代次数;N表示照明阵列中LED的数量,δ是不含任何相位信息的平面光
场在傅立叶频域内的表达,惩罚参数α用来描述光瞳约束中间变量定义的准确程度,参数β
的引入用来稳定光瞳函数的更新过程以防止极端值的出现,参数γ用来度量施加正则化项
的强弱程度,参数η是算法更新的步长。
步骤4,原始误差和对偶误差计算分别如下:
其中分别表示二维傅里叶逆变换,||·||2表示矩阵的2-范数运算;为待测
物体的低分辨强度图像集,q j ,λ j 分别表示光瞳约束值和光瞳约束优化条件;N表示照明阵列
中LED的数量,惩罚参数α用来描述光瞳约束中间变量定义的准确程度,参数η是算法更新的
步长,所有变量所携带的上标表明当前的迭代次数。
步骤5中设置原始误差和对偶误差判断容限条件为:
其中E p 和E d 分别表示原始误差和对偶误差,上标表明当前的迭代次数,ε tol 是预先设定好的误差容限值。
完成步骤2对待测物体的物函数以及成像***的光瞳函数进行初始化后,还包括对待测物体的物函数与光瞳函数防串扰干预,防串扰干预过程如下:
其中“:=”是对“定义为”的符号表示,diag(·)是矩阵相乘运算的数学运算符,
表示二维傅里叶逆变换,||·||2表示矩阵的2-范数运算;p表示光学***的光瞳函数,Q j 表
示对应于第j个LED照明位置的采样矩阵,为待测物体的低分辨强度图像集,s表示待测
物体,N表示照明阵列中LED的数量,δ是不含任何相位信息的平面光场在傅立叶频域内的表
达,γ用来度量施加正则化项的强弱程度。
参见图25,对洋葱鳞片叶表皮细胞切片分别进行20dB和30dB两组不同增益条件下的FPM成像实验。在数据采集过程中人为添加相机增益是为了模拟低成本低性能图像传感器的使用场景,从而获得噪声明显的低分辨率强度图像。实验结果表明,本发明使用的重构方法可以提供更加均匀平滑的背景信息和更好的重构对比度;其他两种重构方法的重建图像在样品的细胞壁周围聚集了较为明显的噪点,而本发明的重构方法则能够有效避免因噪声存在导致的图像重建质量的下降。
5、真彩色图像重构
现有的针对FPM的彩色化技术主要有四种。第一种,单色相机配合三色光源依次照明实现R+G+B三通道合成的全彩色FPM重构。与原始的单波长FPM方法相比,在交叠率和***环境方面没有额外要求,仅是采集和重构时间上翻了三倍且需要进行白平衡工作。第二种方法是直接使用彩色相机和三色光源同时照明,通过拜耳滤波片分离R/G/B三个主要通道。这种方法既省时又不受重叠率和***环境的额外限制,但是彩色相机的像元尺寸通常比黑白相机大,光子效率将非常低,这可能会导致再追求分辨率时采集更高角度信息的乏力,且利用这种方法时需要解决彩色泄露问题。第三种方法是波长复用FPM,称为WMFPM,仍然利用单色相机,但是同时采集三波长照明下的多角度图像,通过复杂的解耦合过程生成彩色高分辨率图像。从表面来看,这种方法所需的采集时间相较于传统R+G+B方法减少了2/3,然而潜在中这种方法所需的交叠率却是R+G+B彩色化方法的三倍,进而导致该方法在采集上并没有减少那么多时间,同时由于复杂的解耦合过程所增加的时间,总体上看效率并没有明显提升,且存在与现有***误差矫正算法难以嵌套的问题。第四种方法是近些年来大为火热的深度学习方法,然而这种方法存在过拟合、难泛化、大型训练集获取困难、缺乏物理原理、缺乏可解释性等一系列问题。
在具体实现彩色传递过程中,技术难点有两个,一是如何确保颜色的真实性及显示过程中的正确性,即保真度;二是如何在效率提升的过程中保证彩色传递的精度。针对第一点,建立了FPM技术的CIE-XYZ彩色空间及不同彩色空间的显示映射关系,使得所恢复的彩色结果在***示器、投影仪上都具有相同的保真度。
通过对不同染色剂染色的生物样本进行CFPM恢复发现CFPM方法具有“通道选择性”,也就是说,在已知染色剂种类的情况下可以找出最适合进行CFPM恢复的颜色通道。
参见图26-27,实现FPM彩色化灵感来源于颜色匹配的思想。所谓颜色匹配,就是将施主图像的色调信息传递到受主图像上,使得受主图像拥有与施主图像相同的“颜色风格”。我们提出的CFPM方案将低分辨率全视场彩色图像与高分辨率的FPM重构图像进行亮度直方图匹配,将低分辨率全彩色图像的颜色传递给高分辨率的FPM重构图像,使得染色后的高分辨率的彩色FPM图像最大程度的接近低分辨率全彩色图像。
通过对不同染色剂染色的生物样本进行CFPM恢复发现CFPM方法具有“通道选择性”,也就是说,在已知染色剂种类的情况下可以找出最适合进行CFPM恢复的颜色通道。
本发明的步骤六:处理器将低分辨率真彩色信息传递给高通量灰度图像,获得高通量真彩色图像,该方法具体包括如下步骤:(1.0版本)
步骤1:建立一个FPM颜色空间,确认sRGB、NTSC常规显示器所用色域和FPM颜色空间的相互映射关系,确保后续在***示设备上的颜色保真;
步骤2:将低分辨率真彩色图像和高通量灰度图像都从RGB颜色空间转换到Lab颜色空间,实现灰度与颜色的分离;
步骤3:将低分辨率真彩色图像的颜色信息迁移到高通量灰度图像上,得到高通量真彩色图像;
步骤4:将颜色迁移后的高通量真彩色图像的Lab颜色空间模式转换成RGB颜色空间,同时根据步骤1映射关系用于最终显示。
进一步地,步骤1:建立一个FPM颜色空间,用以实现sRGB颜色空间和FPM颜色空间的相互转换,具体包括:
1)FPM颜色空间的定量计算;
2)CIE-RGB颜色空间通常转换为CIE-XYZ颜色空间;
3)计算FPM色彩空间中的白平衡系数;
4)计算FPM和sRGB颜色空间的映射关系。
步骤2:将Lab颜色空间与RGB颜色空间的互相转换,将低分辨率彩色图像和高分辨率灰度图像都转换到Lab空间,具体包括:
1)将施主图像的sRGB颜色空间转换为Lab颜色空间;
2)将受体图像的亮度直方图与供体图像进行匹配;
3)将CFPM恢复的Lab彩色模型转换为sRGB彩色空间进行显示。
步骤3:将低分辨率彩色图像的颜色信息转移到高分辨率灰度图像上,具体包括如下步骤:
1)计算出供体图像每个像素的亮度值和邻域亮度方差值;
2)根据受体图像的像素值,结合邻域内的亮度值和亮度方差值,计算得到RA j 值;
3)用受体图像的RA j 值减去供体图像每个像素的RD值,得到最小的供体图像像素,即离受者图像的第j个像素最近的点,然后,将供体图像的这个像素的a和b通道值赋给受体图像的像素值;
对受体图像的每个像素重复上述1)—3)步骤,得到颜色信息转移。
参见图2,遵循上述技术逻辑,结合FPM的原理图和实验结构分别如图26(a)和(b)所示。一个32 × 32可编程的R/G/BLED阵列(Adafruit, 4mm间距,由Arduino控制)放置在样品上方70.0 mm处,而15 × 15中心LED元素依次点亮进行数据采集。通过OceanOptics光谱仪的测试,红色、绿色和蓝色led在波长为630.1、515.0和462.6 nm处有一个主导的窄峰,半峰宽(FWHM)分别为20.8、38.0和34.6 nm(图29(b1)-(b3))。采用紧凑的倒置显微镜,可以方便地进一步结合荧光成像(图14)。所有数据均由4×/0.1NA计划消色差物镜和16位sCMOS相机(Neo 5.5, Andor, 2160 × 2560像素,6.5µmpixelpitch)采集。
色彩空间的建立、白平衡以及对不同色彩空间的显示进行变换对于获得正确、准确的真彩图像具有重要意义。相应地,三个LED光源在1931CIE-XYZ标准图上校准为三基色(图29(a))。然后,得到我们的FPM***的颜色空间,如蓝色三角形区域所示。得到FPM颜色空间的定量计算如下:
其中I r (λ)、I g (λ)和I b (λ)为需要校准的三种原色的光谱曲线,图23(b1)-(b4)中的光谱曲线,r(λ)、g(λ)、b(λ)为CIE-RGB颜色空间的归一化匹配函数。配色函数会出现负值,不易操作;因此,CIE-RGB颜色空间通常转换为CIE-XYZ颜色空间,转换矩阵定义为:
其中CIE-XYZ空间仍为3D空间,不利于显示和计算。因此,XYZ 3D空间被投影到X +Y + Z = 1平面上(图29(a))。因此,我们的FPM-XYZ坐标值为(0.6625,0.2901,0.0474) R ,(0.2410, 0.4521, 0.3069) G , (0.1719, 0.0608, 0.7663) B 。另外,D65是标准白光光源,因其色温为6500 K而命名,是模拟日光的人工光源。通常将D65作为白平衡点,其RGB空间坐标值为(0.95047,1.00000,1.08883)。每个通道的白平衡系数可由下式得出:
W=γ 1 R+γ 2 G+γ 3 B(75)
式中γ 1、γ 2、γ 3为白平衡系数。因此,我们的FPM色彩空间中的白平衡系数为:
我们已经获得了FPM颜色空间及其显示的白平衡系数。然而,显示器或投影机通常使用sRGB颜色空间。因此,对于FPM-XYZ颜色空间中的特定颜色,其对应的sRGB颜色空间坐标值是:
式中γ’ 1、γ’ 2、γ’ 3为sRGB颜色空间的白平衡系数,(x R , y R , z R )为CIE-XYZ中sRGB颜色空间的坐标值。然后,将XYZ转换为sRGB的矩阵为:
矩阵包含负元素,也就是说,从FPM-XYZ空间到sRGB空间的转换可能使某个分量的值为负值。这意味着sRGB颜色空间中的所有颜色都应该位于其三角形内,如果一个特定的颜色落在三角形之外,则显示器或投影仪无法显示它。通常,一个简单的解决方案是强制范围外的值落在范围的边界上,强制负值等于0,强制大于1的值等于1。这种简单的方法相当于直接在三角形的边界上填充颜色,但是会失去一些颜色的准确性。如图29(a)的色域图所示,最终可以显示的颜色是FPM和sRGB颜色空间的重叠区域。对于其他颜色空间(如NTSC颜色空间)的显示的转换也是类似的。
具体来说,颜色转移方案包含三个步骤。首先,将施主图像的sRGB颜色空间转换为Lab颜色空间,其中L是代表亮度的通道,数值范围为0到100作为测量值,a代表红色和绿色通道,b代表黄色和蓝色通道,两者的颜色范围都是[+127, -128]。与sRGB色彩空间相比,Lab色彩空间将亮度和色彩纹理完全分离为两个通道,不会相互干扰,更便于色彩的传递。第二步通过将受体图像的亮度直方图与供体图像进行匹配,将供体的颜色纹理信息转移到Lab颜色模型内的受体图像中,以消除它们之间的整体差异。最后,将CFPM恢复的Lab彩色模型转换为sRGB彩色空间进行显示。
将sRGB颜色模型转换为Lab颜色模型的具体步骤如下。
将sRGB颜色空间转换为人眼三个视锥反应所代表的LMS颜色空间,公式如下:
式中,L为长波长通道,M为中波长通道,S为短波波长通道。
为了使数据在LMS颜色空间中更聚敛,更符合人类对颜色感觉的心理物理学的规律,需要进一步将数据转化为以10为底的对数空间。
对数据进行主成分分析,得到LMS颜色模型到Lab颜色模型的转换参数如下:
从Lab颜色空间到sRGB颜色空间的逆变换原理与上面相同。在对供体和受体图像进行预处理后,对颜色转移进行了如下详细的描述。供体和受体图像的视场是相同的,但分辨率却不相同,这导致了上采样放大倍数的差异。
步骤1 计算出供体图像每个像素的亮度值PD和5×5邻域亮度方差值DD,定义参数RD=0.5PD+ 0.5DD。
步骤2 对于受体图像的像素j,计算5×5邻域内的亮度值PA j 和亮度方差值DA j ,得到RA j ,其中RA j = 0.5PA j +0.5DA j 。
步骤3 用受体图像的RA j 减去供体图像每个像素的RD值,得到最小的供体图像像素,即离受者图像的第j个像素最近的点。然后,将供体图像的这个像素的a和b通道值赋给受体图像的像素j。
步骤4 对受体图像的每个像素重复相同的操作,就可以得到颜色转移结果。
受到双边滤波(bilateralfiltering,BF)和等价局部二值模式(uniformlocalbinarypatterns,ULBP)的启发,本发明进一步提出了一种基于三边滤波的颜色转移方法,即2.0版本,将空间信息引入到FPM颜色转移算法中,在保证图像恢复质量的同时加快了速度,色彩丰富度也得到了显著的提高。
具体将1.0版本的步骤3:将低分辨率彩色图像的颜色信息转移到高分辨率灰度图像上,替换方案如下:
1)将低分辨率彩色图像和高分辨率灰度图像分成若干小图块,分别计算各个块对应的供体图像和受体图像的ULBP特征图;
2)针对各个块对应的供体图像和受体图像的ULBP特征图,确定滤波器大小和滑动步长参数;根据滤波器的和滑动步长的大小,将颜色转移后的若干小图块拼接成一个完整的HR彩色图像;
3)计算HR彩色图像的G通道与高分辨率灰度图像的RMSE值,当RMSE值的相对波动范围小于阈值时,执行步骤4,否则,将HR彩色图像G通道替换成高分辨率灰度图像,然后转化成新的高分辨率灰度图像,将其作为新的受体图像重复执行步骤1)-步骤3);
4)将HR彩色图像的G通道替换成高分辨率灰度图像,得到最终的高分辨率HR彩色图像。
具体来说,本文提出的着色方法是利用滤波器的思想对大尺度图像进行处理,在Lab空间中传递颜色信息,最终生成高分辨率(highresolution,HR)彩色图像。这解决了CFPM计算图像中所有像素的方差的缺陷,避免了生成的HR彩色图像的直方图过于平衡带来的多染色问题。算法中采用高斯滤波因子引入空间信息,使算法更具说服力,符合客观先验信息的事实,提高了图像信息的利用率。颜色转移方案包含三个阶段,整个算法流程如图28所示。CFFPM的具体流程步骤如下:
步骤1 对于大尺寸的低分辨率(lowresolution,LR)彩色图像和HR灰度图像,计算它们各自的ULBP特征图。CFFPM在BF中采用了滤波器和滑动步长的思想,将图像分成若干个具有一定重叠率的小块,便于并行运算减少颜色转移所需的时间。计算供体图像和受体图像的ULBP特征图,选择滤波器大小和滑动步长等参数。在图像中的一个中心像素及其周围的3×3矩形区域上,将该区域内每个像素的灰度值与中心像素进行比较,然后进行二值化。大于或等于中心像素的像素值编码为 1,小于为 0。将得到的二值化值乘以各自位置的权重2 p 并累加,则得到中心像素的LBP编码值。
原始LBP算子计算方法如下:
其中(x c ,y c )为中心像素,i c 为其灰度值,P为除中心像素外的像素总数,第P个像素的灰度值为i p ,s(x)为符号函数。
为了进一步增强LBP的判别能力,将矩形区域修改为圆形区域。圆形区域的半径为r,其上均匀分布有P个像素。通过双线性插值得到圆周上未被准确覆盖的邻域点的灰度值:
从上面可以看出,在二进制计算中,LBP将生成2 p 结果,这将导致计算量随着采样点数量p的增加而几何增加,并且还将导致识别率的降低。有一种模式在图像LBP特征提取中经常出现。这种模式是等价局部二进制模式(ULBP)。这种高频率的LBP值代表了图像的最基本特征。在二进制中,ULBP值的0/1的转换数小于或等于2。其计算公式如下:
在用ULBP分类方法对LBP进行分类时,模式类型从原来的2 p 减少到p(p-1)+2,而不会丢失太多信息。这使得特征向量的维数更小,并减少了高频噪声的影响。
步骤2 双边滤波(BF)是Tomasi等人在1998年提出的基于高斯滤波的改进算法。BF同时考虑了空间信息和灰度信息,以去除噪声并保留大量边缘细节。滤波函数的计算公式如下:
其中I p 是输入图像在p(x,y)点的灰度值,q表示中心像素点,(u,v)表示中心像素点q的坐标,p是领域内像素点,S是领域像素的集合,W p 是归一化因子,G s 是灰度相似因子,G r 是空间相似因子,σ s 和σ r 是基于高斯函数的距离标准差和灰度标准差。
针对1.0算法缺乏空间信息,2.0算法在Lab色彩空间中吸收BF的灰度相似因子和空间相似因子思想,并在其中引入ULBP,得到如下三边滤波函数:
式中,l p 是输入低分辨率彩色图像在Lab颜色空间中点p(x,y)的亮度值,U p 是输入
低分辨率彩色图像的LBP特征图在点p(x,y)的灰度值,l q 是输入高分辨率灰度图像在Lab颜
色空间中点q(u,v)的亮度值,U q 是输入高分辨率灰度图像的LBP特征图在点q(u,v)的灰度
值,是颜色传递概率函数,点(m,n)是概率函数取值最大处。G s 是空间邻近度因子,G r 是
灰度相似因子,G t 是特征相似度因子,σ s 是基于高斯函数的距离标准差,σ r 是基于高斯函数
的灰度标准差,σ t 是基于高斯函数的特征标准差。
接下来,将用滤波器分割后得到的图片从RGB空间转换到Lab空间,然后利用上述描述的颜色转移函数进行颜色转移,最后再将HR图片颜色转移后的Lab空间转化为RGB空间。
步骤3 根据滤波器的和滑动步长的大小,将颜色转移后的几个小块拼接成一个完整的HR彩色图像。计算得到的高分辨率彩色图像的G通道与最初(第一次)的高分辨率灰度图像的RMSE值。当RMSE的相对波动范围小于10%,执行步骤4,否则将高分辨率彩色图像的G通道替换成高分辨率灰度图像,然后转化成新的高分辨率灰度图像,将其作为新的受体图像(替代最初的高分辨率灰度图像接受低分辨率彩色图像的颜色信息)继续执行步骤1-3。
步骤4 将高分辨率彩色图像的G通道替换成高分辨率灰度图像(最初的),得到最终的高分辨率彩色图像。
6、分析器分析
步骤六:处理器将低分辨率真彩色信息传递给高通量灰度图像,获得高通量真彩色图像,之后还包括步骤D:通过分析器分析真彩色图像,判别有无病变,给出初筛置信度;若有病变,则输出病变区域,给出病变阶段及阶段置信度,并通过步骤七显示器输出带分析结果的高通量真彩色图像。
步骤七:显示器输出分析后的高通量真彩色图像之后,显示器还通过通讯器将分析得到有病变的全部高通量真彩色图像实现屏幕共享或远程医疗。
分析器通过人工智能深度学习的方法,进行图像病理分析,标记出图像病理区域与置信度,具体而言,采用MaskRCNN模型网络结构实现图像分类、定位和分割。
基于MaskRCNN模型网络结构数据分析方法具体如下:
MaskRCNN是用于完成目标分类、检测和语义分割等多种任务的一种2阶段算法,其在FasterRCNN的基础上,添加了一个用于预测目标掩膜的分支,即利用全卷积神经网络对RCNN建议的每一个感兴趣的区域完成语义分割,从而实现同时进行分类、定位和分割。
参见图37,MaskRCNN采用了2阶段网络结构。第1阶段:首先利用骨干网络即深度残差网络(deepresidualnetwork,DResNet),从输入图像中提取不同阶段的特征图;其次利用特征金字塔网络自上而下及横向连接结构融合不同尺度的特征,使其同时具有强语义信息和强空间信息;再次对这些特征图上每一点像素设定固定数量的锚框,通过计算每个锚框与该图片上标注的真实框之间的交,获得多个大小不同的候选区域(regionofinterest,ROI);最后利用区域建议网络对候选的ROI进行二值分类(即前景/背景)及边框回归,过滤掉分类分数低的ROI,并将正负样本的比设定为1:3,以缓解类别不平衡问题,同时减少了第2 阶段对不必要信息的计算。第2阶段首先进行2次对齐操作:(1)对第1阶段选出的ROI进行对齐,把原图像中ROI与特征图中的像素对应起来;(2)把特征图上大小不同的ROI转换成统一大小. 其次为了减少由池化过程带来的误差,通过双线性插值法从特征图上相邻网格点计算每个像素值,获取ROI包含的重要特征信息,完成分类、回归和分割任务. 最后通过在全连接层上增加分割分支,对每个ROI上的每个像素进行分类和回归预测,并输出最终的二值掩膜。
1)特征金字塔网络FPN
在卷积神经网络生成的特征图中,低层特征图尺寸大,分辨率高,深层特征图更抽象,利于分类。FasterRCNN模型属于单尺度的目标检测模型,即它只使用了CNN最后一个卷积层的特征图。现假设输入图像是 224 × 224 ,通过VGGNet之后得到的特征图的大小是7 × 7,可以看出图像尺寸缩小了非常大的倍数,所以图像当中的位置信息也会因此损失很多,如果原图当中存在小目标,则可能会因为损失目标信息而无法检测。所以FPN的提出就是用来解决该问题的,它可以融合高层和低层的特征图信息,属于多尺度的目标检测。特征金字塔网络FPN结构参见图30。
原始图像首先经过五个卷积层分别生成C1-C5的特征图,C5 横向连接使用1×1的卷积核来降低通道数,使之与C4通道数相匹配,会得到M5,对M5 进行2 倍的上采样,使之与C4尺寸相同,再将C4 横向连接1×1 卷积核再与上采样之后的M5相加得到M4,接着将M4进行2 倍上采样与C3横向连接并相加得到M3,以此类推,生成图中的M2-M5,对每一个M系列的特征图还需要做一个3×3 的卷积操作得到对应的P系列,3×3 的卷积可以消除相加运算带来的混叠效应,最后FPN真正的输出是P2-P5 特征图。这种做法可以在提取裂缝特征时从多个尺度上进行融合,有利于更好的学习裂缝特征。
FPN输出的多个不同尺度的特征图会送入到RPN网络中产生ROI感兴趣区域,这些ROI会通过公式到指定的特征图上进行切割。其中,k代表第k个P层特征图,w、h代表特征图的宽和高,当k 0=4 且w×h=224×224 时,即k=4,代表在P4层特征图上进行切割。从公式可以看出,特征图尺寸越大,会选择更高的P层,特征图尺寸越小,会选择低层的P层。
2)生成候选区域的RPN
2.1)区域建议网络RPN
区域建议网络RPN的作用是生成高质量的区域候选框。RPN结构参见图31所示,网络的输入是经过CNN输出的特征图,图左侧的红色方框代表的是一个卷积核,经过滑动计算之后则会生成一张新的特征图。每个卷积核的中心都会生成k个anchors(锚点),每个锚点包含k种不同尺寸{128 × 128,256 × 256,512 × 512} 和不同横纵比 {1: 1,1: 2,2:1} 的矩形框,即经过卷积之后特征图上的每个点都会映射回原图上的k个候选区域,每个候选框会有两部分输出,分别是分类输出和位置回归输出。
对生成的anchor进行规定,正样本为:(1)与真实窗口有最高IoU重叠率,(2)与真实窗口的IoU>0.7。负样本为:与真实窗口的IoU<0.3。中性样本为:与所有真实窗口的IoU在0.3 到0.7 之间。舍弃中性样本的anchor,这些anchor将不会用于后续网络的训练。
2.2)边界回归
目标检测中需要确定每个目标所存在的具***置,这个位置信息被称为BoundingBox,也叫BBox。即使是正样本的anchor也不能完全覆盖目标,因此,RPN加入了边界回归应用于锚点的细化,通过移动anchor框使其调整到包含目标的正确边界。边界回归是对区域建议进行修正的一种线性回归算法,它的目的是为了让区域建议提取到的窗口与目标窗口更加吻合。
如图32所示,通常采用(x, y, w, h) 来表示一个窗口,图中的框P代表原始候选
框,即提取出来的区域,框G代表目标的真实所在位置,边界回归要做的就是使框P (P x , P y
, P w , P h )通过某种映射f,移到与框G最为相近的回归窗口框,使得。
具体做法为先平移,再尺度放缩:
先平移(Δx,Δy)
再做尺度缩放(S w ,S h )
输入为(P x , P y , P w , P h ),经过平移变换和尺度缩放d x (P), d y (P), d w (P), d h
(P),最终输出的结果是预测窗口,并非真实窗口G。下述公式表示的是经过框P和框G得到
的平移量(t x , t y )和尺度缩放(t w , t h )。
定义损失函数为下第一个式子,函数优化目标为下第二个式子其中,t ∗ = (t x , t y , t w , t h ),代表真实坐标。这里,是输入候选区域的特征向量,w*是需要学习的
参数,是得到的预测结果。
MaskRCNN模型中的Boundingbox回归公式为式其中,x,y,w和h代表这个框的中心坐标x-y以及宽度和高度。变量x、x a 和x ∗ 分别为预测窗口、anchor窗口和真实窗口。
2.3)非极大值抑制(NMS)
通常,经过RPN之后会生成大量的候选框,这些候选框之间存在很多冗余的部分,非极大值抑制NMS(Non-maxSuppression)算法就是用来解决该问题的,通过删除多余的候选框以减少检测框的数量。NMS具体思路如图33所示:
1)排序图中所有候选框的置信度得分,如图33(b),将置信度最高的框标为红色;
2)遍历图33(b)中剩下的蓝框,保留与图中红色框的重叠面积小于给定阈值的那些蓝框,其余去除,得到图33(c)所示的结果;
3)从未处理的其余蓝框中继续选一个得分最高的,如图33(d),重复上述过程;
4)最后保留下的框即为最终的预测框(如图33(e))。
3)ROIAlign操作
3.1)双线性插值
双线性插值的原理如图34,现有Q 11,Q 12,Q 21,Q 22四个点,需要求出P点坐标,则可以通过做三次线性插值来实现。首先由Q 12和Q22这两个点可以得到点R 2,再由Q 11和Q 21这两个点可以得到点R 1,最后根据R 1和R 2这两个点就可以得到点P,公式如下式所示。
3.2)ROIAlign操作
CNN产生的特征图经过RPN网络之后会得到一系列的候选框,这些候选框尺寸大小不一,ROIAlign的作用就是为了使这些候选框的尺寸进行统一。ROIAlign层的提出是为了解决FasterRCNN中ROIPooling层导致的空间位置错位问题。ROIPooling的具体做法是,如图下所示,原始图像经过CNN之后假设生成了一个7×7 大小尺寸的特征图,现有一个候选区域ROI的尺寸为5×4,希望固定尺寸大小为2×2,则需要对该ROI划分格子,但由于5/2=2.5,并不是整数,所以只能按3 和2 进行格子的划分,然后取每个格子中的平均值或最大值进行池化,使用ROIPooling这种方法会导致格子划分的大小不一样。并且由于网络生成的ROI尺寸并不一定是整数,如图35(a)所示,红色方框ROI并没有和原特征图的网格边界对应,则首先需要把小数部分略去,只取整数,得到图35(b),使得红色框与特征图网格的边界对应,那么这一步就会损失一些位置信息。之后在进行池化变成2×2 大小输出时,由于格子划分的大小不同,也会损失一部分信息。下图很好的说明了经过两次小数取整带来的误差和对不齐的现象。CNN原始生成的边界框是图35(a),经过第一次小数取整变成了图35(b),再经过图35(c)不均匀的池化即第二次小数取整,使得像素与边界框之间无法达到一一对应,即每一次量化操作都会产生区域特征的错位,这种偏差对目标分类也许影响不大,但对目标的分割则会产生非常大的误差影响,从而直接导致精细目标的分割结果不准确现象。
在语义分割这种精细程度很高的任务中,由于需要对每个像素点进行精细的预测,如果使用ROIPooling,则会损失很多信息并且产生错位现象,从而导致分割的不精准,因此提出了ROIAlign方法。通过保留小数信息实现真正像素点之间的对应。
如图35(a)所示,一张 7 × 8 的特征图,绿色框是生成的ROI,假设需要得到一个2×2 大小的输出,首先直接对这个边界框进行均等的划分,划分为4 个小格子,然后对每个小格子再等分为4 个小区域,对每个小区域取它的中心点,这个点的值是由它周边特征图上距离它最近的4 个真实像素点的值经过双线性插值得到的,图35(b)为其中一个中心点采用双线性插值计算的示意图。同理可以得到每个小格子的中心点的值,然后取这4 个点值的最大值,即最大池化。当然也可以不进行二次划分小格子,直接取第一次绿线划分的格子里中心点的值,这个取决于具体的应用,在本文中进行了二次划分格子。从图中可以看出,整个过程均保留了小数信息,避免了ROIPooling层两次整数化带来的空间错位现象。
4)Mask掩膜分支
如上文可知,原始裂缝图像经过ResNet-FPN会输出多个尺度的特征图,这些特征图经过RPN可以得到多个ROI,这些ROI会通过在相应的特征图上切割产生候选区域,然后进行ROIAlign使得这些不同大小的候选区域固定统一尺寸,最后网络分成三路,分别实现目标的分类、检测以及分割任务。与FasterRCNN不同的是,MaskRCNN增加了一个Mask掩膜分支部分,通过连接全卷积神经网络来对特征图中的每个像素点进行分类,从而实现像素级的目标分割工作。
Mask分支的结构没有用到全连接层。现假设预测一个28×28 的掩膜mask,这个mask中的每个点都会对应一个N维的向量(N代表有N个类别),即每个类别都会生成一个28×28 像素的掩膜,由浮点数表示,因此比二进制掩膜可以保存更多细节,这种小尺寸的mask也有助于Mask分支网络的轻量化。传统对像素进行分类的方法采用的都是Softmax,其做法是对每个类别都输出一个打分,属于该类别的分数高一些,其他类别的分数相对低一些,这样就会存在一个不同类别之间的类间竞争问题。Mask分支则采用像素级Sigmoid方法进行预测,每个像素会对应一个预测的向量,向量维度是80,即80 个类别,这个像素属于哪个类,则相应维度就会得到一个比较高的分值。这里需注意的是,不管目标尺寸有多大,在训练期间,都会将真实的mask标签缩小到28×28 大小来计算损失,在测试期间,则会将预测的mask根据ROI的大小进行放大,并给出最终预测的mask,使得不同目标有不同mask的大小。
5)多任务损失Loss
MaskRCNN算法采用多任务的训练方式,在实现目标分类的基础上可以同时进行定位和分割。其对图像的多任务损失函数为下式,它是由五种损失函数相加得到的:
其中,L rc 表示RPN网络分类损失;L rb 表示RPN网络边界回归损失;L cls 表示目标分类损失;L box 表示边界框回归损失;L mask 表示掩膜分割结果的损失。
5.1)分类损失函数L cls 和L rc
对于目标分类任务使用下式计算损失,其中p表示该候选框判定为分类目标的概率;y表示是否属于该类别的二进制参数,为正标签时y=1,为负标签时y=0;M表示类别个数。
5.2)回归损失函数L box 和L rb
对于目标边界框回归任务采用L1 范数,并在smoothL1 函数的基础上加入范围在(-1,+1)的分段函数以保证损失函数在0 点附近的平滑性。L box 和L rb 损失函数的公式均如下式所示,其中,t i = (t x , t y , t w , t h ),表示预测框的坐标,t * i 表示真实窗口的坐标。
5.3)分割损失函数L mask
MaskRCNN算法使用的是逐像素Sigmoid来计算损失,假设有N个类别,对于一个类
别为n的ROI,L mask 代表的是在第n个mask上输出的损失(其余n-1 个mask输出对整个Loss损
失无贡献),如下式所示,其中,y ij 表示a×a掩膜区域内的像素(i,j)标签,代表在相同位
置上第n层mask上像素的预测值。
参见图38,使用基于MaskRCNN模型结构的深度学习技术对宫颈病理切片的成像结果进行分析,并对病变等级分类。该技术能够精确定位与识别病理切片中的病变位置,同时根据组织细胞的形态特征将病变划分为以下四个阶段:0、病变开始(正常或亚正常);1、潜在低恶性(低度鳞状内皮病变);2、潜在高恶性(低度鳞状内皮病变);3、扩散性癌症(扩散性鳞状癌)。
参见图39,本发明中批量样品的检测与分析遵循一定的并行处理时序逻辑。对于一个单独的样品而言,采集19×19张低分辨率图像序列的时间约为4s,包含相机的触发与曝光、LED切换、图像储存等。采集完成后,对该样品依次进行图像重构和图像分析操作,耗时分别为120s和60s左右;与此同时,仪器立刻并行对后续样本进行数据采集,大致在完成第46个样本的采集之后,第一个样本完成分析。当第200个样本完成采集后,可更换样本集开始新一轮的批量样本检测与分析。
Claims (42)
1.一种全划片高通量彩色病理成像分析仪器,其特征在于,该仪器至少包括:自适应强度可变角度照明器及控制器,用于实现不同的时间,以多个入射角度照明样本,将各个入射角度光源到达样本平面亮度达到一致,同时形成高亮度光源;
所述的自适应强度可变角度照明器控制器,其控制方法如下:
1)选择大于阵列式高亮度光源的NA,确保成像探测器拍摄任意角度照明下均为明场图像;
2)固定曝光时间,取灰度平均值作为每一角度抵达样本平面的亮度值;
3)以正入射光源为基准,计算亮度测量系数矩阵;
4)根据亮度测量系数矩阵反馈控制每个角度的到达样本平面亮度,使其均匀一致;所述的反馈控制通过基于亮度测量系数矩阵调节单独光源的电功率;
所述的亮度测量系数矩阵定义为:
光机电***及其控制器,至少包括一个偏振片、一个显微物镜、一个扫描位移台,其中偏振片设置在高亮度光源与样本之间,偏振片结合高亮度光源实现快速曝光和抑制杂散光,进一步通过显微物镜收集经过扫描样本发出的光,再利用扫描位移台实现扫描多个样本;所述的多个入射角和显微物镜的傅立叶频域中的多个重叠区域相对应;
成像探测器,用于获取光机电***及其控制器输出的多个强度图像,每个强度图像分别对应自适应强度可变角度照明器生成的各个入射角度;
处理器,至少用于实现重建比探测器输出的强度图像分辨率更高、图像精度更高的可重复图像,定量恢复样本相位信息。
2.如权利要求1所述的全划片高通量彩色病理成像分析仪器,其特征在于,所述的处理器还用于实现分析仪器的初始化参数标定。
3.如权利要求1或2所述的全划片高通量彩色病理成像分析仪器,其特征在于,所述的处理器还用于实现自动聚焦成像。
4.如权利要求3所述的全划片高通量彩色病理成像分析仪器,其特征在于:所述的处理器在实现自动聚焦成像失败时,还能通过数字重聚焦方法进行景深延拓,所述的景深延拓长度大于0.1mm。
5.如权利要求3所述的全划片高通量彩色病理成像分析仪器,其特征在于,所述的处理器还用于实现快速高精度全彩色成像。
6.如权利要求1所述的全划片高通量彩色病理成像分析仪器,其特征在于,所述的自适应强度可变角度照明器的光源由多个高亮度独立光源单元间隔排列,形成阵列式高亮度光源,从而无须机械扫描实现各个高亮度独立光源单元照明角度切换。
7.如权利要求6所述的全划片高通量彩色病理成像分析仪器,其特征在于,所述的阵列式高亮度光源排布形式为平板矩阵列排布、平板圆环阵列排布、曲面半球体阵列排布或异形变体阵列排布。
8.如权利要求6或7所述的全划片高通量彩色病理成像分析仪器,其特征在于,所述的阵列式高亮度光源为R/G/B三波长光源或者R/G/B/W四波长光源,其中R为红光,中心波长范围为600-740 nm;G为绿光,中心波长范围为492-590 nm;B为蓝光,中心波长范围为400-490nm;W为白光,中心波长范围为540-580 nm;不同波长的峰值半高宽范围为20-40 nm。
9.如权利要求8所述的全划片高通量彩色病理成像分析仪器,其特征在于,所述的自适应强度可变角度照明器的光源为LED、激光二极管、OLED或LCD;所述的自适应强度可变角度照明器的高亮度独立光源单元每个通道刷新帧率为20-1000Hz,电功率为0.2-3W。
10.如权利要求7所述的全划片高通量彩色病理成像分析仪器,其特征在于,所述的阵列式高亮度光源排布形式为平板矩阵列排布或平板圆环阵列排布时,独立光源单元间距d为3 mm~6 mm;阵列式光源距样本的高度h为45-80 mm;高亮度独立光源单元为平板矩阵列排布时,阵列排布数量范围为3×3~64×64;高亮度独立光源单元为平板圆环阵列排布时,阵列排布数量范围为4~1000。
12.如权利要求11所述的全划片高通量彩色病理成像分析仪器,其特征在于,所述的阵列式高亮度光源的傅立叶频域交叠率计算方式如下:
1)当阵列式排布光源的各个高亮度独立光源单元为均匀分布时,该光源的傅立叶频域交叠率计算如下:
其中S t 是NA的增加步长,NA obj 为物镜的数值孔径;
2)当阵列式高亮度光源的各个高亮度独立光源单元为非均匀分布时,该光源的傅立叶频域交叠率通过加权平均计算获得如下:
其中S t,i 是第i圈高亮度独立光源单元NA的增加步长,NA obj 为物镜的数值孔径,m为高亮度独立光源单元总数,n为第i圈高亮度独立光源单元的个数;
所述的阵列式高亮度光源的傅立叶频域交叠率为40%-85%。
14.如权利要求13所述的全划片高通量彩色病理成像分析仪器,其特征在于:所述的自适应强度可变角度照明器,还包括与处理器相连的自适应强度可变角度照明器控制器,用于根据成像探测器的反馈控制每个高亮度独立光源单元,使得阵列式高亮度光源到达样品平面亮度自适应地保持一致。
15.如权利要求14所述的全划片高通量彩色病理成像分析仪器,其特征在于:所述的自适应强度可变角度照明器控制器,还能够作为信号发生器输出上升沿或下降沿信号,用于硬触发成像探测器完成自动序列采集图像,所输出的上升沿或下降沿电压范围为3v-12v。
16.如权利要求1所述的全划片高通量彩色病理成像分析仪器,其特征在于:所述的扫描位移台为单轴或双轴x-y轴电动载物台及其控制器,用于实现大尺寸样本的扫描拼接重构和多个样本扫描;其大尺寸样本为10mm×10mm视场范围以上,多个样本为2个及以上。
17.如权利要求1或16所述的全划片高通量彩色病理成像分析仪器,其特征在于:所述的光机电***及其控制器还包括样本及样本存储器、z轴电动平移台及控制器、筒透镜;其中样本存储器用于存储大批量待测或已测的样本,样本及样本存储器通过z轴电动平移台实现交互;z轴电动平移台及控制器,通过自动聚焦成像方法获取具体离焦量后调整z轴位置至实焦面实现自动聚焦;
自适应强度可变角度照明器照明样本,经过偏振片和高亮度光源抑制杂散光后到达样本平面,显微物镜收集经过扫描样本发出的光,由筒透镜反馈至成像探测器,然后通过z轴电动平移台及控制器实现自动聚焦,借助单轴或双轴x-y轴电动载物台及其控制器实现大尺寸样本的扫描拼接重构和通过与样本存储器交互实现样本批量扫描。
18.如权利要求17所述的全划片高通量彩色病理成像分析仪器,其特征在于:所述的处理器还分别与单轴或双轴x-y轴电动载物台及其控制器、显微物镜控制器、z轴电动平移台控制器相连,显微物镜控制器与显微物镜相连。
19.如权利要求18所述的全划片高通量彩色病理成像分析仪器,其特征在于:所述的光机电***及其控制器内,在样品与探测器之间还设置有一个偏振片,用于实现偏振成像,从而实现多模态成像。
20.如权利要求19所述的全划片高通量彩色病理成像分析仪器,其特征在于:所述的偏振成像包括如下步骤:
21.如权利要求19所述的全划片高通量彩色病理成像分析仪器,其特征在于:所述的筒透镜还通过显微镜架连接有目镜。
22.如权利要求21所述的全划片高通量彩色病理成像分析仪器,其特征在于:所述的显微物镜还通过滤光片连接卤素灯或紫外灯,用于实现非相干成像和荧光多模态成像。
23.如权利要求1所述的全划片高通量彩色病理成像分析仪器,其特征在于:所述的成像探测器为CCD或CMOS,成像探测器与处理器相连。
25.如权利要求5所述的全划片高通量彩色病理成像分析仪器,其特征在于:所述的处理器上还连接有存储器、分析器和显示器,构成通用仪器;所述的存储器用于实现数据存储及压缩处理,分析器利用人工智能识别、分类图像中的病理信息,输出图像中的病理区域与置信度,显示器实现图像显示。
26.如权利要求25所述的全划片高通量彩色病理成像分析仪器,其特征在于:所述的显示器还通过连接通信器实现屏幕共享或远程医疗。
27.一种利用权利要求1的成像分析仪器实现的全划片高通量彩色病理成像分析方法,其特征在于:该方法包括如下步骤:
步骤一:成像分析仪器初始化参数标定;
步骤二:装填待测样本;
步骤三:点亮自适应强度可变角度照明器正中央高亮度光源单元,依次单通道正入射照明样本,成像探测器采集低分辨率图像并存储,合成低分辨率彩色图像;
步骤四:自适应强度可变角度照明器自动采用特定波长扫描样本,通过成像探测器采集序列低分辨率图像并存储,直至完成待测样本采集;
步骤五:处理器重构高通量灰度图像;
步骤六:处理器将低分辨率真彩色信息传递给高通量灰度图像,获得高通量真彩色图像;
步骤七:显示器输出高通量真彩色图像。
28.如权利要求27所述的全划片高通量彩色病理成像分析方法,其特征在于:所述的步骤一:成像分析仪器初始化参数标定,是通过***误差矫正算法矫正分析仪器***光波矢误差和像差影响,完成分析仪器初始化参数标定;***误差矫正算法是在经典高斯牛顿重构算法基础上,依次添加了基于模拟退火、非线性回归波矢矫正及全局像差约束的三个矫正模块予以实现。
29.如权利要求28所述的全划片高通量彩色病理成像分析方法,其特征在于:所述的***误差矫正算法包括如下步骤:
步骤1:计算每个高亮度独立光源单元对应的傅立叶频域位置信息,通过成像探测器捕获一组低分辨率强度图像,带入高斯牛顿重构算法,根据重构结果判断是否存在光波矢误差,若存在则启用***误差矫正算法;若不存在,则认为照明波矢位置准确,记录照明波矢位置,用于设定后续每一待测样本重构算法中波矢位置,同时将高斯牛顿重构算法输出的全视场数字化像差保存,用于后续每一待测样本重构算法中的数字化像差补偿;
步骤2:计算阵列式排布光源的各个高亮度独立光源单元对应的频谱信息;
步骤3:在阵列式高亮度光源照明条件下,向多个方向随机进行傅立叶频域步长搜索,更新待测目标频谱信息;
步骤4:通过频谱信息误差矩阵,更新每个高亮度独立光源单元对应的频谱信息;
步骤5:借助非线性规划计算阵列式高亮度光源整体偏差参数;
步骤6:通过自适应步长多次迭代,实现全局像差约束,得到阵列式高亮度光源全视场像差值。
30.如权利要求28所述的全划片高通量彩色病理成像分析方法,其特征在于:所述的步骤二:装填待测样本之后,还包括步骤A自动聚焦待测样本;如果自动聚焦成像失败,通过数字重聚焦方法进行景深延拓。
31.如权利要求30所述的全划片高通量彩色病理成像分析方法,其特征在于:所述的步骤A自动聚焦待测样本,具体通过聚焦窗口的选择、图像清晰度评价、搜索策略确定实现自动聚焦成像。
32.如权利要求30所述的全划片高通量彩色病理成像分析方法,其特征在于:所述的数字重聚焦方法包括如下步骤:
步骤1:通过光机电***及其控制器获取阵列式高亮度光源的分布情况和待测物体的低分辨强度图像集;
步骤2:针对离焦待测物体计算原始强度图像的偏移量;
步骤3:将低分辨强度图像集的图像进行最大强度归一化处理;
步骤4:根据偏移量的大小分别进行偏移,最终合成一张低分辨率图像。
33.如权利要求27所述的全划片高通量彩色病理成像分析方法,其特征在于:所述的步骤三:点亮自适应强度可变角度照明器正中央高亮度光源单元,依次单通道正入射照明样本,成像探测器采集低分辨率图像并存储,合成低分辨率彩色图像;成像探测器还开启HDR功能获取高动态范围图像,并且依次通过步骤B:判别光学窗口,选取最佳单波长照明,通过步骤C:读取一维或二维条码信息。
34.如权利要求33所述的全划片高通量彩色病理成像分析方法,其特征在于:所述的步骤B:判别光学窗口,选取最佳单波长照明,具体包括如下步骤:
1)获取样本染色的染色剂吸收曲线;
2)分析获取得到的染色剂吸收曲线是否是已知的染色剂吸收曲线;
3)当染色剂吸收曲线为已知的染色剂吸收曲线时,通过其在已知光谱吸收曲线上的位置与红绿蓝三色光的光谱曲线进行对照,确定使用哪种单色光重构高分辨率灰度图像;
4)若染色剂的吸收曲线为未知,拍摄一张低分辨率彩色图片,得到红绿蓝三通道的灰度图片,将低分辨率彩色图片灰度化,通过计算红绿蓝三通道与灰度图的RMSE值,取RMSE值最小时对应的通道即作为高分辨率重构色彩通道。
35.如权利要求27所述的全划片高通量彩色病理成像分析方法,其特征在于:所述的步骤五:处理器重构高通量灰度图像,具体包括如下步骤:
步骤1:读入初始化参数标定的光波矢位置及全视场像差值,构建光瞳约束条件和相位恢复图像重建模型;
步骤2:通过采集待测物体的低分辨强度图像,然后初始化待测物体的物函数,同时使用已标定好成像***的光波矢位置及光瞳函数;
步骤3:将经过步骤2初始化的物函数和标定的光瞳函数输入相位恢复图像重建模型,进一步求解光瞳约束值、图像的复振幅分布、光学成像***光瞳函数、光瞳约束优化条件,实现FPM高分辨率图像重建;其中,光学成像***的光瞳函数需要在每次迭代重构完毕后,用已知标定好的光瞳函数覆盖重建得到的光瞳函数及更新光瞳约束值;
步骤4:计算原始误差和对偶误差,采用原始误差衡量FPM高分辨率图像重建结果与待恢复图像之间的偏差;采用对偶误差评价光瞳约束值、图像的复振幅分布、光学成像***的光瞳函数、光瞳约束优化条件的结合与相位恢复图像重建模型得到的数值之间的偏差;
步骤5:设置原始误差和对偶误差判断容限,当误差容限满足设定要求时,输出图像重建结果,当误差容限不满足设定要求时,返回步骤3,直到输出满足误差设定要求的图像重建为止。
36.如权利要求27所述的全划片高通量彩色病理成像分析方法,其特征在于:所述的步骤六:处理器将低分辨率真彩色信息传递给高通量灰度图像,获得高通量真彩色图像。
37.如权利要求36所述的全划片高通量彩色病理成像分析方法,其特征在于:所述的步骤七显示器输出高通量真彩色图像。
38.如权利要求27所述的全划片高通量彩色病理成像分析方法,其特征在于:所述的步骤六:处理器将低分辨率真彩色信息传递给高通量灰度图像,获得高通量真彩色图像,该方法具体包括如下步骤:
步骤1:建立一个FPM颜色空间,确认sRGB、NTSC常规显示器所用色域和FPM颜色空间的相互映射关系,确保后续在***示设备上的颜色保真;
步骤2:将低分辨率真彩色图像和高通量灰度图像都从RGB颜色空间转换到Lab颜色空间,实现灰度与颜色的分离;
步骤3:将低分辨率真彩色图像的颜色信息迁移到高通量灰度图像上,得到高通量真彩色图像;
步骤4:将颜色迁移后的高通量真彩色图像的Lab颜色空间模式转换成RGB颜色空间,同时根据步骤1映射关系用于最终显示。
39.如权利要求38所述的全划片高通量彩色病理成像分析方法,其特征在于:所述的步骤1:建立一个FPM颜色空间,确认sRGB、NTSC常规显示器所用色域和FPM颜色空间的相互映射关系,确保后续在***示设备上的颜色保真,具体包括:
1)FPM颜色空间的定量计算;
2)CIE-RGB颜色空间通常转换为CIE-XYZ颜色空间;
3)计算FPM色彩空间中的白平衡系数;
4)计算FPM和sRGB颜色空间的映射关系。
40.如权利要求38所述的全划片高通量彩色病理成像分析方法,其特征在于:所述的步骤2:将低分辨率真彩色图像和高通量灰度图像都从RGB颜色空间转换到Lab颜色空间,实现灰度与颜色的分离,具体包括:
1)将施主图像的sRGB颜色空间转换为Lab颜色空间;
2)将受体图像的亮度直方图与供体图像进行匹配;
3)将CFPM恢复的Lab彩色模型转换为sRGB彩色空间进行显示。
41.如权利要求38所述的全划片高通量彩色病理成像分析方法,其特征在于:所述的步骤3:将低分辨率真彩色图像的颜色信息迁移到高通量灰度图像上,得到高通量真彩色图像,具体包括如下步骤:
1)计算出供体图像每个像素的亮度值和邻域亮度方差值;
2)根据受体图像的像素值,结合邻域内的亮度值和亮度方差值,计算得到RAj值;
3)用受体图像的RA j 值减去供体图像每个像素的RD值,得到最小的供体图像像素,即离受者图像的第j个像素最近的点,然后,将供体图像的这个像素的a和b通道值赋给受体图像的像素值;
对受体图像的每个像素重复上述1)—3)步骤,得到颜色信息转移。
42.如权利要求41所述的全划片高通量彩色病理成像分析方法,其特征在于:所述的步骤3:将低分辨率真彩色图像的颜色信息迁移到高通量灰度图像上,得到高通量真彩色图像,替换方案如下:
1)将低分辨率彩色图像和高分辨率灰度图像分成若干小图块,分别计算各个块对应的供体图像和受体图像的ULBP特征图;
2)针对各个块对应的供体图像和受体图像的ULBP特征图,确定滤波器大小和滑动步长参数;根据滤波器和滑动步长的大小,将颜色转移后的若干小图块拼接成一个完整的HR彩色图像;
3)计算HR彩色图像的G通道与高分辨率灰度图像的RMSE值,当RMSE值的相对波动范围小于阈值时,执行步骤4,否则,将HR彩色图像G通道替换成高分辨率灰度图像,然后转化成新的高分辨率灰度图像,将其作为新的受体图像重复执行步骤1)-步骤3);
4)将HR彩色图像的G通道替换成高分辨率灰度图像,得到最终的高分辨率HR彩色图像。
Priority Applications (1)
Application Number | Priority Date | Filing Date | Title |
---|---|---|---|
CN202210961097.7A CN115032196B (zh) | 2022-08-11 | 2022-08-11 | 一种全划片高通量彩色病理成像分析仪器及方法 |
Applications Claiming Priority (1)
Application Number | Priority Date | Filing Date | Title |
---|---|---|---|
CN202210961097.7A CN115032196B (zh) | 2022-08-11 | 2022-08-11 | 一种全划片高通量彩色病理成像分析仪器及方法 |
Publications (2)
Publication Number | Publication Date |
---|---|
CN115032196A CN115032196A (zh) | 2022-09-09 |
CN115032196B true CN115032196B (zh) | 2022-12-13 |
Family
ID=83130881
Family Applications (1)
Application Number | Title | Priority Date | Filing Date |
---|---|---|---|
CN202210961097.7A Active CN115032196B (zh) | 2022-08-11 | 2022-08-11 | 一种全划片高通量彩色病理成像分析仪器及方法 |
Country Status (1)
Country | Link |
---|---|
CN (1) | CN115032196B (zh) |
Families Citing this family (3)
Publication number | Priority date | Publication date | Assignee | Title |
---|---|---|---|---|
CN117058014B (zh) * | 2023-07-14 | 2024-03-29 | 北京透彻未来科技有限公司 | 一种基于lab色彩空间匹配的染色归一化***及方法 |
CN116609942B (zh) * | 2023-07-18 | 2023-09-22 | 长春理工大学 | 一种分孔径压缩感知偏振超分辨率成像方法 |
CN117456544B (zh) * | 2023-12-25 | 2024-03-15 | 苏州可帮基因科技有限公司 | 带笔迹病理图像的修复方法及设备 |
Family Cites Families (9)
Publication number | Priority date | Publication date | Assignee | Title |
---|---|---|---|---|
US20060093809A1 (en) * | 2004-10-29 | 2006-05-04 | Hebrink Timothy J | Optical bodies and methods for making optical bodies |
CN110082900B (zh) * | 2013-08-22 | 2022-05-13 | 加州理工学院 | 可变照明傅立叶重叠关联成像设备、***以及方法 |
CN105976315B (zh) * | 2016-04-26 | 2019-04-09 | 清华大学深圳研究生院 | 基于部分傅里叶空间的显微成像方法 |
CN108388014A (zh) * | 2018-05-15 | 2018-08-10 | 萤欧(上海)汽车科技有限公司 | 一种亮度自适应调节ar hud照明光学*** |
CN110773736B (zh) * | 2018-05-18 | 2022-05-13 | Ii-Vi特拉华有限公司 | 利用光纤阵列激光源和自适应多光束整形的金属中的增材制造 |
CN109581643B (zh) * | 2018-11-28 | 2020-11-17 | 中国科学院西安光学精密机械研究所 | 傅里叶叠层显微成像装置及方法 |
CN110308549B (zh) * | 2019-08-05 | 2024-01-19 | 常丽 | 一种整玻片多光路高通量病理切片扫描仪 |
EP4031923A1 (en) * | 2019-09-19 | 2022-07-27 | Siemens Healthcare Diagnostics, Inc. | Ptychographic imaging system and method for generating images |
CN111929881A (zh) * | 2020-08-06 | 2020-11-13 | 东南大学江北新区创新研究院 | 一种基于色散的相位物体成像设备及方法 |
-
2022
- 2022-08-11 CN CN202210961097.7A patent/CN115032196B/zh active Active
Also Published As
Publication number | Publication date |
---|---|
CN115032196A (zh) | 2022-09-09 |
Similar Documents
Publication | Publication Date | Title |
---|---|---|
CN115032196B (zh) | 一种全划片高通量彩色病理成像分析仪器及方法 | |
JP6915349B2 (ja) | 画像処理装置、画像処理方法、及び画像処理プログラム | |
US9117273B2 (en) | Method and apparatus for single-particle localization using wavelet analysis | |
CA2460801C (en) | Method for quantitative video-microscopy and associated system and computer software program product | |
CN105378538B (zh) | 用于多频谱成像的自动聚焦方法和*** | |
US8743195B2 (en) | Whole slide fluorescence scanner | |
JP6307770B2 (ja) | 撮像素子 | |
Chieco et al. | Image cytometry: protocols for 2D and 3D quantification in microscopic images | |
US20220334371A1 (en) | Intelligent automated imaging system | |
US11776124B1 (en) | Transforming multispectral images to enhanced resolution images enabled by machine learning | |
US8611620B2 (en) | Advanced digital pathology and provisions for remote diagnostics | |
MX2007014016A (es) | Metodos de analisis de imagenes basados en la separacion del cromogen. | |
WO2005114147A1 (ja) | 着色された光透過性の試料を撮像して得られた画像を処理する画像処理装置 | |
JP2000508095A (ja) | 境界マッピングのシステムおよび方法 | |
US9046476B2 (en) | Method and system for the detections of biological objects | |
JP2020169995A (ja) | バイオマーカの同時、多重化、超高感度およびハイスループットの光学的検出のためのバイオセンサプラットフォームおよび方法 | |
WO2021195817A1 (zh) | 一种提取待测物质的光谱信息的方法 | |
CN110785709A (zh) | 从低分辨率图像产生高分辨率图像以用于半导体应用 | |
Ma et al. | Light-field tomographic fluorescence lifetime imaging microscopy | |
Tran et al. | Mobile Fluorescence Imaging and Protein Crystal Recognition | |
WO2023189393A1 (ja) | 生体試料観察システム、情報処理装置及び画像生成方法 | |
KR102562740B1 (ko) | 이미지 정규화를 이용한 기생충 충란 식별 방법 및 장치 | |
Wu | Deep Learning-Enabled Computational Imaging in Optical Microscopy and Air Quality Monitoring | |
Sun et al. | Artificial Intelligence Techniques for Uncovering Resolved Planetary Nebula Candidates from Wide-field VPHAS+ Survey Data | |
CN117546007A (zh) | 信息处理装置、生物样本观察***及图像生成方法 |
Legal Events
Date | Code | Title | Description |
---|---|---|---|
PB01 | Publication | ||
PB01 | Publication | ||
SE01 | Entry into force of request for substantive examination | ||
SE01 | Entry into force of request for substantive examination | ||
GR01 | Patent grant | ||
GR01 | Patent grant |