CUDA Cpp正电子发射断层扫描仪校准和图像重建—蒙特卡洛3D伊辛模型

本文主要是介绍CUDA Cpp正电子发射断层扫描仪校准和图像重建—蒙特卡洛3D伊辛模型,希望对大家解决编程问题提供一定的参考价值,需要的开发者们随着小编来一起学习吧!

要点

  1. GPU对比CPU计算正弦和:使用单CPU、使用OpenMP库和CUDA
  2. CUDA并行计算:3D网格运行内核:线程块,线程线性处理3D数组,并行归约,共享内存,矩阵乘法/平铺矩阵乘法,基本线性代数子程序
  3. 平铺分区,矢量加载,warp级内在函数和子 warp,线程发散和同步,联合组
  4. 使用 2D 和 3D 模板,迭代求解偏微分方程和图像处理
  5. 使用 GPU 纹理硬件执行快速插值,图像配准
  6. 蒙特卡洛模拟3D伊辛模型
  7. CUDA 流
  8. CUDA正电子发射断层扫描仪校准和图像重建
  9. GPU扩展

矩阵乘法示例

假设我们有两个矩阵, A A A B B B。假设 A A A是一个 n × m n \times m n×m矩阵,这意味着它有 n n n行和 m m m列。还假设 B B B m × w m \times w m×w 矩阵。乘法 A ∗ B A * B AB(与 B ∗ A B * A BA不同)的结果是一个 n × w n \times w n×w矩阵,我们称之为 M M M。也就是说,结果矩阵的行数等于第一个矩阵 A A A 的行数和第二个矩阵 B B B​ 的列数。

为什么会发生这种情况以及它是如何运作的?这里两个问题的答案是相同的。我们以 M M M 的单元格 1,1(第一行,第一列)为例。运算 M = A ∗ B M=A * B M=AB 后其中的数字是 A A A 第 1 行中的数字与 B B B 第 1 列中的数字的所有逐元素乘法之和。也就是说,在 M M M 的单元格 i , j i, j i,j 中,我们得到了 A A A 中第 i i i 行和 B B B 中第 j j j​ 列中所有数字的逐元素乘法之和。

下图直观地解释了这个想法:

现在应该很清楚为什么矩阵-矩阵乘法是并行计算的一个很好的例子了。我们必须计算 C C C 中的每个元素,并且每个元素彼此独立,因此我们可以有效地并行化。

我们将看到实现这一目标的不同方法。我们的目标是在本文中添加新概念,最终形成一个 2D 内核,它使用共享内存来有效地优化操作。

网格和块

当我们使用指令 <<< >>> 调用内核时,我们自动定义一个 dim3 类型变量,定义每个网格的块数和每个块的线程数。

事实上,网格和块分别是块和线程的 3D 数组。当我们在调用内核之前定义它们时,这一点很明显,如下所示:

dim3 blocksPerGrid(512, 1, 1)
dim3 threadsPerBlock(512, 1, 1)
kernel<<<blocksPerGrid, threadsPerBlock>>>()

当我们现在处理矩阵时,我们想要指定第二个维度(同样,我们可以省略第三个维度)。这对于使线程正常工作非常有用,有时甚至是必不可少的。

事实上,通过这种方式,我们可以按照前面示例中所遵循的相同方式来引用 x x x y y y​ 轴。 我们看一下代码:

 int row = blockIdx.y * blockDim.y + threadIdx.y;int col = blockIdx.x * blockDim.x + threadIdx.x;

正如您所看到的,它们的代码相似。 在 CUDA 中,blockIdx、blockDim 和 threadIdx 是具有成员 x、y 和 z 的内置函数。 它们在 C 中被索引为法线向量,因此介于 0 和最大数减 1 之间。例如,如果我们的网格维度为blocksPerGrid = (512, 1, 1),则 blockIdx.x 的范围将在 0 到 511 之间。

使用多维块意味着您必须小心地在所有维度之间分配这个数量的线程。 在 1D 块中,x 轴最多可以设置 1024 个线程,但在 2D 块中,如果将 y 的大小设置为 2,则 x 不能超过 512! 例如,允许使用dim3threadsPerBlock(1024, 1, 1),以及dim3threadsPerBlock(512, 2, 1),但不允许使用dim3threadsPerBlock(256, 3, 2)。

内核

