CN113344765B - 一种频域天文图像目标检测方法及系统 - Google Patents

一种频域天文图像目标检测方法及系统 Download PDF

Info

Publication number
CN113344765B
CN113344765B CN202110527356.0A CN202110527356A CN113344765B CN 113344765 B CN113344765 B CN 113344765B CN 202110527356 A CN202110527356 A CN 202110527356A CN 113344765 B CN113344765 B CN 113344765B
Authority
CN
China
Prior art keywords
image
astronomical
sub
frequency domain
convolution kernel
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
Application number
CN202110527356.0A
Other languages
English (en)
Other versions
CN113344765A (zh
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.)
National Space Science Center of CAS
Original Assignee
National Space Science Center of CAS
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 National Space Science Center of CAS filed Critical National Space Science Center of CAS
Priority to CN202110527356.0A priority Critical patent/CN113344765B/zh
Publication of CN113344765A publication Critical patent/CN113344765A/zh
Application granted granted Critical
Publication of CN113344765B publication Critical patent/CN113344765B/zh
Active legal-status Critical Current
Anticipated expiration legal-status Critical

Links

Classifications

    • G—PHYSICS
    • G06—COMPUTING; CALCULATING OR COUNTING
    • G06T—IMAGE DATA PROCESSING OR GENERATION, IN GENERAL
    • G06T1/00—General purpose image data processing
    • G06T1/20—Processor architectures; Processor configuration, e.g. pipelining
    • 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
    • G06F17/14—Fourier, Walsh or analogous domain transformations, e.g. Laplace, Hilbert, Karhunen-Loeve, transforms
    • G06F17/141—Discrete Fourier transforms
    • G06F17/142—Fast Fourier transforms, e.g. using a Cooley-Tukey type algorithm
    • G—PHYSICS
    • G06—COMPUTING; CALCULATING OR COUNTING
    • G06F—ELECTRIC DIGITAL DATA PROCESSING
    • G06F18/00—Pattern recognition
    • G06F18/20—Analysing
    • G06F18/21—Design or setup of recognition systems or techniques; Extraction of features in feature space; Blind source separation
    • G06F18/213—Feature extraction, e.g. by transforming the feature space; Summarisation; Mappings, e.g. subspace methods
    • G—PHYSICS
    • G06—COMPUTING; CALCULATING OR COUNTING
    • G06F—ELECTRIC DIGITAL DATA PROCESSING
    • G06F9/00—Arrangements for program control, e.g. control units
    • G06F9/06—Arrangements for program control, e.g. control units using stored programs, i.e. using an internal store of processing equipment to receive or retain programs
    • G06F9/46—Multiprogramming arrangements
    • G06F9/54—Interprogram communication
    • G06F9/544—Buffers; Shared memory; Pipes
    • G—PHYSICS
    • G06—COMPUTING; CALCULATING OR COUNTING
    • G06T—IMAGE DATA PROCESSING OR GENERATION, IN GENERAL
    • G06T5/00—Image enhancement or restoration
    • G06T5/10—Image enhancement or restoration using non-spatial domain filtering
    • G—PHYSICS
    • G06—COMPUTING; CALCULATING OR COUNTING
    • G06T—IMAGE DATA PROCESSING OR GENERATION, IN GENERAL
    • G06T5/00—Image enhancement or restoration
    • G06T5/70—Denoising; Smoothing
    • G—PHYSICS
    • G06—COMPUTING; CALCULATING OR COUNTING
    • G06T—IMAGE DATA PROCESSING OR GENERATION, IN GENERAL
    • G06T7/00—Image analysis
    • G06T7/30—Determination of transform parameters for the alignment of images, i.e. image registration
    • G06T7/33—Determination of transform parameters for the alignment of images, i.e. image registration using feature-based methods
    • G06T7/337—Determination of transform parameters for the alignment of images, i.e. image registration using feature-based methods involving reference images or patches
    • 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/10004—Still image; Photographic image

Landscapes

  • Engineering & Computer Science (AREA)
  • Physics & Mathematics (AREA)
  • General Physics & Mathematics (AREA)
  • Theoretical Computer Science (AREA)
  • Mathematical Physics (AREA)
  • Data Mining & Analysis (AREA)
  • General Engineering & Computer Science (AREA)
  • Software Systems (AREA)
  • Computer Vision & Pattern Recognition (AREA)
  • Mathematical Analysis (AREA)
  • Mathematical Optimization (AREA)
  • Pure & Applied Mathematics (AREA)
  • Computational Mathematics (AREA)
  • Algebra (AREA)
  • Databases & Information Systems (AREA)
  • Discrete Mathematics (AREA)
  • Life Sciences & Earth Sciences (AREA)
  • Artificial Intelligence (AREA)
  • Bioinformatics & Cheminformatics (AREA)
  • Bioinformatics & Computational Biology (AREA)
  • Evolutionary Biology (AREA)
  • Evolutionary Computation (AREA)
  • Image Processing (AREA)

Abstract

本发明公开了一种频域天文图像目标检测方法及系统，该方法基于CPU‑GPU异构处理器实现，该方法包括：基于预先得到的参考图像，由CPU对采集的原始天文图像进行预处理；采用重叠保存的方法对参考图像和预处理后的天文图像分别进行图像分块，采用高斯基函数和多项式相乘得到卷积核的n组基向量并输入GPU；由GPU根据参考子图像和天文子图像，结合n组基向量进行拟合，得到对应每个参考子图像的卷积核，使用卷积核对每个参考子图像进行频域滤波降晰处理，得到模板图像并输入CPU；由CPU将模板图像的边缘丢弃，并将剩余部分连接并与原始天文图像作差，得到差异图像，进而实现对天文图像的目标检测。

Description

