CN110717882A - 图像处理方法、计算机程序以及记录介质 - Google Patents

图像处理方法、计算机程序以及记录介质 Download PDF

Info

Publication number
CN110717882A
CN110717882A CN201910568602.XA CN201910568602A CN110717882A CN 110717882 A CN110717882 A CN 110717882A CN 201910568602 A CN201910568602 A CN 201910568602A CN 110717882 A CN110717882 A CN 110717882A
Authority
CN
China
Prior art keywords
image
block
processing method
original image
image processing
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.)
Pending
Application number
CN201910568602.XA
Other languages
English (en)
Inventor
藤本博己
Current Assignee (The listed assignees may be inaccurate. Google has not performed a legal analysis and makes no representation or warranty as to the accuracy of the list.)
Screen Holdings Co Ltd
Original Assignee
Screen Holdings Co Ltd
Priority date (The priority date 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 date listed.)
Filing date
Publication date
Application filed by Screen Holdings Co Ltd filed Critical Screen Holdings Co Ltd
Publication of CN110717882A publication Critical patent/CN110717882A/zh
Pending legal-status Critical Current

Links

Images

Classifications

    • 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
    • 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
    • G06T7/00—Image analysis
    • G06T7/10—Segmentation; Edge detection
    • G06T7/11—Region-based segmentation
    • G—PHYSICS
    • G06—COMPUTING; CALCULATING OR COUNTING
    • G06T—IMAGE DATA PROCESSING OR GENERATION, IN GENERAL
    • G06T7/00—Image analysis
    • G06T7/10—Segmentation; Edge detection
    • G06T7/194—Segmentation; Edge detection involving foreground-background segmentation
    • 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/62—Analysis of geometric attributes of area, perimeter, diameter or volume
    • H—ELECTRICITY
    • H04—ELECTRIC COMMUNICATION TECHNIQUE
    • H04N—PICTORIAL COMMUNICATION, e.g. TELEVISION
    • H04N19/00—Methods or arrangements for coding, decoding, compressing or decompressing digital video signals
    • H04N19/10—Methods or arrangements for coding, decoding, compressing or decompressing digital video signals using adaptive coding
    • H04N19/169—Methods or arrangements for coding, decoding, compressing or decompressing digital video signals using adaptive coding characterised by the coding unit, i.e. the structural portion or semantic portion of the video signal being the object or the subject of the adaptive coding
    • H04N19/17—Methods or arrangements for coding, decoding, compressing or decompressing digital video signals using adaptive coding characterised by the coding unit, i.e. the structural portion or semantic portion of the video signal being the object or the subject of the adaptive coding the unit being an image region, e.g. an object
    • H04N19/176—Methods or arrangements for coding, decoding, compressing or decompressing digital video signals using adaptive coding characterised by the coding unit, i.e. the structural portion or semantic portion of the video signal being the object or the subject of the adaptive coding the unit being an image region, e.g. an object the region being a block, e.g. a macroblock
    • 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/10—Image acquisition modality
    • G06T2207/10056—Microscopic image
    • 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/20—Special algorithmic details
    • G06T2207/20021—Dividing image into blocks, subimages or windows
    • 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/20—Special algorithmic details
    • G06T2207/20092—Interactive image processing based on input by user
    • G06T2207/20104—Interactive definition of region of interest [ROI]
    • 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/30004—Biomedical image processing
    • G06T2207/30024—Cell structures in vitro; Tissue sections in vitro

Landscapes

  • Engineering & Computer Science (AREA)
  • Physics & Mathematics (AREA)
  • General Physics & Mathematics (AREA)
  • Theoretical Computer Science (AREA)
  • Computer Vision & Pattern Recognition (AREA)
  • Data Mining & Analysis (AREA)
  • Mathematical Physics (AREA)
  • Databases & Information Systems (AREA)
  • General Engineering & Computer Science (AREA)
  • Software Systems (AREA)
  • Multimedia (AREA)
  • Pure & Applied Mathematics (AREA)
  • Signal Processing (AREA)
  • Mathematical Optimization (AREA)
  • Algebra (AREA)
  • Computational Mathematics (AREA)
  • Mathematical Analysis (AREA)
  • Geometry (AREA)
  • Quality & Reliability (AREA)
  • Radiology & Medical Imaging (AREA)
  • Nuclear Medicine, Radiotherapy & Molecular Imaging (AREA)
  • Medical Informatics (AREA)
  • General Health & Medical Sciences (AREA)
  • Health & Medical Sciences (AREA)
  • Image Analysis (AREA)
  • Image Processing (AREA)
  • Apparatus Associated With Microorganisms And Enzymes (AREA)
  • Measuring Or Testing Involving Enzymes Or Micro-Organisms (AREA)

Abstract

根据图像的特征将包含细胞的图像适当地分割为多种区域。为了达到该目的，本发明的图像处理方法包括如下工序：获取包含有培养的细胞的原始图像(步骤S101)；将原始图像分割为由规定数量的像素构成的块，求出各块内的图像所具有的空间频率分量(步骤S104～S106)；将每个块分类为多维特征量空间中的多个集群，该多维特征量空间是将空间频率分量中的直流分量的强度和空间频率彼此不同的多个交流分量的各自的强度作为特征量(步骤S107、S108)；以及基于分类结果将原始图像分割为多个区域(步骤S109)。

Description

