KD-tree比赛题目CUDA

Nvidia CUDA全国程序设计大赛命题2：光线和三角形求交算法的并行化在光线跟踪或光线投射算法中，涉及到大量的光线和三角形面片求交计算，具体来说，对于光线投射或一遍光线跟踪来说，我们需要从视点出发向屏幕每个像素引出一条光线，然后计算出这条光线跟物体的交点，然后把交点的颜色作为屏幕上对应像素的颜色值，这样得到整个三维场景的渲染效果图。

显然，每根光线跟三角形面片的计算是相互独立的，而且某一根光线跟每个三角形之间也是相互独立的。

这里存在着大量平行的可能性。

问题：利用CUDA计算出每一条光线跟三角形的交点，注意有可能其中一些光线跟所有的三角形都不相交，如果某根光线跟多个三角形相交，则只返回离视点最近的那个交点的坐标。

输入：视点位置，屏幕坐标点，三角形面片集合其中三角形集合的格式类似标准的Wave front的obj格式，这里用更简化的形式来描述，即只包含顶点序列和三角形序列。

文件格式如下：顶点数x,y,z（即顶点1三维坐标值）x,y,z（即顶点2三维坐标值）…三角形数v1 v2 v3（即三角形1顶点序号）v1 v2 v3（即三角形2顶点序号）…一些提示：1.暴力算法不是好算法；2.光线由视点与屏幕坐标点连线构成；3.场景可以利用二叉树划分，包围盒等技术来加速；4.某些光线可以按簇来处理参考文献：1.Timothy J. Purcell, Ian Buck, William R. Mark, and Pat Hanrahan, Ray Tracing onProgrammable Graphics Hardware,ACM Transactions on Graphics. 21 (3), pp. 703-712, 2002.2.Min Shih, Yung-Feng Chiu, Ying-Chieh Chen and Chun-Fa Chang, Real-Time Ray Tracingwith CUDA,Lecture Notes in Computer Science, V olume 5574/2009.4./object/cuda_home_new.html5./object/nvision08-IRT.html注：下面是计算机图形学中整个光线跟踪的伪代码（参考彭群生教授编写的《计算机图形学应用基础》），但本竞赛题只需求出一遍光线与三角形面片的交点，即黑体部分。

KdTree算法详解kd树（k-dimensional树的简称），是⼀种分割k维数据空间的数据结构，主要应⽤于多维空间关键数据的近邻查找(Nearest Neighbor)和近似最近邻查找(Approximate Nearest Neighbor)。

⼀、Kd-tree其实KDTree就是⼆叉查找树（Binary Search Tree，BST）的变种。

⼆叉查找树的性质如下：1）若它的左⼦树不为空，则左⼦树上所有结点的值均⼩于它的根结点的值；2）若它的右⼦树不为空，则右⼦树上所有结点的值均⼤于它的根结点的值；3）它的左、右⼦树也分别为⼆叉排序树；例如：如果我们要处理的对象集合是⼀个K维空间中的数据集，我们⾸先需要确定是：怎样将⼀个K维数据划分到左⼦树或右⼦树？在构造1维BST树类似，只不过对于Kd树，在当前节点的⽐较并不是通过对K维数据进⾏整体的⽐较，⽽是选择某⼀个维度d，然后⽐较两个K维数据在该维度 d上的⼤⼩关系，即每次选择⼀个维度d来对K维数据进⾏划分，相当于⽤⼀个垂直于该维度d的超平⾯将K维数据空间⼀分为⼆，平⾯⼀边的所有K维数据在d维度上的值⼩于平⾯另⼀边的所有K维数据对应维度上的值。

也就是说，我们每选择⼀个维度进⾏如上的划分，就会将K维数据空间划分为两个部分，如果我们继续分别对这两个⼦K维空间进⾏如上的划分，⼜会得到新的⼦空间，对新的⼦空间⼜继续划分，重复以上过程直到每个⼦空间都不能再划分为⽌。

以上就是构造 Kd-Tree的过程，上述过程中涉及到两个重要的问题：1. 每次对⼦空间的划分时，怎样确定在哪个维度上进⾏划分；2. 在某个维度上进⾏划分时，怎样确保建⽴的树尽量地平衡，树越平衡代表着分割得越平均，搜索的时间也就是越少。

1、在哪个维度上进⾏划分？⼀种选取轴点的策略是median of the most spread dimension pivoting strategy，统计样本在每个维度上的数据⽅差，挑选出对应⽅差最⼤值的那个维度。

cuda实践5---cudakdtree、octreecuda kdtree前⾔：将kdtree 查询部分移植到GPU端，在很多应⽤中对提⾼算法的执⾏效率很有帮助，本⽂使⽤英伟达GPU语⾔cuda，完成了kdtree GPU端的移植。