一种频域天文图像目标检测方法及系统
技术领域
本发明属于目标检测领域，涉及一种频域天文图像目标检测方法及系统，可用于星载嵌入式系统。
背景技术
在我国空间发展计划公布的空间科学和探测任务中，许多任务对航天器的人工智能控制、自主导航与控制、实时在轨处理速度、存储数据吞吐率、海量数据处理、复杂科学计算、图像处理等方面提出了迫切的要求。尤其是深空任务，由于通信延迟和带宽的限制，比如通常到火星的通信延迟20分钟、带宽250kbps，到木星的通信延迟2小时、带宽120kbps，有些探测任务来不及依赖地面控制，而科学探测的原始数据不能全部回传到地面，因此对航天器的智能自治和实时数据处理能力提出了更高的要求。而我国的卫星数据处理流程基本采用的是“地面任务规划-星上采集-星地数传-地面处理-分发应用”的传统应用模式。该数据处理方法不仅实时性低，而且还需要占用大量的通信资源。数据获取的时效性和智能化水平难以满足新形势下对卫星数据获取的实时性要求。星上实时数据处理是解决这些问题的有效途径，已成为未来深空自适应探测技术发展的一个重要方向。因此需要在星上进行自适应的数据处理，对目标进行直接的识别和判断。星上自适应科学探测系统一般有两种需求：首先，上述系统的数据处理必须实时。其次，用于执行识别的算法不应影响精度。
对瞬变、移动或亮度变化的天体和暂现源的观测是天文学领域研究的一项重要任务。Tomaney和Crotts在1996年首次提出频域的图像相减算法，首先使用数据图像匹配中应用最广泛的PSF作为参考来进行差分成像。该方法通过简单地获取每个图像上光点的傅立叶变换的比值来确定卷积核，由于这种反卷积方法会降低图像质量，因此它只适用于信噪比较高的数据。Kochanski,Tyson&Fischer(1996)提出的非线性PSF拟合则更为准确，但在数值计算上比较耗时。Alard和Lupton(1998)开发了一种最优图像减影算法(OpticalImage Subtraction，OIS)，随后，Alard(2000)引入了空间变换卷积核的概念，该算法使用空间变换卷积核对两个图像进行卷积。因此，该方法不仅适用于有天体密度高的区域，也适用于天体密度较低的区域。Israel、Bramich和Zackay分别于2007年、2008年和2016年在Alard研究的基础上提出了改进算法,降低了由配准误差引起的残差。目前，解决问题的最新技术是训练一种机器学习算法来过滤大部分残影，并将伪目标的数量减少到最小(e.g.,Cabrera et al.2017；Masci et al.2017；Sedaghat et al.2018)。但上述的所有算法都是关注精度，而不是运行速度。在处理速度的研究上，李继良等人在2010年使用Open MP对图像差异算法进行并行优化加速，在Dell Power Edge T300 System服务器上，加速比可达1.754。Yan Zhao等人在2013年使用GPU(图形处理单元)对空间变化的内核进行图像卷积，在Intel i7 CPU和NVIDIAGTX580 GPU的台式机上算法的运行速度提升了4倍，基本满足了地面实时处理的需要。虽然GPU在时域天文学的应用大大提高了数据的处理能力，但是这并不能从根本上解决数据获取的实时性问题，那是因为目前的天文观测仍然采用在终端收集数据，通过星地数传将数据下传给地面，地面再进行处理的数据处理模式。面对越来越多的天文数据，这种传统的数据处理模式费时又费力，但如果能在卫星上进行一部分的数据处理和筛选工作，过滤冗余信息，提取出有价值的观测信息下传给地面，这不仅能大大提高信息获取的实时性，还能提高星上智能信息处理水平，降低星地之间的数传压力。
发明内容
本发明的目的在于克服现有技术缺陷，提出了一种频域天文图像目标检测方法及系统。
为了实现上述目的，本发明提出了一种频域天文图像目标检测方法，基于CPU-GPU异构处理器实现，所述方法包括：
基于预先得到的参考图像，由CPU对采集的原始天文图像进行预处理；采用重叠保存的方法对参考图像和预处理后的天文图像分别进行图像分块，得到尺寸和个数均相同的参考图像子图像和天文图像子图像；采用高斯基函数和多项式相乘得到卷积核的n组基向量并输入GPU；
由GPU根据参考子图像和天文子图像，结合n组基向量进行拟合，得到对应每个参考子图像的卷积核，使用卷积核对每个参考子图像进行频域滤波降晰处理，得到模板图像并输入CPU；
由CPU将模板图像的边缘丢弃，并将剩余部分连接并与原始天文图像作差，得到差异图像，进而实现对天文图像的目标检测。
作为上述方法的一种改进，所述基于预先得到的参考图像，由CPU对采集的原始天文图像进行预处理；具体包括：
由CPU从预先得到的参考图像中选取参考坐标系，通过分别在参考图像和输入图像上选取相同数量的天体以确定二维多项式形式的变换，利用该变换，在参考坐标系上对原始天文图像进行重新采样，实现与参考图像的对准；所述重新采样采用双三次样条法进行插值处理。
作为上述方法的一种改进，所述预处理前后的天文图像和参考图像的长宽尺寸均为N*N。
作为上述方法的一种改进，所述采用重叠保存的方法对参考图像和预处理后的天文图像分别进行图像分块，得到尺寸和个数均相同的参考图像子图像和天文图像子图像；具体包括：
设置卷积模板的长和宽均为正整数h；
根据空间可变卷积核的概念，对参考图像和预处理后的天文图像均进行以下相同操作：
按照尺寸为n*n进行划分，得到个不重叠的子图像，/>表示向上取整，不足的区域补0，随后根据重叠保存方法对子图像进行补0，将图像填充至大小。
作为上述方法的一种改进，所述采用高斯基函数和多项式相乘得到卷积核的n组基向量；具体包括：
使用三个高斯函数作为卷积核K的基函数，σk表示第k个基函数的方差，k＝1,2,3，(x,y)表示每个子图像的坐标，令
σ1＝0.7；i+j≤6
σ2＝1.5；i+j≤4
σ3＝3.0；i+j≤2
其中，i和j表示基函数的空间自由度，根据i和j的组合，得到n组基向量Kn(x,y)为：
作为上述方法的一种改进，所述由GPU根据参考子图像和天文子图像，结合n组基向量进行拟合，得到对应每个参考子图像的卷积核，具体包括：
令∑(T(x,y)*K(x,y)-I(x,y))2最小，以找到图像减法的内核，其中，T(x,y)表示参考子图像，I(x,y)表示天文子图像，K(x,y)为待求参考子图像的卷积核，(x,y)为每个子图像的坐标，表示内核的空间变化特征，基于每个子图像内点扩散函数恒定的设定条件，用子图像的中心坐标表示(x,y)；将最优化问题转换为线性最小二乘问题，求解以下线性系统，依据参考子图像和天文子图像得到对应每个参考子图像的卷积核的系数
其中，M为最小二乘的元素构成的矩阵，为向量Wm1和向量Wm2的标量积，系数的第n个元素为an(x,y)，/>为向量Wm1与I(x,y)的标量积，矩阵M的第m1m2个元素Mm1,m2满足下式：
Mm1,m2＝Wm1(x,y)Wm2(x,y)dxdy
其中，m1,m2＝0,1,2,…,n；
的第m1个元素Bm1满足下式：
采用n组基向量Kn(x,y)和系数an(x,y)进行拟合，得到卷积核K(x,y)。
作为上述方法的一种改进，所述使用卷积核对每个参考子图像进行频域滤波降晰处理，得到模板图像并输入CPU；具体包括：
按照重叠保存方法，并行的对每个参考子图像和对应的卷积核进行傅里叶变换，然后对计算结果相乘，并对矩阵相乘结果进行逆傅里叶变换，实现频域滤波降晰处理，得到模板图像并输入CPU。
一种频域天文图像目标检测系统，其特征在于，所述系统包括：部署在CPU的预处理模块、基函数计算模块和目标检测模块，以及部署在GPU的卷积核计算模块和频域滤波降晰模块，其中，
所述预处理模块，用于基于预先得到的参考图像，对采集的原始天文图像进行预处理；还用于采用重叠保存的方法对参考图像和预处理后的天文图像分别进行图像分块，得到尺寸和个数均相同的参考图像子图像和天文图像子图像；
所述基函数计算模块，用于采用高斯基函数和多项式相乘得到卷积核的n组基向量并输入卷积核计算模块；
所述卷积核计算模块，用于根据参考子图像和天文子图像，结合n组基向量进行拟合，得到对应每个参考子图像的卷积核；
所述频域滤波降晰模块，用于使用卷积核对每个参考子图像进行频域滤波降晰处理，得到模板图像并输入目标检测模块；
所述目标检测模块，用于将模板图像的边缘丢弃，并将剩余部分连接并与原始天文图像作差，得到差异图像，进而实现对天文图像的目标检测。
与现有技术相比，本发明的优势在于：
1、本发明的方法将原算法中最耗时的图像卷积部分转换到频域中，利用频域滤波来代替时域的卷积运算，优化了天文图像差异算法中的性能瓶颈，降低算法的计算复杂度，缩短了算法的运行时间；
2、本发明的方法通过重叠保存设计保留分块图像一定位数的原图像数据来延拓图像，在频域滤波完成后再将这些重叠部分舍弃，再按位相加，最大限度地利用图像数据，该设计既解决了卷积的循环特性产生频谱混叠问题，又避免了传统补零方法导致的内存占用过大的问题，保证了图像的连续处理，降低了算法的内存占用，提高了算法的计算效率；
3、本发明的方法采用重叠保存设计将每个子块的重叠部分分布在图像边缘，通过丢弃重叠部分消除了由扩充图像引入高频分量导致的干扰，保证了算法的精度；
4、本发明的方法充分利用了CPU和GPU的计算资源，将GPU用于卷积核的确定以及频域滤波计算，将图像卷积工作的高斯基函数Kn的计算分配给多核CPU并行处理，让CPU真正地参与到计算中，实现了异构平台的并行计算的效率提升；
5、本发明通过对方法的设计和优化，能够满足天文图像在星载设备上进行实时处理的约束，在此CPU-GPU异构并行计算系统内，基于重叠保存的频域天文图像差异算法既可以满足深空探测任务对于快速、实时的要求，也可以满足星载计算机低内存、低功耗的目标，为实现航天器的在轨自适应科学探测提供了可能。
附图说明
图1是基于CPU-GPU异构嵌入式平台处理天文图像差异算法的系统框图；
图2是基于重叠保存的频域天文图像差异算法在CPU-GPU异构平台上的映射；
图3是基于重叠保存的频域天文图像差异算法流程图；
图4是采用重叠保存法对输入矩阵进行重叠分割，其中，图4(a)为第一步示意，图4(b)为第二步示意，图4(c)为第三步示意，图4(d)为第四步示意；
图5是利用重叠保存法，通过FFT计算卷积结果示意图；
图6是基于重叠保存的频域天文图像差异算法的伪代码；
图7(a)是本发明的方法与时域卷积的计算复杂度对比；
图7(b)是本发明的方法与时域卷积的内存占用对比。
具体实施方式
本发明的目的在于针对上述现有星载计算资源的不足，以天文图像差异算法为出发点，提出一种基于重叠保存的频域天文图像差异算法，并基于CPU和GPU设计了一种基于CPU-GPU异构嵌入式计算平台的频域天文图像差异算法，整个处理器系统以CPU-GPU为基本处理单元，以天文图像为背景，利用CPU/GPU的并行技术，优化策略来对传统的差异算法进行优化加速，提高了天文图像差异算法的处理速度，基本满足在轨实时性处理的约束，为天文图像在星载嵌入式设备上进行实时处理提供可能。
本发明对天文图像差异算法的加速主要体现在以下三个方面：
1)首先针对信息处理的运算速度、实时性的需求，详尽分析了天文图像差异算法的处理流程以及性能瓶颈，并基于此设计了基于重叠保存的频域优化方案，减少了算法的计算复杂度和内存占用。
2)从计算平台出发，充分利用了CPU和GPU的计算资源，根据对基于重叠保存的频域加速算法的性能分析和流程分析，确定算法中可以并行化的部分，并根据OpenCL计算模型将并行化的部分映射到NDRange中的多个维度的多线程，完成异构平台任务的划分与映射。
3)从对关键算子优化的角度出发，通过对算法的分析，从算法中提取其主要算子，并在异构平台上根据算法本身和GPU异构计算的性能优化方法对设计进行深度优化和并行加速。
参照图1，本实例基于CPU-GPU的频域天文图像处理方法。嵌入式天文图像在轨处理系统，包括图像采集模块和图像计算模块，所述图像计算模块采用CPU与GPU异构处理器架构的星上在轨处理系统，所述图像采集模块对CCD相机信号进行实时采集并存储至LPDDR3，然后LPDDR3内的图像数据将传送至图像计算模块。
本实施方式中，所述的图像计算模块包括数据接收和数据发送。所述的数据接收是采用图像采集模块完成对CCD相机信号的实时采集，CPU作为控制器对图像采集模块进行逻辑控制并将接收数据写入LPDDR3缓存。
本实施方式中，所述图像计算模块包括图像对准、图像分块、确定卷积核、频域滤波以及输出五个步骤，实现了对空间中瞬变、移动或亮度变化的天体和暂现源的观测。为了能够完成实时的图像差异计算，提高观测效率，图像计算功能采用基于CPU-GPU异构处理器架构。图像处理部分以GPU作为协处理器，充分利用GPU片上数量众多的计算单元来对海量数据进行高性能通用的并行计算，提高图像处理的速度。另外，相比于FPGA，GPU具有通用性高，算法可移植性较高的特点，能够大幅降低系统设计难度，提高了星载系统的通用性和功能升级性。
随后本文所述方法将根据算法的特点将算法的不同部分划分给不同的处理器进行计算，将改进后的基于重叠保存的频域天文图像差异算法技术上述的计算不相关且密集型步骤并行化，并移植到GPU上运行。对于传统的CPU-GPU异构计算模式下，通常CPU是处理分支预测和逻辑性较强的运算，而GPU是处理数据相对比较密集且适合并行进行的计算。为了能够更加有效的利用异构系统的计算能力，本发明不仅将GPU用于两个数据并行计算组件，还将图像卷积工作的一部分分配给多核CPU，以充分利用整个异构计算平台的计算资源，让CPU真正的参与到计算中。在本发明中计算模块中的CPU的计算平台主要负责数据的预处理，基函数的计算以及接收结果，计算模块中的GPU计算平台是将所述方法中计算密集且独立的卷积核的计算以及卷积降晰这两个步骤，按照每个子图像进行任务划分，计算任务被分割成数百万个互不相关的微任务进行并行计算。基于重叠保存的频域天文图像差异算法在CPU-GPU异构平台上的实现方案如图2所示。
下面结合附图和实施例对本发明的技术方案进行详细的说明。
实施例1
本发明的实施例1提出了一种频域天文图像目标检测方法，包括以下步骤：
基于预先得到的参考图像，由CPU对采集的原始天文图像进行预处理；采用重叠保存的方法对参考图像和预处理后的天文图像分别进行图像分块，得到尺寸和个数均相同的参考图像子图像和天文图像子图像；采用高斯基函数和多项式相乘得到卷积核的n组基向量并输入GPU；
由GPU根据参考子图像和天文子图像，结合n组基向量进行拟合，得到对应每个参考子图像的卷积核，使用卷积核对每个参考子图像进行频域滤波降晰处理，得到模板图像并输入CPU；
由CPU将模板图像的边缘丢弃，并将剩余部分连接并与原始天文图像作差，得到差异图像，进而实现对天文图像的目标检测。
本方法使用CPU-GPU异构并行架构，并在移动中央处理单元(CPU)以及在移动图形处理单元(GPU)上使用OpenCL来加速基于重叠保存的频域天文图像差异算法。首先该算法通过将天文图像差异算法中计算量最大、最耗时的时域图像卷积降晰部分转换到频域上进行处理，降低了算法的计算复杂度；并通过重叠保存设计最大限度的利用原图像数据，解决了卷积的循环特性导致的混叠问题，提高了该算法的计算效率，降低了内存占用。随后通过利用算法的并行性，示例实施例优化基于重叠保存的频域天文图像差异算法以提高处理速度，并进一步利用现代移动处理器的异构架构，对工作负载进行划分并将它们分配到CPU和GPU二者上。在示例实施例中，优化使得能够在嵌入式设备上进行天文图像差异计算，并且可以应用于星载平台，实现航天器的在轨自适应天文探测。
对瞬变、移动或亮度变化的天体和暂现源的观测是空间巡天任务的至关重要一环，也是整个天文图像处理流程中最耗时的部分。它主要通过天文图像相减技术将图像与在天空的同一块上拍摄的参考图像进行比较，产生差异图像，并从海量数据中分析并识别瞬变天体、具有时变亮度的物体以及一些天文学现象。图3给出了基于重叠保存的频域天文图像差异算法的流程图。
所述的基于重叠保存的频域天文图像差异算法如下：
(1)对图像进行预处理：通过执行一个天体测量转换来匹配参考图像坐标系的坐标。本文通过选取参考坐标系和拍摄图像上的150颗恒星来拟合一个二维多项式，从而确定这个变换。利用这个变换，可以在参考图像定义的坐标系上对拍摄到的图像进行重新采样。这种重新采样是通过使用双三次样条法来进行插值，能够得到很高的精度。随后根据图4对拍摄到的天文图像按照重叠保存设计进行补0和分块。其中模板图像定义为T(x,y)，输入图像为I(x,y)，空间变卷积核函数为K(x,y)，(x,y)是图像的位置。
(2)其中，图4对重叠保存进行了图形描述。假设输入图像的长和宽相等且都为正整数N，每个子图像的长和宽相等且都为正整数n，每个子图像对应的卷积模板的长和宽也相等且都为正整数h,并令N′＝n+h-1。如图4(a)所示，根据空间可变卷积核的概念，先将输入图像划分为小的、不重叠的子图像n*n，因此可将图像划分成个子图像(/>表示向上取整，不足的补0)，随后根据重叠保存设计，对输入图像进行补0，将图像填充至/>大小。参数k1和k2可以被认为是移位参数，因为改变k1和k2都会选择不同的子矩阵。第一个输入子矩阵(参数为k1＝0,k2＝0)如图4(b)所示，其中N′*N′大小的矩形代表的是循环卷积，即补0后需要处理的部分。n*n大小的左斜线矩阵则表示所述的线性卷积，即每个子图像最后输出的部分。第二个输入子矩阵(参数为k1＝1,k2＝0)如图4(c)所示，每一行之间每个子块都从与前一个块重叠N′*(h-1)的区域开始，其中前/>的区域(大网格部分)表示利用前一段子矩阵数据替代补零的部分。第三个输入子矩阵(参数为k1＝0,k2＝1)如图4(d)所示，每一列之间每个子块都从与上一个块重叠(h-1)*N′的区域开始，前/>的区域(宽竖线部分)表示利用前一段子矩阵数据替代补零的部分。最后，在完成计算后，可以将每个子块的污染(重叠部分)丢弃，并将剩余的输出部分连接起来，这既确保沿信号操作的连续性，最大程度利用了图像自身的数据，并为避免混叠效应提供了必要的填充，提高了算法的计算效率。另外，扩充图像会引入高频分量，从而带来干扰，通常图像的频域的高频分量都集中在图像的四周，该重叠保存法将每个子块的重叠部分集中在四周，通过丢弃掉重叠部分从而消去了由高频分量引入的干扰，保证了图像的精度。这种方法的另一个优点是，它允许选择最佳子图像大小以实现FFT的最佳性能。图5给出了基于重叠保存法用FFT计算卷积结果的示意图。
(3)随后，依据I，T构建卷积核K，使用K对T卷积降晰，使T和I清晰度相同。
T(x,y)*K(x,y)＝I(x,y) (1)
卷积核可以被分解成一系列具有位置加权基函数，那么寻找最优解问题可以转换为线性最小二乘问题。由于卷积核可用高斯基函数同多项式相乘之后来精确地拟合。本发明采用高斯函数作为基函数：
其中an(x,y)和Kn分别是卷积核的系数和基函数。其下索引n由(i,j,k)共同表示，k是用于拟合的基函数的个数，σk表示高斯基函数的方差。公式(2)中的卷积核考虑空间位置对卷积核的影响，通过将积分域划分为矩形子图像。在每个小区域内，可以假定内核是恒定的。并用子图像的中心坐标代替(x,y)代入得到卷积核。因此，该方法即使在图像尺寸较大的情况下也能获得最优的结果。
(4)通过求解以下线性系统，可以计算出卷积核:
其中：
Mm1,m2＝∫Wm1(x,y)Wm2(x,y)dxdy (5)
(5)根据公式(1)可知，该运算是一个匹配滤波过程，根据时域卷积定理，可以用大规模图像的傅里叶变换与若干短的卷积模板的傅里叶变换的乘积来替代图像与卷积模板的时域卷积，因此公式(1)可以替换为：
Iij(x,y)＝IFFT(Tij(u,v)*Kij(u,v)) (7)
其中：Tij(u,v)＝FFT(Tij(x,y))，Kij(u,v)＝FFT(Kij(x,y))，Tij(x,y)表示按照重叠保存设计分割的子图像，Kij(u,v)表示子图像对应的卷积核。
(6)对拍摄到的图像和降晰后的模板图像作差，得到差异图像。
该算法的伪代码如图6所示，其中
1)初始化(第1行)：对输入图像和参考图像进行对准，并根据程序配置参数。
2)图像分块(第2-9行)：根据重叠保存的方法对输入图像
3)计算卷积核(第10-14行)：根据公式(4)(5)(6)求解随空和模板图像进行分块，不足的地方补0；间位置变化的卷积核的参数，并计算每个子区域对应的卷积核，为匹配两幅天文图像的PSF做准备；
4)频域滤波替代时域卷积(第20-26行)：根据改进后的重叠保存的方法对图像补0、选取子图像，并通过应用前向FFT、与若干预先计算的模板进行复数乘法，以及对结果进行逆FFT来执行匹配滤波；
5)结果输出：将输出中的重叠部分丢弃，并将剩下的输出部分按位相加形成最终输出。
基于重叠保存的频域天文图像差异算法在CPU-GPU异构平台上的实现方案如下：
(1)将CCD相机拍摄到的指定天区的天文图像数据拷贝到主机端中，对拍摄到的天文图像进行图像配准的预处理，随后根据重叠保存设计对图像进行补0输入图像和模板图像从磁盘存储到内存中。
(2)从主机端将数据复制到设备以存储设备中的数据：在这一步中，应该在设备中分配常量内存来存储模板图像、卷积图像、核系数和核向量。共享内存被分配用于存储卷积核及其系数，然后将数据从主机复制到设备。
(3)工作项分配：将数据从主机复制到设备后，模板图像需要分成多个子区域，每个子区域的大小由算法本身定义。每个子区域都被分派到一个工作组中。工作组中的每个工作项负责计算该子区域中的一个像素的卷积，每个工作组的大小为L*L，其中L＝n+h-1，n为子区域的大小,h为卷积核的尺寸。
(4)计算内核并存储到共享存储器中：在计算傅里叶变换之前，必须计算出子区域的空间变化卷积核。工作组负责计算分配给它的子区域的卷积核，并将卷积核存储到共享内存中。所以在计算子区域和子区域对应的卷积核的傅里叶变换时，工作组中的每个工作项都可以访问卷积核。
(5)通过barrier(CLK_LOCAL_MEM_FENCE)命令同步：同步工作组中的工作项，确保卷积核的计算已经完成。
(6)频域降晰：按照重叠保存设计，对每个子区域和卷积核进行傅里叶变换，然后对计算结果相乘，并对矩阵相乘结果进行逆傅里叶变换。这一部分的数据是相互独立的，可以并行计算，工作组中的每个线程都负责计算子区域内和卷积核的一个蝶形计算单元。
(7)同步：确保所有的子区域都已经完成了频域降晰计算。
(8)将降晰后图像从设备复制到主机：计算完参考图像的所有子区域后，可以将频域傅里叶变换降晰后的图像从设备复制到主机。
(9)释放设备中分配的空闲内存：释放计算过程中所需要的所有内存，避免由于占用过多的显存而导致启动新的内核函数失败。
(10)图像相减：在最后一步中，输入图像与频域降晰后的模板图像相减得到结果图像。
从基于重叠保存的频域天文图像差异算法中可以看出，其主要算子为矩阵乘法和FFT，其中FFT是整个算法的重点。本发明主要介绍快速傅里叶变换以及在OpenCL上的实现和优化。
(1)确定二维快速傅里叶变换的分解方式：
通常计算二维FFT的方法是通过一维FFT变换计算二维图像所有的行，然后通过一维FFT变换计算二维图像所有的列。但是，对于很多内存系统，在对所有的列进行FFT的过程中可能会出现性能问题。这是因为大多数内存系统的存储方式都是行优先存储的方式，这样就导致在读取列数据时造成了延迟。从行时间与列时间之间的巨大差异可以发现另外一种计算二维FFT的方法：
1.使用1维FFT变换计算二维图像所有的行。
2.对二维图像进行转置。
3.使用1维FFT变换计算二维图像所有的行。
4.对二维图像进行转置。
(2)减少数据的交互
在计算中，N点FFT一般可以分成M级，其中，考虑到FFT的多通结构，通常可以使用局部内存将一部分需要交换的数据存储在局部内存中，而Mali T860并没有提供专用的局部存储器和独立的常量存储器，那么每一通都需要与全局内存进行数据交换，大大降低了算法的性能。因此必须通过使用多基FFT算法，采用长度更长的FFT来尽早地停止递归，从而减少数据与全局内存的交换。
(3)旋转因子
FFT中使用的旋转因子是将多个较小的傅里叶变换合并成一个较大的傅里叶变换的关键，其中旋转因子的形式为(i是虚数单位，N是FFT的长度)。为了优化其性能，一般可通过将旋转因子存储到常量内存中，在计算时直接调用常量内存来减少计算量和数据读取的延迟。由于ARM Mali系列GPU并没有提供常用内存，但Mali系列GPU提供了图像(纹理)缓存，可以通过将旋转因子数据存放在纹理缓存中以供计算调用，此时必须使用图像对象来访问旋转因子。
实施例2
本发明的实施例2提出了一种频域天文图像目标检测系统，所述系统包括：部署在CPU的预处理模块、基函数计算模块和目标检测模块，以及部署在GPU的卷积核计算模块和频域滤波降晰模块，其中，
所述预处理模块，用于基于预先得到的参考图像，对采集的原始天文图像进行预处理；还用于采用重叠保存的方法对参考图像和预处理后的天文图像分别进行图像分块，得到尺寸和个数均相同的参考图像子图像和天文图像子图像；
所述基函数计算模块，用于采用高斯基函数和多项式相乘得到卷积核的n组基向量并输入卷积核计算模块；
所述卷积核计算模块，用于根据参考子图像和天文子图像，结合n组基向量进行拟合，得到对应每个参考子图像的卷积核；
所述频域滤波降晰模块，用于使用卷积核对每个参考子图像进行频域滤波降晰处理，得到模板图像并输入目标检测模块；
所述目标检测模块，用于将模板图像的边缘丢弃，并将剩余部分连接并与原始天文图像作差，得到差异图像，进而实现对天文图像的目标检测。
本发明的效果可通过以下仿真进一步阐述。
1.仿真条件：
测试仿真使用的CPU为一个基于big.little大小核架构的6核64位处理器(双核Cortex-A72，最高频率1.8GHz，四核Cortex-A53，最高频率1.4GHz)，操作系统采用Fedora28。测试仿真使用的GPU主要采用ARM Mali-T860MP4，该GPU支持OpenGL ES1.1/2.0/3.0/3.1，OpenVG1.1，OpenCL，DX11，AFBC(帧缓冲压缩)。其具有4个渲染核心，每个核心有4个计算单元、一个标量单元和一个纹理单元。每个向量计算单元是一个128位的乘加单元。单精度浮点计算能力＝4*(4*4+1)*2*主频。在存储器方面，所有渲染核共享一个统一的256KB的二级缓存，二级缓存分成两块，一共可提供128位访存带宽。在内存方面，采用双通道内存LPDDR3 4GB，每个内存通道带宽为64位，其内存带宽＝LPDDR等效频率*64/8*2。没有独立的常量存储器缓存，没有专用的局部存储器硬件。
本发明所使用的天文图像处理软件是最新版本的HOTPANTS v5.1.11。其中包括：(1)使用高斯基函数；(2)在每个stamp中选择三个子图像；(3)使用卷积核全部元素的和(kernel sum)作为Sigma剔除异常值的阈值；(4)分别对输入图像和模板生成卷积核，然后从两者中选择平均像素残差较低的一个。另外，在HOTPANTS中，本文还使用以下预定义参数及参数值，如表1所示。GCC版本4.4.5和NVCC版本V0.2.1221，采用OpenCL1.2异构并行计算标准。
表1预定义的参数表
2.仿真内容：
仿真1：在上述仿真条件下，分别对本发明和常规天文图像差异方法进行了算法复杂度和内存占用的仿真。
在计算复杂度方面，为了通过频域傅里叶变换计算时域卷积，需要对每一个子区域进行3次FFT：一次是子图像的傅里叶变换，一次是子图像对应的卷积核的傅里叶变换，一次是两者乘积的傅里叶逆变换。而二维图像的傅里叶变换的复杂度为Offt＝MNlog2(MN)，其中M为图像的长，N为图像的宽。此外，我们还必须考虑卷积核和图像的FFT的逐元乘法。在内存占用方面，改进后的算法以及时域卷积的内存消耗主要包括卷积内存和存储内存两个部分。由于基于重叠保留的频域天文图像差异算法的输出重叠部分直接丢弃，且传统的时域卷积没有重叠部分，因此两算法的存储内存均为零。因此，对于传统卷积和频域滤波算法，最占用内存的步骤都是对每个卷积核、子图像和最终计算结果的存储。假设整个流程使用双精度64bit格式表示数据，对于卷积核随空间变换的情况，传统时域卷积和改进后的频域滤波的计算复杂度和内存占用如下：
不同规模图像下，不同卷积核大小对基于重叠保存的天文图像差异算法与OIS算法的计算复杂度和内存占用的影响如图7(a)和图7(b)所示。
仿真2：在上述仿真条件下，用本发明和常规天文图像差异算法分别对四幅真实的天文图像进行处理，其运行时间对比结果如表2所示。
表2天文图像差异算法在不同计算平台上的处理时间
从表2可以看出，在CPU计算平台处理天文图像(1Kx1K到4Kx4K)时，相比于OIS算法，改进后的算法的运行速度明显加快。对于4Kx4K的图像，基于重叠保存的频域天文图像差异算法的处理时间是31.5秒，而原来的OIS算法则需要46秒。与原天文图像差异算法相比，改进后的频域天文图像差异算法的平均性能提高了17.6％到31.5％。而在CPU-GPU异构嵌入式计算平台上，改进后的频域天文图像差异算法的处理时间只需要14.7秒，改进后的算法的平均性能提高了1.77到2.14倍，验证了本发明对天文图像差异算法的加速处理效果。另外随着图像尺寸的增加，改进后算法的加速比随之增大，这是由于随着输入图像尺寸的增加，并行部分耗时占总耗时的比例也随之增加，给算法提供了更多的优化空间。
本发明主要创新点如下：
1)在复杂度方面：卷积降晰部分是整个系统的性能瓶颈，也是制约整个系统满足实时性的关键因素。因此本文将原算法中最耗时的图像卷积部分转换到频域中，利用频域滤波来代替时域的卷积运算，优化了天文图像差异算法中的性能瓶颈，降低算法的计算复杂度，缩短了算法的运行时间。
2)在内存占用、计算效率方面：通过重叠保存设计保留分块图像一定位数的原图像数据来延拓图像，在频域滤波完成后再将这些重叠部分舍弃，再按位相加，最大限度地利用图像数据。该设计既解决了卷积的循环特性产生频谱混叠问题，又避免了传统补零方法导致的内存占用过大的问题，保证了图像的连续处理，降低了算法的内存占用，提高了算法的计算效率。
3)该算法的重叠保存设计将每个子块的重叠部分分布在图像边缘，通过丢弃重叠部分消除了由扩充图像引入高频分量导致的干扰，保证了算法的精度。
在计算资源分配方面，本文的算法加速方案充分利用了CPU和GPU的计算资源，将GPU用于卷积核的确定以及频域滤波计算，将图像卷积工作的高斯基函数Kn的计算分配给多核CPU并行处理，让CPU真正地参与到计算中，实现了异构平台的并行计算的效率提升。
通过对算法的设计和优化，使算法能够满足天文图像在星载设备上进行实时处理的约束。在此CPU-GPU异构并行计算系统内，基于重叠保存的频域天文图像差异算法既可以满足深空探测任务对于快速、实时的要求，也可以满足星载计算机低内存、低功耗的目标，为实现航天器的在轨自适应科学探测提供了可能。
最后所应说明的是，以上实施例仅用以说明本发明的技术方案而非限制。尽管参照实施例对本发明进行了详细说明，本领域的普通技术人员应当理解，对本发明的技术方案进行修改或者等同替换，都不脱离本发明技术方案的精神和范围，其均应涵盖在本发明的权利要求范围当中。