图像处理方法、计算机程序以及记录介质
技术领域
本发明涉及一种图像处理技术，用于从以包括培养的细胞及其背景的方式拍摄的图像中分割细胞区域和除该区域之外的区域。
背景技术
在以包括培养的细胞的方式拍摄的图像中，有时在图像内将作为本来的拍摄对象物的细胞的区域与作为其背景的区域进行区别处理。例如为了求出细胞的大小，需要划定细胞与背景的边界。以往，例如通过使用规定的药剂对细胞进行染色或者投入荧光试剂来进行细胞区域的确定。然而，例如以当观察细胞随时间的变化等为目的时，有时需要不使用像改变作为拍摄对象的细胞这样的处理而可以确定细胞区域。
培养基中的细胞接近透明，折射率也与水几乎没有差别。因此，在拍摄的图像中，细胞与背景的边界不清楚的情况也很多。为了解决该问题，本申请的发明人目前为止提出了一种在明亮视野下拍摄的技术，尤其是能够获取清晰地呈现细胞轮廓的图像(例如，参考日本特开2016-071588号公报)。为了区分图像中的细胞与背景，能够可以使用例如使用了细胞与背景之间的亮度差的二值化处理、用于提取细胞轮廓的边缘检测等图像处理。
发明内容
在随着时间变化观察培养环境下的细胞的情况下，在培养的初始阶段，在拍摄的图像中包含的是细胞和作为大致相同背景的培养基。因此，对细胞与背景之间的识别比较容易。然而，随着培养的行进，从细胞排出的废物或死亡的细胞残骸等(碎片；debris)出现在作为拍摄对象的细胞的周围，并包含在图像中。因此，在二值化或轮廓提取等简单的处理中，难以适当地识别细胞与碎片。
另外，例如在使用了ES细胞、IPS细胞等干细胞的培养中，有时干细胞与饲养细胞(feeder cell)一起培养。另外，在伴随着培养干细胞分化为体细胞的过程中，出现具有与干细胞不同形态特征的细胞。因此，在拍摄的图像中混合存在多个种类的细胞，并且需要在图像内识别他们。
由此，需要一种能够根据例如在图像中出现的纹理(texture)等图像特征的差异，从而将图像分割成多种、尤其是三种以上的区域的技术。
本发明是鉴于上述问题而提出的，其目的在于提供一种能够根据图像的特征将包含细胞的图像适当地分割为多种区域的技术。
在本发明的图像处理方法的一个方面，为了达到上述的目的，具有如下的工序：获取包含有培养的细胞的原始图像；将所述原始图像分割为由规定数量的像素构成的块，求出各块内的图像所具有的空间频率分量；将所述各块分类为多维特征量空间中的多个集群(cluster)，所述多维特征量空间是将所述空间频率分量中的直流分量的强度和空间频率彼此不同的多个交流分量中的每一个的强度作为特征量；以及基于所述分类的结果将所述原始图像分割为多个区域。
在这种结构的发明中，基于图像内的局部空间频率分量的分布分割图像。即，将图像分割为块，求出每块的空间频率分量。空间频率分量中的直流分量表示该块中的平均亮度。另一方面，交流分量表示块内的亮度的变化，比较低的频率分量表示缓慢的亮度变化，较高的频率分量表示更细微的亮度变化的状态。因此，在这样根据空间频率分量的强度来表现块的情况下，分别由直流分量反映该块的亮度，或者由交流分量反映块内的纹理。尤其是，频率越高，表示越精细的构造纹理。
如果将这些各频率分量的强度认为是该块的特征量，则能够通过这些特征量以定量的方式表现块内的图像内容的亮度或纹理的复杂度。在图像中，例如即使在包含多种细胞的情况下，其亮度或纹理的差异也由特征量的差异表现。
利用上述情况，能够根据图像内容将原始图像分割为多个区域。具体来说，认为在将各块所具有的直流分量以及多个交流分量的强度作为特征量的多维特征量空间中，各块占据的位置。这样，在符合细胞的块与除其以外的其他的块之间、或者对应于彼此不同种类的细胞的块之间，其分布出现显著的差别。因此，如果基于该特征量对作为原始图像的一部分的各块进行分类，则具有相同的纹理的块在特征量空间中形成一个集群。认为在特征量空间接近的多个块彼此包含相同种类的细胞。另一方面，当两个块在特征量空间中位于较远的位置时，认为这些块包含不同种类的细胞。这样，利用基于特征量的分类结果，能够根据该图像内容将原始图像分割为多个区域。
本发明的图像处理方法能够将计算机作为其执行主体。在该意义上来说，本发明能够作为用于使计算机执行上述处理的计算机程序来实现。另外，也能够作为非临时性地存储该计算机程序的存储介质来实现。
如上所述，根据本发明的图像处理方法，针对将原始图像进行分割而得到的每个块，基于该块内的图像所具有的空间频率分量的特征，能将原始图像分类为多个区域。通过分别将直流分量的强度和频率不同的多个交流分量中的每一个的强度作为特征量，从而能够根据块整体的亮度或纹理的差异进行适当的分类。
附图说明
图1是示出本实施方式的区域分割处理的流程图。
图2是执行该实施方式的区域分割处理的计算机装置的结构例。
图3A至图3C是示出原始图像的例子的图。
图4是示出DCT运算的概念的图。
图5A至图5C是示出DCT系数的分组的例子的图。
图6是示出用于制定判定基准的处理的流程图。
图7A以及图7B是例示处理过程中的指示输入方式的图。
图8A以及图8B是示意性地示出特征量空间中的集群分布的例子的图。
附图标记说明
1 计算机装置
2 存储光盘(存储介质)
B 背景
BL 块
C1～C5 细胞
Ia 原始图像
具体实施方式
下面，对本发明的图像处理方法的一实施方式进行说明。该实施方式中的图像处理是区域分割处理，该区域分割处理用于对在培养基中培养的细胞或者细胞的集合体进行拍摄而得到图像。具体来说，将细胞的区域与其他的区域进行分割，另外，在包含多种细胞的情况下将与这些对应的区域彼此分割。该图像处理可用于处理由拍摄装置拍摄的未加工的原始图像，该拍摄装置具有对包含细胞的试料进行拍摄的拍摄功能。拍摄装置可以执行该图像处理作为对拍摄后的图像数据的后续处理之一。另外，从拍摄装置或者适当的存储单元接收到图像数据的计算机装置也可以执行该图像处理。
其中，以通用计算机装置对已经执行拍摄生成的原始图像数据执行本实施方式的图像处理的方式为例进行说明。拍摄装置的结构并没有特别的限定，只要具有对将包含培养的细胞的试料与培养基一起进行拍摄，并输出为数字图像数据的功能即可。另外，优选的是，图像是在明亮视野下拍摄的图像。另外，对于计算机装置，例如能够使用具有产品化为个人计算机的通常的硬件结构的计算机。以下，省略对这些硬件的详细说明。
图1是示出本实施方式的区域分割处理的流程图。首先，参考图1对本实施方式中的区域分割处理的概要进行说明。然后，对处理的具体的内容进行详细地说明。计算机装置首先获取原始图像(步骤S101)。原始图像是通过用适当的拍摄装置以及拍摄条件在明亮视野下对包含培养基中培养的细胞的试料进行拍摄而获得。原始图像是包含至少一个细胞整体的图像。因此，原始图像可以包含细胞和细胞以外的碎片、背景等。
细胞可以包含形态特征不同的多个种类。例如在干细胞的培养中共同培养干细胞和饲养细胞的环境下，这些可能均包含在原始图像中。另外，在干细胞的一部分分化为体细胞的情况下，可能处于干细胞与体细胞混合存在的状态。
针对获取的原始图像，判断是否需要进行分辨率变换(步骤S102)。分辨率变换是将拍摄时的原始图像的像素尺寸变换成适合于区域分割处理的像素尺寸的处理。在判断需要的情况下，执行分辨率变换(步骤S103)。以下，将没有进行分辨率变换的拍摄时的图像、变换后的图像均称为“原始图像”。
在下面的描述中，当简称为“像素”而无特别地描述时，是指在步骤S104以后的处理中使用的图像的像素。执行分辨率变换的变换后的图像的像素或者如果没有执行分辨率变换的拍摄时的原始图像的像素在这里称为“像素”。需要说明的是，当需要特别区分时，有时将变换前的(即拍摄时的)像素称为“原始像素”。即，“原始像素”与是否执行分辨率变换无关，是指拍摄时的像素。
原始图像被分割为多个块(步骤S104)。在该实施方式中，将原始图像中的4个像素×4个像素的正方形视为一个块。当进行分辨率变换时，如上所述，这里所称的像素为变换后的像素。在一行中的像素数量不能被一块的一边的像素数量4整除的情况下，在原始图像的端部可以产生未被计入处理的像素。另外，也可以预先去除原始图像的端部以免出现这样的小数。
接下来，从各像素的像素值减去原始图像整体中的平均像素值作为偏移值(步骤S105)。为了便于以后的运算，该处理是用于电平位移(level shift)的处理。当像素值表示为例如8位数据时，各像素可以获取从0到255的256级的值。然而，通过减去平均像素值，也会出现获取负值作为像素值。
针对这样的方式分割的每个块，求出空间频率分量。为了达到该目的，在本实施方式中使用二维离散余弦变换(Discrete Cosine Transform；以下简称为“DCT”)运算(步骤S106)。例如，如在JPEG(Joint Photographic Experts Group)图像压缩技术中使用的那样，通过二维的DCT运算，能将各块的图像内容分解为二维空间频率分量。
该块的图像内容能够通过对各频率分量进行合成而再现。根据图像的内容，每个频率分量的强度都不同。如在JPEG方式的图像压缩理论中已知的那样，求出一个块中的各频率分量的强度作为通过DCT运算获得的DCT系数。
空间频率分量中的低频分量表示块内的图像中的缓慢的亮度变化。其中的直流分量相当于该块的亮度的平均值。另一方面，高频分量表示更加细微的亮度变化，并具有与块内的图像的纹理相关的信息。频率越高表示越细微的纹理。因此，基于求出的DCT系数，计算与各频率分量的强度对应的特征量(步骤S107)。具体来说，将与直流分量对应的DCT系数作为一个特征量，并且根据与DCT系数对应的频率的高度将交流分量分为几个组，将组内的DCT系数的合计作为一个特征量。分组的具体的方法稍后将描述。
这样，能够通过表示亮度的直流分量特征量和表示纹理的多个交流分量特征量来表现一个块。根据这些特征量所表示的块是较强地具有作为任一种类的细胞的特征，或者具有细胞以外的非细胞的特征，能够判断该块对应于哪个种类的细胞，或者对应于除此以外的非细胞(步骤S108)。基于该结果，将原始图像分割为多个区域，例如分割为细胞的区域与非细胞的区域，或者根据细胞种类分割区域(步骤S109)。以上是本实施方式的区域分割处理的概要。
图2是执行本实施方式的区域分割处理的计算机装置的结构例。计算机装置1例如具有作为个人计算机的通常的结构，包括CPU(Central Processing Unit)10、内存14、存储器15、输入设备16、显示部17、接口18以及磁盘驱动器(Disk Driver)19等。
CPU10执行预先准备的控制程序。这样，以软件的方式实现用于执行区域分割处理的各功能块，即，图像加工部11、DCT运算部12以及分类部13等功能块。需要说明的是，这些可以由专用硬件构成。内存14临时性地存储在CPU10的运算过程中生成的各种数据。除了CPU10应该执行的控制程序之外，存储器15还长期地存储原始图像的图像数据或处理后的图像数据等。
输入设备16是用于接受来自操作人员的指示输入的设备，例如包含鼠标、键盘等。另外，显示部17具有显示图像的功能，例如为液晶显示器，显示原始图像或处理后的图像，并向操作人员显示消息等各种信息。需要说明的是，输入设备和显示部可以设置为一体化的触摸面板。
接口18经由电通信线路与外部装置之间进行各种数据的交换。磁盘驱动器19接收外部的记录光盘2，该外部的记录光盘2记录有图像数据或控制程序等各种数据。存储在记录光盘2中的图像数据或控制程序等通过磁盘驱动器19读取，并存储在存储器15中。磁盘驱动器19可以具有将在计算机装置1内生成的数据写入到记录光盘2中的功能。
对于用于使计算机装置1执行本实施方式的区域分割处理的控制程序，磁盘驱动器19也可以是访问并读取非暂时性地记录有该控制程序的记录光盘2的方式。另外，也可以是经由接口18从外部装置提供控制程序的方式。对于原始图像数据也相同。
图像加工部11执行原始图像的分辨率变换、块分割以及最终的图像的区域分割等。另外，DCT运算部12对每个块执行上述的DCT运算并求出空间频率分量。分类部13基于针对每个块求出的特征量，来对各块进行分类。
图3A至图3C是示出原始图像的例子的图。如图3A所示，在步骤S101中获取的原始图像Ia作为图像对象包含细胞C1～C5以及与培养基对应的背景B。其中，细胞C1与细胞C4是相同的种类，其内部具有相同的纹理。同样地，细胞C2与细胞C3是相同的种类。细胞C1、C2、C5是彼此为不同的种类。在图中，细胞内部的浓度的差异表示种类的差异。
在这些图中，将图像的横向坐标设为X坐标，并将纵向坐标设为Y坐标。这些图像对象的形状以及浓淡仅表示一个例子，在实际获得的图像中可以获取各种形状以及浓淡。需要说明的是，除了这些细胞C1～C5与背景B之外，碎片也可以包含在原始图像中。另外，也可以是图像被细胞C1～C5以及围绕他们周围的碎片所占据而不包含与培养基对应的背景。另外，由于培养基的量或拍摄时的照明条件等，背景B的亮度也并不总是相同的。
因此，应该认为图像中的细胞以及背景的配置或浓淡等的关系因每个原始图像不同而不同。这里描述的图像处理无论这些关系如何均成立。在该实施方式中的区域分割处理的基本想法是：将原始图像Ia分割为较小的块，根据各块具有的空间频率分量，来判断该块是较强地具有任一种类的细胞的图像特征，还是较强地具有除此以外的特征。
在该区域分割处理中，将原始图像Ia分割为具有规定尺寸的正方形的块而进行处理。如图3B所示，一个块BL是由在X方向以及Y方向上分别4个像素，总计16个像素P组成的正方形的区域。优选的是，块BL的大小是根据作为图像中的提取对象的细胞的大小而设定。
如果块尺寸过小，则一个块只包含细胞的极小一部分。因此，无法能使图像中的细胞所具有的特征性质充分地反映在处理中，导致区域分割的精度降低。另外，如果块尺寸过大，则就会在一个块中混合存在细胞的特征与除此以外的对象的特征，还会成为精度降低的原因。根据本申请发明人的见解，优选的是，块BL的大小选择为在作为图像中的提取对象的细胞的内部放入大致一个块左右的大小。
当细胞是例如人体的体细胞时，其直径大致为10μm～20μm左右，但干细胞稍小于此。因此，优选的是，对于块BL的尺寸，当换算为试料尺寸时的一边的长度等于或稍小于此。因此，优选的是，一个像素P的尺寸为一边是2.5μm～5μm或者稍小于此。
如果将拍摄中的分辨率设为例如6400dpi(dots per inch)，则像素尺寸为约4μm。另外，如果将分辨率设为16000dpi，则像素尺寸为约1.6μm。因此，作为用于获得上述的像素P的尺寸所需的分辨率，优选是6400dpi至16000dpi。但是，这些是适用于进行本实施方式的区域分割处理的值，在拍摄时，不受此限制，基本上以尽可能高的分辨率进行拍摄。
当能够以例如32000dpi(像素尺寸0.8μm)的分辨率拍摄原始图像时，为了便于区域分割处理并不一定使分辨率降低。在该情况下，如图3C所示，将拍摄时的原始图像Ia中的4个像素以伪方式(pseudo manner)视为一个像素P，将这些平均像素值设为该像素P的像素值。这样，不需要使拍摄时的分辨率降低，就能够获得适合于区域分割处理的像素尺寸(1.6μm)。在本实施方式中将这样的处理称为“分辨率变换”。
是否需要进行分辨率变换是通过拍摄时的像素尺寸与处理对象即细胞的尺寸之间的关系而确定。在预先已知这些关系的情况下，能够根据这些已知的关系判断是否需要进行分辨率变换。另外，是否进行分辨率变换也可以基于来自用户的指示输入而确定。在判断需要的情况下，执行分辨率变换(步骤S103)。
在步骤S104中，当进行分辨率变换时，将变换后的图像中的4个像素×4个像素的正方形作为1个块BL。如果不进行分辨率变换，则将在拍摄原始图像Ia时的4个像素×4个像素的正方形作为一个块BL。在X方向以及Y方向的至少一个方向的一行中的像素数量不能被一块的一边的像素数量4整除的情况下，在原始图像的端部可以产生未被计入处理的像素。另外，可以预先去除原始图像的端部以免出现这样的小数
接下来，将从各像素的像素值减去原始图像Ia整体中的平均像素值作为偏移值(步骤S105)。当像素值表示为例如8位数据时，各像素可以获取从0到255的256级的值，但通过减去平均像素值，也会出现获取负值作为像素值。根据本申请发明人的见解，在图像中的最大像素值与最小像素值之间的差比较小的细胞图像中，通过从各像素减去图像整体的平均像素值，从而各像素值的正负平衡变得良好。当将各像素P的像素值设为p(X，Y)、将平均像素值设为pav时，进行减法后的像素P的像素值g(X，Y)通过下式(数学式1)来表示。
g(X,Y)＝p(X,Y)-pav…(数学式1)
对于通过这样方式分割的每个块，执行DCT运算(步骤S106)。通过DCT运算，各块的图像内容分别在X方向以及Y方向上分解为二维空间频率分量。众所周知，DCT运算处理也用于JPEG方式的图像压缩处理。因此，下面，将省略对能够应用与JPEG图像压缩处理相同的想法的事项的说明。
图4是示出DCT运算的概念的图。图左侧示出了原始图像的一部分即由(4×4)像素组成的一个块BL。使用小写字符是为与原始图像的坐标位置相区分以表示各像素在块内的位置，将横向坐标设为x坐标，将纵向坐标设为y坐标。DCT运算能够认为是根据块BL内的各像素的像素值g(x，y)来求出图右侧所示的4行×4列的DCT系数表T的运算。
DCT系数表T中的行位置u(u＝0、1、2、3)表示在x方向上的空间频率分量，具体来说，u＝0表示直流分量、u＝1表示基本频率分量、u＝2、3分别表示基本频率的2倍、3倍的频率分量。同样地，列位置v(v＝0、1、2、3)表示在y方向上的空间频率分量，v＝0表示直流分量、v＝1表示基本频率分量、v＝2、3分别表示基本频率的2倍、3倍的频率分量。
表示各频率分量的强度的DCT系数C(u，v)通过下式(数学式2)表示。
Figure BDA0002110193120000091
其中，
Figure BDA0002110193120000092
其中，定数N是块尺寸。在JPEG方式中，将N＝8，即，(8×8)像素作为1块，如上所述，但在该实施方式中由于将(4×4)像素作为1块，因此N＝4。根据本申请发明人的见解，在将DCT系数利用为特征量而不以图像的压缩/恢复为目的本处理方法中，将块尺寸设为(4×4)像素的方式更能得到良好的结果。另外，这样，也能够降低运算量并适合于高速化。
在构成DCT系数表T的DCT系数C(u，v)中，表示直流分量的系数C(0，0)为直流分量特征量。另外，表示交流分量的各系数被分为多组，并在组内进行合计。将该合计值、或者合计值除以合计的系数的数量来归一化的值作为与该组对应的一个特征量。需要说明的是，DCT系数也可以取正负任意值。为了避免在合计时抵消这些值，例如DCT系数的合计能够是他们各自的绝对值的和。
图5A至图5C是示出DCT系数的分组的例子的图。基本的想法是：如图5A所示，将最大频率为基本频率(u＝1、v＝1)的三个DCT系数，如果用符号表示则将C(1，0)、C(0，1)、以及C(1，1)设为第一组G1。另外，将最大频率是基本频率的2倍(u＝2、v＝2)的五个DCT系数C(2，0)、C(2，1)、C(0，2)、C(1，2)、以及C(2，2)设为第二组G2。进一步地，将除去上述以及直流分量C(0，0)的七个DCT系数C(3，0)、C(3，1)、C(3，2)、C(0，3)、C(1，3)、C(2，3)以及C(3，3)设为第三组G3。这基于对在图像中表示相同程度的纹理的细微度的频率分量的强度进行合计的想法。
能够将分别对各组G1、G2、G3中包含的DCT系数进行合计的值，或者合计值除以该组的DCT系数的个数(例如在组G1中为3)来归一化的值作为与该频率对应的交流分量特征量。代表各组的频率能够通过例如该组中的最大频率来表示。
进行归一化的理由是因为在后面的分类工序中使用相同的权重来评价各频率分量。由此，能避免仅与确定的频率分量对应的纹理对分类结果造成很大的影响。需要说明的是，直流分量特征量与纹理无关，表示块整体的亮度。因此，与表示纹理的交流分量特征量之间的归一化不一定是必须的。
另外，能够适当地执行分组而不限于上述原则。如图5B所示，例如能够将除了直流分量的系数C(0，0)以及第一组G1之外的12个系数全部包含在第二组G2中。
但是，如图5A中虚线所示，优选的是，在DCT系数表T中，相对于将变量u、v视为坐标时的直线v＝u的彼此对称的位置处的DCT系数包含在相同的组中。这是因为，变量u与图像的x坐标对应，变量v与图像的y坐标对应，在作为拍摄对象物的培养细胞等图像中，在原理上由于x方向与y方向之间没有特别的异向性，因此并不应该使在u方向与v方向上的处理不同。
另外，在图5C所示的分组的例子中，系数C(1，0)、C(2，0)、C(0，1)、C(1，1)、C(0，2)的五个包含在第一组G1中。另外，与这些相邻的系数C(3，0)、C(2，1)、C(1，2)、C(2，2)、C(0，3)的五个包含第二组G2中。进一步地，其他的五个系数C(3，1)、C(3，2)、C(1，3)、C(2，3)、C(3，3)包含在第3组G3中。如果通过如上的操作，则值在各组中进行合计的DCT系数的个数全部相同，由于在各组求出的交流分量特征量之间权重相同，因此不需要进行归一化。
虽然构成的组的数量越多就越多地保存表示图像的特征的信息，但由于特征量空间的维度变高，因此在后面的处理中的运算量增大。通过将频率分量比较接近的系数彼此进行分组，从而既能够维持图像的信息，又能够降低运算负载。
由通过上述方式求出的各频率分量的特征量构成表示各块BL的亮度以及纹理的特征量矢量。即，表示该块BL的亮度的直流分量特征量和表示块内的纹理的多个交流分量特征量构成特征量矢量。
基于该特征量矢量，能判定各块相当于多种“细胞”、“非细胞”中的哪一种类别(步骤S108)。该判定是基于在由各频率分量的特征量的组合构成的多维多维特征量空间中与该块对应的特征量矢量指示的点的位置而进行的。为了进行判定，需要预先确定该判定基准，具体来说确定在多维特征量空间中各类别的特征占据的范围。
在多个拍摄之间，在固定拍摄的细胞的种类、培养基的种类以及量、拍摄条件(拍摄倍率或照明强度等)等的情况下，收集多个以相同条件拍摄的典型的细胞的图像作为事例，并使用适当的机器学习算法。通过这样操作，在多维特征量空间中，能够预先确定与各类别对应的集群占据的范围。在该情况下，根据从未知的原始图像中选出的块的特征量属于哪个集群，能够判定该块与哪个种类的细胞或者非细胞对应。
但是，在通常情况下，拍摄条件等因每个原始图像不同而不同。因此，为了进行更高精度的判定，如果能够根据从作为处理对象的原始图像读取的信息来制定判定基准，则是方便的。以下，对该方法的一例进行说明。
图6是示出用于制定判定基准的处理的流程图。另外，图7A以及图7B是例示处理过程中的指示输入方式的图。该处理接收来自操作人员针对在获取的原始图像Ia中较强地具有多种细胞的各自的特征的典型例的教示，并根据该信息确定多维特征量空间中的类别之间的边界。该处理能够在计算机装置1获取到原始图像Ia(步骤S101)以后到开始对块进行分类(步骤S108)为止的任意时刻执行。处理通过由CPU10执行存储在存储器15中的控制程序来实现。
最初，将获取的原始图像Ia显示在显示部17上提示给操作人员(步骤S201)。然后，在显示的原始图像Ia的中，接受来自操作人员的教示输入，该教示输入指定属于细胞的区域的代表点和与非细胞的区域对应的代表点(步骤S202)。在存在多种细胞并需要识别他们的情况下，接受每个种类的教示输入。需要说明的是，在细胞的区域与非细胞的区域在外观上有明显的差异的情况下，可以仅指定属于细胞的区域的代表点。
为了提高精度，优选指定多个代表点。因此，如图7A所示，例如，通过鼠标操作或者触摸笔输入等，使操作人员执行在细胞C2的区域内连续移动指针的操作或者选择几点的操作。这样，能够获取并收集在该轨迹上的多个点作为细胞C2的代表点。图中的白色空心箭头表示与由用户操作鼠标等连动地在画面上移动的指针。对于种类不同的细胞C4、C5或背景B，也同样地接受代表点的指定。
另外，如图7B的虚线所示，可以将属于细胞内的矩形区域、细胞外的矩形区域等指定为二维区域。优选的是，当图像内存在碎片时，使指定的细胞外的区域中包括碎片的区域与背景的区域两者。这样，在从细胞外的区域提取出的代表点中均包含具有碎片的特征的点和具有背景的特征的点。因此，当确定特征量空间中的细胞的区域与非细胞的区域的边界时，能够更强烈地反映与非细胞相关的信息。当然，也可以对碎片和背景分别指定代表点。
如果通过这样的方式指定若干个代表点，则与上述区域分割处理同样地，从各像素的像素值减去原始图像Ia整体中的平均像素值(步骤S203)，并对包含代表点的(4×4)像素尺寸的虚拟的块执行DCT运算(步骤S204)，根据DCT系数来计算直流分量特征量以及多个交流分量特征量(步骤S205)。由此，将包含代表点的块作为典型例的多维特征量空间中的分布变清晰，根据其结果能确定多维特征量空间中的各种“细胞”的集群占据的区域和“非细胞”的集群占据的区域(步骤S206)。
图8A以及图8B是示意性地示出特征量空间中的集群分布的例子的图。虽然特征量空间的维度能取各种值，但在这里图示了三维的例子。使用图7A或者图7B所示的方法，如果由用户指定多种细胞区域以及非细胞区域，则针对包含这些区域内的代表点的若干个块能求出直流分量特征量以及交流分量特征量。如图8A所示，认为如果在多维特征量空间绘制每个块，则由于各种细胞以及非细胞各自的图像特征的差异，因此形成与细胞C2对应的集群CL1、与细胞C4对应的集群CL2、与细胞C5对应的集群CL3、与除这些以外的非细胞即背景B以及碎片对应的集群CL4。
然后，当对构成原始图像Ia的各块BL进行分类时，通过表示该块的特征量矢量在多维特征量空间中占据的位置应该属于哪个集群，从而将该块分类为任一集群。例如，如果特征量矢量指示的位置处于任一集群区域的内部，则能够判定该块属于该集群。
另外，如图8B所示，根据应该分类的块的特征量矢量分别求出示出位置的点Pa与各集群之间的距离，可以将该距离分类为最近的集群。对于点Pa与集群之间的距离，例如能够通过如下的方式求出。对于各集群，预先确定每个距离基准点P1～P4。其中，作为距离基准点P1～P4可以是例如各集群在特征量空间中占据的空间区域的重心点。另外，可以将与从用户指示的代表点的平均位置对应的点设为距离基准点。然后，根据点Pa与距离基准点之间的距离，能够表示点Pa与各集群之间的距离。
如上所述，在该实施方式的区域分割处理中，将包含细胞的原始图像分割为多个小块，并求出每个块的空间频率分量。然后，构成特征量矢量，该特征量矢量包含基于空间频率分量中的直流分量的强度的直流分量特征量和基于交流分量中的多个频率分量各自的强度的多个交流分量特征量。然后，根据该特征量矢量在多维特征量空间所示的位置属于哪个集群，能判定该块的图像是与细胞对应的图像还是与除此以外的非细胞对应的图像，或者与哪个种类的细胞相对应。
直流分量特征量表示块的亮度。另外，在交流分量特征量中与频率比较低的分量对应的特征量是表示结构比较简单的特征量，与高频分量对应的特征量是表示块内的细微的纹理的特征量。这样，通过对亮度和纹理进行综合评价，从而能够通过比基于简单的二值化处理或轮廓提取处理的区域分割更优的精度来区分图像内的细胞和非细胞。因此，当多个种类的细胞包含在图像中时，能够区分它们。
需要说明的是，本发明并不限于上述的实施方式，能够在不脱离其主旨的情况下进行除上述实施形式之外的各种变更。例如，在上述实施方式中的原始图像的区域分割(步骤S109)中，将原始图像以块为单位分割为细胞或非细胞。然而，块边界是自动制定的，并不限于与原始图像中的细胞和其周围的边界一致。鉴于此，例如可以与公知的轮廓提取处理组合地进行原始图像的分割。即，通过上述的方式，可以通过根据轮廓提取处理的结果来修正以块为单位确定的区域的边界的位置。
另外，在上述实施方式中，确定在多维特征量空间的集群的处理仅是示例性的，除此以外能够应用各种分群方法。另外，对于收集细胞以及非细胞的典型例的方法，并不限定于上述方法。
另外，在上述实施方式中，由通过离散余弦变换处理求出的DCT系数表示各频率分量的强度。然而，除了DCT运算以外，如果存在能求出适当的空间频率分量的方法，也可以利用该方法。
另外，在上述实施方式中，虽然由通用的计算机装置1执行区域分割处理，但如上所述例如可以在拍摄装置装入该处理功能。另外，还能够在现有的拍摄装置中追加安装用于执行本实施方式的区域分割处理的程序来进行功能扩展。
另外，虽然在上述描述中例示了像素尺寸或块尺寸的具体数值，但是即使在这些数值范围之外的条件下也可以执行本实施方式的区域分割处理。另外，图像中包含的细胞的种类是单一的，并且当分割细胞与其他区域时也能够应用本发明的处理方法。
以上，如所例示并说明具体的实施方式那样，在本发明中，对每个块可以执行基于该块中包括的像素的像素值的二维离散余弦变换，并且可以基于各频率分量的离散余弦变换系数来求出特征量。这种处理也用于JPEG的图像压缩处理，并且由于已经建立了运算算法，所以这种处理也能够应用于本发明。
另外，可以将根据空间频率的高度对由二维离散余弦变换求出的多个离散余弦变换系数进行分组而得到的多组中的每一组的离散余弦变换系数的合计，以及根据合计的离散余弦变换系数的数量来归一化合计的值中的任意一个作为特征量。根据这样的结构，通过将例如比较接近的频率分量的离散余弦变换系数彼此结合为一个，就能降低特征量空间的维度而不会丢失与图像的特征有关的信息，从而降低在后面工序中的运算负载。
在该情况下，可以基于从像素各自的像素值减去原始图像整体的平均像素值而得到的值执行二维离散余弦变换。由于在对接近透明的细胞进行拍摄而得到的图像中像素值的分布宽度较小，因此优选的是，使用原始图像整体的平均像素值来进行用于离散余弦变换运算的电平位移。
另外，可以求出多个频率分量的强度的合计来作为各频率分量的离散余弦变换系数的绝对值的和。由于离散余弦变换系数可以取正负任意值，因此在简单的合计中有可能抵消各频率分量而不能获得良好的结果。通过获取绝对值的和，能消除该问题。
另外，各块可以是一边相当于4个像素的块尺寸的正方形。在JPEG方式的图像压缩中，块尺寸为8个像素，但根据本申请发明人的见解，在不以图像的压缩/恢复为目的的本发明中，将块尺寸设为4个像素的方式更能得到良好的结果。
另外，在将原始图像分割为块之前，还可以设置有如下工序：可以将原始图像的像素尺寸从获取时的原始像素尺寸变换为与原始像素尺寸不同的块分割用的像素尺寸。在利用细胞和除此以外的非细胞之间的亮度和纹理的差异来进行区域分割的本发明中，运算时优选的像素尺寸取决于例如细胞等图像中的对象的尺寸。因此，虽然像素尺寸不一定与拍摄时的优选的像素尺寸一致，但通过进行像素尺寸的变换，能够谋求分别单独优化像素尺寸。
另外，在分类之前，还可以设置有如下工序：在原始图像中，接受来自用户的关于与多个集群各自对应的代表点的教示输入；对包含代表点的块，基于教示输入的结果求出特征量，并基于该结果，确定多维特征量空间中的集群各自的范围。根据这样的结构，即使是细胞的种类或拍摄条件不恒定的情况，也能够稳定地进行细胞或除其以外的区域的分割。
另外，原始图像可以是通过在明亮视野下拍摄细胞而得到的图像。通常细胞接近透明，另外，由于与培养基之间的折射率的差也较小，在明亮视野下拍摄的图像中很难通过目视来区分细胞与非细胞。通过将本发明应用于这样的图像，能够良好且稳定地分割图像中的细胞的区域与除其以外的区域。
工业实用性
本发明能够适用于拍摄细胞来评价图像的生物化学或医疗领域，特别适用于需要不进行像损伤培养的细胞这样的加工的拍摄的技术领域。

