CN108765382B - 基于gpu的丰度估计并行计算方法 - Google Patents

基于gpu的丰度估计并行计算方法 Download PDF

Info

Publication number
CN108765382B
CN108765382B CN201810461103.6A CN201810461103A CN108765382B CN 108765382 B CN108765382 B CN 108765382B CN 201810461103 A CN201810461103 A CN 201810461103A CN 108765382 B CN108765382 B CN 108765382B
Authority
CN
China
Prior art keywords
matrix
gpu
memory
kernel function
cpu
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
CN201810461103.6A
Other languages
English (en)
Other versions
CN108765382A (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.)
Dalian Maritime University
Original Assignee
Dalian Maritime University
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 Dalian Maritime University filed Critical Dalian Maritime University
Priority to CN201810461103.6A priority Critical patent/CN108765382B/zh
Publication of CN108765382A publication Critical patent/CN108765382A/zh
Application granted granted Critical
Publication of CN108765382B publication Critical patent/CN108765382B/zh
Active legal-status Critical Current
Anticipated expiration legal-status Critical

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
    • 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
    • G06T2207/00—Indexing scheme for image analysis or image enhancement
    • G06T2207/10—Image acquisition modality
    • G06T2207/10032—Satellite or aerial image; Remote sensing
    • 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

Landscapes

  • Engineering & Computer Science (AREA)
  • Physics & Mathematics (AREA)
  • General Physics & Mathematics (AREA)
  • Theoretical Computer Science (AREA)
  • Quality & Reliability (AREA)
  • Computer Vision & Pattern Recognition (AREA)
  • Image Processing (AREA)

Abstract

本发明公开了一种基于GPU的丰度估计并行计算方法，包括：通过CPU载入原始高光谱像元数据H以及端元数据M，并将所述原始高光谱像元数据H以及端元数据M发送至GPU显存；CPU调用GPU核函数，基于ROVP算法计算各端元mi的丰度αi；所述GPU将计算结果返回至CPU端，并由CPU端输出。通过上述方案实现了基于CUDA库的ROVP‑C算法和基于CUBLAS库的ROVP‑L算法，与传统的串行算法比较分析可知，本发明提出的算法计算速度快，有效提高了丰度估计运行效率。

Description