步骤⽐较简单：1、cpu端创建kdtree； 2、迁移kdtree node 节点到GPU端；3、GPU端实现临近检索（注：⾥⾯会有很多处理⼩技巧，望相互学习）核⼼代码：cuda_kdtree.cu//main.cu#include <stdio.h>#include <iostream>#include "cuda_runtime.h"#include "device_launch_parameters.h"#include <math.h>#include "cuda_kdtree.h"void CheckCUDAError(const char *msg){cudaError_t err = cudaGetLastError();if (cudaSuccess != err) {fprintf(stderr, "Cuda error: %s: %s.\n", msg, cudaGetErrorString(err));exit(EXIT_FAILURE);}}/////////////////////////////////////////////////////////////////////////////////////////////////////////////////////////////////////////////////////////__device__ float distance3d(float* pt1, float* pt2){float temp = sqrtf((pt1[0] - pt2[0])*(pt1[0] - pt2[0]) +(pt1[1] - pt2[1])*(pt1[1] - pt2[1]) +(pt1[2] - pt2[2])*(pt1[2] - pt2[2]));return temp;}/////////////////////////////////////////////////////////////////////////////////////////////////////////////////////////////////////////////////////////__device__ int findLeafNode_GPU(float* searchPoint,int current_node ,int node_num, CUDA_KDNode* kdnode_vector) {if (node_num == 0){//std::cout << "CUDA_KDNode is empty ..." << std::endl;return -1;}while (true){if (kdnode_vector[current_node].left == -1 && kdnode_vector[current_node].right == -1) //leaf node{break;}int slipnum0 = 1;slipnum0 = kdnode_vector[current_node].dim;switch (slipnum0){case1:if (searchPoint[0] < kdnode_vector[current_node].split_value){current_node = kdnode_vector[current_node].left;}else{current_node = kdnode_vector[current_node].right;}break;case2:if (searchPoint[1] < kdnode_vector[current_node].split_value){current_node = kdnode_vector[current_node].left;}else{current_node = kdnode_vector[current_node].right;}break;case3:if (searchPoint[2] < kdnode_vector[current_node].split_value){current_node = kdnode_vector[current_node].left;}else{current_node = kdnode_vector[current_node].right;}break;}}return current_node;}/////////////////////////////////////////////////////////////////////////////////////////////////////////////////////////////////////////////////////////__device__ int findNode_GPU(float* searchPoint, float searchRadius, int current_node, int node_num, CUDA_KDNode* kdnode_vector) { if (node_num == 0){//std::cout << "CUDA_KDNode is empty ..." << std::endl;return -1;}while (true){if (kdnode_vector[current_node].left == -1 || kdnode_vector[current_node].right == -1) //leaf node{break;}int slipnum0 = 1;slipnum0 = kdnode_vector[current_node].dim;float dis_temp = fabs(searchPoint[slipnum0 - 1] - kdnode_vector[current_node].split_value);if (dis_temp < searchRadius){break;}if (searchPoint[slipnum0-1] < kdnode_vector[current_node].split_value){current_node = kdnode_vector[current_node].left;}else{current_node = kdnode_vector[current_node].right;}}return current_node;}/////////////////////////////////////////////////////////////////////////////////////////////////////////////////////////////////////////////////////////__device__ int searchRadius_GPU(float* pointCloud,int* pt_indexs, int node_num, CUDA_KDNode* kdnode_vector,float* searchPoint,float searchRadius,int* searchIndex, float* searchDistance){if (node_num == NULL){//std::cout << "CUDA_KDNode is empty ..." << std::endl;return -1;}int current_node = 0; //根节点int templeafnode = 0;// findNode_GPU(searchPoint, searchRadius, current_node, node_num, kdnode_vector);float distanceTemp = 0.0;int temp = 0;for (int j = 0; j < kdnode_vector[templeafnode].numOfData; j++){float pt2[3];pt2[0] = pointCloud[3 * pt_indexs[kdnode_vector[templeafnode].start_index + j]];pt2[1] = pointCloud[3 * pt_indexs[kdnode_vector[templeafnode].start_index + j] + 1];pt2[2] = pointCloud[3 * pt_indexs[kdnode_vector[templeafnode].start_index + j] + 2];distanceTemp = distance3d(searchPoint, pt2);if (searchRadius > distanceTemp && temp < maxN){searchIndex[temp] = pt_indexs[kdnode_vector[templeafnode].start_index + j];searchDistance[temp] = distanceTemp;temp++;}}while (true){int indexTemp = -1;float tempdis;int changetimes = 0;if (temp == 0){break;}for (int i = 0; i < temp-1; i++){if (searchDistance[i] > searchDistance[i + 1]){tempdis = searchDistance[i];searchDistance[i] = searchDistance[i + 1];searchDistance[i + 1] = tempdis;indexTemp = searchIndex[i];searchIndex[i] = searchIndex[i + 1];searchIndex[i + 1] = indexTemp;changetimes++;}}if (changetimes == 0){break;}}return temp;}/////////////////////////////////////////////////////////////////////////////////////////////////////////////////////////////////////////////////////////template <int BLOCK_SIZE> __global__ void runKNNSearchRadius_GPU(float* pointCloud, int* pt_indexs, int node_num, CUDA_KDNode* kdnode_vector,float searchRadius, int* searchIndex_d, float* searchDistance_d){int searchIndex[maxN];float searchDistance[maxN];for (int i = 0; i < maxN; i++){searchIndex[i] = -999;searchDistance[i] = 999.;}int Row = blockIdx.y * BLOCK_SIZE + threadIdx.y;int Col = blockIdx.x * BLOCK_SIZE + threadIdx.x;float searchPoint[2];searchPoint[0] = Col * 1.0 / imgWidth_d;searchPoint[1] = (imgHeight_d - Row) * 1.0 / imgHeight_d;int searchNum = searchRadius_GPU(pointCloud, pt_indexs, node_num, kdnode_vector,searchPoint, searchRadius,searchIndex, searchDistance);int threadId = Row * imgWidth_d + Col;for (int i = 0; i < maxN; i++){searchIndex_d[maxN*threadId + i] = searchIndex[i];searchDistance_d[maxN*threadId + i] = searchDistance[i];}}/////////////////////////////////////////////////////////////////////////////////////////////////////////////////////////////////////////////////////////CUDA_KDTREE::~CUDA_KDTREE(){cudaFree(m_gpu_nodes);cudaFree(m_gpu_indexes);cudaFree(m_gpu_points);}/////////////////////////////////////////////////////////////////////////////////////////////////////////////////////////////////////////////////////////void CUDA_KDTREE::createKDTree(std::vector<Point3d>& pointCloud, vector<CUDA_KDNode>& cpu_nodes, vector<int>& indexes){kd::KDTree<Point3d> tree;tree.setInputPointCloud(pointCloud);tree.setNumOfLeafData(100);tree.buildKDTree();std::queue<kd::kdnode*> kdnodePtrQue;kdnodePtrQue.push(tree.ptree->root);std::vector<kd::kdnode*> kdnode_vect;//go through kdtree nodeint node_t = 0;while (!kdnodePtrQue.empty()){kd::kdnode* tempkdnode = kdnodePtrQue.front();if (tempkdnode==NULL){kdnodePtrQue.pop();continue;}tempkdnode->node_index = node_t;kdnode_vect.push_back(tempkdnode);CUDA_KDNode gpu_node_temp;gpu_node_temp.dim = tempkdnode->dim;gpu_node_temp.nodelr = tempkdnode->nodelr;gpu_node_temp.numOfData = tempkdnode->numOfData;gpu_node_temp.split_value = tempkdnode->split_value;cpu_nodes.push_back(gpu_node_temp);kdnodePtrQue.pop();if (tempkdnode->left!=NULL){kdnodePtrQue.push(tempkdnode->left);}if (tempkdnode->right != NULL){kdnodePtrQue.push(tempkdnode->right);}node_t++;}int index_t = 0;for (int i = 0; i < node_t; i++){if (kdnode_vect[i]->father==NULL)cpu_nodes[i].father = -1;elsecpu_nodes[i].father = kdnode_vect[i]->father->node_index;if (kdnode_vect[i]->left == NULL)cpu_nodes[i].left = -1;elsecpu_nodes[i].left = kdnode_vect[i]->left->node_index;if (kdnode_vect[i]->right == NULL)cpu_nodes[i].right = -1;elsecpu_nodes[i].right = kdnode_vect[i]->right->node_index;cpu_nodes[i].start_index = index_t;index_t += cpu_nodes[i].numOfData;for (int j = 0; j < cpu_nodes[i].numOfData; j++){indexes.push_back(kdnode_vect[i]->data[j]);}}}/////////////////////////////////////////////////////////////////////////////////////////////////////////////////////////////////////////////////////////int CUDA_KDTREE::search(std::vector<Point3d>& pointCloud, vector<CUDA_KDNode>& cpu_nodes, vector<int>& indexes, int* queries_indexes, float *queries_dists) {// Create the nodes again on the CPU, laid out nicely for the GPU transferint m_num_index = indexes.size();int num_nodes = cpu_nodes.size();int point_num = pointCloud.size();//int* queries_indexes = NULL;//float *queries_dists = NULL;//queries_indexes = (int*)malloc(sizeof(int) * imgHeight_d * imgWidth_d* maxN);//queries_dists = (float*)malloc(sizeof(int) * imgHeight_d * imgWidth_d* maxN);int* searchIndex_d = NULL;float* searchDistance_d = NULL;cudaMalloc((void**)&m_gpu_nodes, sizeof(CUDA_KDNode)*num_nodes);cudaMalloc((void**)&m_gpu_indexes, sizeof(int)*m_num_index);cudaMalloc((void**)&m_gpu_points, sizeof(float) * 3 * point_num);cudaMalloc((void**)&searchIndex_d, sizeof(int) * imgHeight_d * imgWidth_d* maxN);cudaMalloc((void**)&searchDistance_d, sizeof(float) * imgHeight_d * imgWidth_d* maxN);cudaMemcpy(m_gpu_nodes, &cpu_nodes[0], sizeof(CUDA_KDNode) * num_nodes, cudaMemcpyHostToDevice);cudaMemcpy(m_gpu_indexes, &indexes[0], sizeof(int) * m_num_index, cudaMemcpyHostToDevice);cudaMemcpy(m_gpu_points, &pointCloud[0], sizeof(float) * 3 * point_num, cudaMemcpyHostToDevice);int BLOCK_SIZE = 16;// Setup execution parametersdim3 block(BLOCK_SIZE, BLOCK_SIZE); //16*16 或者 32*32 block⼤⼩dim3 grid(imgWidth_d / block.x, imgHeight_d / block.y); //计算grid ⼤⼩runKNNSearchRadius_GPU<16><<<grid, block>>> (m_gpu_points, m_gpu_indexes, num_nodes, m_gpu_nodes, 0.02, searchIndex_d, searchDistance_d); cudaError_t cudaStatus;cudaStatus = cudaGetLastError();// cudaDeviceSynchronize waits for the kernel to finish, and returns// any errors encountered during the launch.cudaStatus = cudaDeviceSynchronize();cudaMemcpy(&queries_indexes[0], searchIndex_d, sizeof(int)*imgHeight_d * imgWidth_d* maxN, cudaMemcpyDeviceToHost);cudaMemcpy(&queries_dists[0], searchDistance_d, sizeof(float)*imgHeight_d * imgWidth_d* maxN, cudaMemcpyDeviceToHost);CheckCUDAError("CreateKDTree");// freecudaFree(searchIndex_d);cudaFree(searchDistance_d);return0;}参考：。