Claims (12)

1.一种图像处理方法，其中，所述图像处理方法包括如下工序：
获取包含有培养的细胞的原始图像；
将所述原始图像分割为由规定数量的像素构成的块，求出各块内的图像所具有的空间频率分量；
将所述各块分类为多维特征量空间中的多个集群，所述多维特征量空间将所述空间频率分量中的直流分量的强度和空间频率彼此不同的多个交流分量的每一个交流分量的强度作为特征量；以及
基于所述分类的结果将所述原始图像分割为多个区域。
2.根据权利要求1所述的图像处理方法，其中，
针对每个所述块，执行基于该块中包含的所述像素的像素值的二维离散余弦变换，并基于各频率分量的离散余弦变换系数来求出所述特征量。
3.根据权利要求2所述的图像处理方法，其中，
将根据空间频率的高度对由所述二维离散余弦变换求出的多个所述离散余弦变换系数进行分组而得到的多个组中的每一组的离散余弦变换系数的合计、以及根据合计的所述离散余弦变换系数的数量将所述合计进行归一化而得到的值中的任意一个作为特征量。
4.根据权利要求2所述的图像处理方法，其中，
基于从所述像素的各个的像素值减去所述原始图像整体的平均像素值而得到的值来执行所述二维离散余弦变换。
5.根据权利要求3所述的图像处理方法，其中，
基于从所述像素的各个的像素值减去所述原始图像整体的平均像素值而得到的值来执行所述二维离散余弦变换。
6.根据权利要求3所述的图像处理方法，其中，
求出所述合计，作为所述离散余弦变换系数的绝对值的和。
7.根据权利要求2所述的图像处理方法，其中，
所述块是一边相当于4个像素的块尺寸的正方形。
8.根据权利要求1～7中任一项所述的图像处理方法，其中，
在将所述原始图像分割为所述块之前，具有如下工序：将所述原始图像的像素尺寸从获取时的原始像素尺寸变换为与所述原始像素尺寸不同的块分割用的像素尺寸。
9.根据权利要求1～7中任一项所述的图像处理方法，其中，
在所述分类之前，执行如下工序：
接受来自用户的关于与所述原始图像中所述多个集群的各集群对应的代表点的教示输入；以及
针对包含所述代表点的所述块，基于所述教示输入的结果求出所述特征量，并基于该结果，确定所述多维特征量空间中的所述集群的各自的范围。
10.根据权利要求1～7中任一项所述的图像处理方法，其中，
所述原始图像是通过在明亮视野下拍摄所述细胞而得到的图像。
11.一种计算机程序，使计算机执行权利要求1～7中任一项所述的图像处理方法的各工序。
12.一种计算机可读的存储介质，存储有权利要求11中所述的计算机程序。
CN201910568602.XA 2018-06-27 2019-06-27 图像处理方法、计算机程序以及记录介质 Pending CN110717882A (zh)