Claims (5)

1.一种频域天文图像目标检测方法，基于CPU-GPU异构处理器实现，所述方法包括：
基于预先得到的参考图像，由CPU对采集的原始天文图像进行预处理；采用重叠保存的方法对参考图像和预处理后的天文图像分别进行图像分块，得到尺寸和个数均相同的参考图像子图像和天文图像子图像；采用高斯基函数和多项式相乘得到卷积核的n组基向量并输入GPU；
由GPU根据参考子图像和天文子图像，结合n组基向量进行拟合，得到对应每个参考子图像的卷积核，使用卷积核对每个参考子图像进行频域滤波降晰处理，得到模板图像并输入CPU；
由CPU将模板图像的边缘丢弃，并将剩余部分连接并与原始天文图像作差，得到差异图像，进而实现对天文图像的目标检测；
所述预处理前后的天文图像和参考图像的长宽尺寸均为N*N；
所述采用重叠保存的方法对参考图像和预处理后的天文图像分别进行图像分块，得到尺寸和个数均相同的参考图像子图像和天文图像子图像；具体包括：
设置卷积模板的长和宽均为正整数h；
根据空间可变卷积核的概念，对参考图像和预处理后的天文图像均进行以下相同操作：
按照尺寸为n*n进行划分，得到个不重叠的子图像，/>表示向上取整，不足的区域补0，随后根据重叠保存方法对子图像进行补0，将图像填充至大小；
所述使用卷积核对每个参考子图像进行频域滤波降晰处理，得到模板图像并输入CPU；具体包括：
按照重叠保存方法，并行的对每个参考子图像和对应的卷积核进行傅里叶变换，然后对计算结果相乘，并对矩阵相乘结果进行逆傅里叶变换，实现频域滤波降晰处理，得到模板图像并输入CPU。
2.根据权利要求1所述的频域天文图像目标检测方法，其特征在于，所述基于预先得到的参考图像，由CPU对采集的原始天文图像进行预处理；具体包括：
由CPU从预先得到的参考图像中选取参考坐标系，通过分别在参考图像和输入图像上选取相同数量的天体以确定二维多项式形式的变换，利用该变换，在参考坐标系上对原始天文图像进行重新采样，实现与参考图像的对准；所述重新采样采用双三次样条法进行插值处理。
3.根据权利要求2所述的频域天文图像目标检测方法，其特征在于，所述采用高斯基函数和多项式相乘得到卷积核的n组基向量；具体包括：
使用三个高斯函数作为卷积核K的基函数，σk表示第k个基函数的方差，k＝1,2,3，(x,y)表示每个子图像的坐标，令
σ1＝0.7；i+j≤6
σ2＝1.5；i+j≤4
σ3＝3.0；i+j≤2
其中，i和j表示基函数的空间自由度，根据i和j的组合，得到n组基向量Kn(x,y)为：
4.根据权利要求3所述的频域天文图像目标检测方法，其特征在于，所述由GPU根据参考子图像和天文子图像，结合n组基向量进行拟合，得到对应每个参考子图像的卷积核，具体包括：
令∑(T(x,y)*K(x,y)-I(x,y))2最小，以找到图像减法的内核，其中，T(x,y)表示参考子图像，I(x,y)表示天文子图像，K(x,y)为待求参考子图像的卷积核，(x,y)为每个子图像的坐标，表示内核的空间变化特征，基于每个子图像内点扩散函数恒定的设定条件，用子图像的中心坐标表示(x,y)；将最优化问题转换为线性最小二乘问题，求解以下线性系统，依据参考子图像和天文子图像得到对应每个参考子图像的卷积核的系数
其中，M为最小二乘的元素构成的矩阵，为向量Wm1和向量Wm2的标量积，系数的第n个元素为an(x,y)，/>为向量Wm1与I(x,y)的标量积，矩阵M的第m1m2个元素Mm1,m2满足下式：
Mm1,m2＝∫Wm1(x,y)Wm2(x,y)dxdy
其中，m1,m2＝0,1,2,…,n；
的第m1个元素Bm1满足下式：
采用n组基向量Kn(x,y)和系数an(x,y)进行拟合，得到卷积核K(x,y)。
5.一种基于权利要求1的频域天文图像目标检测方法的检测系统，其特征在于，所述系统包括：部署在CPU的预处理模块、基函数计算模块和目标检测模块，以及部署在GPU的卷积核计算模块和频域滤波降晰模块，其中，
所述预处理模块，用于基于预先得到的参考图像，对采集的原始天文图像进行预处理；还用于采用重叠保存的方法对参考图像和预处理后的天文图像分别进行图像分块，得到尺寸和个数均相同的参考图像子图像和天文图像子图像；
所述基函数计算模块，用于采用高斯基函数和多项式相乘得到卷积核的n组基向量并输入卷积核计算模块；
所述卷积核计算模块，用于根据参考子图像和天文子图像，结合n组基向量进行拟合，得到对应每个参考子图像的卷积核；
所述频域滤波降晰模块，用于使用卷积核对每个参考子图像进行频域滤波降晰处理，得到模板图像并输入目标检测模块；
所述目标检测模块，用于将模板图像的边缘丢弃，并将剩余部分连接并与原始天文图像作差，得到差异图像，进而实现对天文图像的目标检测。
CN202110527356.0A 2021-05-14 2021-05-14 一种频域天文图像目标检测方法及系统 Active CN113344765B (zh)

