CN109410136B - 基于最短传递路径的匀色方法及处理装置 - Google Patents
基于最短传递路径的匀色方法及处理装置 Download PDFInfo
- Publication number
- CN109410136B CN109410136B CN201811178782.2A CN201811178782A CN109410136B CN 109410136 B CN109410136 B CN 109410136B CN 201811178782 A CN201811178782 A CN 201811178782A CN 109410136 B CN109410136 B CN 109410136B
- Authority
- CN
- China
- Prior art keywords
- image
- standard deviation
- color
- processed
- gray
- Prior art date
- Legal status (The legal status is an assumption and is not a legal conclusion. Google has not performed a legal analysis and makes no representation as to the accuracy of the status listed.)
- Active
Links
- 238000000034 method Methods 0.000 title claims abstract description 70
- 238000012545 processing Methods 0.000 title claims abstract description 35
- 230000005540 biological transmission Effects 0.000 title claims abstract description 23
- 238000010586 diagram Methods 0.000 claims abstract description 15
- 230000015654 memory Effects 0.000 claims description 35
- 230000009466 transformation Effects 0.000 claims description 20
- 238000004364 calculation method Methods 0.000 claims description 19
- 230000009467 reduction Effects 0.000 claims description 13
- 230000000903 blocking effect Effects 0.000 claims description 7
- 238000004590 computer program Methods 0.000 claims description 3
- 238000012546 transfer Methods 0.000 claims description 3
- 238000004422 calculation algorithm Methods 0.000 description 16
- 230000001133 acceleration Effects 0.000 description 12
- 238000002474 experimental method Methods 0.000 description 9
- 238000004458 analytical method Methods 0.000 description 6
- 230000008569 process Effects 0.000 description 6
- 230000000694 effects Effects 0.000 description 4
- 230000006870 function Effects 0.000 description 4
- 238000000265 homogenisation Methods 0.000 description 4
- 238000005457 optimization Methods 0.000 description 4
- 230000000007 visual effect Effects 0.000 description 4
- 238000012937 correction Methods 0.000 description 3
- 238000011426 transformation method Methods 0.000 description 3
- 238000012442 analytical experiment Methods 0.000 description 2
- 238000011161 development Methods 0.000 description 2
- 238000011156 evaluation Methods 0.000 description 2
- 238000013178 mathematical model Methods 0.000 description 2
- 230000008520 organization Effects 0.000 description 2
- 238000013441 quality evaluation Methods 0.000 description 2
- 238000011158 quantitative evaluation Methods 0.000 description 2
- 239000000654 additive Substances 0.000 description 1
- 230000000996 additive effect Effects 0.000 description 1
- 230000015556 catabolic process Effects 0.000 description 1
- 230000008859 change Effects 0.000 description 1
- 239000003086 colorant Substances 0.000 description 1
- 230000000052 comparative effect Effects 0.000 description 1
- 230000008878 coupling Effects 0.000 description 1
- 238000010168 coupling process Methods 0.000 description 1
- 238000005859 coupling reaction Methods 0.000 description 1
- 238000006731 degradation reaction Methods 0.000 description 1
- 238000013461 design Methods 0.000 description 1
- 238000001514 detection method Methods 0.000 description 1
- 239000006185 dispersion Substances 0.000 description 1
- 238000005286 illumination Methods 0.000 description 1
- 238000010191 image analysis Methods 0.000 description 1
- 238000013507 mapping Methods 0.000 description 1
- 238000012847 principal component analysis method Methods 0.000 description 1
- 238000003672 processing method Methods 0.000 description 1
- 230000005855 radiation Effects 0.000 description 1
- 238000010187 selection method Methods 0.000 description 1
Images
Classifications
-
- G—PHYSICS
- G06—COMPUTING; CALCULATING OR COUNTING
- G06T—IMAGE DATA PROCESSING OR GENERATION, IN GENERAL
- G06T5/00—Image enhancement or restoration
- G06T5/90—Dynamic range modification of images or parts thereof
-
- 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/28—Indexing scheme for image data processing or generation, in general involving image processing hardware
-
- 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
Landscapes
- Physics & Mathematics (AREA)
- General Physics & Mathematics (AREA)
- Engineering & Computer Science (AREA)
- Theoretical Computer Science (AREA)
- Image Processing (AREA)
Abstract
本发明涉及一种基于最短传递路径的匀色方法及处理装置。通过获取设定区域内的各影像，并计算各影像的清晰度，以清晰度最大的影像作为初始参考影像；根据设定区域内的各影像的中心点构建Voronoi图，计算初始参考影像的中心点到其余各影像的中心点的最短路径；根据初始参考影像和待处理影像对应最短路径上的前一影像的匀色结果对该待处理影像进行Wallis匀色，实现了一定区域范围内多幅影像匀色，解决了区域范围内的多幅影像匀色时仅采用单一的参考影像导致匀色结果出现偏色或退化的问题。
Description
技术领域
本发明涉及图像处理技术领域，特别是基于最短传递路径的匀色方法及处理装置。
背景技术
在遥感领域中，大范围的镶嵌影像是影像分析和应用中的重要数据源。但是，受获取时刻光照、天气以及时相等条件的影响，待镶嵌影像之间存在着不同程度的色彩差异和对比度差异，使影像镶嵌更为复杂，也更为困难。因此，对区域范围内的多幅影像进行匀色，也称色彩一致性处理，消除影像之间色彩差异，使测区内影像在色调上保持一致，降低影像镶嵌的难度，具有重要的现实意义。
现有匀色方法可以分为两类，一类是非线性方法，如直方图匹配法和Gamma校正法；另一类是线性变换法，也是最受关注的方法，这类方法倾向于从影像的重叠区域选取不变像素作为样本参与统计，然后利用线性模型进行相对辐射校正。线性变换法中像素样本一般采用人工选取的方法，迭代加权的多元变化检测法，迭代慢特征分析法或者加权主成分分析法等进行筛选。上述的匀色方法虽然可以不同程度的解决色彩不一致问题，但是通常不能很好地消除对比度不一致现象。因此，一种基于均值和方差的特殊线性变换即Wallis变换，引入到色彩一致性处理中。Wallis变换是基于参考影像的处理方法，然而当一个区域范围内存在多幅影像，各幅影像之间的内容往往变化较大，单一的参考影像不适合区域区域内多幅图像之间的匀色处理，如果待处理影像与参考影像的影像内容差距较大，处理结果容易偏色或退化。
发明内容
本发明的目的是提供基于最短传递路径的匀色方法及处理装置，用以解决区域范围内的多幅影像匀色时仅采用单一的参考影像导致匀色结果出现偏色或退化的问题。
为了实现一定区域范围内多幅影像匀色，解决区域范围内的多幅影像匀色时仅采用单一的参考影像导致匀色结果出现偏色或退化的问题。本发明提供一种基于最短传递路径的匀色方法，包括以下步骤：
1)获取设定区域内的各影像，并计算各影像的清晰度，以清晰度最大的影像作为初始参考影像；
2)根据所述设定区域内的各影像的中心点构建Voronoi图，计算初始参考影像的中心点到其余各影像的中心点的最短路径；
3)根据初始参考影像和待处理影像对应最短路径上的前一影像的匀色结果对该待处理影像进行Wallis匀色。
进一步地，为了避免多次传递导致的色彩偏移，保持区域内影像整体的色彩一致性，步骤3)中当待处理影像为i，待处理影像i对应最短路径上的前一影像为j时，若影像j的匀色结果为j′，则在对待处理影像i处理时，对应的参考影像的灰度均值mf和标准偏差sf的计算公式为：
mf＝wmj′+(1-w)m1
sf＝wsj′+(1-w)s1
式中，mj′和sj′分别为影像j′的灰度均值和标准偏差，m1和s1分别为初始参考影像的灰度均值和标准偏差，w∈[0,1]为权值常数。
进一步地，为了实现影像的匀色以便于大范围的影像镶嵌，解决现有的Wallis变换在对大范围影像进行匀色时影像之间存在色彩和对比度都不一致的问题，所述Wallis匀色为分块处理Wallis匀色，该分块处理Wallis匀色包括以下步骤：
(1)对所述待处理影像进行互不重叠分块并统计各影像块的灰度均值和标准偏差；
(2)根据各影像块的灰度均值和标准偏差计算各影像块的四个角点对应的灰度均值和标准偏差；
(3)根据每个像素到所在影像块边缘的距离和该像素所在影像块的四个角点的灰度均值和标准偏差，计算每个像素的灰度均值和标准偏差；
(4)根据每个像素的灰度均值和标准偏差，并结合参考影像的灰度均值和标准偏差，对每个像素进行Wallis变换处理。
进一步地，为了角点赋值更加精确，步骤(2)中角点若只属于一个影像块，则将该影像块的灰度均值和标准偏差赋给该角点；角点若为多个相邻影像块之间的公共角点，则将所属多个影像块的灰度均值和标准偏差的平均值赋给该角点。
进一步地，为了获得较好地匀色效果，步骤(1)中对待处理的影像进行互不重叠分块，分块的个数为W×H，其中
W＝r×w,H＝r×h
式中，CV为影像的变异系数，CVRef为参考影像的变异系数，w和h分别为预设的行、列方向的参考分块数。
进一步地，为了实现不同影像间的色彩均衡，所述Wallis变换的公式为：
f(x,y)＝[g(x,y)-m(x,y)]·(sf/s(x,y))+mf
式中，g(x,y)为待处理影像的灰度值，f(x,y)为处理后的结果影像的灰度值，mf为参考影像的灰度均值，sf为参考影像的标准偏差，m(x,y)和s(x,y)分别为像素的灰度均值和标准偏差。
进一步地，为了提高像素灰度均值和标准偏差的计算速度，步骤(3)中每一个像素的灰度均值和标准偏差的计算为并行计算。
进一步地，为了实现对影像块均值和标准偏差的并行计算，步骤(1)中各影像块的灰度均值和标准偏差为通过并行归约求和法求得。
为了便于实现上述匀色方法，本发明还提供一种基于最短传递路径的匀色处理装置，包括存储器、处理器以及存储在存储器中并可在处理器上运行的计算机程序，所述处理器执行所述程序时实现以下步骤：
1)获取设定区域内的各影像，并计算各影像的清晰度，以清晰度最大的影像作为初始参考影像；
2)根据所述设定区域内的各影像的中心点构建Voronoi图，计算初始参考影像的中心点到其余各影像的中心点的最短路径；
3)根据初始参考影像和待处理影像对应最短路径上的前一影像的匀色结果对该待处理影像进行Wallis匀色。
进一步地，为了实现影像的匀色以便于大范围的影像镶嵌，所述Wallis匀色为分块处理Wallis匀色，该分块处理Wallis匀色包括以下步骤：
(1)对所述待处理影像进行互不重叠分块并统计各影像块的灰度均值和标准偏差；
(2)根据各影像块的灰度均值和标准偏差计算各影像块的四个角点对应的灰度均值和标准偏差；
(3)根据每个像素到所在影像块边缘的距离和该像素所在影像块的四个角点的灰度均值和标准偏差，计算每个像素的灰度均值和标准偏差；
(4)根据每个像素的灰度均值和标准偏差，并结合参考影像的灰度均值和标准偏差，对每个像素进行Wallis变换处理。
进一步地，为了满负荷利用计算资源，可以大幅度缩短计算时间，所述处理器包括GPU端，所述分块处理Wallis匀色的步骤在GPU端运行，所述GPU端包括至少两个线程块，一个影像块对应一个线程块，影像块中的像素与对应该影像块的线程块中的线程一一对应；步骤(3)中一个像素的灰度均值和标准偏差的计算对应一个线程，所有线程并行运算。
进一步地，为了达到并行加速的目的，避免存储片冲突并保持线程块的相邻线程处于活跃状态，所述GPU端在对各影像块的灰度均值和标准偏差求解时通过并行归约求和法进行运算。
进一步地，为了提高访问存储器的速度，所述存储器包括常量存储器，所述常量存储器用于存储四个角点的均值和标准偏差以及参考影像的灰度均值和标准偏差。
附图说明
图1是本发明的Voronoi示意图；
图2是本发明的一种分块处理Wallis匀色方法中任意像素的计算原理图；
图3是GPU线程组织方式的原理图；
图4是本发明的归约求和法的原理图；
图5是本发明的并行归约求和的原理图；
图6是本发明的基于最短传递路径的匀色方法的流程图；
图7(a)是分析实验中第一张原始影像图；
图7(b)是第一张原始影像图经过全局Wallis法处理后的影像图；
图7(c)是第一张原始影像图经过Inpho法处理后的影像图；
图7(d)是第一张原始影像图经过本发明的匀色方法处理后的影像图；
图8(a)是分析实验中第二张原始影像图；
图8(b)是第二张原始影像图经过全局Wallis法处理后的影像图；
图8(c)是第二张原始影像图经过Inpho法处理后的影像图；
图8(d)是第二张原始影像图经过本发明的匀色方法处理后的影像图；
图9(a)是分析实验中第二张原始影像中第一区域A的放大图；
图9(b)是第二张原始影像中第一区域A的放大图经过全局Wallis法处理后的影像图；
图9(c)是第二张原始影像中第一区域A的放大图经过Inpho法处理后的影像图；
图9(d)是第二张原始影像中第一区域A的放大图经过本发明的匀色方法处理后的影像图；
图10(a)是分析实验中第二张原始影像中第二区域B的放大图；
图10(b)是第二张原始影像中第二区域B的放大图经过全局Wallis法处理后的影像图；
图10(c)是第二张原始影像中第二区域B的放大图经过Inpho法处理后的影像图；
图10(d)是第二张原始影像中第二区域B的放大图经过本发明的匀色方法处理后的影像图；
图11是不同尺寸影像的并行加速对比图。
图中，A为第一区域；B为第二区域。
具体实施方式
下面结合附图对本发明做进一步详细的说明。
本发明提供一种基于最短传递路径的匀色方法，采用Voronoi图与Dijkstra算法相结合的方法确定匀色处理顺序；根据影像中心点位置构造Voronoi图，对区域范围内的影像有效的组织，便于邻近关系的查询；利用Dijkstra算法计算最短传递路径，对应处理的相邻影像具有较大重叠区域的相邻影像，传递次数较少，可以减少传递过程中的误差，提高区域内影像之间的色彩一致性。具体包括以下步骤：
1)获取设定区域内的各影像，并计算各影像的清晰度，以清晰度最大的影像作为初始参考影像。
初始参考影像的质量会影响区域内影像的整体质量，优选选取具有最大清晰度的影像作为初始参考影像。清晰度是利用影像的平均梯度来衡量影像在纹理和细节方面表达能力的质量评价参数，其值越大表示影像越清晰，影像质量越高，计算公式为：
Δx＝f(x+1,y)-f(x,y)
Δy＝f(x,y+1)-f(x,y)
上述公式中，M和N分别为影像的宽和高，f(x,y)为影像在(x,y)处的灰度值。
2)根据设定区域内的各影像的中心点构建Voronoi图，计算初始参考影像的中心点到其余各影像的中心点的最短路径。
如图1所示，为根据区域内的各影像的中心点构建Voronoi图，此图1为一种示意图，其中的v1、v2等为各影像的中心点，各中心点连接线上的数字代表了对应连接线的长度即两个中心点的距离；图1中虚线所示，每个中心点对应一个Voronoi多边形，根据Voronoi多边形之间是否存在公共边，即可判断对应的影像是否相邻。由于Voronoi图是以距离中心点最近为原则划分的，可以合理地判断多度重叠影像之间的邻近关系，只有距离较近的影像才会判断为相邻，对应影像的重叠度较大，可以提高匀色效果。
3)根据初始参考影像和待处理影像对应最短路径上的前一影像的匀色结果对该待处理影像进行Wallis匀色。
利用Dijkstra算法计算初始参考影像中心点到其他影像中心点的最短路径，即图1中v1到其余各点vi的最短路径，并记录vi的前一点vj，影像j的匀色结果为影像i的参考影像。
为了避免多次传递导致的色彩偏移，保持区域内影像整体的色彩一致性，步骤3)为根据初始参考影像的灰度均值和待处理影像对应最短路径上的前一影像匀色后的灰度均值以第一设定权重关系计算得到待处理影像的参考灰度均值；并且根据初始参考影像的标准偏差和待处理影像对应最短路径上的前一影像匀色后的标准偏差以第二设定权重关系计算得到待处理影像的参考标准偏差，根据待处理影像的参考灰度均值和参考标准偏差对待处理影像进行匀色处理。具体公式，利用初始参考影像的灰度均值和标准偏差对新的参考影像的灰度均值和标准偏差进行约束，设影像j的匀色结果为影像j'，在对影像i处理时，参考灰度均值mf和标准偏差sf的计算公式为：
mf＝wmj′+(1-w)m1
sf＝wsj′+(1-w)s1
式中，mj′和sj′分别为影像j'的灰度均值和标准偏差，m1和s1分别为初始参考影像的灰度均值和标准偏差，w∈[0,1]为权值常数，一般取0.5。
根据参考影像的灰度均值和标准偏差以及待处理影像的灰度均值和标准偏差进行Wallis变换，Wallis是一种特殊的线性变换法，通过将待处理影像的均值和标准偏差映射到参考影像的均值和标准偏差，实现不同影像间的色彩均衡。线性数学模型如下：
f(x,y)＝g(x,y)r1+r0
式中，g(x,y)为待处理影像的灰度值，f(x,y)为处理后的结果影像的灰度值。r0为加性系数，r1为乘性系数，可以表示为：
r0＝bmf+(1-b-r1)mg
式中，mg和mf分别为待处理影像和参考影像的灰度均值，sg和sf分别为待处理影像和参考影像的标准偏差。b为影像亮度系数，b∈[0,1]。c为影像方差扩展系数，c∈[0,1]。
Wallis变换的目的是将处理影像的灰度均值和标准偏差分别被强制到mf和sf，通常取b＝1、c＝1，此时线性数学模型变为：
f(x,y)＝[g(x,y)-mg]·(sf/sg)+mf。
具体的Wallis匀色法，包括以下步骤：
1)对待处理的影像进行互不重叠分块并统计各影像块的灰度均值和标准偏差。
Wallis匀色法是根据影像整体的均值和标准偏差，利用同一个线性关系对影像的每一个像素进行处理。但是影像中的地物复杂多样，颜色信息也各不相同，影像整体的均值和标准偏差无法准确地反映局部地物色彩特征，采用同一个线性关系显然是不合理的。为了解决上述问题，本发明对待处理的影像进行互不重叠分块，分块的个数为W×H，其中
W＝r×w,H＝r×h
式中，CV为待处理影像的变异系数，CVRef为参考影像的变异系数，w和h分别为预设的行、列方向的参考分块数。
采用分块策略时，匀色处理的质量受分块数目影响。若分块数目太多，即各影像块太小时，容易过度校正造成地物失真偏色并且计算量较大；若分块数目太少，即各影像块太大时，统计的均值方差不能准确反映地物分布，不能很好地消除影像之间的色彩差异。变异系数是标准偏差与均值之比，也称离散系数，可以描述影像内地物的丰富程度。当变异系数越大时，影像内地物种类越丰富，对应的分块数目应该越大即影像块越小，才能获得较好地匀色效果。
2)根据各影像块的灰度均值和标准偏差计算各影像块的四个角点对应的灰度均值和标准偏差。
其中，角点若只属于一个影像块，则将该影像块的灰度均值和标准偏差赋给该角点；角点若为多个相邻影像块之间的公共角点，则将所属多个影像块的灰度均值和标准偏差的平均值赋给该角点，最终计算得到各影像块的四个角点对应的灰度均值和标准偏差。
3)根据每个像素到所在影像块边缘的距离和该像素所在影像块的四个角点的灰度均值和标准偏差，计算每个像素的灰度均值和标准偏差。
如图2所示，对于影像块B(w,h)中的一点p(x,y)像素，其均值和标准偏差由角点Pw,h、Pw+1,h、Pw,h+1、Pw+1,h+1的均值和标准偏差以及该点到影像块边缘的距离Δx、Δy确定，利用双线性插值进行计算的公式如下：
式中，Δx为像素到m(w,h)对应的角点Pw,h的横向距离，Δy为像素到m(w,h)对应的角点Pw,h的纵向距离，m(x,y)和s(x,y)分别为像素的灰度均值和标准偏差，m(w,h)、m(w+1,h)、m(w,h+1)、m(w+1,h+1)分别为角点Pw,h、Pw+1,h、Pw,h+1、Pw+1,h+1的灰度均值，s(w,h)、s(w+1,h)、s(w,h+1)、s(w+1,h+1)分别为角点Pw,h、Pw+1,h、Pw,h+1、Pw+1,h+1的标准偏差，X和Y分别为该像素所在影像块B(w,h)的宽和高。
4)根据每个像素的灰度均值和标准偏差，并结合参考影像的灰度均值和标准偏差，对每个像素进行Wallis变换处理。
Wallis变换的方程为：
f(x,y)＝[g(x,y)-m(x,y)]·(sf/s(x,y))+mf
式中，g(x,y)为待处理影像的灰度值，f(x,y)为处理后的结果影像的灰度值，mf为参考影像的灰度均值，sf为参考影像的标准偏差。
利用双线性插值计算每个像素的线性变换参数，可以保证相邻影像块之间的平滑性。此外，利用影像的各角点而非中心点参与计算，可以避免影像边缘的分块出现锯齿现象。
本发明还提供一种基于最短传递路径的匀色处理装置，包括存储器、处理器以及存储在存储器中并可在处理器上运行的计算机程序，处理器执行程序实现上述的最短传递路径的运算，该最短传递路径的运算可以在处理器中的CPU端计算。
另外，该处理器还包括GPU端，该GPU端执行程序实现Wallis匀色方法；该Wallis匀色方法的主要计算任务包括三个步骤：计算各影像块的灰度均值和标准偏差、逐像素双线性插值和线性变换计算新的灰度值。分析可知：
(1)双线性插值和线性变换需要逐像素计算，是计算量最大的部分，各像素的计算是相互独立的，非常适合GPU并行处理，一个线程执行一个像素点的计算任务。
GPU端包括至少两个线程块，一个影像块对应一个线程块，影像块中的像素与对应该影像块的线程块中的线程一一对应，双线性插值和线性变换为重复的密集计算任务，并行化比较简单，可以直接分配给各个线程同时计算，满负荷利用计算资源，可以大幅度缩短计算时间。线程是GPU的最小执行单元，具体执行时，按照“线程格网—线程块—线程”的层次结构进行组织，如图3所示。假设待处理影像大小为M×N，分配相同大小的线程格网，设置线程块大小为l×k，则线程块数目为
一个影像块对应一个线程块，一个线程对应一个像素。各线程同时进行运算，并将结果按索引号赋给对应的像素。
(2)影像块灰度均值和标准偏差的本质均为累计求和计算，耦合度低，不能直接进行并行设计，因此处理器在对各影像块的灰度均值和标准偏差求解时通过并行归约求和法进行运算。
归约求和是一种缩减计算方法，基于对数步长交替两两求和，如图4所示，可以将求和的时间复杂度由O(N)(N为数据个数)降为O(log2N)。在每一个线程块内利用共享内存进行归约求和，可以达到并行加速的目的，并且交替策略可以避免存储片冲突并保持线程块的相邻线程处于活跃状态。
针对GPU的层次结构，本发明采用两遍归约求和策略。第一阶段内核执行n个并行归约，其中n是指线程块数，得到一个中间结果数组；第二个阶段通过调用一个线程块对这个中间数组进行归约，从而得到最终结果，如图5所示。具体步骤如下：
1、对输入数组落入每个线程中的数据进行求和。每个线程把它得到的累计值写入共享内存，交替因子为n*m，并在执行对数步长的归约前进行同步操作。
2、对共享内存中的值进行对数步长的归约操作。共享内存中后半部分的值被加到前半部分，即a[i]＝a[i]+a[m/2],(0≤i＜m/2)，参与的线程依次减半。在此操作执行log2m次后，共享内存中第一个线程对应的值a[0]即为该线程块的和。共享内存的大小等于线程块的线程数量m，并且m必须是2的幂次。
3、线程块的和写入全局内存。
利用并行归约求和方法计算影像像素值之和以及像素个数，即可求出影像均值，同理可求标准偏差，本发明的方法的流程如图6所示，CPU端运行最短路径运算，GPU端通过并行计算的方式进行逐像素的Wallis变换。
另外，为了达到尽可能高的并行加速比，本发明结合算法的自身特点，从配置划分、存储器带宽和指令吞吐量等方面进行优化。具体如下：
(1)合理组织线程
GPU中线程块独立被调度到流式多处理器中，来自同一个线程块的线程在同一个流式多处理器中执行。为了使流式多处理器的性能达到最优，流式多处理器中线程块的数量要小于8、线程数等于1536，线程块中的线程数要小于1024。分析可知，当线程块大小设置为256或512时，流式多处理器的性能最优，可以提高并行算法的运算速度。
(2)存储器优化
由于常量存储器的访问速度明显优于全局存储器，合理利用常量存储器代替全局存储器可以有效地减少内存流量和内存带宽。双线性插值中的四个角点的均值和标准偏差和线性变换中的参考均值和参考标准偏差数据量不大且每个线程都会读取，将这些参数分配为常量内存，可以大大加速访问数据速度和减少内存带宽，有效地提升程序的运行效率。
针对共享内存中的归约求和，为了减少不必要的线程同步，使用线程束同步优化。由于每个线程块中的线程束是按照锁步方式执行每条指令的，当线程块中的活动线程数低于硬件线程束的大小32时，无须再调用_syncthreads()内置函数。线程束同步可以避免每个Warp中出现分支导致效率低下，减小线程的闲置，提高并行度。
(3)指令优化
为了使指令吞吐量最大化，在满足精度的前提下，尽可能使用单精度float类型代替双精度double类型，使用硬件函数代替常规函数，用最少的指令完成相同的运算。
为了验证本发明的有效性和高效性，本发明设计两组对比实验分别从匀色质量和匀色效率两个方面进行验证。实验硬件环境为Intel(R)Core(TM)i5-6300HQ CPU 2.30GHz、8GB内存、NVIDIA GeForce GTX960M显卡和2GB显存的便携式计算机。软件环境为Windows 7操作系统、Microsoft Visual Studio 2010开发环境、C++编程语言和Cuda 7.5并行开发包。实验包括匀色质量对比实验和匀色效率对比实验，匀色质量对比实验如下：
为了验证匀色质量，分别利用全局Wallis匀色法、商业软件Inpho5.6以及本发明的方法对两组航空正射影像进行对比实验，结果如图7(a～d)和图8(a～d)所示，图9(a～d)和图10(a～d)为图8(a)中第一区域A和第二区域B的局部放大图。第一组实验数据(数据Ⅰ)，如图7(a)所示，为4幅存在明显的色彩不一致和对比度不一致的彩色影像，两两重叠度约为40％。第二组实验数据(数据Ⅱ)，如图8(a)所示，为5条航带的234张影像，航向重叠度约为70-80％，旁向重叠度约为40-50％。为了直观对比三种方法的匀色质量，图7(a～d)～图10(a～d)均为几何镶嵌的结果，未进行镶嵌缝羽化处理。
可以看出，全局Walllis匀色法只能在一定程度上减小影像间的色彩差异，影像间仍存在明显的色调差异，如图7(b)、图8(b)和图9(b)所示，这是由于影像整体的统计信息无法代表影像局部的地物特征，并且采用单一参考影像会产生过增强或偏色现象，如图10(b)所示。Inpho能够较好地消除影像间的色彩差异，使区域影像整体色调一致，但是影像整体对比度下降并且色调与原始影像产生了一定的偏差，如图8(c)和图9(c)所示，同时，当影像的色调与区域整体色调差异较大时，Inpho的处理结果不理想，如图7(c)和图10(c)所示。相比而言，本发明的方法对两组影像的处理结果都比较理想，在使区域影像整体色彩和对比度一致的同时消除或减小了相邻影像重叠区域的局部差异，并且没有产生偏色现象，同时，本发明自动选取了清晰度最大的影像作为参考影像，区域内影像整体清晰度较高，具有较好的目视效果。
为了对匀色影像进行定量评价，统计区域范围内相邻影像重叠区域的灰度均值之差的平均值
和标准偏差之差的平均值
计算公式如下：
式中，mij和mji分别为影像i和影像j在重叠区域中灰度均值，sij和sji分别为影像i和影像j在重叠区域中标准偏差，n为重叠区域的个数。
由表1可知，数据Ⅰ的原始影像的
和
都较大，说明原始影像间的色调差异和对比度差异都较大。三种方法匀色影像的
都小于原始影像，说明三种方法不同程度地减小了影像间的色彩差异。本发明的方法中
最小，Inpho次之，全局Wallis方法最大，说明本发明的方法消除色彩差异的能力最强。同理，本发明的方法中
最小，说明本发明的方法消除对比度不一致的能力最好。Inpho的
虽然较小，但是
要大于原始影像，说明Inpho消除影像间对比度差异的能力较差。而全局Wallis则相反，可以较好地消除影像间对比度不一致，但是消除影像间色彩差异的能力较弱。对比表2中统计结果，可得到类似于表1的结论，与目视评价结果一致。
由上可知，本发明的方法在目视效果评价和定量评价两方面均优于全局Wallis匀色法和Inpho，可以同时有效地消除影像间的色彩差异和对比度差异，体现了明显的优势。
表1
表2
匀色效率对比实验如下：
为了验证本发明并行加速策略的有效性，分别统计CPU串行算法和GPU并行算法处理不同尺寸影像的运行时间和加速比，分别如表3和图11所示。实验影像为数据Ⅱ的缩放影像，加速比为CPU串行算法运行时间与GPU并行算法运行时间的比值。统计的运行时间不包括硬盘与内存之间数据传输消耗的时间，统计均值和标准偏差的运行时间为统计影像整体均值和标准偏差以及统计分块均值和标准偏差的运行时间之和，分块数目均设置为8×8。
表3
分析表3和图11可知：
(1)当影像尺寸较小如512×512时，统计影像均值和标准偏差的串行算法比并行算法快。这是由于统计每一个影像块的均值和标准偏差时都存在内存与显存之间的数据传输，而GPU并行的加速时间小于数据传输消耗的时间。但是，双线性插值和线性变换的并行加速效果较好，因此，整个过程GPU并行算法的运行时间仍优于CPU串行算法的运行时间。
(2)随着影像尺寸的增加，GPU的计算资源得到更充分的利用，GPU并行算法的优势越来越明显，加速比越来越大。
(3)当影像尺寸增加到一定程度时，GPU并行算法的加速比趋于稳定不再明显增加，最高加速比可以达到60倍以上。
以上给出了本发明涉及的具体实施方式，但本发明不局限于所描述的实施方式。在本发明给出的思路下，采用对本领域技术人员而言容易想到的方式对上述实施例中的技术手段进行变换、替换、修改，并且起到的作用与本发明中的相应技术手段基本相同、实现的发明目的也基本相同，这样形成的技术方案是对上述实施例进行微调形成的，这种技术方案仍落入本发明的保护范围内。
Claims (12)
1.一种基于最短传递路径的匀色方法，其特征在于，包括以下步骤：
1)获取设定区域内的各影像，并计算各影像的清晰度，以清晰度最大的影像作为初始参考影像；
2)根据所述设定区域内的各影像的中心点构建Voronoi图，计算初始参考影像的中心点到其余各影像的中心点的最短路径；
3)根据待处理影像对应最短路径上的前一影像的匀色结果和初始参考影像对该待处理影像进行Wallis匀色；
步骤3)中当待处理影像为i，待处理影像i对应最短路径上的前一影像为j时，若影像j的匀色结果为j′，则在对待处理影像i处理时，对应的参考影像的灰度均值mf和标准偏差sf的计算公式为：
mf＝wmj′+(1-w)m1
sf＝wsj′+(1-w)s1
式中，mj′和sj′分别为影像j′的灰度均值和标准偏差，m1和s1分别为初始参考影像的灰度均值和标准偏差，w∈[0,1]为权值常数。
2.根据权利要求1所述的基于最短传递路径的匀色方法，其特征在于，所述Wallis匀色为分块处理Wallis匀色，该分块处理Wallis匀色包括以下步骤：
(1)对所述待处理影像进行互不重叠分块并统计各影像块的灰度均值和标准偏差；
(2)根据各影像块的灰度均值和标准偏差计算各影像块的四个角点对应的灰度均值和标准偏差；
(3)根据每个像素到所在影像块边缘的距离和该像素所在影像块的四个角点的灰度均值和标准偏差，计算每个像素的灰度均值和标准偏差；
(4)根据每个像素的灰度均值和标准偏差，并结合参考影像的灰度均值和标准偏差，对每个像素进行Wallis变换处理。
3.根据权利要求2所述的基于最短传递路径的匀色方法，其特征在于，步骤(2)中角点若只属于一个影像块，则将该影像块的灰度均值和标准偏差赋给该角点；角点若为多个相邻影像块之间的公共角点，则将所属多个影像块的灰度均值和标准偏差的平均值赋给该角点。
4.根据权利要求2或3所述的基于最短传递路径的匀色方法，其特征在于，步骤(1)中对待处理的影像进行互不重叠分块，分块的个数为W×H，其中
W＝r×w,H＝r×h
式中，CV为影像的变异系数，CVRef为参考影像的变异系数，w和h分别为预设的行、列方向的参考分块数。
5.根据权利要求4所述的基于最短传递路径的匀色方法，其特征在于，所述Wallis变换的公式为：
f(x,y)＝[g(x,y)-m(x,y)]·(sf/s(x,y))+mf
式中，g(x,y)为待处理影像的灰度值，f(x,y)为处理后的结果影像的灰度值，mf为参考影像的灰度均值，sf为参考影像的标准偏差，m(x,y)和s(x,y)分别为像素的灰度均值和标准偏差。
6.根据权利要求2所述的基于最短传递路径的匀色方法，其特征在于，步骤(3)中每一个像素的灰度均值和标准偏差的计算为并行计算。
7.根据权利要求6所述的基于最短传递路径的匀色方法，其特征在于，步骤(1)中各影像块的灰度均值和标准偏差为通过并行归约求和法求得。
8.一种基于最短传递路径的匀色处理装置，包括存储器、处理器以及存储在存储器中并可在处理器上运行的计算机程序，其特征在于，所述处理器执行所述程序时实现以下步骤：
1)获取设定区域内的各影像，并计算各影像的清晰度，以清晰度最大的影像作为初始参考影像；
2)根据所述设定区域内的各影像的中心点构建Voronoi图，计算初始参考影像的中心点到其余各影像的中心点的最短路径；
3)根据待处理影像对应最短路径上的前一影像的匀色结果和初始参考影像对该待处理影像进行Wallis匀色；
步骤3)中当待处理影像为i，待处理影像i对应最短路径上的前一影像为j时，若影像j的匀色结果为j′，则在对待处理影像i处理时，对应的参考影像的灰度均值mf和标准偏差sf的计算公式为：
mf＝wmj′+(1-w)m1
sf＝wsj′+(1-w)s1
式中，mj′和sj′分别为影像j′的灰度均值和标准偏差，m1和s1分别为初始参考影像的灰度均值和标准偏差，w∈[0,1]为权值常数。
9.根据权利要求8所述的基于最短传递路径的匀色处理装置，其特征在于，所述Wallis匀色为分块处理Wallis匀色，该分块处理Wallis匀色包括以下步骤：
(1)对所述待处理影像进行互不重叠分块并统计各影像块的灰度均值和标准偏差；
(2)根据各影像块的灰度均值和标准偏差计算各影像块的四个角点对应的灰度均值和标准偏差；
(3)根据每个像素到所在影像块边缘的距离和该像素所在影像块的四个角点的灰度均值和标准偏差，计算每个像素的灰度均值和标准偏差；
(4)根据每个像素的灰度均值和标准偏差，并结合参考影像的灰度均值和标准偏差，对每个像素进行Wallis变换处理。
10.根据权利要求9所述的基于最短传递路径的匀色处理装置，其特征在于，所述处理器包括GPU端，所述分块处理Wallis匀色的步骤在GPU端运行，所述GPU端包括至少两个线程块，一个影像块对应一个线程块，影像块中的像素与对应该影像块的线程块中的线程一一对应；步骤(3)中一个像素的灰度均值和标准偏差的计算对应一个线程，所有线程并行运算。
11.根据权利要求10所述的基于最短传递路径的匀色处理装置，其特征在于，所述GPU端在对各影像块的灰度均值和标准偏差求解时通过并行归约求和法进行运算。
12.根据权利要求10或11所述的基于最短传递路径的匀色处理装置，其特征在于，所述存储器包括常量存储器，所述常量存储器用于存储四个角点的均值和标准偏差以及参考影像的灰度均值和标准偏差。
Priority Applications (1)
| Application Number | Priority Date | Filing Date | Title |
|---|---|---|---|
| CN201811178782.2A CN109410136B (zh) | 2018-10-10 | 2018-10-10 | 基于最短传递路径的匀色方法及处理装置 |
Applications Claiming Priority (1)
| Application Number | Priority Date | Filing Date | Title |
|---|---|---|---|
| CN201811178782.2A CN109410136B (zh) | 2018-10-10 | 2018-10-10 | 基于最短传递路径的匀色方法及处理装置 |
Publications (2)
| Publication Number | Publication Date |
|---|---|
| CN109410136A CN109410136A (zh) | 2019-03-01 |
| CN109410136B true CN109410136B (zh) | 2020-10-27 |
Family
ID=65466927
Family Applications (1)
| Application Number | Title | Priority Date | Filing Date |
|---|---|---|---|
| CN201811178782.2A Active CN109410136B (zh) | 2018-10-10 | 2018-10-10 | 基于最短传递路径的匀色方法及处理装置 |
Country Status (1)
| Country | Link |
|---|---|
| CN (1) | CN109410136B (zh) |
Families Citing this family (3)
| Publication number | Priority date | Publication date | Assignee | Title |
|---|---|---|---|---|
| CN111754590B (zh) * | 2020-05-14 | 2024-04-02 | 北京吉威空间信息股份有限公司 | 基于全球色彩特征库遥感影像自动匀色的方法 |
| CN111652826B (zh) * | 2020-05-18 | 2023-04-25 | 哈尔滨工业大学 | 基于Wallis滤波+直方图匹配的有重多/高光谱遥感图像匀色方法 |
| CN112164006B (zh) * | 2020-09-25 | 2023-11-03 | 航天宏图信息技术股份有限公司 | 一种影像匀色方法、装置、电子设备及存储介质 |
Citations (5)
| Publication number | Priority date | Publication date | Assignee | Title |
|---|---|---|---|---|
| CN101770639A (zh) * | 2010-01-14 | 2010-07-07 | 北京航空航天大学 | 一种低照度图像增强方法 |
| CN105528797A (zh) * | 2015-12-02 | 2016-04-27 | 深圳飞马机器人科技有限公司 | 一种光学影像色彩一致性自适应处理与快速镶嵌方法 |
| CN105957111A (zh) * | 2016-04-27 | 2016-09-21 | 武汉大学 | 序列遥感影像的色调一致性校正方法及系统 |
| CN106485664A (zh) * | 2015-08-27 | 2017-03-08 | 上海沃韦信息科技有限公司 | 一种基于小波变换和Wallis变换的卫星图像色彩平衡方法 |
| CN108550129A (zh) * | 2018-04-20 | 2018-09-18 | 北京航天宏图信息技术股份有限公司 | 基于地理模板的匀色方法及装置 |
Family Cites Families (1)
| Publication number | Priority date | Publication date | Assignee | Title |
|---|---|---|---|---|
| US9594955B2 (en) * | 2014-05-08 | 2017-03-14 | Rolta India Ltd | Modified wallis filter for improving the local contrast of GIS related images |
-
2018
- 2018-10-10 CN CN201811178782.2A patent/CN109410136B/zh active Active
Patent Citations (5)
| Publication number | Priority date | Publication date | Assignee | Title |
|---|---|---|---|---|
| CN101770639A (zh) * | 2010-01-14 | 2010-07-07 | 北京航空航天大学 | 一种低照度图像增强方法 |
| CN106485664A (zh) * | 2015-08-27 | 2017-03-08 | 上海沃韦信息科技有限公司 | 一种基于小波变换和Wallis变换的卫星图像色彩平衡方法 |
| CN105528797A (zh) * | 2015-12-02 | 2016-04-27 | 深圳飞马机器人科技有限公司 | 一种光学影像色彩一致性自适应处理与快速镶嵌方法 |
| CN105957111A (zh) * | 2016-04-27 | 2016-09-21 | 武汉大学 | 序列遥感影像的色调一致性校正方法及系统 |
| CN108550129A (zh) * | 2018-04-20 | 2018-09-18 | 北京航天宏图信息技术股份有限公司 | 基于地理模板的匀色方法及装置 |
Non-Patent Citations (3)
| Title |
|---|
| "An auto-adapting global-to-local color balancing method for optical imagery mosaic";Lei Yu 等;《ISPRS Journal of Photogrammetry and Remote Sensing》;20171031;第132卷;1-19 * |
| "最小传递路径的正射影像色彩一致性处理方法";李烁等;《测绘科学技术学报》;20161231;第33卷(第6期);第3节 * |
| "遥感影像色彩一致性处理技术研究";易磊;《中国优秀硕士学位论文全文数据库-信息科技辑》;20160715;第2016年卷(第7期);I140-349 * |
Also Published As
| Publication number | Publication date |
|---|---|
| CN109410136A (zh) | 2019-03-01 |
Similar Documents
| Publication | Publication Date | Title |
|---|---|---|
| CN109300083B (zh) | 一种分块处理Wallis匀色方法及装置 | |
| Wehrens et al. | Flexible self-organizing maps in kohonen 3.0 | |
| CN109410136B (zh) | 基于最短传递路径的匀色方法及处理装置 | |
| WO2018040463A1 (zh) | DeMura表的数据压缩、解压缩方法及Mura补偿方法 | |
| US7538779B2 (en) | Method of rendering pixel images from abstract datasets | |
| US8330763B2 (en) | Apparatus and method for volume rendering on multiple graphics processing units (GPUs) | |
| CN111028327B (zh) | 一种三维点云的处理方法、装置及设备 | |
| CN107657599B (zh) | 基于混合粒度划分和动态负载分配的遥感图像融合系统并行实现方法 | |
| US8872826B2 (en) | System and method for decoupled ray marching for production ray tracking in inhomogeneous participating media | |
| Weiss et al. | Differentiable direct volume rendering | |
| CN111798467A (zh) | 一种图像分割方法、装置、设备及存储介质 | |
| WO2018113502A1 (zh) | 一种自动生成网格与着色器多层次细节的方法 | |
| Happ et al. | A region-growing segmentation algorithm for GPUs | |
| CN113744136A (zh) | 基于通道约束多特征融合的图像超分辨率重建方法和系统 | |
| Cui et al. | Real-time stereo vision implementation on Nvidia Jetson TX2 | |
| CN114745532A (zh) | 混合色温场景白平衡处理方法、装置、存储介质及终端 | |
| US20090322781A1 (en) | Anti-aliasing techniques for image processing | |
| DE102013021044A1 (de) | Erzeugung fehlerbefreiter Voxel-Daten | |
| US20110052059A1 (en) | Generating image histogram by parallel processing | |
| CN116228753B (zh) | 肿瘤预后评估方法、装置、计算机设备和存储介质 | |
| CN115564812A (zh) | 基于高精度视差细化的立体匹配方法、系统、设备及介质 | |
| Carabaño et al. | Efficient implementation of a fast viewshed algorithm on SIMD architectures | |
| CN115994935A (zh) | 应用于激光通讯的光斑质心确定方法及系统 | |
| CN112598663B (zh) | 基于视觉显著性的粮食害虫检测方法和装置 | |
| Gao et al. | Single image dehazing based on single pixel energy minimization |
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 |