算法比赛考试题及答案一、单项选择题（每题5分，共20分）1. 在算法中，用于描述算法性能的指标是：A. 时间复杂度B. 空间复杂度C. 算法正确性D. 算法稳定性答案：A2. 快速排序算法的平均时间复杂度为：A. O(n)B. O(n^2)C. O(nlogn)D. O(logn)答案：C3. 以下哪个数据结构不适合用链表实现：A. 队列B. 栈C. 哈希表D. 双向队列答案：C4. 动态规划和贪心算法的主要区别在于：A. 是否使用递归B. 是否有最优子结构性质C. 是否有重叠子问题D. 是否需要回溯答案：B二、填空题（每题5分，共20分）1. 算法的时间复杂度为O(n^2)的排序算法是______。

答案：冒泡排序2. 在图的遍历中，深度优先搜索使用的栈是______。

答案：非必需的3. 哈希表解决冲突的方法之一是______。

答案：链地址法4. 二叉树的前序遍历顺序是______。

答案：根左右三、简答题（每题10分，共30分）1. 请简述分治法的基本思想。

答案：分治法的基本思想是将一个难以直接解决的大问题分解成若干个规模较小的相同问题，递归地解决这些子问题，然后将子问题的解合并得到原问题的解。

2. 什么是递归？请给出一个递归算法的例子。

答案：递归是一种算法设计技巧，它将问题分解成更小的、与原问题相似的子问题，然后递归地解决这些子问题。

例如，计算n的阶乘可以通过递归实现：n! = n * (n-1)!，其中递归的基本情况是0! = 1。

3. 请解释什么是动态规划，并给出一个动态规划解决问题的例子。

答案：动态规划是一种通过将复杂问题分解成更简单的子问题来求解的方法，它通常用于求解具有重叠子问题和最优子结构性质的问题。

例如，斐波那契数列的计算可以通过动态规划实现，避免重复计算相同的子问题。

四、编程题（共30分）1. 给定一个整数数组，请编写一个函数，找出数组中第二大的数。

答案：以下是使用Python语言实现的代码示例：```pythondef find_second_max(nums):if len(nums) < 2:return Nonefirst_max = second_max = float('-inf')for num in nums:if num > first_max:second_max = first_maxfirst_max = numelif num > second_max and num != first_max:second_max = numreturn second_max if second_max != float('-inf') else None```请注意，以上代码仅为示例，实际编程题的解答可能需要根据具体的编程语言和题目要求进行调整。