基于GPU的丰度估计并行计算方法
技术领域
本发明涉及高光谱遥感技术领域，具体说涉及一种基于GPU的丰度估计并行计算方法。
背景技术
高光谱遥感技术是近年来迅速发展起来的一种全新遥感技术。高光谱遥感数据涵盖了自然界中所有地物，由于数据复杂性以及空间分辨率的限制使其每个端元包含了众多物质信息，其中包括大量混合像元的存在，从而增加了数据分析的难度。丰度估计是高光谱混合像元分解技术中最重要的主题之一，其目的是精确分析混合像素的比重。目前常用的高光谱遥感混合像元分解算法均采用串行处理思路，受限于高光谱遥感图像自身空间分辨率、光谱分辨率等因素造成的大数据、冗余多等问题，算法计算过程复杂度难以降低。采用传统的串行处理方法，执行过程中高达数百亿浮点运算严重影响计算机执行速度，时间消耗巨大，无法满足混合像元分解的实时处理需求。因此更适于进行数据密集型和计算密集型计算GPU在高光谱遥感领域应用前景广阔。
正交向量投影算法(Orthogonal Vector Projection，OVP)采用Gram-Schmidt正交化估计混合像元中端元的丰度，不涉及任何矩阵求逆过程，仅存在类似于最小二乘误差(Linear Square Estimation，LSE)和正交子空间投影(Orthogonal SubspaceProjection，OSP)算法的重新计算问题，当有一个新的端元被添加到端元矩阵M中时，OVP算法还需要重新计算新的端元mp+1,当p很大时，会使计算时间大大增加。近期有研究者提出了一种名为递归的正交向量投影(Recursive Orthogonal Vector Projection,ROVP)的新算法，该算法是OVP算法的延伸，可以将一些已经计算过的重要的结果用到下一次迭代中，则计算成本将明显降低，并且通过实验证明当估计一个或者所有端元的丰度时，ROVP算法都是最快的，优于OVP算法，并且该算法适合于并行计算，本发明正以此为契机，对ROVP算法的并行设计进行深入分析。
发明内容
鉴于现有技术的不足，本发明的目的是要提供一种基于GPU的丰度估计并行计算方法，采用ROVP算法对高光谱图像进行丰度估计，以提高运算效率。
本发明的技术方案如下：
一种基于GPU的丰度估计并行计算方法，其特征在于，步骤包括：
通过CPU载入原始高光谱像元数据H以及端元数据M，并将所述原始高光谱像元数据H以及端元数据M发送至GPU显存；
CPU调用GPU核函数，基于ROVP算法计算各端元mi的丰度αi；
所述GPU将计算结果返回至CPU端，并由CPU端输出。
根据本发明实施例的另一方面，还提供了一种基于GPU的丰度估计并行计算系统，其特征在于包括：图像数据载入单元，通过CPU载入原始高光谱像元数据H以及端元数据M，并将所述原始高光谱像元数据H以及端元数据M发送至GPU显存；计算单元，通过CPU调用GPU核函数，基于ROVP算法计算各端元mi的丰度αi；输出单元，用以将所述GPU计算结果返回至CPU端，并由CPU端输出。
根据本发明实施例的另一方面，还提供了一种存储介质，所述存储介质包括存储的程序，其中，所述程序执行上述任意一项所述的方法。
根据本发明实施例的另一方面，还提供了一种处理器，所述处理器用于运行程序，其中，所述程序运行时执行上述任意一项所述的方法。
本发明是在ROVP算法基础上提出的基于GPU的丰度估计并行计算方法，实现了基于CUDA的ROVP-C算法和基于CUBLAS库的ROVP-L算法，与传统的串行算法比较分析可知，本发明提出的算法运行速度快，有效提高了丰度估计的运行效率。
附图说明
为了更清楚的说明本发明的实施例或现有技术的技术方案，下面将对实施例或现有技术描述中所需要使用的附图做一简单地介绍，显而易见地，下面描述中的附图仅仅是本发明的一些实施例，对于本领域普通技术人员来讲，在不付出创造性劳动的前提下，还可以根据这些附图获得其他的附图。
图1为本发明丰度估计计算方法流程图；
图2为基于CUDA的ROVP-C算法流程图；
图3为基于CUBLAS的ROVP-L算法流程图；
图4为实验数据图像；
图5为CPU串行、ROVP-C及ROVP-L三种算法的时间对比情况；
图6为ROVP-C算法和ROVP-L算法的加速比情况；
图7为OVP-GPU、ROVP-C及ROVP-L三种并行算法的时间对比情况；
图8(a)为CPU串行算法在端元个数为5的情况下解混第5个端元所对应的丰度估计结果图；
图8(b)为CPU串行算法在端元个数为10的情况下解混第5个端元所对应的丰度估计结果图；
图8(c)为CPU串行算法在端元个数为15的情况下解混第5个端元所对应的丰度估计结果图；
图9(a)为ROVP-C算法在端元个数为5个的情况下解混第5个端元所对应的丰度估计结果图；
图9(b)为ROVP-C算法在端元个数为10个的情况下解混第5个端元所对应的丰度估计结果图；
图9(c)为ROVP-C算法在端元个数为15个的情况下解混第5个端元所对应的丰度估计结果图；
图10(a)表示ROVP-L算法在端元个数为5个的情况下解混第5个端元所对应的丰度估计结果图；
图10(b)表示ROVP-L算法在端元个数为10个的情况下解混第5个端元所对应的丰度估计结果图；
图10(c)表示ROVP-L算法在端元个数为15个的情况下解混第5个端元所对应的丰度估计结果图。
具体实施方式
为使本发明的实施例的目的、技术方案和优点更加清楚，下面结合本发明实施例中的附图，对本发明实施例中的技术方案进行清楚完整的描述：
如图1所示：一种基于GPU的丰度估计并行计算方法，其特征在于，步骤包括：
A、通过CPU载入原始高光谱像元数据H以及端元数据M，并将所述原始高光谱像元数据H以及端元数据M发送至GPU显存。
B、CPU调用GPU核函数，基于ROVP算法计算各端元mi的丰度αi。具体包括：
S101、由GPU显存中调取高光谱像元数据H，其中H＝{B1,B2,...,Bi,...,Bl}，
Figure GDA0003510801590000041
N为像元数目，l为波段数目，r为高光谱影像中各混合像元的光谱向量，以及端元数据M，其中M＝[m1,m2,...mp-1,mp]是大小为l×p的端元矩阵，mi＝[mi1,mi2,...mip-1,mip]，
并设置初始迭代条件：
令
Figure GDA0003510801590000042
其中
Figure GDA0003510801590000043
表示第一个端元的正交向量，
Figure GDA0003510801590000044
表示向量
Figure GDA0003510801590000045
的模，定义Kernel矩阵相乘核函数求P，其中
Figure GDA00035108015900000413
具体包括：
a.利用cudaMalloc函数在GPU中取得一块显卡内存给矩阵v1分配空间，再利用cudaMemcpy函数将矩阵v1复制到显卡内存中。
b.设置调用内核函数时的线程数目，将线程数目设定为波段数的平方，即l×l。
c.调用内核函数，每个线程都要执行的是计算矩阵v1和矩阵v1的行号和列号对应下标的值的乘积运算
Figure GDA0003510801590000046
m和n分别为矩阵的行号和列号，将结果存储到矩阵P中。
d.再利用cudaMemcpy函数将计算结果从显存中复制回内存。
e.最后释放设备存储器中的空间。
S102、根据矩阵乘法和矩阵减法核函数对mi进行递归运算，具体包括：
(1)在第j次递归中，将新添加的第j个端元mj置于端元集M的最后，即[m1,m2…mp,mj]，其中2≤j≤P，调用GPU矩阵相乘和矩阵减法核函数根据以下公式计算正交于向量空间[m1,m2…mp]的向量
Figure GDA0003510801590000047
Figure GDA0003510801590000048
(2)对于第i个端元，其中1≤i<j，调用GPU矩阵相乘和矩阵减法核函数根据以下公式计算正交于向量空间[m1,m2…mi-1]的向量
Figure GDA0003510801590000049
Figure GDA00035108015900000410
(3)调用GPU矩阵相乘和矩阵减法核函数根据以下公式计算正交于向量空间[m1,m2…,mj-1,mj]的向量
Figure GDA00035108015900000411
Figure GDA00035108015900000412
其中矩阵减法的具体步骤如下：
a.利用cudaMalloc函数在GPU中取得一块显卡内存给矩阵分配空间，再利用cudaMemcpy函数将矩阵复制到显卡内存中。
b.设置调用内核函数时的线程数目，将线程数目设定为波段数的平方，即l×l。
c.调用内核函数，每个线程都要执行的是计算矩阵
Figure GDA0003510801590000051
和矩阵B的行号和列号对应的下标的减法运算，将计算结果返回到矩阵C。
d.再利用cudaMemcpy函数将计算结果从显存中复制回内存，此时结果矩阵C中存放着波段数的平方个数据。
e.最后释放设备存储器中的空间。
(4)调用GPU矩阵加法核函数根据以下公式更新P：
Figure GDA0003510801590000052
其中矩阵加法具体步骤如下：
a.利用cudaMalloc函数在GPU中取得一块显卡内存给矩阵分配空间，再利用cudaMemcpy函数将矩阵复制到显卡内存中。
b.设置调用内核函数时的线程数目，将线程数目设定为波段数的平方，即l×l。
c.调用内核函数，每个线程都要执行的是计算矩阵
Figure GDA0003510801590000053
和矩阵
Figure GDA0003510801590000054
的行号和列号对应的下标的加法运算，将计算结果返回到矩阵P。
d.再利用cudaMemcpy函数将计算结果从显存中复制回内存，此时结果矩阵P中存放着波段数的平方个数据。
e.最后释放设备存储器中的空间。
S103、判断迭代次数j是否与P相等，如相等则执行步骤S104，否则执行步骤S102；
S104、提取满足停止迭代要求的
Figure GDA0003510801590000055
并计算端元mi的丰度αi：
Figure GDA0003510801590000056
其中
Figure GDA0003510801590000057
表示第i个端元正交于向量空间[mi1,mi2,...mij-1,mij]的向量，
Figure GDA0003510801590000058
表示向量
Figure GDA0003510801590000059
的转置，
Figure GDA00035108015900000510
表示像元r在
Figure GDA00035108015900000511
方向的投影。
C、所述GPU将计算结果返回至CPU端，并由CPU端输出。
如图2所示为本发明基于CUDA的ROVP-C算法流程图。
作为本发明的优选，本实施例还提供了基于CUBLAS库的ROVP-L算法的丰度估计并行计算方法，算法步骤与前述ROVP-C算法基本一样，区别在于该算法不需要建立矩阵相乘核函数，也不需调用矩阵相乘核函数，而是调用CUBLAS库函数来实现矩阵相乘。具体步骤包括：
A、通过CPU载入原始高光谱像元数据H以及端元数据M，并将所述原始高光谱像元数据H以及端元数据M发送至GPU显存。
B、CPU调用GPU核函数，基于ROVP算法计算各端元mi的丰度αi。具体包括：
S201、由GPU显存中调取高光谱像元数据H，其中H＝{B1,B2,...,Bi,...,Bl}，
Figure GDA0003510801590000061
N为像元数目，l为波段数目，r为高光谱影像中各混合像元的光谱向量，以及端元数据M，其中M＝[m1,m2,...mp-1,mp]是大小为l×p的端元矩阵，mi＝[mi1,mi2,...mip-1,mip]，
并设置初始迭代条件：
令
Figure GDA0003510801590000062
其中
Figure GDA0003510801590000063
表示第一个端元的正交向量，
Figure GDA0003510801590000064
表示向量
Figure GDA0003510801590000065
的模，定义Kernel矩阵相乘核函数求P，其中
Figure GDA00035108015900000613
S202、根据CUBLAS库函数cublasDgemm以及矩阵减法核函数对P进行递归运算，包括：
(1)在第j次递归中，将新添加的第j个端元mj置于端元集M的最后，即[m1,m2…mp,mj]，其中2≤j≤P，调用CUBLAS库函数cublasDgemm和矩阵减法核函数根据以下公式计算正交于向量空间[m1,m2…mp]的向量
Figure GDA0003510801590000066
Figure GDA0003510801590000067
(2)对于第i个端元，其中1≤i<j，调用CUBLAS库函数cublasDgemm和矩阵减法核函数根据以下公式计算正交于向量空间[m1,m2…mj-1]的向量
Figure GDA0003510801590000068
Figure GDA0003510801590000069
(3)调用CUBLAS库函数cublasDgemm和矩阵减法核函数根据以下公式计算正交于向量空间[m1,m2…,mj-1,mj]的向量
Figure GDA00035108015900000610
Figure GDA00035108015900000611
(4)调用GPU矩阵加法核函数根据以下公式更新P：
Figure GDA00035108015900000612
上述内容中基于CUBLAS库实现矩阵乘法的步骤包括：
a.使用cublasCreateHandle创建一个CUBLAS句柄。
b.使用cudaMalloc可以分配用于输入输出的设备内存。
c.使用cublasSetVector向分配好的设备内存填充输入数据。
d.调用cublasDgemm库来让GPU执行矩阵乘法操作。
e.使用cublasGetVector从设备内存中取出结果。
f.使用cudaFree和cublasDestroy来释放CUDA和CUBLAS资源。
S203、判断迭代次数j是否与P相等，如相等则执行步骤S204，否则执行步骤S202；
S204、提取满足停止迭代要求的
Figure GDA0003510801590000071
并计算端元mi的丰度αi：
Figure GDA0003510801590000072
其中
Figure GDA0003510801590000073
表示第i个端元正交于向量空间[mi1,mi2,...mij-1,mij]的向量，
Figure GDA0003510801590000074
表示向量
Figure GDA0003510801590000075
的转置，
Figure GDA0003510801590000076
表示像元r在
Figure GDA0003510801590000077
方向的投影。
C、所述GPU将计算结果返回至CPU端，并由CPU端输出。
图3所示为本发明基于CUBLAS的ROVP-L算法流程图。
下面通过具体实施例对本发明的技术方案及效果做进一步说明和验证：
本发明在模拟图像以及真实图像上的验证
实验平台及实验数据
实验平台搭建：实验机器为HP-PC Z240，硬件配置处理器为Intel(R)Core(TM)i7-6700 CPU@3.40GHz四核、内存64GB；显卡信息：NvidiaQuadro M2000(4GB/惠普)，显卡内存为4GB。软件平台中操作系统为Windows 7，开发环境Visual Studio 2013以及CUDA7.5。
实验数据采用的是拍摄于1997年美国内华达州的一个赤铜矿区的赤铜矿图像，大小为350×350，有189个波段。图4显示的是图像中第100个波段的图像。本文实验提取不定个数的端元进行实验操作，端元个数从3到80，记录串行ROVP，ROVP-C以及ROVP-L的结果以及并行OVP-GPU的结果。
实验结果比较
当选择端元个数依次为10、40、80时，测试串行的ROVP算法，ROVP-C以及ROVP-L时间的实验结果如表1所示，为了更明显的对比每一个算法的加速比情况，把实验结果以折线图的形式表示，由图5可清楚的看出各版本的执行时间。
表1各版本执行时间对比(单位：ms)
Figure GDA0003510801590000081
ROVP-C和ROVP-L算法的加速比情况如表2所示，其中加速比是通过用CPU串行算法所用的时间分别除以ROVP-C和ROVP-L算法所用时间得到的比值。
表2 ROVP-C与ROVP-L的加速比情况
Figure GDA0003510801590000082
表2中数据显示，ROVP-C算法可加速3.2～10.1倍不等；ROVP-L算法加速1.8～13.8倍不等。为了更明显的对比每一个算法的加速比情况，把实验结果以折线图的形式表示，图6为ROVP-C算法和ROVP-L算法的加速比情况。
当选择端元个数依次为10、40、80时，并行的OVP-GPU算法，ROVP-C以及ROVP-L三种并行算法时间比较如表3所示：
表3 OVP-GPU，ROVP-C以及ROVP-L三种并行算法时间比较
Figure GDA0003510801590000083
Figure GDA0003510801590000091
从表3可以看出，OVP-GPU和ROVP-C算法执行时间相当，随着端元个数的增加，算法所用时间也随着增加，其中ROVP-C算法时间优势稍有体现，而ROVP-L算法时间很稳定，随着端元个数增加时间几乎不变，为了更明显的对比每一个算法的时间走势，把实验结果以折线图的形式表示，图7表示三种并行算法的时间随着端元个数增加的变化，从图中可以清楚的看出当端元个数达到50个的时候，ROVP-L算法所用时间与ROVP-C所用时间几乎相等，但当端元个数大于50时，ROVP-L算法所用时间明显少于其他两种算法。
图8、图9、图10中的(a)、(b)、(c)三张图分别表示CPU串行算法、ROVP-C算法和ROVP-L算法在端元个数为5、15、30个的情况下解混第5个端元所对应的丰度情况的灰度图像，从图中可以看出，随着端元个数的增加解混的效果越好，并且ROVP-C和ROVP-L算法的解混效果都与CPU一致，从而证明了结果的准确性。
本发明实施例还提供了一种基于GPU的丰度估计并行计算系统，其特征在于包括：图像数据载入单元，通过CPU载入原始高光谱像元数据H以及端元数据M，并将所述原始高光谱像元数据H以及端元数据M发送至GPU显存；计算单元，通过CPU调用GPU核函数，基于ROVP算法计算各端元mi的丰度αi；输出单元，用以将所述GPU计算结果返回至CPU端，并由CPU端输出。需要说明的是，本发明实施例的基于GPU的丰度估计并行计算系统可以用于执行本发明实施例所提供的基于GPU的丰度估计并行计算方法；本发明实施例的基于GPU的丰度估计并行计算方法也可以通过本发明实施例所提供的基于GPU的丰度估计并行计算系统来执行。
在本发明的上述实施例中，对各个实施例的描述都各有侧重，某个实施例中没有详述的部分，可以参见其他实施例的相关描述。
在本申请所提供的几个实施例中，应该理解到，所揭露的技术内容，可通过其它的方式实现。其中，以上所描述的装置实施例仅仅是示意性的，例如所述单元的划分，可以为一种逻辑功能划分，实际实现时可以有另外的划分方式，例如多个单元或组件可以结合或者可以集成到另一个系统，或一些特征可以忽略，或不执行。另一点，所显示或讨论的相互之间的耦合或直接耦合或通信连接可以是通过一些接口，单元或模块的间接耦合或通信连接，可以是电性或其它的形式。
所述作为分离部件说明的单元可以是或者也可以不是物理上分开的，作为单元显示的部件可以是或者也可以不是物理单元，即可以位于一个地方，或者也可以分布到多个单元上。可以根据实际的需要选择其中的部分或者全部单元来实现本实施例方案的目的。
另外，在本发明各个实施例中的各功能单元可以集成在一个处理单元中，也可以是各个单元单独物理存在，也可以两个或两个以上单元集成在一个单元中。上述集成的单元既可以采用硬件的形式实现，也可以采用软件功能单元的形式实现。
所述集成的单元如果以软件功能单元的形式实现并作为独立的产品销售或使用时，可以存储在一个计算机可读取存储介质中。基于这样的理解，本发明的技术方案本质上或者说对现有技术做出贡献的部分或者该技术方案的全部或部分可以以软件产品的形式体现出来，该计算机软件产品存储在一个存储介质中，包括若干指令用以使得一台计算机设备(可为个人计算机、服务器或者网络设备等)执行本发明各个实施例所述方法的全部或部分步骤。而前述的存储介质包括：U盘、只读存储器(ROM，Read-Only Memory)、随机存取存储器(RAM，Random Access Memory)、移动硬盘、磁碟或者光盘等各种可以存储程序代码的介质。以上所述，仅为本发明较佳的具体实施方式，但本发明的保护范围并不局限于此，任何熟悉本技术领域的技术人员在本发明揭露的技术范围内，根据本发明的技术方案及其发明构思加以等同替换或改变，都应涵盖在本发明的保护范围之内。