Applications Claiming Priority (2)

Application Number Priority Date Filing Date Title
JP2018122123A JP2020004033A (ja) 2018-06-27 2018-06-27 画像処理方法、コンピュータプログラムおよび記録媒体
JP2018-122123 2018-06-27

Publications (1)

Publication Number Publication Date
CN110717882A true CN110717882A (zh) 2020-01-21

Family

ID=66998174

Family Applications (1)

Application Number Title Priority Date Filing Date
CN201910568602.XA Pending CN110717882A (zh) 2018-06-27 2019-06-27 图像处理方法、计算机程序以及记录介质

Country Status (4)

Country Link
US (1) US11025924B2 (zh)
EP (1) EP3588435B1 (zh)
JP (1) JP2020004033A (zh)
CN (1) CN110717882A (zh)

Citations (5)

* Cited by examiner, † Cited by third party
Publication number Priority date Publication date Assignee Title
JPH0216682A (ja) * 1988-07-04 1990-01-19 Sumitomo Electric Ind Ltd 細胞増殖状態のモニタ方法及びそのモニタ装置
CN1622137A (zh) * 2005-01-07 2005-06-01 北京北大方正电子有限公司 一种对图像进行快速压缩和解压缩的方法
CN1925550A (zh) * 2005-08-30 2007-03-07 夏普株式会社 图像处理方法、图像处理装置、图像形成装置及记录介质
JP2007275346A (ja) * 2006-04-07 2007-10-25 Toshiba Corp 放射線画像処理装置、放射線画像処理方法及び放射線画像処理プログラム
WO2008139825A1 (ja) * 2007-05-08 2008-11-20 Olympus Corporation 画像処理装置および画像処理プログラム