Priority Applications (1)

Application Number Priority Date Filing Date Title
CN202110527356.0A CN113344765B (zh) 2021-05-14 2021-05-14 一种频域天文图像目标检测方法及系统

Applications Claiming Priority (1)

Application Number Priority Date Filing Date Title
CN202110527356.0A CN113344765B (zh) 2021-05-14 2021-05-14 一种频域天文图像目标检测方法及系统

Publications (2)

Publication Number Publication Date
CN113344765A CN113344765A (zh) 2021-09-03
CN113344765B true CN113344765B (zh) 2023-11-03

Family

ID=77469995

Family Applications (1)

Application Number Title Priority Date Filing Date
CN202110527356.0A Active CN113344765B (zh) 2021-05-14 2021-05-14 一种频域天文图像目标检测方法及系统

Country Status (1)

Country Link
CN (1) CN113344765B (zh)

Families Citing this family (1)

* Cited by examiner, † Cited by third party
Publication number Priority date Publication date Assignee Title
CN115774830B (zh) * 2022-10-28 2023-09-29 贵州大学 一种用于稀疏干涉阵列的快速成图方法

Citations (2)

* Cited by examiner, † Cited by third party
Publication number Priority date Publication date Assignee Title
CN107451955A (zh) * 2017-06-20 2017-12-08 昆明理工大学 一种k‑t算法重建天文图像中斑点图的并行化实现方法
CN112101440A (zh) * 2020-09-09 2020-12-18 昆明理工大学 一种基于改进U-Net网络的日冕物质抛射检测方法