Claims (4)

1.一种基于GPU的丰度估计并行计算方法，其特征在于，步骤包括：
通过CPU载入原始高光谱像元数据H以及端元数据M，并将所述原始高光谱像元数据H以及端元数据M发送至GPU显存；
CPU调用GPU核函数，采用基于CUDA的ROVP-C算法计算各端元
Figure DEST_PATH_IMAGE001
的丰度
Figure 632933DEST_PATH_IMAGE002
，具体包括：
S101、由GPU显存中调取高光谱像元数据H，其中
Figure DEST_PATH_IMAGE003
，
Figure 953187DEST_PATH_IMAGE004
，
Figure DEST_PATH_IMAGE005
为像元数目，l为波段数目，
Figure 210336DEST_PATH_IMAGE006
为高光谱影像中各混合像元的光谱向量，以及端元数据M，其中
Figure DEST_PATH_IMAGE007
是大小为
Figure 198014DEST_PATH_IMAGE008
的端元矩阵，
Figure DEST_PATH_IMAGE009
，
并设置初始迭代条件：
令
Figure 736443DEST_PATH_IMAGE010
，
Figure DEST_PATH_IMAGE011
，其中
Figure 167162DEST_PATH_IMAGE012
表示第一个端元的正交向量，
Figure DEST_PATH_IMAGE013
表示向量
Figure 449239DEST_PATH_IMAGE014
的模，定义Kernel矩阵相乘核函数求
Figure DEST_PATH_IMAGE015
，其中
Figure 884900DEST_PATH_IMAGE016
，具体包括：
a.利用cudaMalloc函数在GPU中取得一块显卡内存给矩阵
Figure DEST_PATH_IMAGE017
分配空间，再利用cudaMemcpy函数将矩阵
Figure 236640DEST_PATH_IMAGE017
复制到显卡内存中，
b.设置调用内核函数时的线程数目，将线程数目设定为波段数的平方，即
Figure 718437DEST_PATH_IMAGE018
，
c.调用内核函数，每个线程都要执行的是计算矩阵
Figure 273046DEST_PATH_IMAGE017
和矩阵
Figure 156688DEST_PATH_IMAGE017
的行号和列号对应下标的值的乘积运算
Figure DEST_PATH_IMAGE019
，
Figure 177865DEST_PATH_IMAGE020
和
Figure DEST_PATH_IMAGE021
分别为矩阵的行号和列号，将结果存储到矩阵
Figure 52018DEST_PATH_IMAGE015
中，
d.再利用cudaMemcpy函数将计算结果从显存中复制回内存，
e.最后释放设备存储器中的空间，
S102、根据矩阵乘法和矩阵减法核函数对
Figure 269373DEST_PATH_IMAGE022
进行递归运算，具体包括：
（1）在第j次递归中，将新添加的第j个端元
Figure DEST_PATH_IMAGE023
置于端元集
Figure 679625DEST_PATH_IMAGE024
的最后，即
Figure DEST_PATH_IMAGE025
，其中
Figure 871704DEST_PATH_IMAGE026
，调用GPU矩阵相乘和矩阵减法核函数根据以下公式计算正交于向量空间
Figure 328093DEST_PATH_IMAGE025
的向量
Figure DEST_PATH_IMAGE027
：
Figure 257128DEST_PATH_IMAGE028
（1）
（2）对于第
Figure DEST_PATH_IMAGE029
个端元，其中
Figure 256308DEST_PATH_IMAGE030
，调用GPU矩阵相乘和矩阵减法核函数根据以下公式计算正交于向量空间
Figure DEST_PATH_IMAGE031
的向量
Figure 619287DEST_PATH_IMAGE032
：
Figure DEST_PATH_IMAGE033
（2）
（3）调用GPU矩阵相乘和矩阵减法核函数根据以下公式计算正交于向量空间
Figure 468032DEST_PATH_IMAGE034
的向量
Figure DEST_PATH_IMAGE035
：
Figure 699293DEST_PATH_IMAGE036
（3）
其中矩阵减法的具体步骤如下：
a.利用cudaMalloc函数在GPU中取得一块显卡内存给矩阵分配空间，再利用cudaMemcpy函数将矩阵复制到显卡内存中，
b.设置调用内核函数时的线程数目，将线程数目设定为波段数的平方，即
Figure 412034DEST_PATH_IMAGE018
，
c.调用内核函数，每个线程都要执行的是计算矩阵
Figure DEST_PATH_IMAGE037
和矩阵
Figure 680336DEST_PATH_IMAGE038
的行号和列号对应的下标的减法运算，将计算结果返回到矩阵
Figure DEST_PATH_IMAGE039
，
d.再利用cudaMemcpy函数将计算结果从显存中复制回内存，此时结果矩阵
Figure 19307DEST_PATH_IMAGE039
中存放着波段数的平方个数据，
e.最后释放设备存储器中的空间，
（4）调用GPU矩阵加法核函数根据以下公式更新
Figure 116576DEST_PATH_IMAGE015
：
Figure 90348DEST_PATH_IMAGE040
（4）
其中矩阵加法具体步骤如下：
a.利用cudaMalloc函数在GPU中取得一块显卡内存给矩阵分配空间，再利用cudaMemcpy函数将矩阵复制到显卡内存中，
b.设置调用内核函数时的线程数目，将线程数目设定为波段数的平方，即
Figure 982081DEST_PATH_IMAGE018
，
c.调用内核函数，每个线程都要执行的是计算矩阵
Figure DEST_PATH_IMAGE041
和矩阵
Figure 447828DEST_PATH_IMAGE042
的行号和列号对应的下标的加法运算，将计算结果返回到矩阵
Figure 348788DEST_PATH_IMAGE015
，
d.再利用cudaMemcpy函数将计算结果从显存中复制回内存，此时结果矩阵
Figure DEST_PATH_IMAGE043
中存放着波段数的平方个数据，
e.最后释放设备存储器中的空间，
S103、判断迭代次数j是否与
Figure 941181DEST_PATH_IMAGE015
相等，如相等则执行步骤S104，否则执行步骤S102，
S104、提取满足停止迭代要求的
Figure 144760DEST_PATH_IMAGE044
，并计算端元
Figure 284755DEST_PATH_IMAGE001
的丰度
Figure DEST_PATH_IMAGE045
，
Figure 536876DEST_PATH_IMAGE046
其中
Figure DEST_PATH_IMAGE047
表示第
Figure 974986DEST_PATH_IMAGE048
个端元正交于向量空间
Figure DEST_PATH_IMAGE049
的向量，
Figure 146205DEST_PATH_IMAGE050
表示向量
Figure 773495DEST_PATH_IMAGE047
的转置，
Figure DEST_PATH_IMAGE051
表示像元
Figure 563728DEST_PATH_IMAGE052
在
Figure DEST_PATH_IMAGE053
方向的投影；
所述GPU将计算结果返回至CPU端，并由CPU端输出。
2.一种基于GPU的丰度估计并行计算系统，其特征在于，用于执行权利要求1所述的方法，包括：
图像数据载入单元，通过CPU载入原始高光谱像元数据H以及端元数据M，并将所述原始高光谱像元数据H以及端元数据M发送至GPU显存；
计算单元，通过CPU调用GPU核函数，基于ROVP算法计算各端元
Figure 130713DEST_PATH_IMAGE054
的丰度
Figure DEST_PATH_IMAGE055
；
输出单元，用以将所述GPU计算结果返回至CPU端，并由CPU端输出。
3.一种存储介质，其特征在于，所述存储介质包括存储的程序，其中，所述程序执行权利要求1中所述的方法。
4.一种处理器，其特征在于，所述处理器用于运行程序，其中，所述程序运行时执行权利要求1中所述的方法。
CN201810461103.6A 2018-05-15 2018-05-15 基于gpu的丰度估计并行计算方法 Active CN108765382B (zh)