Family Cites Families (3)

* Cited by examiner, † Cited by third party
Publication number Priority date Publication date Assignee Title
US8199989B2 (en) * 2009-05-15 2012-06-12 General Electric Company Automatic fly through review mechanism
JP5568970B2 (ja) * 2009-11-30 2014-08-13 ソニー株式会社 情報処理装置、情報処理方法及びそのプログラム
JP6333145B2 (ja) 2014-09-30 2018-05-30 株式会社Ｓｃｒｅｅｎホールディングス 画像処理方法および画像処理装置

Patent Citations (5)

* Cited by examiner, † Cited by third party
Publication number Priority date Publication date Assignee Title
JPH0216682A (ja) * 1988-07-04 1990-01-19 Sumitomo Electric Ind Ltd 細胞増殖状態のモニタ方法及びそのモニタ装置
CN1622137A (zh) * 2005-01-07 2005-06-01 北京北大方正电子有限公司 一种对图像进行快速压缩和解压缩的方法
CN1925550A (zh) * 2005-08-30 2007-03-07 夏普株式会社 图像处理方法、图像处理装置、图像形成装置及记录介质
JP2007275346A (ja) * 2006-04-07 2007-10-25 Toshiba Corp 放射線画像処理装置、放射線画像処理方法及び放射線画像処理プログラム
WO2008139825A1 (ja) * 2007-05-08 2008-11-20 Olympus Corporation 画像処理装置および画像処理プログラム