现在我们已经有了执行任务所需的所有信息,让我们看一下内核代码。为了简单起见,我们将在示例中使用方阵 N × N N \times N N×N 矩阵。

正如我们已经看到的,要做的第一件事是确定 x x x y y y 轴索引(即行号和列号):

__global__ void multiplication(float *A, float* B, float *C, int N){int ROW = blockIdx.y*blockDim.y+threadIdx.y;int COL = blockIdx.x*blockDim.x+threadIdx.x;

代码

矩阵类

#include <stdexcept>
#include <algorithm>
#include <cuda_runtime.h>template <class T>
class dev_array
{
// public functions
public:explicit dev_array(): start_(0),end_(0){}// constructorexplicit dev_array(size_t size){allocate(size);}// destructor~dev_array(){free();}// resize the vectorvoid resize(size_t size){free();allocate(size);}// get the size of the arraysize_t getSize() const{return end_ - start_;}// get dataconst T* getData() const{return start_;}T* getData(){return start_;}// setvoid set(const T* src, size_t size){size_t min = std::min(size, getSize());cudaError_t result = cudaMemcpy(start_, src, min * sizeof(T), cudaMemcpyHostToDevice);if (result != cudaSuccess){throw std::runtime_error("failed to copy to device memory");}}// getvoid get(T* dest, size_t size){size_t min = std::min(size, getSize());cudaError_t result = cudaMemcpy(dest, start_, min * sizeof(T), cudaMemcpyDeviceToHost);if (result != cudaSuccess){throw std::runtime_error("failed to copy to host memory");}}// private functions
private:// allocate memory on the devicevoid allocate(size_t size){cudaError_t result = cudaMalloc((void**)&start_, size * sizeof(T));if (result != cudaSuccess){start_ = end_ = 0;throw std::runtime_error("failed to allocate device memory");}end_ = start_ + size;}// free memory on the devicevoid free(){if (start_ != 0){cudaFree(start_);start_ = end_ = 0;}}T* start_;T* end_;
};#endif

矩阵乘法处理

#include <stdexcept>
#include <algorithm>
#include <cuda_runtime.h>template <class T>
class dev_array
{
// public functions
public:explicit dev_array(): start_(0),end_(0){}// constructorexplicit dev_array(size_t size){allocate(size);}// destructor~dev_array(){free();}// resize the vectorvoid resize(size_t size){free();allocate(size);}// get the size of the arraysize_t getSize() const{return end_ - start_;}// get dataconst T* getData() const{return start_;}T* getData(){return start_;}// setvoid set(const T* src, size_t size){size_t min = std::min(size, getSize());cudaError_t result = cudaMemcpy(start_, src, min * sizeof(T), cudaMemcpyHostToDevice);if (result != cudaSuccess){throw std::runtime_error("failed to copy to device memory");}}// getvoid get(T* dest, size_t size){size_t min = std::min(size, getSize());cudaError_t result = cudaMemcpy(dest, start_, min * sizeof(T), cudaMemcpyDeviceToHost);if (result != cudaSuccess){throw std::runtime_error("failed to copy to host memory");}}// private functions
private:// allocate memory on the devicevoid allocate(size_t size){cudaError_t result = cudaMalloc((void**)&start_, size * sizeof(T));if (result != cudaSuccess){start_ = end_ = 0;throw std::runtime_error("failed to allocate device memory");}end_ = start_ + size;}// free memory on the devicevoid free(){if (start_ != 0){cudaFree(start_);start_ = end_ = 0;}}T* start_;T* end_;
};#endif

内核头文件

内核CUDA处理

正电子发射断层扫描

临床正电子发射断层扫描的工作原理是,首先向受试者注射含有短寿命放射性同位素(例如 18F)的药物,该同位素会因发射正电子而衰变。 发射的正电子将与受试者组织中附近的电子湮灭,产生一对背靠背的伽马射线,这可以在周围的伽马探测器阵列中作为巧合事件进行检测。

PET 扫描仪中使用的典型探测器由光电倍增管观察的闪烁晶体组成。 晶体吸收伽马射线的能量并将其转换为可见光范围内的光子,进而在光电倍增管中产生信号。 重合事件中的两个探测器定义了一条响应线,衰变事件沿着这条响应线发生。 较旧的 PET 探测器的时间分辨率为几纳秒,足以找到真正湮没事件之间的重合,同时拒绝大多数其他背景来源。 然而,这个时间并不足以提供有关特定衰变事件发生在响应线上的位置的任何有用信息。 临床 PET 扫描可能涉及采集数十亿个 LOR,因此可以通过统计分析来重建受试者体内最可能的活动分布。 重建的活动通常显示为 3D 笛卡尔体素网格,其中 x − y x-y xy 平面位于扫描仪的横向平面中,z 轴与探测器形成的圆柱体的轴对齐。

一种流行的图像重建算法是最大似然期望最大化方法,如果每个响应线中观察到的衰变数量遵循基础放射性衰变过程的泊松分布,则该算法在理论上是最佳的。

最大似然期望最大化方法是迭代的,并使用如下方程式中所示的看似简单的公式:
a v n + 1 = a v n ∑ l ′ = 1 N l S v l ′ ∑ l = 1 N l m l S v l ∑ v ′ = 1 N v a v ′ n S v ′ l a_v^{n+1}=\frac{a_v^n}{\sum_{l^{\prime}=1}^{N_l} S_{v l^{\prime}}} \sum_{l=1}^{N_l} \quad \frac{m_l S_{v l}}{\sum_{v^{\prime}=1}^{N_v} a_{v^{\prime}}^n S_{v^{\prime} l}} avn+1=l=1NlSvlavnl=1Nlv=1NvavnSvlmlSvl
其中 a v n a_v^n avn 是迭代 $n 时体素 v v v 的估计活动,m_l$ 是响应线 l l l 中测量的衰减数, S v l S_{v l} Svl 是 PET 系统矩阵 (SM)。 SM 元素 S v l S_{v l} Svl 定义为在响应线 l l l 中检测到体素 v v v 衰减的概率。 迭代通常是通过为所有 a v 0 a_v^0 av0 分配相等的值来开始的。 总和限值 N l N_l Nl N v N_v Nv 是涉及的响应线和体素的总数。

伊辛模型

该模型在固态物理学中众所周知,涉及在任意维数的等间隔网格上排列的自旋二分之一粒子的演化。 在这里我们将讨论适用于实体的 3D 网格。 在每个网格点,我们放置一个自旋为 S 的对象,其中 S 可以为 ± 1 \pm 1 ±1。 每个自旋都受到其邻居的作用,使得平行自旋的能量低于反平行自旋,并且自旋网格位置 ( x , y , z ) (x, y, z) (x,y,z)​ 的最终能量可以取为:
E = J ( S x − 1 , y , z + S x + 1 , y , z + S x , y − 1 , z + S x , y + 1 , z + S x , y , z − 1 + S x , y , z + 1 ) S x , y , z E=J\left(S_{x-1, y, z}+S_{x+1, y, z}+S_{x, y-1, z}+S_{x, y+1, z}+S_{x, y, z-1}+S_{x, y, z+1}\right) S_{x, y, z} E=J(Sx1,y,z+Sx+1,y,z+Sx,y1,z+Sx,y+1,z+Sx,y,z1+Sx,y,z+1)Sx,y,z
此式是最简单的伊辛模型,我们只考虑最近邻居之间的相互作用。 常数 J J J 测量自旋-自旋相互作用的强度,在我们的模拟中我们将其设置为 1。 请注意,如果大多数相邻自旋与 S x , y , z S_{x, y, z} Sx,y,z 平行,则 E E E 将为正;如果它们大多反平行,则 E E E 将为负。 在模拟中,我们测试单个任意旋转并根据框中显示的标准翻转它。

蒙特卡洛

蒙特卡洛方法的科学应用一直很重要,并且随着计算能力的增长和可用性,其应用的重要性不断增长。 如今,它们已成为许多科学领域的重要工具。 简单来说,这些应用程序可以分为两类: (a) 积分,在其域中的随机点上对某些函数进行采样;(b) 模拟,使用随机数来模拟随机过程来研究某些物理系统或实验设备的行为。 好的随机数生成器是计算机在科学领域的许多应用的重要工具。

CUDA Cpp并行计算编程

参阅一:亚图跨际
参阅二:亚图跨际

这篇关于CUDA Cpp正电子发射断层扫描仪校准和图像重建—蒙特卡洛3D伊辛模型的文章就介绍到这儿,希望我们推荐的文章对编程师们有所帮助!



http://www.chinasem.cn/article/676780

相关文章

Java的IO模型、Netty原理解析

《Java的IO模型、Netty原理解析》Java的I/O是以流的方式进行数据输入输出的,Java的类库涉及很多领域的IO内容:标准的输入输出,文件的操作、网络上的数据传输流、字符串流、对象流等,这篇... 目录1.什么是IO2.同步与异步、阻塞与非阻塞3.三种IO模型BIO(blocking I/O)NI

基于Flask框架添加多个AI模型的API并进行交互

《基于Flask框架添加多个AI模型的API并进行交互》:本文主要介绍如何基于Flask框架开发AI模型API管理系统,允许用户添加、删除不同AI模型的API密钥,感兴趣的可以了解下... 目录1. 概述2. 后端代码说明2.1 依赖库导入2.2 应用初始化2.3 API 存储字典2.4 路由函数2.5 应

使用Python开发一个图像标注与OCR识别工具

《使用Python开发一个图像标注与OCR识别工具》:本文主要介绍一个使用Python开发的工具,允许用户在图像上进行矩形标注,使用OCR对标注区域进行文本识别,并将结果保存为Excel文件,感兴... 目录项目简介1. 图像加载与显示2. 矩形标注3. OCR识别4. 标注的保存与加载5. 裁剪与重置图像

C#集成DeepSeek模型实现AI私有化的流程步骤(本地部署与API调用教程)

《C#集成DeepSeek模型实现AI私有化的流程步骤(本地部署与API调用教程)》本文主要介绍了C#集成DeepSeek模型实现AI私有化的方法,包括搭建基础环境,如安装Ollama和下载DeepS... 目录前言搭建基础环境1、安装 Ollama2、下载 DeepSeek R1 模型客户端 ChatBo

SpringBoot快速接入OpenAI大模型的方法(JDK8)

《SpringBoot快速接入OpenAI大模型的方法(JDK8)》本文介绍了如何使用AI4J快速接入OpenAI大模型,并展示了如何实现流式与非流式的输出,以及对函数调用的使用,AI4J支持JDK8... 目录使用AI4J快速接入OpenAI大模型介绍AI4J-github快速使用创建SpringBoot

0基础租个硬件玩deepseek,蓝耘元生代智算云|本地部署DeepSeek R1模型的操作流程

《0基础租个硬件玩deepseek,蓝耘元生代智算云|本地部署DeepSeekR1模型的操作流程》DeepSeekR1模型凭借其强大的自然语言处理能力,在未来具有广阔的应用前景,有望在多个领域发... 目录0基础租个硬件玩deepseek,蓝耘元生代智算云|本地部署DeepSeek R1模型,3步搞定一个应

Deepseek R1模型本地化部署+API接口调用详细教程(释放AI生产力)

《DeepseekR1模型本地化部署+API接口调用详细教程(释放AI生产力)》本文介绍了本地部署DeepSeekR1模型和通过API调用将其集成到VSCode中的过程,作者详细步骤展示了如何下载和... 目录前言一、deepseek R1模型与chatGPT o1系列模型对比二、本地部署步骤1.安装oll

Spring AI Alibaba接入大模型时的依赖问题小结

《SpringAIAlibaba接入大模型时的依赖问题小结》文章介绍了如何在pom.xml文件中配置SpringAIAlibaba依赖,并提供了一个示例pom.xml文件,同时,建议将Maven仓... 目录(一)pom.XML文件:(二)application.yml配置文件(一)pom.xml文件:首

如何在本地部署 DeepSeek Janus Pro 文生图大模型

《如何在本地部署DeepSeekJanusPro文生图大模型》DeepSeekJanusPro模型在本地成功部署,支持图片理解和文生图功能,通过Gradio界面进行交互,展示了其强大的多模态处... 目录什么是 Janus Pro1. 安装 conda2. 创建 python 虚拟环境3. 克隆 janus

本地私有化部署DeepSeek模型的详细教程

《本地私有化部署DeepSeek模型的详细教程》DeepSeek模型是一种强大的语言模型,本地私有化部署可以让用户在自己的环境中安全、高效地使用该模型,避免数据传输到外部带来的安全风险,同时也能根据自... 目录一、引言二、环境准备(一)硬件要求(二)软件要求(三)创建虚拟环境三、安装依赖库四、获取 Dee