Priority Applications (1)

Application Number Priority Date Filing Date Title
CN201810461103.6A CN108765382B (zh) 2018-05-15 2018-05-15 基于gpu的丰度估计并行计算方法

Applications Claiming Priority (1)

Application Number Priority Date Filing Date Title
CN201810461103.6A CN108765382B (zh) 2018-05-15 2018-05-15 基于gpu的丰度估计并行计算方法

Publications (2)

Publication Number Publication Date
CN108765382A CN108765382A (zh) 2018-11-06
CN108765382B true CN108765382B (zh) 2022-06-24

Family

ID=64006774

Family Applications (1)

Application Number Title Priority Date Filing Date
CN201810461103.6A Active CN108765382B (zh) 2018-05-15 2018-05-15 基于gpu的丰度估计并行计算方法

Country Status (1)

Country Link
CN (1) CN108765382B (zh)

Citations (2)

* Cited by examiner, † Cited by third party
Publication number Priority date Publication date Assignee Title
CN106788714A (zh) * 2016-12-05 2017-05-31 重庆工商大学 一种基于光学计算的稀疏解混方法
CN107644393A (zh) * 2017-09-28 2018-01-30 大连海事大学 一种基于gpu的丰度估计算法的并行实现方法

Family Cites Families (1)

* Cited by examiner, † Cited by third party
Publication number Priority date Publication date Assignee Title
US9449244B2 (en) * 2013-12-11 2016-09-20 Her Majesty The Queen In Right Of Canada, As Represented By The Minister Of National Defense Methods for in-scene atmospheric compensation by endmember matching