Non-Patent Citations (2)

* Cited by examiner, † Cited by third party
Title
PRAMOD K. SINGH ET AL.: "Unsupervised Segmentation of Medical Images using DCT Coefficients" *
刘凌霞: "基于分块DCT的视频文字检测算法" *

Also Published As

Publication number Publication date
US11025924B2 (en) 2021-06-01
EP3588435B1 (en) 2020-11-18
EP3588435A1 (en) 2020-01-01
JP2020004033A (ja) 2020-01-09
US20200007871A1 (en) 2020-01-02

Similar Documents

Publication Publication Date Title
CN111583227B (zh) 一种荧光细胞自动计数方法、装置、设备和介质
JP6660313B2 (ja) 画像解析を用いた核のエッジの検出
CN109978839B (zh) 晶圆低纹理缺陷的检测方法
CN108416766B (zh) 双侧入光式导光板缺陷视觉检测方法
US10168526B2 (en) Cell contour formation apparatus and method of the same, and non-transitory computer readable storage medium storing a cell contour formation program
WO2019181072A1 (ja) 画像処理方法、コンピュータプログラムおよび記録媒体
CN110956612A (zh) 一种细胞快速计数方法、计算机设备及可读存储介质
CN109949294A (zh) 一种基于OpenCV的断口形貌图裂纹缺陷提取方法
CN110648312A (zh) 一种基于鳞片形态特征分析的羊毛与羊绒纤维识别方法
CN114862765A (zh) 一种细胞分层图像处理方法
Soni et al. CT scan based brain tumor recognition and extraction using Prewitt and morphological dilation
CN108717699B (zh) 一种基于连续最小割的超声波图像分割方法
CN112184696B (zh) 一种细胞核和细胞器计数及其面积计算方法与系统
CN117765330A (zh) 基于mri图像的数据标注方法及系统
CN110458042B (zh) 一种荧光ctc中的探针数目检测方法
CN109448010B (zh) 一种基于内容特征的四方连续纹样自动生成方法
CN110717882A (zh) 图像处理方法、计算机程序以及记录介质
Khan et al. Segmentation of single and overlapping leaves by extracting appropriate contours
Guatemala-Sanchez et al. Nuclei segmentation on histopathology images of breast carcinoma
CN113470028A (zh) 染色体核型图像质量评估方法、染色体分析仪和存储介质
RU118774U1 (ru) Автоматизированная система анализа биомедицинских микроизображений для обнаружения и характеризации информативных объектов заданного вида на неоднородном фоне
Singh et al. Non Uniform Background Removal using Morphology based Structuring Element for Particle Analysis
RU136612U1 (ru) Устройство диагностики области новообразования типа "рак молочной железы" на маммографических снимках
CN113344810B (zh) 基于动态数据分布的图像增强方法
CN114240988B (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
WD01 Invention patent application deemed withdrawn after publication

Application publication date: 20200121

WD01 Invention patent application deemed withdrawn after publication