CN107451955A - 一种k‑t算法重建天文图像中斑点图的并行化实现方法 - Google Patents
一种k‑t算法重建天文图像中斑点图的并行化实现方法 Download PDFInfo
- Publication number
- CN107451955A CN107451955A CN201710466806.3A CN201710466806A CN107451955A CN 107451955 A CN107451955 A CN 107451955A CN 201710466806 A CN201710466806 A CN 201710466806A CN 107451955 A CN107451955 A CN 107451955A
- Authority
- CN
- China
- Prior art keywords
- sub
- block
- spot
- image
- gpu
- 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
Links
- 238000000034 method Methods 0.000 title claims abstract description 95
- 230000008569 process Effects 0.000 claims abstract description 43
- HPTJABJPZMULFH-UHFFFAOYSA-N 12-[(Cyclohexylcarbamoyl)amino]dodecanoic acid Chemical compound OC(=O)CCCCCCCCCCCNC(=O)NC1CCCCC1 HPTJABJPZMULFH-UHFFFAOYSA-N 0.000 claims abstract description 18
- 238000012545 processing Methods 0.000 claims abstract description 9
- 238000001228 spectrum Methods 0.000 claims description 40
- 230000006870 function Effects 0.000 claims description 38
- 230000015654 memory Effects 0.000 claims description 14
- 230000005540 biological transmission Effects 0.000 claims description 12
- 230000000694 effects Effects 0.000 claims description 10
- 238000005305 interferometry Methods 0.000 claims description 8
- 230000000007 visual effect Effects 0.000 claims description 6
- 230000007246 mechanism Effects 0.000 abstract description 4
- 238000005516 engineering process Methods 0.000 abstract description 3
- 230000009466 transformation Effects 0.000 description 8
- 238000004364 calculation method Methods 0.000 description 6
- 238000010586 diagram Methods 0.000 description 4
- 238000003860 storage Methods 0.000 description 4
- 230000008901 benefit Effects 0.000 description 3
- 238000009826 distribution Methods 0.000 description 3
- 208000002173 dizziness Diseases 0.000 description 3
- 239000000203 mixture Substances 0.000 description 3
- 238000005457 optimization Methods 0.000 description 3
- 230000001360 synchronised effect Effects 0.000 description 3
- 238000012546 transfer Methods 0.000 description 3
- 230000001133 acceleration Effects 0.000 description 2
- 230000009286 beneficial effect Effects 0.000 description 2
- 230000004087 circulation Effects 0.000 description 2
- 238000004891 communication Methods 0.000 description 2
- 235000013399 edible fruits Nutrition 0.000 description 2
- 238000001914 filtration Methods 0.000 description 2
- 239000011159 matrix material Substances 0.000 description 2
- 230000011218 segmentation Effects 0.000 description 2
- 238000013519 translation Methods 0.000 description 2
- 241000208340 Araliaceae Species 0.000 description 1
- 240000007594 Oryza sativa Species 0.000 description 1
- 235000007164 Oryza sativa Nutrition 0.000 description 1
- 235000005035 Panax pseudoginseng ssp. pseudoginseng Nutrition 0.000 description 1
- 235000003140 Panax quinquefolius Nutrition 0.000 description 1
- 235000013339 cereals Nutrition 0.000 description 1
- 238000006243 chemical reaction Methods 0.000 description 1
- 230000021615 conjugation Effects 0.000 description 1
- 230000001186 cumulative effect Effects 0.000 description 1
- 238000005520 cutting process Methods 0.000 description 1
- 238000013461 design Methods 0.000 description 1
- 206010016256 fatigue Diseases 0.000 description 1
- 238000009432 framing Methods 0.000 description 1
- 230000008014 freezing Effects 0.000 description 1
- 238000007710 freezing Methods 0.000 description 1
- 235000008434 ginseng Nutrition 0.000 description 1
- 238000003384 imaging method Methods 0.000 description 1
- 230000006872 improvement Effects 0.000 description 1
- 238000011835 investigation Methods 0.000 description 1
- 238000007781 pre-processing Methods 0.000 description 1
- 230000002035 prolonged effect Effects 0.000 description 1
- 230000011514 reflex Effects 0.000 description 1
- 230000010076 replication Effects 0.000 description 1
- 235000009566 rice Nutrition 0.000 description 1
- 239000011435 rock Substances 0.000 description 1
- 230000003595 spectral effect Effects 0.000 description 1
- 230000003313 weakening effect Effects 0.000 description 1
Classifications
-
- 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/4038—Image mosaicing, e.g. composing plane images from plane sub-images
-
- 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
- G06T—IMAGE DATA PROCESSING OR GENERATION, IN GENERAL
- G06T2200/00—Indexing scheme for image data processing or generation, in general
- G06T2200/32—Indexing scheme for image data processing or generation, in general involving image mosaicing
-
- 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/30181—Earth observation
- G06T2207/30192—Weather; Meteorology
Landscapes
- Physics & Mathematics (AREA)
- General Physics & Mathematics (AREA)
- Engineering & Computer Science (AREA)
- Theoretical Computer Science (AREA)
- Image Processing (AREA)
Abstract
本发明涉及一种K‑T算法高分辨重建天文图像中斑点图的并行化实现方法，属天文技术图像处理计算领域。本发明首先通过CPU主进程读取图像及参数文件并且计算相关物理参数进行图像的预处理和数据对齐，接着把全部帧数的整个视场图像依据等晕区大小分割成很多的子块进行图像分块，然后这些子块帧被分发到独立的子进程中采用CUDA并行编程模型进行并行计算处理，进行傅里叶的振幅和相位重建，最后子进程对每一个子块通过联合傅里叶振幅和相位进行傅里叶逆变换获得重建后图像，主进程对所有的子块拼接处理完成重建。本发明将CPU+GPU混合异构并行模式并结合GPU‑aware MPI机制应用于高分辨重建天文图像中斑点图问题上，解决了对斑点图高分辨率重建的问题，具备较高的扩展性。
Description
技术领域
本发明涉及一种K-T算法高分辨重建天文图像中斑点图的并行化实现方法，属天文技术、图像处理和分布式并行计算领域。
背景技术
天文图像高分辨率频域重建的开创性贡献是拉贝里（Labeyrie）作出的。Labeyrie法也常被称为斑点干涉术(speckle interferometry)在此之前，人们已经认识到天文图像的长曝光像几乎失去了所有的高频信息，短曝光图像也因为波前畸变而破碎变形。拉贝里仔细研究了短曝光图像的斑点结构后指出，短曝光天文图像中含有可以利用的高频信息。
天文目标发出的光通过大气层时，由于空气的温度和湿度不均匀产生的折射率不均匀导致了光波的振幅和相位的随机扰动，严重降低了望远镜的成像质量。要使成像分辨力达到望远镜的衍射极限就必须克服大气湍流的影响。“天文图像高分辨率重建”希望解决的是由地球大气造成的天文图像像质衰减。频域重建所依据的理论基础均是湍流大气模型以及线性空间不变综合系统模型，重建法使用的数学工具是频域变换和统计矩，统计的样本就是被称为“斑点图”的短曝光天文图像。
当曝光时间很短时，由于大气的“冻结”效应，波前的相位起伏会保持一个相对稳定的空间分布，因此，虽然高频信噪比很低，但斑点图中还是保留了很多可以提取的高频信息。如果延长曝光时间，大气将不能保持恒定的相位起伏，从而导致各频率条纹的初相位随机变化，最终丢失高频信息。
当重建对象为一个扩展目标（如太阳表面）或是在一个空间范围成一定分布的一组目标时（典型如密集星团），必须关心大气-望远镜综合系统传递函数的有效空间范围。当考察相隔α（角距离）的两个天文点源目标A和B的瞬间波前起伏几乎同步，甚至结构相同。此时综合系统对目标有相同的瞬时传递函数，它使在描述一个瞬间样本时也成立，就称这样的综合系统为线性空间不变系统。而能保证满足线性空不变系统基本成立的α就称为等晕区或等晕角。
为了重建目标相位，Knox和Thompson于1974年提出了一种基于计算频谱的互相关（又称为互谱）来获取目标相位的方法，称为Knox-Thompson (简称K-T)法。太阳斑点成像术的执行和实现过程分为以下几个步骤：图像预处理、数据对齐、图像分块、振幅和相位重建以及图像拼接。其中振幅重建由斑点干涉术完成，而相位重建可由频域算法K-T法或斑点掩模法完成，亦可在空域中由选帧位移叠加法完成。
新型异构混合体系结构对大规模并行算法研究提出了新的挑战，迫切需要深入研究与该体系结构相适应的并行算法。针对MPI+GPU异构混合体系结构的高性能计算平台，研究相应的协同并行计算技术，设计并实现大型科学及工程计算问题的新型并行算法，具有重大的理论和实际意义。
发明内容
本发明提供了一种K-T算法高分辨重建天文图像中斑点图的并行化实现方法，将CPU+GPU混合异构并行模式结合GPU-aware MPI机制应用于高分辨重建天文图像中斑点图。同时，具备较高的扩展性，可根据硬件节点数的增加提高重建效率。
本发明K-T算法高分辨重建天文图像中斑点图的并行化实现方法是这样实现的：
首先通过CPU主进程（master）读取图像及参数文件并且计算相关物理参数进行图像的预处理和数据对齐；接着把全部帧数的整个视场图像依据等晕区大小分割成很多的子块进行图像分块；然后这些子块帧被分发到独立的子进程（slave）中去采用CUDA并行编程模型进行并行计算处理，进行傅里叶的振幅和相位重建；最后子进程对每一个子块通过联合傅里叶振幅和相位进行傅里叶逆变换获得各个子块的重建结果，主进程对所有子块的重建结果进行拼接得到一帧整个视场的重建图像，保存相关数据后完成重建。
所述K-T算法高分辨重建天文图像中斑点图的并行化实现方法的具体方法如下：
主进程中对全视场图像分块，对待处理的全部帧数张的斑点图像依据等晕区大小进行整体分块，为了减少边缘效应对最后重建结果的影响，采取对各个子块50%重叠的方式进行分割，即分块的时候子块之间有交叠，同时为了减少重建结果中切趾效应的影响采用海明窗和汉宁窗对图像进行加窗处理。
子块的振幅重建在频域中由斑点干涉术退卷积传递函数完成，在子进程中由GPU多线程并行计算，首先通过对系列短曝光图像进行谱比计算，然后将计算所得谱比和标准谱比进行比对得到大气相干长度r0（Fried参数，又称视宁度参数）。斑点传递函数STF（Speckle Transfer Function）用Korff的理论模型代替，在Korff模型中STF只是r0的函数，由r0得出STF。最后退卷积STF后得出傅里叶振幅。
子块的相位重建由扩展Knox-Thompson法计算得出，将频谱中相隔较小频域距离的相位相减来消除由大气带来的相位起伏，留下的部分是这对频率的相位差。有了整个频域面上的相位差，通过由低频向高频的递推实现相位重建。在子进程中由GPU多线程并行计算K-T交叉谱、振幅加权整合、误差最小化迭代。
将CPU+GPU混合异构并行模式结合GPU-aware MPI机制进行并行优化，MPI实现任务分配及数据的发送和接受，使CPU主机到GPU设备的直接通信。CPU与GPU协同工作，并执行各自的计算程序。GPU加速计算密集且耗时的代码，CPU则执行串行代码。主要以计算密集度高、循环迭代次数多、需要多次调用的函数进行并行化为kernel核函数。对于算法中傅里叶变换使用CUFFT库函数，矩阵和向量的运算使用CUBLAS库函数。
本发明的有益效果是：
K-T算法高分辨重建天文图像中斑点图的并行化实现方法，将CPU+GPU混合异构并行模式并结合GPU-aware MPI机制应用于高分辨重建天文图像中斑点图问题上是前所未有的、有效的，该方法在很大程度上解决了对斑点图高分辨率重建的问题。同时，具备较高的扩展性，可根据硬件节点数的增加提高重建效率。
附图说明
图1是本发明K-T算法高分辨重建天文图像中斑点图的并行化实现方法的总体流程图；
图2是本发明通过主进程对100帧待重建图像分块示意图；
图3是本发明通过主进程对子块的重建结果进行拼接示意图；
图4是本发明中并行策略下主进程与子进程执行任务示意图；
图5是本发明中待重建的一帧短曝光的太阳斑点图；
图6是本发明中的方法重建后的结果图。
具体实施方式
实施例：如图1-6所示，一种K-T算法高分辨重建天文图像中太阳斑点图的并行化实现方法，图像重建过程以此为开始，即首先通过主进程（master）读取待重建的云南澄江抚仙湖ONSET望远镜在白光波段的一组太阳局部像观测数据原始100帧短曝光斑点图及参数文件并且计算相关物理参数，例如衍射极限、图像尺寸等。接着继续把整个视场图像分割成很多的子块以克服非等晕性问题，其中子块的大小依据等晕区划分，这个大小可以在参数文件里面进行修改调整。每个子块的频域拟合变换通过调用快速傅里库CUFFT库函数处理。这些子块帧接着被分发到独立的子进程（slave）中去计算，所使用的子进程数量也可以通过MPI的运行程序命令来调整（例如mpirun –np 16 ./task）。在这些子进程中，各子块斑点图通过和帧数平均的互相关来计算他们各自的子块移位。计算所得的信息发送回主进程，被用于对整个视场的图像进行相关跟踪帧数的分割。下一步是计算径向谱比并且发送给主进程。主进程可以通过谱比估算出Fried参数并且确定出合适的模型库以外的斑点传递函数（STF）。在子进程中计算傅里叶振幅的平均值后同样也发回主进程。接着，为了估算傅里叶相位，将频谱中相隔较小频域距离的相位相减来消除由大气带来的相位起伏，留下的部分是这对频率的相位差。有了整个频域面上的相位差，通过由低频向高频的递推实现相位重建。在子进程中由GPU多线程并行计算K-T交叉谱、振幅加权整合、误差最小化迭代。估算出的傅里叶相位和计算出的维纳滤波（Wiener filter）也要发回主进程。这个过程一直重复直到所有的子块都被处理完。主进程等待直到所有的子块完成，然后用之前计算的Fried参数来估算，接着对傅里叶振幅进行噪声过滤并且用STF来校正。对每一个子块通过联合傅里叶振幅和相位进行傅里叶逆变换获得各个子块的重建结果。主进程对所有子块的重建结果进行拼接得到一帧整个视场的重建图像，保存相关数据后完成重建。
结合图1，所述K-T算法高分辨重建太阳斑点图的并行化实现方法的步骤如下：
步骤1、通过CPU主进程（master）读取100帧图像及参数文件并且计算相关物理参数进行图像的预处理数据对齐和图像分块；
步骤2、通过CPU主进程把子块帧分发到独立的子进程（slave）中去采用CUDA并行编程模型进行并行计算处理，进行傅里叶的振幅和相位重建，把相关信息发送回主进程；
步骤3、主进程对所有子块的重建结果进行拼接得到一帧整个视场的重建图像，保存相关数据后完成重建。
在步骤1中，重建信息包括重建参数信息和100帧图像合成的数据本身。具体步骤如下：
首先读取配置信息包括文件名、保存文件名、图像尺寸、角秒/像素、望远镜直径、波长等，根据读取配置信息时获取的图像数据路径读取图像数据，根据读取的配置信息确定子块的大小sfs为5，确定子块,行列大小ssizex、ssizey为5/0.042。为了消除望远镜跟踪误差或者风载而引起的图像晃动，并减弱近地面湍流导致的波前倾斜而引起的图像偏移进行数据初对齐。频域重建算法在等晕假设下才成立，所以在完成对齐后全视场的图像需要分成若干个子块，每个子块大小与等晕区大小相当。对ssizex循环除以2直到等于0，得到循环次数pot，重新确定子块大小为2的pot次幂（取整为2的幂大小过程）。
接着根据子块大小确定每张图划分的子块个数nfrx，nfry为15。为了减少边缘效应对最后重建结果的影响，采取对各个子块重叠50%的方式进行分割，即分块的时候子块之间有交叠，总共划分225个子块。
然后计算抵消边界offx，offy，此处计算为0。主进程根据子块大小个数和抵消边界将图像分割分配给子进程，切割方式为根据进程号确定一张图的子块位置p，然后将100张图相对应的子块分配给该进程。其中，子块大小除以2是因为分配给子进程重建的子块有一半重叠。
最后将重建相关参数(info)分配给子进程，对主进程发过来的数据进行计算100张子块中相应每个像素的平均值和相关移位数据。为了减少重建结果中切趾效应的影响采用海明窗和汉宁窗对图像进行加窗处理。同时计算平均值的傅里叶变换temp1->fref，这里采用CUFFT库的双精度的实数到复数的傅里叶变换；计算原始图像的傅里叶变换temp2->fimage，这里采用CUFFT库的双精度的实数到复数的傅里叶变换；计算以上两者的共同傅里叶变换fccr->ccr，这里采用CUFFT库的双精度的复数到实数的傅里叶变换。将计算的结果发给主进程，主进程根据子进程发送过来的数据，进行对齐分块完成后，还需将子块图像乘上窗函数，以保证边缘部分强度是缓变的，这样能有效抑制边缘效应。之后重新发给相应子进程进行后续的重建过程，子进程根据分配到的图像数据(data)和参数(sl_info)进行重建。
在步骤2中，采用CUDA并行编程模型并行计算重建傅里叶的振幅和相位。具体步骤如下：
振幅重建基于斑点干涉术（Labeyrie法）的方法。在满足等晕性假设下，每帧短曝光图像i(x) 可以看作是目标o(x)和系统点扩展函数h(x)的卷积：
， （1）
这里，x 是二维空间坐标变量，*表示卷积算符，在频域上式满足：
（2）
I(q)、O(q)和H(q)分别表示i(x)、o(x)和h(x)的傅里叶变换，q是二维空间频率变量。斑点干涉术，是通过统计斑点图的平均能谱来复原目标傅里叶振幅。斑点图的能谱统计结果为：
（3）
其中括号 表示算术平均， 被称为斑点干涉术传递函数(STF)。
1983年由O. von der Luhe提出的谱比法由于其计算简单、精度高被广泛用于斑点图像重建中大气相干长度（r0）的计算。这里使用谱比法计算谱比，公式为：
（4）
在子进程中由GPU多线程并行计算，通过对系列短曝光图像进行上述公式的谱比计算，然后将计算所得谱比和标准谱比进行比对得到大气相干长度r0。斑点传递函数STF用Korff的理论模型代替，在Korff模型中STF只是r0的函数，由r0得出STF。最后由上述公式（3）退卷积STF后得出子块图像的傅里叶振幅。通过谱比获取r0并计算STF，就可以退卷积得到目标振幅，但一般短曝光图像含有多种随机噪声。为抑制高频噪声在做除法时被放大，退卷积STF时平均功率谱还需乘上一个维纳滤波器。
相位重建基于扩展Knox-Thompson算法的方法。1973年Knox和Thompson对Labeyrie法进行改进提出了一种基于计算频谱的互相关（又称为互谱）来获取目标傅里叶相位的方法，称为Knox-Thompson（K-T）法，又称互谱法或交叉谱算法。该方法在统计二阶矩的时候加上一个频率的平移。 函数 i(x) 的互谱表达式如下：
， （5）
其中，I(q)表示i(x)的傅里叶变换，q是二维空间频率变量；式中*表示共轭，∆q 是空间频率平移量，当∆q = 0时，互谱就等同于能谱。对于∆q大小，Knox和Thompson在取其为频谱图像大小的1/6。随后Goodman在其著作里对K-T法进行了解释，并认为|∆q|应该小于r0/λf，这里r0是大气相干长度，λ是观测波长，f是系统焦距。由上式，斑点图的互谱统计结果为：
（6）
其中是K-T法的传递函数，von der Luhe曾推导了其系综平均，并给出了实函数形式的简化公式。Roggemann在其著作里对K-T传递函数也进行了推导，并证明当|∆q| > r0/λf时K-T法系综平均传递函数近似为零，进一步说明了计算互谱时|∆q|大小不应超过r0/λf。对式（6）退卷积K-T传递函数可得目标的互谱，对式（5）两边取相位有：
（7）
式中Φc(q,∆q)是互谱相位，Φo(q)是目标相位。由上式可看出，目标频谱相位与互谱相位存在递推关系，互相关谱的相位值就代表着目标各相邻频率处的相位差。通过从零频处的迭代计算其他频率处的相位值，为了提高相位重建的精度，减少递推过程中的误差累积，在K-T法的递推过程中加了权重，当迭代得到的相位与原始各频率处的相位差足够小的时候认为迭代终止。最终可由互谱相位实现目标相位从低频到高频依次重建。
在步骤3中，在子进程中，若当前子进程重建完分配的子块后就判断是否还有子块没有处理，有的话则分别分配给各个已完成工作的子进程处理，处理流程和以上描述一致，在各个子进程将所有子块重建完成后，主进程对所有子块的重建结果进行拼接得到一帧整个视场的重建图像，保存相关数据后完成重建。
在K-T算法高分辨重建天文图像中太阳斑点图的并行化实现方法中，使用GPU并行加速的总体思想为：
1、在CUDA模型下，先由CPU将运算所需要的数据复制到GPU内存中；
2、再由GPU内部核函数进行处理；
3、最后将结果数据从GPU内存复制回CPU内存。
具体来说，运行在GPU上的程序称为Kernel（内核）函数，Kernel函数以Grid（线程网格）的形式组织，每个Grid由若干Block（线程块）组成，每个Block又包含若干Thread（线程）。在程序执行时，Grid被加载到SPA（流式处理器矩阵级别）上，Block则被分发到各个SM（流多处理器），Thread以Warp（线程束）为单位运行。Block可以是1维、2维或者3维。一个Block内的线程可彼此协作，通过共享存储器来共享数据，并同步其执行来协同存储器访问。通过Block的粗粒度分解和Thread的细粒度并行执行就可以完成任务的并行计算，达到加速的目的。
GPU并行运算的耗时主要分为两块，裸数据从内存传输至显卡的时间及并行运算的时间。在本方法中，尽量减少主机和设备的数据传输，减少两者的数据拷贝。在子进程拿到子块数据后，在求振幅和相位之前，一次性发送给GPU。
根据各个子任务的特点将其调度到不同的硬件上执行，子任务分为串行子任务和并行子任务，其中，CPU负责子任务的调度和串行子任务的执行，执行主机代码；GPU负责并行子任务的执行，执行设备代码例如__global__和__device__中的函数。CPU与GPU之间通过PCI-E总线进行通信。
尽管GPU在解决计算密集型问题上比CPU具有较大优势，然而不可能仅仅依靠GPU来完成计算任务，很多串行指令和控制指令任然需要CPU来完成。因此，基于CUDA编程模型的GPU通用计算实际上是充分发挥CPU与GPU在解决特定问题中各自的优势，构建CPU和GPU协作的CPU-GPU异构计算模式来更好地解决计算问题。CPU-GPU异构计算模式有两个特定：一是CPU与GPU的执行是异步的，即CPU将任务分配给GPU后，它将继续执行下面的任务而不会等待GPU端的运算完成，直到执行同步指令__syncthreads()。这种模式虽然需要人为地对CPU与GPU的运算进行同步，但却保证了运算的准确性。二是通信开销量较大。在执行任务的过程中，CPU将数据通过PCI-E总线发送给GPU，GPU在完成运算之后再将结果回传到内存中，这种模式将增加大数据量应用场合的执行时间。
本方法中，对全局存储器采用合并访问优化。在满足合并访问条件的情况下，只需要一次传输便可以处理来自一个Half-Warp的线程对全局内存的访问请求。在重建中，利用CUDA编程模型中自带的CUFFT库进行频域傅里叶变换操作。通过使用常数存储器存储部分计算参数的中间量，减少GPU的重复计算量。通过优化一个内核函数中多重循环迭代，减少访问全局存储器的次数。
在求相位迭代终止时，求累加和的并行优化策略如下：首先获得处理像素对应的线程映射；然后由多个并行线程在满足不超过总线程数的情况下；并行执行部分；不同线程号对应不同的vartmp数组下标，存放当前线程运算结果；同步以上并行过程的线程，等待所有线程结束计算，避免影响后续累加运算。
本方法在CUDA的线程数和线程块的策略上采用：根据GPU的特性设置每个线程块（Block）的尺寸，根据待重建图像的尺寸设置线程块的个数。
if(nx >= ny) maxnxny = nx; else maxnxny = ny;
frachamming_kernel1<<<(maxnxny+15)/16, 16>>>(nx, ny, Hx, Hy, nxfh, nyfh);
即线程块个数是(maxnxny+15)/16，每个块16个线程。
分配CUDA核函数的线程块数以及每个线程块所包含的线程数，将每个线程与每个像素点一一对应，进行所有像素点的并行计算。从像素点出发，每次处理该像素点数据的运算。在CUDA架构下提供内建变量threadIdx.x和blockIdx.x来进行寻址，使用二维的线程块直接索引像素点的坐标：
int x = threadIdx.x + blockIdx.x * blockDim.x;表示x方向像素点坐标，
int y = threadIdx.y + blockIdx.y * blockDim.y;表示y方向像素点坐标，
int offset = x + y * blockDim.x * gridDim.x;转换为图形中的唯一索引。
本方法在GPU的SM占用率和每个线程完成的计算量间的分配策略采用：在减少全局存储器读写次数的同时，算法会增加内核函数（Kernel）中每个线程的计算负担，若一味的增大每个线程块的线程数量，势必会降低整个GPU中活动Block和活动Warp的数量，这样又会反过来影响GPU隐藏长延时操作的优势。因此，按输出对GPU的任务进行划分。在线程分配时考虑一个很重要的原则，即流多处理器（SM）的占用率不能太低，SM占用率指的是每个SM上的活动Warp数量与允许的最大活动Warp数量之比；这主要是因为GPU是通过线程间的切换来隐藏长延时操作（访问全局存储器、线程间同步等）的，当一个Warp中的线程进行长延时操作时，另外一个活动Warp中的线程会自动进入运算状态，这样便可以隐藏一部分延迟；但这并不代表SM占用率越高越好，SM占用率越高则每个线程占用的GPU资源越少，即每个线程完成的计算量越少，而每个SM上的最大活动Warp数量是一定的，因此就可能出现即时SM中活动的Warp数量达到了最大，由于Warp中每个线程计算量太少，从而使所有活动Warp中线程同时进入长延时操作，不能充分的隐藏延迟；因此，在GPU的SM占用率和每个线程完成的计算量间选择一个平衡点，使GPU的性能发挥到最佳。
同时，本方法中利用CUDA流技术，使应用程序实现任务级的并行化，即GPU可以并行执行两个或多个不同的任务，使GPU在执行核函数的同时，能在主机与设备之间执行复制操作。两个或两个以上的流处理数据互不影响。
通过CPU控制多个CUDA流的运行顺序，同步运行当前CUDA流所进行的数据复制运算和相邻CUDA流所进行的核函数运算。具体来说，CUDA流表示一个GPU操作队列，并且该队列中的操作将以指定的顺序执行。在一个流中可以添加操作，例如核函数启动，内存复制，以及事件的启动和结束。GPU的硬件基础有两个基本的引擎，内存复制引擎和核函数执行引擎，这两个引擎在核函数访问的数据已经复制好的前提下是可以并行运作的，从而都有益于本方法中重建效率的提高。
完成本发明所述K-T算法高分辨重建太阳斑点图的三个步骤后，针对待重建的一组图像中的一帧（如图5），对比观察通过本发明中方法重建后的图像（如图6），可以明显看到清晰度有很大程度的提高，米粒细节可辨识度提升明显。
上面结合附图对本发明的具体实施方式作了详细说明，但是本发明并不限于上述实施方式，在本领域普通技术人员所具备的知识范围内，还可以在不脱离本发明宗旨的前提下作出各种变化。
Claims (6)
1.一种K-T算法高分辨重建天文图像中斑点图的并行化实现方法，其特征在于：首先通过CPU主进程读取图像及参数文件并且计算相关物理参数进行图像的预处理和数据对齐；接着把全部帧数的整个视场图像依据等晕区大小分割成很多的子块进行图像分块；然后这些子块被分发到独立的子进程中去采用CUDA并行编程模型进行并行计算处理，进行傅里叶的振幅和相位重建；最后子进程对每一个子块通过联合傅里叶振幅和相位进行傅里叶逆变换获得各个子块的重建结果，主进程对所有子块的重建结果进行拼接得到一帧整个视场的重建图像，保存相关数据后完成重建。
2.根据权利要求1所述的K-T算法高分辨重建天文图像中斑点图的并行化实现方法，其特征在于：对待处理的全部帧数的整个视场图像依据等晕区大小进行整体分块，为了减少边缘效应对最后重建结果的影响，采取对各个子块50%重叠的方式进行分割，即分块的时候子块之间有交叠，同时采用海明窗和汉宁窗对图像进行加窗处理。
3.根据权利要求1所述的K-T算法高分辨重建天文图像中斑点图的并行化实现方法，其特征在于：所述子块的振幅重建在频域中由斑点干涉术退卷积传递函数完成，相位重建由扩展Knox-Thompson法计算得出。
4.根据权利要求3所述的K-T算法高分辨重建天文图像中斑点图的并行化实现方法，其特征在于：所述子块的振幅重建具体过程为，在子进程中由GPU多线程并行计算，首先通过对系列短曝光图像进行谱比计算，然后将计算所得谱比和标准谱比进行比对得到大气相干长度r0；斑点传递函数STF用Korff的理论模型代替，在Korff模型中STF只是r0的函数，由r0得出STF；最后退卷积STF后得出傅里叶振幅。
5.根据权利要求3所述的K-T算法高分辨重建天文图像中斑点图的并行化实现方法，其特征在于：所述子块的相位重建具体过程为，将频谱中相隔较小频域距离的相位相减来消除由大气带来的相位起伏，留下的部分是这对频率的相位差；有了整个频域面上的相位差，通过由低频向高频的递推实现相位重建；在子进程中由GPU多线程并行计算K-T交叉谱、振幅加权整合、误差最小化迭代。
6.根据权利要求1所述的子进程中采用CUDA并行模型进行并行计算处理，其特征在于：所述CUDA并行编程模型进行并行计算处理过程为在CUDA模型下，先由CPU将运算所需要的数据复制到GPU内存中；再由GPU内部核函数进行处理；最后将结果数据从GPU内存复制回CPU内存。
Priority Applications (1)
| Application Number | Priority Date | Filing Date | Title |
|---|---|---|---|
| CN201710466806.3A CN107451955A (zh) | 2017-06-20 | 2017-06-20 | 一种k‑t算法重建天文图像中斑点图的并行化实现方法 |
Applications Claiming Priority (1)
| Application Number | Priority Date | Filing Date | Title |
|---|---|---|---|
| CN201710466806.3A CN107451955A (zh) | 2017-06-20 | 2017-06-20 | 一种k‑t算法重建天文图像中斑点图的并行化实现方法 |
Publications (1)
| Publication Number | Publication Date |
|---|---|
| CN107451955A true CN107451955A (zh) | 2017-12-08 |
Family
ID=60486567
Family Applications (1)
| Application Number | Title | Priority Date | Filing Date |
|---|---|---|---|
| CN201710466806.3A Pending CN107451955A (zh) | 2017-06-20 | 2017-06-20 | 一种k‑t算法重建天文图像中斑点图的并行化实现方法 |
Country Status (1)
| Country | Link |
|---|---|
| CN (1) | CN107451955A (zh) |
Cited By (13)
| Publication number | Priority date | Publication date | Assignee | Title |
|---|---|---|---|---|
| CN106650925A (zh) * | 2016-11-29 | 2017-05-10 | 郑州云海信息技术有限公司 | 一种基于MIC集群的深度学习框架Caffe系统及算法 |
| CN108510492A (zh) * | 2018-04-04 | 2018-09-07 | 中国科学院云南天文台 | 一种基于太阳色球幸运图像的多核并行实时重建方法 |
| CN108710753A (zh) * | 2018-05-18 | 2018-10-26 | 上海精密计量测试研究所 | 临近空间中子环境仿真方法 |
| CN109683018A (zh) * | 2018-12-24 | 2019-04-26 | 电子科技大学 | 一种实时多帧频域数据的并行处理方法 |
| CN109949256A (zh) * | 2019-01-14 | 2019-06-28 | 昆明理工大学 | 一种基于傅里叶变换的天文图像融合方法 |
| CN109961397A (zh) * | 2018-04-12 | 2019-07-02 | 华为技术有限公司 | 图像重建方法及设备 |
| CN111539997A (zh) * | 2020-04-23 | 2020-08-14 | 中国科学院自动化研究所 | 基于gpu计算平台的图像并行配准方法、系统、装置 |
| CN111754393A (zh) * | 2020-06-28 | 2020-10-09 | 展讯通信（上海）有限公司 | 图像处理方法、系统、电子设备和介质 |
| CN112446001A (zh) * | 2020-11-05 | 2021-03-05 | 西安中科星图空间数据技术有限公司 | 一种适用于分布式架构的波段计算方法、装置及存储介质 |
| CN112991141A (zh) * | 2021-02-23 | 2021-06-18 | 昆明理工大学 | 一种基于gpu并行加速的频域幸运成像方法 |
| CN113344765A (zh) * | 2021-05-14 | 2021-09-03 | 中国科学院国家空间科学中心 | 一种频域天文图像目标检测方法及系统 |
| CN113409183A (zh) * | 2021-08-02 | 2021-09-17 | 广州汇图计算机信息技术有限公司 | 一种基于gpu的快速重建成像方法 |
| CN117061759A (zh) * | 2023-10-11 | 2023-11-14 | 苏州元脑智能科技有限公司 | 图像压缩方法、装置、计算机设备及存储介质 |
Citations (2)
| Publication number | Priority date | Publication date | Assignee | Title |
|---|---|---|---|---|
| CN106570890A (zh) * | 2016-11-08 | 2017-04-19 | 昆明理工大学 | 一种提取太阳高分辨率序列图像中不同速度区间内动态信息的方法 |
| CN106791859A (zh) * | 2015-11-24 | 2017-05-31 | 龙芯中科技术有限公司 | 视频编码方法和视频编码器 |
-
2017
- 2017-06-20 CN CN201710466806.3A patent/CN107451955A/zh active Pending
Patent Citations (2)
| Publication number | Priority date | Publication date | Assignee | Title |
|---|---|---|---|---|
| CN106791859A (zh) * | 2015-11-24 | 2017-05-31 | 龙芯中科技术有限公司 | 视频编码方法和视频编码器 |
| CN106570890A (zh) * | 2016-11-08 | 2017-04-19 | 昆明理工大学 | 一种提取太阳高分辨率序列图像中不同速度区间内动态信息的方法 |
Non-Patent Citations (3)
| Title |
|---|
| 向永源 等: "高分辨率太阳图像重建方法", 《天文学进展》 * |
| 李雪宝: "太阳望远镜海量数据并行处理技术研究", 《中国博士学位论文全文数据库 基础科学辑》 * |
| 钟立波: "太阳高分辨力图像重建技术研究", 《中国博士学位论文全文数据库 信息科技辑》 * |
Cited By (21)
| Publication number | Priority date | Publication date | Assignee | Title |
|---|---|---|---|---|
| CN106650925A (zh) * | 2016-11-29 | 2017-05-10 | 郑州云海信息技术有限公司 | 一种基于MIC集群的深度学习框架Caffe系统及算法 |
| CN108510492A (zh) * | 2018-04-04 | 2018-09-07 | 中国科学院云南天文台 | 一种基于太阳色球幸运图像的多核并行实时重建方法 |
| CN108510492B (zh) * | 2018-04-04 | 2023-08-25 | 中国科学院云南天文台 | 一种基于太阳色球幸运图像的多核并行实时重建方法 |
| CN109961397A (zh) * | 2018-04-12 | 2019-07-02 | 华为技术有限公司 | 图像重建方法及设备 |
| CN109961397B (zh) * | 2018-04-12 | 2023-07-28 | 华为技术有限公司 | 图像重建方法及设备 |
| CN108710753A (zh) * | 2018-05-18 | 2018-10-26 | 上海精密计量测试研究所 | 临近空间中子环境仿真方法 |
| CN109683018B (zh) * | 2018-12-24 | 2020-12-01 | 电子科技大学 | 一种实时多帧频域数据的并行处理方法 |
| CN109683018A (zh) * | 2018-12-24 | 2019-04-26 | 电子科技大学 | 一种实时多帧频域数据的并行处理方法 |
| CN109949256B (zh) * | 2019-01-14 | 2023-04-07 | 昆明理工大学 | 一种基于傅里叶变换的天文图像融合方法 |
| CN109949256A (zh) * | 2019-01-14 | 2019-06-28 | 昆明理工大学 | 一种基于傅里叶变换的天文图像融合方法 |
| CN111539997B (zh) * | 2020-04-23 | 2022-06-10 | 中国科学院自动化研究所 | 基于gpu计算平台的图像并行配准方法、系统、装置 |
| CN111539997A (zh) * | 2020-04-23 | 2020-08-14 | 中国科学院自动化研究所 | 基于gpu计算平台的图像并行配准方法、系统、装置 |
| CN111754393A (zh) * | 2020-06-28 | 2020-10-09 | 展讯通信（上海）有限公司 | 图像处理方法、系统、电子设备和介质 |
| CN112446001A (zh) * | 2020-11-05 | 2021-03-05 | 西安中科星图空间数据技术有限公司 | 一种适用于分布式架构的波段计算方法、装置及存储介质 |
| CN112991141A (zh) * | 2021-02-23 | 2021-06-18 | 昆明理工大学 | 一种基于gpu并行加速的频域幸运成像方法 |
| CN112991141B (zh) * | 2021-02-23 | 2022-05-20 | 昆明理工大学 | 一种基于gpu并行加速的频域幸运成像方法 |
| CN113344765A (zh) * | 2021-05-14 | 2021-09-03 | 中国科学院国家空间科学中心 | 一种频域天文图像目标检测方法及系统 |
| CN113344765B (zh) * | 2021-05-14 | 2023-11-03 | 中国科学院国家空间科学中心 | 一种频域天文图像目标检测方法及系统 |
| CN113409183A (zh) * | 2021-08-02 | 2021-09-17 | 广州汇图计算机信息技术有限公司 | 一种基于gpu的快速重建成像方法 |
| CN117061759A (zh) * | 2023-10-11 | 2023-11-14 | 苏州元脑智能科技有限公司 | 图像压缩方法、装置、计算机设备及存储介质 |
| CN117061759B (zh) * | 2023-10-11 | 2024-02-06 | 苏州元脑智能科技有限公司 | 图像压缩方法、装置、计算机设备及存储介质 |
Similar Documents
| Publication | Publication Date | Title |
|---|---|---|
| CN107451955A (zh) | 一种k‑t算法重建天文图像中斑点图的并行化实现方法 | |
| DE102018114929B4 (de) | System und verfahren für ein rendering eines lichtfeldes | |
| US20200090383A1 (en) | Multi-level image reconstruction using one or more neural networks | |
| US8330763B2 (en) | Apparatus and method for volume rendering on multiple graphics processing units (GPUs) | |
| DE102018133555A1 (de) | Berechnungsoptimierungsmechanismus für tiefe neuronale Netze | |
| DE102018110371A1 (de) | Intelligente speicherhandhabung und datenmanagement für maschinenlernnetzwerke | |
| DE102021120605A1 (de) | Effiziente softmax-berechnung | |
| DE102020107080A1 (de) | Grafiksysteme und Verfahren zum Beschleunigen von Synchronisation mittels feinkörniger Abhängigkeitsprüfung und Planungsoptimierungen basierend auf verfügbarem gemeinsam genutztem Speicherplatz | |
| DE102020132272A1 (de) | Verfahren und vorrichtung zum codieren basierend auf schattierungsraten | |
| DE102019115130A1 (de) | Vorrichtung und Verfahren für konservatives morphologisches Anti-Aliasing mit Mehrfachabtastung | |
| CN109300083A (zh) | 一种分块处理Wallis匀色方法及装置 | |
| DE102020131852A1 (de) | Vorrichtung und verfahren zum ausführen eines stabilen sortiervorgangs mit kurzer latenz | |
| DE102018128699A1 (de) | Einstellen einer Winkelabtastrate während eines Renderings unter Verwendung von Blickinformationen | |
| DE102020129409A1 (de) | Dynamisches unterteilen von aktivierungen und kernels zum verbessern von speichereffizienz | |
| DE102020129432A1 (de) | System und Verfahren zum Anpassen eines ausführbaren Objekts an eine Verarbeitungseinheit | |
| DE102020108526A1 (de) | Adaptive pixelabtastreihenfolge für zeitlich dichtes rendern | |
| DE102021104310A1 (de) | Reservoir-basiertes räumlich-zeitliches resampling nach wichtigkeit unter verwendung einer globalen beleuchtungsdatenstruktur | |
| CN106484532B (zh) | 面向sph流体模拟的gpgpu并行计算方法 | |
| CN106971369A (zh) | 一种基于gpu的地形可视域分析的数据调度与分发方法 | |
| DE102021116364A1 (de) | Vorrichtung und verfahren für strahlverfolgungs-detaillierungsgradübergänge von hoher qualität | |
| Hansen et al. | Synthetic aperture beamformation using the GPU | |
| Heigl et al. | High-speed reconstruction for C-arm computed tomography | |
| DE102020130081A1 (de) | Erweiterte prozessorfunktionen für berechnungen | |
| Melo et al. | Real-time HD image distortion correction in heterogeneous parallel computing systems using efficient memory access patterns | |
| DE112018007596T5 (de) | Vorrichtung und verfahren für merkmalspunktverfolgung unter verwendung von inter-frame-voraussage |
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 | ||
| RJ01 | Rejection of invention patent application after publication | ||
| RJ01 | Rejection of invention patent application after publication |
Application publication date: 20171208 |