Patent Citations (2)

* Cited by examiner, † Cited by third party
Publication number Priority date Publication date Assignee Title
CN106788714A (zh) * 2016-12-05 2017-05-31 重庆工商大学 一种基于光学计算的稀疏解混方法
CN107644393A (zh) * 2017-09-28 2018-01-30 大连海事大学 一种基于gpu的丰度估计算法的并行实现方法

Non-Patent Citations (2)

* Cited by examiner, † Cited by third party
Title
Recursive orthogonal vector projection algorithm for linear spectral unmixing;Meiping Song et al.;《2014 6th Workshop on Hyperspectral Image and Signal Processing Evolution in Remote Sensing (WHISPERS)》;20171026;1-4 *
用于光谱解混的正交向量投影算法;宋梅萍等;《光谱学与光谱分析》;20151231;3465-3470 *

Also Published As

Publication number Publication date
CN108765382A (zh) 2018-11-06

Similar Documents

Publication Publication Date Title
US11847550B2 (en) Sparse convolutional neural network accelerator
Gu et al. Self-guided network for fast image denoising
US20220284638A1 (en) Method for image processing, computer device, and storage medium
US9508181B2 (en) Ordering and rendering buffers for complex scenes with cyclic dependency
US20200257902A1 (en) Extraction of spatial-temporal feature representation
CN111382867A (zh) 神经网络压缩的方法、数据处理的方法及相关装置
CN108765282B (zh) 基于fpga的实时超分辨方法及系统
Mazumdar et al. A hardware-friendly bilateral solver for real-time virtual reality video
JP2022173321A (ja) オブジェクトの検出方法、装置、デバイス、媒体及びプログラム
Mahmoudi et al. Towards a smart selection of resources in the cloud for low‐energy multimedia processing
CN107644393B (zh) 一种基于gpu的丰度估计算法的并行实现方法
CN108520532B (zh) 识别视频中物体运动方向的方法及装置
CN114049491A (zh) 指纹分割模型训练、指纹分割方法、装置、设备及介质
CN108765382B (zh) 基于gpu的丰度估计并行计算方法
CN112614108A (zh) 基于深度学习检测甲状腺超声图像中结节的方法和装置
CN115146226B (zh) 基于张量压缩方法的流数据处理方法、装置及设备
CN107622498B (zh) 基于场景分割的图像穿越处理方法、装置及计算设备
CN116309158A (zh) 网络模型的训练方法、三维重建方法、装置、设备和介质
CN111860557A (zh) 图像处理方法及装置、电子设备及计算机存储介质
CN108765259B (zh) 一种基于gpu的高光谱图像ratgp和rosp并行优化方法
US11636569B1 (en) Matrix transpose hardware acceleration
Kisačanin et al. Algorithmic and software techniques for embedded vision on programmable processors
CN111027670B (zh) 特征图处理方法、装置、电子设备及存储介质
CN113762173A (zh) 人脸光流估计及光流值预测模型的训练方法和装置
CN114004751A (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