Patent Citations (2)

* Cited by examiner, † Cited by third party
Publication number Priority date Publication date Assignee Title
CN107451955A (zh) * 2017-06-20 2017-12-08 昆明理工大学 一种k‑t算法重建天文图像中斑点图的并行化实现方法
CN112101440A (zh) * 2020-09-09 2020-12-18 昆明理工大学 一种基于改进U-Net网络的日冕物质抛射检测方法

Non-Patent Citations (2)

* Cited by examiner, † Cited by third party
Title
基于OpenMP的并行图像相减算法实现与分析;李继良;于策;孙济洲;商朝晖;陈锦言;曹玮;张旭明;;天文研究与技术(第02期);56-62 *
天文图像序列中弱目标的实时检测算法;韩建涛;张月;陈曾平;;光电工程(第12期);4-7 *

Also Published As

Publication number Publication date
CN113344765A (zh) 2021-09-03

Similar Documents

Publication Publication Date Title
US10467763B1 (en) System and method for optical flow estimation
US20230410375A1 (en) Temporally stable data reconstruction with an external recurrent neural network
US10922793B2 (en) Guided hallucination for missing image content using a neural network
CN109472858B (zh) 用于逆向图形的可微分渲染管线
CN106991665B (zh) 基于cuda图像融合并行计算的方法
Ma et al. Parallel programing templates for remote sensing image processing on GPU architectures: design and implementation
US20230419611A1 (en) Reducing level of detail of a polygon mesh to decrease a complexity of rendered geometry within a scene
US20230237313A1 (en) Layout Parasitics and Device Parameter Prediction using Graph Neural Networks
CN109859178B (zh) 一种基于fpga的红外遥感图像实时目标检测方法
DE102018114929A1 (de) System und Verfahren für augennahes Lichtfeldrendering für eine interaktive dreidimensionale Computergraphik mit weitem Blickfeld
US11494879B2 (en) Convolutional blind-spot architectures and bayesian image restoration
CN103019671A (zh) 面向数据密集型遥感图像处理的泛型编程的框架编程方法
CN111539997B (zh) 基于gpu计算平台的图像并行配准方法、系统、装置
EP3678037A1 (en) Neural network generator
Fang et al. CPU/GPU near real-time preprocessing for ZY-3 satellite images: Relative radiometric correction, MTF compensation, and geocorrection
CN113344765B (zh) 一种频域天文图像目标检测方法及系统
Poostchi et al. Efficient GPU implementation of the integral histogram
Rahman et al. Parallel implementation of a spatio-temporal visual saliency model
CN113822975B (zh) 用于对图像进行有效采样的技术
Jensen et al. A two-level real-time vision machine combining coarse-and fine-grained parallelism
Lin et al. The parallelization of video processing
Lei et al. Stream model-based orthorectification in a GPU cluster environment
DE102019101845B4 (de) Photorealistische Bildstilisierung unter Verwendung eines Neuronalen-Netzwerkmodells
Kyo et al. An integrated memory array processor for embedded image recognition systems
CN116127685A (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