Now Loading ...
-
How to compute integral image on GPU?
Related Work
GPU Gems 3 (2005). Chapter: “Summed-Area Tables on the GPU” by Dave Johnson (NVIDIA).
https://developer.nvidia.com/gpugems/gpugems3/part-ii-light-and-shadows/chapter-8-summed-area-variance-shadow-maps
Parallel Prefix Sums
https://developer.nvidia.com/gpugems/gpugems3/part-vi-gpu-computing/chapter-39-parallel-prefix-sum-scan-cuda
Scan primitives for GPU computing
GPU-efficient recursive filtering and summed-area tables
https://dl.acm.org/doi/10.1145/2070781.2024210
技术重点：
并行前缀和（Scan）：大多数GPU实现将积分图像计算分解为两个主要的并行前缀和传递，一个是水平方向的，另一个是垂直方向的。
内存合并与冲突：高效的GPU积分图像算法注意内存布局，确保线程访问连续的内存段以减少延迟并提高吞吐量。
工作负载划分（Tile）：大型图像被划分为独立且并行计算的小块。在每个小块上计算出部分总和后，额外进行一次合并操作以得到全局积分图像。
层次化方法：层次化求和策略通过将求和问题拆解成多个阶段来降低处理大型图像时所需处理复杂度。
使用共享内存与寄存器：高性能实现将中间总和保存在在快速GPU共享内存中，从而减少全局内存流量。
Summed-Area Variance Shadow Maps(GPU gem3)
https://developer.nvidia.com/gpugems/gpugems3/part-ii-light-and-shadows/chapter-8-summed-area-variance-shadow-maps
该章节指出SAT的计算本质上是两个维度的前缀和过程。Hensley et al. 2005.提出了一种更并行的方法。
本章节主要关注积分图像计算时的精度问题，Hensley et al. 2005有很多关于精度问题的解决方案， Donnelly and Lauritzen 2006也给出了一种解决方案。
Parallel Prefix Sum (Scan) with CUDA（GPU Gem3）
Naive Method
CPU的前缀和计算的复杂度为O(n),GPU算法的复杂度如果不超过该复杂度，我们称之为efficient work。
最原始的GPU版本就是一个不高效的实现，该算法需要两个Buffer来保证正确性。同时，该算法还假设处理器单元的数量和数组大小一致。
Work-Efficient Method
使用平衡树来优化复杂度。
分成两个过程，Red1uce(归约过程，也被称作上扫过程)；和下扫过程。
同一个位置上的操作在共享内存中完成。
Reduce
该过程结束后，根结点上保存了整个数组的和。
down-sweep
在上扫描阶段，我们已经获得了总和，但前缀和还未完全计算出来。下扫描阶段通过逆向操作，逐步填充前缀和数组。
上扫描结束后，位置0，1，3，7的值是前缀和数组中会出现的值（为inclusive的，但通常需要计算exclusive的）。例如，25的右子树14需要11的值来更新子节点。
因此，下扫描过程可以理解为一个传递到左子树、累加右子树的过程。
初始化：设置根节点为0.
遍历步长数组[4,2,1]
设置左子树为根节点的值，右子树为左子树和右子树的和（左子树使用设置后的值）。
该过程一个线程操作两个元素，需要在同一个work group中完成。
使用上扫阶段的结果：
树结构：
25
/ \
11 14
/ \ / \
4 7 5 9
/ \ / \ / \ / \
3 1 7 0 4 1 6 3
初始化根节点为 0：
0
/ \
11 14
/ \ / \
4 7 5 9
/ \ / \ / \ / \
3 1 7 0 4 1 6 3
传递并更新：
根节点 0 分配给左子树 left = 0 和右子树 right = 11（左子树的部分和）。
左子树 11 分配给其左子树 left = 0 和右子树 right = 4。
右子树 14 分配给其左子树 left = 11 和右子树 right = 5。
继续递归，最终得到前缀和数组。
最终前缀和数组 S = [0, 3, 4, 11, 11, 15, 16, 22]。
Fast Summed‐Area Table Generation and its Applications[Hensley et al. 2005]
提到论文：Simple Blurry Reflections with Environment Maps 似乎是使用miomap近似模糊的lerp blur。
方法
使用图形管线实现，整个过程分解为两个phase，每个phase包括log(n)个pass。
使用两张纹理图像，互相作为输入输出。
积分图像的精度问题
累加时机器误差的积累。
使用积分图像时通常是两个接近的值相减，尤其是两个值很大时，误差会更加明显。
积分图像的最大值大小为$w \times h \times 255$
Using Signed-Offset Pixel Representation
将像素值进行偏移，从[0,1]偏移到[-0.5,0.5],这样可以使得：
数值不总是单调的。
有两种方式实现这一点：
所有像素值-0.5.
所有像素值-平均值。
Using Origin-Centered Image Representation
以图片的中心点作为原点，相当于计算四张更小的积分图像，从而避免了极大值的出现。
但是使用时增加了额外的计算量。
Scan Primitives for GPU Computing 2007
GPU并行模型
符合流水线结构的算法天然适合GPU实现。即，每个kernel处理单独的输入，产生单独的输出。
某些问题，例如前缀求和问题，需要输入数据的全局知识。
相关工作
原始方法
horn 2005
[[#Fast Summed‐Area Table Generation and its Applications[Hensley et al. 2005]]]
这些方法的复杂度为O(nlogn)，是non-work-efficient的。
reduce and down-sweep
提出：Blelloch in 1990[Vector Models for Data-Parallel Computing]
GPU实现：Sengupta(本文作者) et al. and Greß et al. in 2006
CUDA实现：本文
方法
本文方法的主要贡献是引入了分段扫描的概念.
分段扫描Segmented Scan，是一种将输入序列分成多个互不影响的字序列并进行后续操作的过程。
使用flag数组标记不同的分段，例如：
Data: [a, b, c, d, e]
Flags: [1, 0, 1, 0, 0]
该算法用于解决输入向量超过线程块大小的情况。
在归约过程和下扫描之间插入部分和的计算。
（影响下扫描的初始化）
Efficient Integral Image Computation on the GPU 2010
使用work-efficient的前缀和计算，利用了分段扫描解决大数组问题。积分图像的计算转化为：前缀和-转置-前缀和过程。
使用tanspose操作的优点在于，可以使用同一个kernel完成两次扫描操作。
GPU-Efficient Recursive Filtering and Summed-Area Tables 2011
Efficient Algorithms for the Summed Area Tables Primitive on GPUs 2018
Related work
Scan-scan algorithms:
Compute prefix sums directly
Limited by memory access patterns
Scan-transpose-scan algorithms:
Use matrix transposition between steps
Rely heavily on scratchpad memory
Have expensive matrix transpose operations
Face memory bandwidth limitations
该方法来自于Efficient Integral Image Computation on the GPU。本文指出，该方法需要对global memory的聚合访问和昂贵的transpose操作。
[!聚合访问]
对全局内容的访问较慢，只有每一个线程以相同stride的方式访问全局内存，才能实现峰值吞吐：
Thread 0: accesses address N
Thread 1: accesses address N+1
Thread 2: accesses address N+2
…and so on
[!bank conflict]
现代GPU的共享内容被分成一系列的bank。当同一个warp中的不同线程访问同一个bank时，会发生冲突。
例如，一个warp中存在12个线程，每个线程都需要访问共享内存对应行的32个数据，即第i个线程访问sharedMem[i][0…31]由于bank的存在，每个线程访问数据j时都对应到同一个bank中，因此存在冲突。
本文为了避免这种冲突，定义数组大小为[32][33]
本文方法的动机
SAT 计算的瓶颈在于数据移动。以聚合模式高效访问全局内存对算法性能至关重要。此外，设计减少为了减少数据移动，我们采用了寄存器缓存方法。问题在于如何在有效使用寄存器的同时避免争用。
线程间通信寄存器的技术（即shuffle指令）只能在单个 warp（CUDA 中一起执行的一组线程）内工作。因此，我们调整了算法，在warp级别进行扫描，同时避免warp内部通信。
最直接的SAT计算
Prefix Sum
Basic method
Kogge-Stone
http://lumetta.web.engr.illinois.edu/408-S20/slide-copies/ece408-lecture16-S20.pdf
Step1得到长度小于等于2的累加和
Step2得到长度小于等于4的累加和
Step3得到长度小于等于8的累加和
LF-scan
warp shuffle
In theory, the LFscan achieves the highest computing efficiency with logN stages and $\frac{NlogN}{2}$
addition operations as Fig. 2c shows,
‘reduce (up-sweep) and down-sweep’
https://developer.nvidia.com/gpugems/gpugems3/part-vi-gpu-computing/chapter-39-parallel-prefix-sum-scan-cuda
问题
算法流程决定加法运算的次数和算法阶段数。
数据存储、传递的方式(GPU显存访问模式)也会显著影响算法性能。
本文贡献
使用寄存器存储中间数据。
提出BRLT
方法
Caching Data Using Register Files
SAT的计算问题是内存限制的。
每个线程都先使用寄存器缓存数据。
并且warp中线程的通信泗洪shuffle操作完成。
GPU聚合访问
用于降低对Global Memory的访问成本。具体做法是，在同一个warp中的所有线程访问临近的内存位置。GPU可以将这些操作打包成一次数据传输，从而降低延迟、提升带宽。
data向16 bits or 32 bits对齐。
globalID = (x + width * y + width * height * z),对内存位置的访问和GlobalID一致（尤其是2D图像，y连续访问不是聚合访问）。
Google
Ch
问题回顾
计算前缀和在两个方向上进行。
但是：
两个方向需要使用不同的kernel函数（不同的着色器程序，或者引入不同的分支）。
列方向处理的数据不连续，无法实现聚合访问。
因此，SAT计算更常见的做法是Scan-transpose-scan。
接下来需要考虑矩阵转置是否是聚合操作。
矩阵转置
warp level
element level
在每个warp内部：
Method：直接做法
写入操作不是聚合操作。
本文提出了BRLT，在Shared memory中对warp中的数据进行转置。
BRLT
Block-Register-Local-Transpose Method
这是一种将数据从寄存器复制到共享内容再复制回寄存器的矩阵transpose方法。
这种方法的仅使用共享内存作为缓冲区，transpose操作是在寄存器-共享内存之间完成的；传统算法的操作是在共享内存中完成，并且需要从主存中加载数据。
Bank conflict
google
GPU共享内存以bank的形式被组织。N卡通常32banks，bank中数据单元大小为4或8。当一个warp中的不同线程访问到同一个bank时，会发生conflict，导致IO操作无法并行。
例如，bank number=4.该任务是一个矩阵转置任务，Clock0时，四个线程写入第一列：
此时发生bank conflict。
但是如果修改共享内存大小为4*5:
部分和计算
将warp计算出的部分和存储在共享内存中。
在共享内存中计算数据的前缀和。
将数据从共享内存加载回对应的warp。
方法一：ScanRow-BRLT
Efficient Integral Image Computation on the GPU 2010中方法需要将行扫描的结果存储在全局内存中，transpose操作也在全局内存中完成。
本文中的方法:
将输入2D图像的一块直接加载到block中,
然后执行行扫描-BRLT-写回主存。因此转置操作在共享内存中完成。
如上图所示，每个block负责32行的数据，每个warp处理32*32的块。
每个线程读取32个数据，线程之间使用LF-scan（利用shuffle操作）。
warp之间的通信：
在一个线程中串行操作，复杂度O(n^2)：
思考
warp是被调度的基本单元，但并不意味着不同warp之间有先后顺序，也不意味着warp被完整调度。
为什么转置以后再写回主存，直接rigister-shared memory-rigister的模式完成两次算法IO更少。
如何处理更大的图像。
方法二：BRLT-ScanRow
与方法1类似，但是先转置再执行scan操作。重复两次。但是论文指出，该方法使用串行扫描算法，并且这种方法效率更高。
warp之间使用共享内存通信，但是该方法的每个线程只需要一个位置的preSum，复杂度降低为O(n)
方法三：Register-based ScanRowColumn Method
每个block处理32行的数据，每个warp处理一行的数据。
例如warp0，扫描0-31号元素，使用shuffle操作传递到warp0的第二次扫描的第一个线程上。
-
Vulkan Compute Shader
Reference
P.A. Minerva’s Vulkan Tutorial on Compute Shaders
vulkan tutorial
AMD opengpu
GPU Program Blog
Example for Compute Shader
Intro
核心组成：GPU由数千个处理核心组成，这些核心专门用于执行由许多核心同时处理的单独任务。GPU通过并行执行相同的指令序列（描述特定任务）来实现高速处理大量数据。
核心组织：这些核心被组织成多处理器（multiprocessors），能够在多个线程上并行执行相同的指令。这种硬件模型被称为SIMT架构（Single-Instruction, Multiple-Thread）。
多处理器内存：每个多处理器包括以下几种类型的片上内存：
32位寄存器，分配给核心。
共享内存，由所有核心共享。
只读常量缓存，加速从设备内存的只读区域读取。
只读纹理缓存，加速从纹理内存空间读取。
线程组和执行：多处理器创建、管理、调度并执行称为warp（或wavefront）的32或64个并行线程组。每个warp在同一程序地址开始，但具有自己的指令地址计数器和寄存器状态，因此可以独立分支和执行。
线程块和warp调度：多处理器接收一个或多个线程块进行执行时，将它们划分为warps，并由warp调度器管理和调度。线程块到warp的划分是一致的，每个warp包含连续递增的线程ID。
执行上下文和调度：每个warp的执行上下文在其整个生命周期内都保持在芯片上，因此从一个执行上下文切换到另一个没有成本。每个指令发出时，warp调度器可以选择一个准备好执行下一指令的warp，并向这些线程发出指令。
线程块（Thread Blocks）
定义：线程块是一组在GPU上同时执行的线程。它是由程序员定义的，用于组织和执行并行任务。
大小和形状：线程块的大小（即包含的线程数）和形状（如1D、2D或3D）可以根据特定的计算任务进行调整。
资源共享：线程块内的线程可以共享一定量的快速访问内存（称为共享内存），并且可以进行同步操作。
Warp（线程束）
定义：warp是线程块中的一小部分线程，这些线程在GPU上以单一的指令流同时执行相同的操作。在NVIDIA的GPU中，一个warp通常包含32个线程。
硬件调度单位：warp是GPU硬件调度和执行的基本单位。GPU的warp调度器负责管理这些warp的执行。
线程块与Warp之间的关系
线程块划分为Warp：当线程块被提交到GPU执行时，它被划分为多个warp。这个划分是自动进行的，基于warp的大小（如32个线程）。
连续线程ID：每个warp包含具有连续线程ID的线程。例如，在32线程的warp中，第一个warp包含线程ID 0-31，第二个warp包含线程ID 32-63，依此类推。
并行执行：线程块内的所有warp可以在GPU上并行执行，但每个warp内的线程同时执行相同的指令。
执行效率：合理地组织线程块和warp对于实现高效的GPU并行计算至关重要。线程块的大小应该是warp大小的整数倍，以最大化GPU核心的利用率并减少空闲线程。
SIMD vs SIMT
GPU并不提供对SIMD的支持，因此考虑并行应该是SIMT，多个线程执行相同的指令。
GPU可以理解为一系列Multiprocessor的集合，如上图。Multiprocessor内部的processor可以共享一部分数据。
因此线程组的数量应该和multiprocessor数量对应(线程组可以多于multiprocessor数量)。
而warp是一组线程，是最小的调度单位。warp内的线程并行执行，但不同的warp之间是否并行取决于资源的分配（实际上这是由于某些操作会被挂起等待，转而执行其他的warp）。
CS
Thread Blocks：一个三维逻辑结构
Coding
Check maxComputeWorkGroupInvocations and maxComputeWorkGroupSize in VkPhysicalDeviceLimits
查询线程组数量和线程组大小限制
线程组的数量由API制定，线程组的大小在GLSL代码中被指定
Sync
共享内存
// Array allocated in shared memory to be shared by all invocations in a work group.
// Check VkPhysicalDeviceLimits::maxComputeSharedMemorySize for the maximum
// total storage size, in bytes, available for variables shared by all invocations in a work group.
shared vec4 gCache_0[256];
在一个线程组中的线程的执行很难确定先后顺序，如果多个线程访问到相同的数据，可能会造成线程冲突。
void main()
{
ivec2 textureLocation = ivec2(gl_GlobalInvocationID.xy);
// Read texel from input image at mipmap level 0 and
// save the result in the shared memory
gCache_0[gl_LocalInvocationID.x] = imageLoad(InputTex, textureLocation);
// Wait for all invocations in the work group
barrier();
// OK!
vec4 left = gCache_0[gl_LocalInvocationID.x - 1];
vec4 right = gCache_0[gl_LocalInvocationID.x + 1];
// ...
}
使用屏障保证数据写入完成。
barrier()函数的作用
同步工作组内的着色器调用：在计算着色器中，barrier()函数确保在一个工作组（workgroup）内的所有着色器调用都到达这个屏障点之前，任何一个着色器调用都不会继续执行超过这个点的代码。
在曲面细分控制着色器中的应用：barrier()函数也可以在曲面细分控制着色器中使用，以同步单个输入补丁的所有着色器调用。
控制流的一致性
控制流必须是一致的：在计算着色器中使用barrier()时，控制流必须是一致的。这意味着如果任何一个着色器调用进入了一个条件语句，那么所有的调用都必须进入它。
barrier()与内存同步
控制流和共享变量的同步：barrier()函数影响控制流，并同步对共享变量的内存访问（以及曲面细分控制输出变量）。
其他内存访问：对于非共享变量的内存访问，barrier()函数并不能保证一个着色器调用在barrier()之前写入的值可以被其他调用在barrier()之后安全地读取。
使用内存屏障函数
确保内存访问的顺序：为了确保一个着色器调用写入的值可以被其他调用安全地读取，需要同时使用barrier()和特定的内存屏障函数。
内存屏障函数的作用：内存屏障函数用于对可被其他着色器调用访问的内存中的变量进行读写操作的排序。这些函数在被调用时，会等待调用者之前执行的所有读写操作的完成，然后再返回。
不同类型的内存屏障函数
特定类型变量的内存屏障：如memoryBarrierAtomicCounter()、memoryBarrierBuffer()、memoryBarrierImage()和memoryBarrierShared()等函数，分别用于等待对原子计数器、缓冲区、图像和共享变量的访问的完成。
全局内存屏障：memoryBarrier()和groupMemoryBarrier()函数用于等待对所有上述变量类型的访问的完成。
着色器类型的可用性：memoryBarrierShared()和groupMemoryBarrier()只在计算着色器中可用，而其他函数在所有着色器类型中都可用。
例子：粒子系统
使用GPU更新顶点位置，避免了受到总线带宽的限制。
SSBO
需要一个能读能写的buffer，SSBO可以满足这一点。
flag为：
VK_BUFFER_USAGE_VERTEX_BUFFER_BIT | VK_BUFFER_USAGE_STORAGE_BUFFER_BIT | VK_BUFFER_USAGE_TRANSFER_DST_BIT
VK_MEMORY_PROPERTY_DEVICE_LOCAL_BIT
Buffer type
. VK_BUFFER_USAGE_STORAGE_BUFFER_BIT
• Description: This flag indicates that the buffer can be used as a storage buffer. Storage buffers are used to store large blocks of data that shaders (especially compute shaders) can read from and write to.
• Common Use Case:
• It is commonly used in compute shaders or fragment shaders where you need random-access read and write operations.
• It’s suitable for large datasets that may change frequently, like results of computations or intermediate data.
• Access in Shaders:
• In GLSL (Vulkan’s shading language), a buffer marked with VK_BUFFER_USAGE_STORAGE_BUFFER_BIT is typically accessed using the layout(std430) storage qualifier, like this:
• Both read and write operations are possible in shaders.
VK_BUFFER_USAGE_TRANSFER_SRC_BIT
• Description: This flag indicates that the buffer can be used as a source for data transfer operations. Specifically, this buffer can be used in a memory transfer operation, where data from this buffer will be copied to another buffer or image (e.g., via vkCmdCopyBuffer or vkCmdCopyBufferToImage).
• Common Use Case:
• When you want to copy data from one buffer to another buffer or image, the source buffer should be created with this flag.
• It’s useful when you’re doing staging operations: you might upload data to a buffer that’s visible to the CPU (with VK_MEMORY_PROPERTY_HOST_VISIBLE_BIT), then transfer it to a GPU-only buffer with VK_MEMORY_PROPERTY_DEVICE_LOCAL_BIT.
VK_BUFFER_USAGE_VERTEX_BUFFER_BIT用于顶点着色器，VK_BUFFER_USAGE_STORAGE_BUFFER_BIT用于计算着色器的写入和读取。
VK_BUFFER_USAGE_TRANSFER_DST_BIT表明该buffer是数据复制的目的buffer。例如从staging buffer到该buffer。
VK_BUFFER_USAGE_TRANSFER_SRC_BIT表明该buffer是数据复制的起源。例如staging buffer。
Storage Image
借助SI可以完成对图片的操作，例如后处理、生成mip-maps等。
信号量和栅栏
计算cmd需要栅栏避免冗余提交（确保之前的
指令已经被执行再提交，cpu和gpu之间），计算管线和图形管线之间需要信号量，保证图形管线开始时计算任务已经完成（gpu内部）。
[!note]
解释了为什么信号量以数组形式给出：
VkSemaphore waitSemaphores[] = { computeFinishedSemaphores[currentFrame], imageAvailableSemaphores[currentFrame] };
延伸：
Asynchronous compute
[[TBDR]]
Touch background to close
Calendar's Blog
