一种GPU加速的参数曲线弧长计算及弧长参数化方法与流程

本发明涉及计算机辅助设计领域，尤其是涉及一种gpu加速的参数曲线弧长计算及弧长参数化方法。

背景技术：

参数曲线的弧长计算和弧长参数化广泛地应用在运动控制、路径规划等问题之中。

参数曲线c(u)的弧长函数定义为从最小参数值u0起到参数t为止的，曲线逐点切线即一阶导矢模长之积分，即弧长公式其反方向的任务，给出弧长a，由弧长函数的逆函数求出对应的参数t，即t＝a^-1(a)，称之为弧长参数化。

由于参数曲线的形式较复杂，其弧长公式往往无法化简为只用初等函数表达的解析形式，因此无论求解曲线弧长还是弧长参数化都难以直接计算。

传统的cpu方法[1]，使用递归算法求解弧长，若其不满足用户指定的精度条件，则将曲线以二分法进行细化，然后求解每一个细化后区间的弧长，直到满足精度条件为止。在递归计算的过程中，按次序记录下每个满足精度条件的区间的弧长及端点，称作弧长表。最终的弧长为所有满足精度条件的弧长相加得到。在求解弧长参数化问题时，首先在弧长表中查找弧长对应的参数区间，然后在区间中递归地查找弧长所对应的参数，即二分区间、保留其中弧长所在的区间、在该区间中继续搜索直到满足用户指定的误差条件为止。

在弧长计算和弧长参数化中，将区间进行二分并计算其弧长本质上是二叉树的遍历任务。在每个区间中计算弧长时，有多种求积分近似解的方法，例如高斯-勒让德求积法(gauss-legendrequadrature)([1]guenterb,parentr.computingthearclengthofparametriccurves[j].ieeecomputergraphicsandapplications,1990,10(3):0-78)、折线逼近法(polylinemethod)、辛普森法(simpsonmethod)和文森特-福西法(vincent-forseyrule)([2]vincents,forseyd.fastandaccurateparametriccurvelengthcomputation[j].2001；[3]floaterms,rasmussenaf.point-basedmethodsforestimatingthelengthofaparametriccurve[j].journalofcomputationalandappliedmathematics,2006,196(2):512-522)。

传统的cpu方法采用递归方法计算弧长、求解弧长参数化，这种方法由于计算量高，因此在需要较长的计算时间，并且其方法相比于gpu的单指令多线程(simt)架构，具有函数调用开销大、指令分歧率高、占用局部存储存储空间多等多种缺点。并且方法中的弧长表动态构建，需要在运行时对其存储空间进行动态分配和释放，在gpu中开销较大，因此并不适合直接移植到gpu中。

技术实现要素：

本发明的目的在于针对现有技术存在的上述不足，提供数据的存储方式与gpu的数据存取方式友好，指令分歧率小，减少时间消耗的一种gpu加速的参数曲线弧长计算及弧长参数化方法。

本发明包括以下步骤：

1)输入参数曲线的数量、每条参数曲线区间段的个数以及曲线参数；

2)将每一个区间段等分为n份首尾相接的参数“小区间段”；

3)构建“小区间弧长表”和“区间弧长表”；

4)每条曲线对应的区间弧长表中的最后一项为弧长，将其输出；

5)输入待求解弧长参数化的弧长值及其所在的曲线编号i；

6)确定该弧长所在的“搜索区间”，初始化搜索任务，每个弧长的搜索区间的初始化由1个gpu线程完成；

7)使用“基于二进制位操作加速二叉树遍历回溯的深度优先搜索”，将搜索区间的长度缩小到不大于用户指定的误差范围ξ，每个弧长的弧长参数化由1个gpu线程完成；

8)将最终搜索区间中点作为弧长所对应的参数输出。

在步骤1)中，所述输入参数曲线应满足如下要求：每条曲线与其他曲线彼此独立；每条曲线包括一个或若干个由用户指定的参数“区间段”，每个区间段上的曲线上点及切线的计算与其他区间段彼此独立，在同一区间段内曲线上不同参数的点及切线的计算彼此独立。

在步骤2)中，所述小区间段的划分方法可为：n为16、32、64、128或256，所有的区间段的划分时所取的n均相同；划分时，令每一个小区间段均保存该小区间段对应的参数区间端点值[ucij,uci(j+1)]、曲线编号i、区间段编号j以及曲线在该区间段上的参数值p(cij)，使得每一个小区间段均可被独立计算。

在步骤3)中，所述包括构建“小区间弧长表”和“区间弧长表”的具体方法可为：

步骤31.采用“双阶段的广度优先搜索”计算满足用户指定误差ε的曲线在该小区间段上的弧长值；小区间段弧长的计算结果按顺序存储，称为“小区间弧长表”，其排序优先级为，曲线编号>区间段编号>小区间段首端点值；

步骤32.计算小区间弧长表中每个区间段的弧长的前缀和，存储，称之为小区间弧长表；

步骤33.计算每条曲线各个区间段所对应的小区间弧长表中每个区间的最后一项的前缀和，按顺序存储，称为“区间弧长表”，排序优先级为，曲线编号>区间段编号。

在步骤31中，所述“双阶段的广度优先搜索”在搜索时所使用的节点包括待求弧长的区间的首尾端点值α和β；所在曲线编号i；所在区间段编号j；区间段所对应曲线的参数信息p(cij)；上一层搜索中，计算该节点的父节点的线程的编号pid；节点的父节点的弧长pal；

所述“双阶段的广度优先搜索”的具体步骤可为：

步骤311.从小区间段所包含的参数中构建出搜索树的根节点，其中首尾端点值、曲线编号、区间段编号、参数信息均与小区间段中的数据相同，pid无须初始化，pal为本区间的弧长值，即所有搜索树的根节点按照其对应的小区间段的顺序存储在第0层的任务队列fq0中，初始化层数level＝0；

步骤312.线程从当前任务队列fqlevel中读取线程编号tid所对应的搜索节点，计算其区间中点计算两个子区间[α,γ]和区间[γ,β]的弧长alleft和alright；

步骤313.对比|alleft+alright-pal|与用户指定误差ε的大小，若小于ε，则接受alleft+alright作为该节点的弧长，并将本节点的pid和弧长alleft+alright存储到反向队列bqlevel的tid位置；若不小于ε，则创建左、右两个子节点，两个子节点的i、j、p(cij)和本节点相同，它们的pid赋值为本节点的线程编号tid，左、右节点的pal分别赋值为alleft和alright，左、右节点的(α,β)分别赋值为(α,γ)和(γ,β)，并将这两个子节点连续存储到下一层的任务队列fqlevel+1中；

步骤314.等待当前任务队列中所有节点均被处理且子节点都已经生成并存储到下一层的任务队列之后，若fqlevel为空，则转入步骤315；否则level＝level+1，转入步骤312；

步骤315.对于反向队列bqlevel中的每一项，读取出其中保存的弧长和pid，将其中保存的弧长加回到bqlevel-1的pid位置；

步骤316.当bqlevel-1每一项都被读取并计算完成后，修改level＝level-1；若level＞0，转入步骤315；否则，转入步骤317；

步骤317.将bq0中所保存的每个小区间段的弧长用于步骤32的输入。

在步骤6中，所述确定该弧长所在的“搜索区间”的具体步骤可为：

(1)对于每一个弧长值，通过二分查找在其所在曲线对应的区间弧长表中查找到其所对应的参数段区间编号j；

(2)通过二分查找在其参数段所对应的小区间弧长表中搜索到该弧长对应的小区间段的首尾端点[ucij,uci(j+1)]；

(3)根据弧长对应的曲线编号i和查找到的区间段的编号j，加载步骤2中所述区间段所应具备曲线参数信息p(cij)，根据查找到的小区间段的首尾端点[ucij,uci(j+1)]确定该弧长对应参数所在的区间，称为“搜索区间”。

在步骤7中，所述使用“基于二进制位操作加速二叉树遍历回溯的深度优先搜索”，包括以下步骤：

步骤71.对于给定待搜索弧长al和搜索区间[α,β]及其对应的区间弧长表和小区间弧长表的位置，求出弧长al在这个小区间中从小区间首端点起的相对弧长，初始化节点的二进制编号为bid＝01b，初始化pal为曲线在区间[α,β]的弧长，初始化递归深度r＝0，初始化累积弧长值alacc为0；

步骤72.求搜索区间的中点计算区间、[α,γ]和[γ,β]对应的弧长值,alleft和alright；若不满足精度条件|alleft+alright-pal|＜ε，则转入步骤73；否则转入步骤74；

步骤73.保存alright，修改pal＝alleft，修改(α,β)＝(α,γ)，修改bid为bid＝bid＜＜1，即将bid向左移一位，这相当于进入了二叉树节点的左子树进行搜索，转入步骤72；

步骤74.针对al与alacc+alleft、alacc+alleft+alright的大小关系，会出现三种情况：

①若al≤alacc+alleft，则说明弧长al对应参数值在左半区间，若满足精度条件|γ-α|＜ξ，则转入步骤75,否则，转入步骤73；

②若alacc+alleft＜al＜alacc+alleft+alright，这说明弧长al对应的参数在右半区间，若满足精度条件|β-γ|＜ξ，则转入步骤75；否则修改α＝γ，pal为当前节点右兄弟的弧长值，即当前节点的弧长作为alleft时对应的alright，修改bid＝(bid＜＜1)|01b，即将bid左移一位，且将最右端新增加的一位修改为1，修改alacc＝alacc+alleft这相当于进入二叉树节点的右子树进行搜索，转入步骤72；

③若al＞alacc+alleft+alright，这说明弧长al对应的参数值在区间外的右侧，此节点到根节点路径上(含自身)必存在一个最低的、是其父节点左节点的节点，这个节点的右兄弟节点尚未被搜索，它即是应回溯到的节点，此时需要对二叉树的搜索进行回溯；上述是其父节点左孩子的祖先节点，对应着二进制编号bid中最右侧的0，这个0右边跟随的1是连续的，其数量n1等于当前节点与应回溯到的节点的高度之差；这些连续的1，转化为对应的十进制数后再加一，即正是两个节点区间长度倍数之比；当前节点的尾端点是应回溯到的节点的右端点，利用这一特征将当前节点的端点α和β快速回溯到应回溯到的节点；修改(α,β)＝(β,β+(β-α)×2ⁿ¹)，修改bid＝(bid＞＞n1)|01b，修改pal为此前保存的计算回溯后节点的父节点时对应的alright，修改alacc＝alacc+alleft+alright，完成回溯，转入步骤72；

步骤75.结束树的搜索，保存最终的搜索区间首尾端点。

在步骤311、步骤312、步骤71、步骤72中，所述计算一个区间对应的曲线弧长时，即求解积分时，可采用的方法包括并不限于高斯勒让德求积法(要求曲线可微)、折线逼近法、文森特法或文森特-福西法，其他可用于求解曲线弧长的方法均可被采用，但是在上述4个步骤中，需采用同样的求解方法。

本发明重建设计了弧长表的构建过程，将任务进行了分解，增大了任务的可并行性，避免了弧长表的动态构建。步骤1)-步骤4)可实现参数曲线弧长计算，执行步骤1)-步骤8)可实现参数曲线弧长参数化。

本发明提出一种新的适用于任务中弧长计算时的二叉树遍历算法，称作“双阶段的广度优先搜索”，它可用于处理搜索时树中更低层的计算结果需要汇聚到更高层的情况，并且数据的存储方式与gpu的数据存取方式友好，指令分歧率小，减少了因线程块中少数线程递归深度高致使其他线程空等的时间浪费，它相比于朴素递归算法没有精度损失，并且计算时间大幅减少。

本发明中提出了一种适用于任务中弧长参数化时的二叉树遍历算法，称作“基于二进制位操作加速二叉树遍历回溯的深度优先搜索”。它将树中节点的左右孩子进行二进制编码，利用先序搜索中的节点访问次序、区间长度与树中层次的倍数关系，实现了弧长参数化时所搜索的二叉树的快速回溯，相比朴素方法具有较小精度误差，使用更少的gpu资源，特别是寄存器和局部存储空间，大幅加速二叉树的遍历速度。

与现有技术相比，本发明具有如下优点：

1.本发明提出的算法框架在构建弧长表时，就采用了并行化的方法，将区间的弧长计算分解为了小区间的弧长计算，减少了递归的深度，同时避免了弧长表长度的动态变化，更适合gpu的计算模型，计算速度大幅提升；

2.本发明提出的“双阶段的广度优先搜索”将递归搜索任务分层执行，每次只执行一层，避免了由于不同线程具有不同递归深度造成的忙等，降低了线程块内线程的指令分歧，大幅减少了弧长计算阶段中gpu中核函数的执行时间；

3.本发明提出的“基于二进制位操作加速二叉树遍历回溯的深度优先搜索”将二进制的01表达与二叉树的左右孩子相结合，将二进制表达对应的十进制数与区间二分搜索过程中区间长度的倍数关系相结合，以更少的寄存器、更少局部存储空间的消耗，实现了更快速的回溯。由于资源使用的降低，gpu执行算法时的占用率显著提升，弧长参数化任务的计算速度大幅加快。

附图说明

图1为本发明实施例的流程图。

图2为本发明实施例中的输入参数曲线中的一条。其中菱形为曲线的控制顶点，虚线段为控制顶点的连线，实线段为曲线。

图3为图2中的参数曲线，通过本发明的均匀弧长参数化采点得到的散点图。

图4为实施例中本发明提出的方法与cpu方法的计算耗时对比图。其中，黑色柱为本发明耗时，白色柱为cpu方法耗时，1为弧长计算任务耗时，2为完成弧长表构建后的弧长参数化耗时，3为弧长计算及弧长参数化总的耗时。

具体实施方式

以下实施例将结合附图对要发明作进一步的说明。应当理解，优选实施例仅为了说明本发明，而不是为了限制本发明的保护范围。

本发明提出的一种gpu加速的参数曲线弧长计算及弧长参数化方法其特征在于，具有以下主要步骤，其中执行步骤1-步骤4可实现参数曲线弧长计算，执行步骤1-步骤8可实现参数曲线弧长参数化：

步骤1.输入参数曲线的数量、每条参数曲线区间段的个数、以及曲线参数；

步骤2.将每一个区间段等分为n份首尾相接的参数“小区间段”；

步骤3.构建“小区间弧长表”和“区间弧长表”；

步骤4.每条曲线对应的区间弧长表中的最后一项为弧长，将其输出；

步骤5.输入待求解弧长参数化的弧长值及其所在的曲线编号i；

步骤6.确定该弧长所在的“搜索区间”，初始化搜索任务，每个弧长的搜索区间的初始化由1个gpu线程完成；

步骤7.使用本发明提出的“基于二进制位操作加速二叉树遍历回溯的深度优先搜索”，将搜索区间的长度缩小到不大于用户指定的误差范围ξ，每个弧长的弧长参数化由1个gpu线程完成；

步骤8.将最终搜索区间中点作为弧长所对应的参数输出。

步骤2中的小区间段划分具有以下特征：步骤2中的n为16、32、64、128或256，所有的区间段的划分时所取的n均相同。划分时，令每一个小区间段均保存该小区间段对应的参数区间端点值[ucij,uci(j+1)]、曲线编号i、区间段编号j以及曲线在该区间段上的参数值p(cij)，使得每一个小区间段均可被独立计算。

步骤3包括以下步骤：

步骤31.采用方法为本发明所提出的“双阶段的广度优先搜索”计算满足用户指定误差ε的曲线在该小区间段上的弧长值。小区间段弧长的计算结果按顺序存储，称为“小区间弧长表”，其排序优先级为，曲线编号>区间段编号>小区间段首端点值；

步骤32.计算小区间弧长表中每个区间段的弧长的前缀和，存储，称之为小区间弧长表；

步骤33.计算每条曲线各个区间段所对应的小区间弧长表中每个区间的最后一项的前缀和，按顺序存储，称为“区间弧长表”，排序优先级为，曲线编号>区间段编号。

步骤31中所使用的“双阶段的广度优先搜索”其搜索时所使用的节点具有如下属性：

1)待求弧长的区间的首尾端点值α和β；

2)所在曲线编号i；

3)所在区间段编号j；

4)区间段所对应曲线的参数信息p(cij)；

5)上一层搜索中，计算该节点的父节点的线程的编号pid；

6)节点的父节点的弧长pal。

步骤31中所使用的“双阶段的广度优先搜索”算法步骤如下，其中步骤311-步骤314称作“前向阶段”，步骤315-316称作“后向阶段”：

步骤311.从小区间段所包含的参数中构建出搜索树的根节点，其中首尾端点值、曲线编号、区间段编号、参数信息均与小区间段中的数据相同，所有搜索树的根节点按照其对应的小区间段的顺序存储在第0层的任务队列fq0中，初始化层数level＝0。

步骤312.线程从当前任务队列fqlevel中读取线程编号tid所对应的搜索节点，计算其区间中点计算两个子区间[α,γ]和区间[γ,β]的弧长alleft和alright。

步骤313.对比|alleft+alright-pal|与用户指定误差ε的大小，若小于ε，则接受alleft+alright作为该节点的弧长，并将本节点的pid和弧长alleft+alright存储到反向队列bqlevel的tid的位置；若不小于ε，则创建左、右两个子节点，两个子节点的i、j、p(cij)和本节点相同，它们的pid赋值为本节点的线程编号tid，左、右节点的pal分别赋值为alleft和alright，左、右节点的(α,β)分别赋值为(α,γ)和(γ,β)，并将这两个子节点连续存储到下一层的任务队列fqlevel+1中；

步骤314.等待当前任务队列中所有节点均被处理且子节点都已经生成并存储到下一层的任务队列之后，若fqlevel为空，则转入步骤315；否则level＝level+1，转入步骤312；

步骤315.对于反向队列bqlevel中的每一项，读取出其中保存的弧长和pid，将其中保存的弧长加回到bqlevel-1的pid位置。

步骤316.当bqlevel-1每一项都被读取并计算完成后，修改level＝level-1。若level＞0，转入步骤315；否则，转入步骤317；

步骤317.将bq0中所保存的每个小区间段的弧长用于步骤32的输入。

步骤6包括：

步骤61.对于每一个弧长值，通过二分查找在其所在曲线对应的区间弧长表中查找到其所对应的参数段区间编号j；

步骤62.通过二分查找在其参数段所对应的小区间弧长表中搜索到该弧长对应的小区间段的首尾端点[ucij,uci(j+1)]；

步骤63.根据弧长对应的曲线编号i和查找到的区间段的编号j，加载步骤2中所述区间段所应具备曲线参数信息p(cij)，根据查找到的小区间段的首尾端点[ucij,uci(j+1)]确定该弧长对应参数所在的区间，称为“搜索区间”。

步骤7中所述“基于二进制位操作加速二叉树遍历回溯的深度优先搜索”包括：

步骤71.对于给定待搜索弧长al和搜索区间[α,β]及其对应的区间弧长表和小区间弧长表的位置，求出弧长al在这个小区间中从小区间首端点起的相对弧长，初始化节点的二进制编号为bid＝01b，初始化pal为曲线在区间[α,β]的弧长，初始化递归深度r＝0，初始化累积弧长值alacc为0；

步骤72.求搜索区间的中点计算区间、[α,γ]和[γ,β]对应的弧长值,alleft和alright。若不满足精度条件|alleft+alright-pal|＜ε，则转入步骤73；否则转入步骤74；

步骤73.保存alright，修改pal＝alleft，修改(α,β)＝(α,γ)，修改bid为bid＝bid＜＜1，即将bid向左移一位，这相当于进入了二叉树节点的左子树进行搜索，转入步骤72；

步骤74.针对al与alacc+alleft、alacc+alleft+alright的大小关系，会出现三种情况：

1)若al≤alacc+alleft，则说明弧长al对应参数值在左半区间，若满足精度条件|γ-α|＜ξ，则转入步骤75,否则，转入步骤73；

2)若alacc+alleft＜al＜alacc+alleft+alright，这说明弧长al对应的参数在右半区间，若满足精度条件|β-γ|＜ξ，则转入步骤75；否则修改α＝γ，pal为当前节点右兄弟的弧长值，即当前节点的弧长作为alleft时对应的alright，修改bid＝(bid＜＜1)|01b，即将bid左移一位，且将最右端新增加的一位修改为1，修改alacc＝alacc+alleft这相当于进入二叉树节点的右子树进行搜索，转入步骤72；

3)若al＞alacc+alleft+alright，记bid小端连续的1的数量为n1，修改(α,β)＝(β,β+(β-α)×2ⁿ¹)，修改bid＝(bid＞＞n1)|01b，修改pal为此前保存的计算回溯后节点的父节点时对应的alright，修改alacc＝alacc+alleft+alright，完成回溯，转入步骤72；

步骤75.结束树的搜索，保存最终的搜索区间首尾端点。

以下给出具体实施例。

如图1所示，本实施例提供的一种gpu加速的参数曲线弧长计算及弧长参数化方法，包括以下步骤：

步骤1.输入参数曲线的数量、每条参数曲线区间段的个数、以及曲线参数

从用户指定的文件中读入数据。本实施例使用的是3次(4阶)非有理b样条曲线，每条曲线的控制顶点数不等，至少为4，每条曲线的区间段为控制顶点数数减1，每条曲线的参数区间均为[0,1]。曲线数量为1000条，所有曲线的区间段数总和为5014段，曲线参数为每条b样条曲线的节点区间和控制顶点坐标。其中图2为这1000条曲线中的一条，其曲线编号为9。本例中曲线参数中，所有的实数类型以单精度浮点数表达。

步骤2.将每一个区间段等分为n份首尾相接的参数“小区间段”

本实施例中中使用n＝64，即将b样条曲线的每个区间段分为64个首尾相接的小区间段。

为了方便叙述，仅取这1000条曲线中的一条作为说明，即以如图2所示的曲线为例，它有7个区间段，依次是[0.00,0.15]，[0.15,0.31]，[0.31,0.37]，[0.37,0.44]，[0.44,0.63]，[0.63,0.78]，[0.78,1.00]。

每个区间段被划分为64份，共448个小区间段，它们首尾相接，即[0.00,1/64×0.15]，[1/64×0.15,2/64×0.15]等。这1000条曲线总共的5014个区间段，按照这样的方式被划分为320896个小区间段。

每个小区间段除了区间端点外，还保存着对应的曲线编号、区间段编号和对应的参数信息，以小区间段[0.00,1/64×0.15]为例，它除了保存首尾端点0和1/64×0.15外，还保存着曲线编号9、区间段编号0和这一段对应的节点区间以及计算这一段时涉及到的4个控制顶点。

步骤3.构建“小区间弧长表”和“区间弧长表”

采用方法为本发明所提出的“双阶段的广度优先搜索”，对每个小区间段的弧长进行计算，然后构建小区间弧长表和区间弧长表。

步骤31.采用方法为本发明所提出的“双阶段的广度优先搜索”计算满足用户指定误差ε的曲线在该小区间段上的弧长值。小区间段弧长的计算结果按顺序存储，得到“小区间弧长表”，其排序优先级为，曲线编号>区间段编号>小区间段首端点。

本实施例中，采用的误差值为ε＝1.0×10^-7。“双阶段的广度优先搜索”中的使用节点具有如下属性：待求弧长的区间的首尾端点值α和β；所在曲线编号i；所在区间段编号j；区间段所对应曲线的参数信息p(cij)，即计算该小区间涉及的b样条曲线的节点区间和4个控制顶点；上一层搜索中，计算该节点的父节点的线程的编号pid；节点的父节点的弧长pal。

步骤311.从小区间段所包含的参数中构建出搜索树的根节点，其中首尾端点值、曲线编号、区间段编号、参数信息均与小区间段中的数据相同，所有搜索树的根节点按照其对应的小区间段的顺序存储在第0层的任务队列fq0中，初始化层数level＝0。

步骤312.线程从当前任务队列fqlevel中读取线程编号tid所对应的搜索节点，计算其区间中点计算两个子区间[α,γ]和区间[γ,β]的弧长alleft和alright。

步骤313.对比|alleft+alright-pal|与用户指定误差ε的大小，若小于ε，则接受alleft+alright作为该节点的弧长，并将本节点的pid和弧长alleft+alright存储到反向队列bqlevel的tid位置；若不小于ε，则创建左、右两个子节点，两个子节点的i、j、p(cij)和本节点相同，它们的pid赋值为本节点的线程编号tid，左、右节点的pal分别赋值为alleft和alright，左、右节点的(α,β)分别赋值为(α,γ)和(γ,β)，并将这两个子节点连续存储到下一层的任务队列fqlevel+1中。

步骤314.等待当前任务队列中所有节点均被处理且子节点都已经生成并存储到下一层的任务队列之后，若fqlevel为空，则转入步骤315；否则level＝level+1，转入步骤312。

步骤315.对于反向队列bqlevel中的每一项，读取出其中保存的弧长和pid，将其中保存的弧长加回到bqlevel-1的pid位置。

步骤316.当bqlevel-1每一项都被读取并计算完成后，修改level＝level-1。若level＞0，转入步骤315；否则，转入步骤317。

步骤317.将bq0中所保存的每个小区间段的弧长用于步骤32的输入。

在本实施例中，每个节点的每一层的前向阶段和反向阶段，都由一个gpu线程执行，每个gpu线程在每层内，完成一个节点的计算任务后，如还有未被计算的节点，可以被重复使用。

步骤32.计算小区间弧长表中每个区间段的弧长的前缀和，存储，得到小区间弧长表。

为了方便叙述，仅取这1000条曲线中的一条作为说明，即以如图2所示的曲线为例，它的第0段所对应的经步骤31计算出的小区间弧长为0.0310486145，0.0297047738，…，0.00301694591，共64个值，计算前缀和得到0.0310486145，0.0607533902，…，0.696337581，这就是区间段0对应的小区间弧长表。

步骤33.计算每条曲线各个区间段所对应的小区间弧长表中每个区间的最后一项的前缀和，按顺序存储，得到为“区间弧长表”，排序优先级为，曲线编号>区间段编号。

为了方便叙述，仅取这1000条曲线中的一条作为说明，即以图2所示的曲线为例，经过步骤33的计算小区间弧长表最后一项的前缀后，它对应区间弧长表为0.696337581，1.21599603，…，2.57353115，共有7项。

步骤4.读取每条曲线对应的区间弧长表中的最后一项，即弧长，将其输出

为了方便叙述，仅取这1000条曲线中的一条作为说明，即以图2所示曲线为例，从它对应区间弧长表中可以读出最后一项为2.57353115，即它所对应的弧长。

步骤5.输入待求解弧长参数化的弧长值及其所在的曲线编号i

本实施例中，对步骤1中输入的每一条曲线求解其20个点的均匀弧长采样，即每条曲线等弧长地采样20个点，包括首尾端点，每两点间在曲线上的距离均为该曲线弧长的1/19。这1000条曲线以这样的方法，共采样20000个点，步骤的输入即为这些采样点的弧长值和该点对应的曲线的编号。

为了方便叙述，仅取这1000条曲线中的一条作为说明，即以如图2所示的曲线为例，它上面均匀采20个点，对应的弧长值为0，0.135449，…，2.438082,2.573531等共20个值，输入还包括了这20个值对应的该曲线的编号9，9，…，9共20个编号。

步骤6.确定该弧长所在的“搜索区间”，初始化搜索任务，每个弧长的搜索区间的初始化由1个gpu线程完成

步骤61.对于每一个弧长值，通过二分查找在其所在曲线对应的区间弧长表中查找到其所对应的参数段区间编号j。

步骤62.通过二分查找在其参数段所对应的小区间弧长表中搜索到该弧长对应的小区间段的首尾端点[ucij,uci(j+1)]。

步骤63.根据弧长对应的曲线编号i和查找到的区间段的编号j，加载步骤2中所述区间段所应具备曲线参数信息p(cij)，根据查找到的小区间段的首尾端点[ucij,uci(j+1)]确定该弧长对应参数所在的区间，作为“搜索区间”。

为了方便叙述，仅取这1000条曲线中的一条作为说明，即以如图2所示的曲线为例，它上面所采的即弧长为0.135449对应的点，它的弧长值所在的搜索区间由一个gpu线程进行搜索。首先在编号为9的曲线的区间弧长表中，搜索到它的位于首个区间段。然后再在它首个区间段对应的小区间弧长表中二分搜索得到其弧长位于0.116266891和0.142144069之间，对应的参数小区间段的首尾端点4/64×0.15和5/64×0.15，确定这一段为搜索区间。

步骤7.使用本发明提出的“基于二进制位操作加速二叉树遍历回溯的深度优先搜索”，将搜索区间的长度缩小到不大于用户指定的误差范围ξ，每个弧长的弧长参数化由1个gpu线程完成。

在本实施例中，误差范围取值为ξ＝1×10^-7。

步骤71.对于给定待搜索弧长al和搜索区间[α,β]及其对应的区间弧长表和小区间弧长表的位置，求出弧长al在这个小区间中从小区间首端点起的相对弧长，初始化节点的二进制编号为bid＝01b，初始化pal为曲线在区间[α,β]的弧长，初始化递归深度r＝0，初始化累积弧长值alacc为0。

步骤72.求搜索区间的中点计算区间、[α,γ]和[γ,β]对应的弧长值,alleft和alright。若不满足精度条件|alleft+alright-pal|＜ε，则转入步骤73；否则转入步骤74。

步骤73.保存alright，修改pal＝alleft，修改(α,β)＝(α,γ)，修改bid为bid＝bid＜＜1，即将bid向左移一位，这相当于进入了二叉树节点的左子树进行搜索，转入步骤72。

步骤74.针对al与alacc+alleft、alacc+alleft+alright的大小关系，会出现三种情况：

1)若al≤alacc+alleft，则说明弧长al对应参数值在左半区间，若满足精度条件|γ-α|＜ξ，则转入步骤75,否则，转入步骤73；

2)若alacc+alleft＜al＜alacc+alleft+alright，这说明弧长al对应的参数在右半区间，若满足精度条件|β-γ|＜ξ，则转入步骤75；否则修改α＝γ，pal为当前节点右兄弟的弧长值，即当前节点的弧长作为alleft时对应的alright，修改bid＝(bid＜＜1)|01b，即将bid左移一位，且将最右端新增加的一位修改为1，修改alacc＝alacc+alleft这相当于进入二叉树节点的右子树进行搜索，转入步骤72；

3)若al＞alacc+alleft+alright，记bid小端连续的1的数量为n1，修改(α,β)＝(β,β+(β-α)×2ⁿ¹)，修改bid＝(bid＞＞n1)|01b，修改pal为此前保存的计算回溯后节点的父节点时对应的alright，修改alacc＝alacc+alleft+alright，完成回溯，转入步骤72。

步骤75.结束树的搜索，保存最终的搜索区间首尾端点。

为了方便叙述，仅取这1000条曲线中的一条作为说明，即以如图2所示的曲线为例，它上面所采的即弧长为0.135449对应的点，它相对于区间首端点处的相对弧长为0.019182，它首先进入节点的右子树进行搜索，然后再到这个右子树的左子树中进行搜索。经过如步骤72-步骤75之后，最终的搜索区间首尾端点为0.0110914441829664和0.0110937310782901。

步骤8.将最终搜索区间中点作为弧长所对应的参数输出。

为了方便叙述，仅取这1000条曲线中的一条作为说明，即以如图2所示的曲线为例，它上面所采的即弧长为0.135449对应的点，最终的参数值为首尾端点的中点，即0.01109259。

为了方便叙述，仅取这1000条曲线中的一条作为说明，即以如图2所示的曲线为例，它所求出的20个均匀弧长参数化对应的曲线上的点如图3所示。

在本实施例中，在步骤311、步骤312、步骤71、步骤72中所述计算一个区间对应的曲线弧长时，即求解积分时，所采用的方法为3个点的高斯勒让德求积法。

本实施例所采用的gpu型号为nvidiagt1030，cpu为ryzen1700x。本实施例中，本实施例运算耗时与cpu方法的运算耗时对比如图4所示。在弧长计算任务上，本发明耗时为9.7515毫秒，相比于cpu方法的78.4943毫秒，加速比8.05倍；在弧长参数化任务上，在完成弧长表构建后，本发明耗时为1.4065毫秒，相比于cpu方法的56.8079毫秒，加速比40.39倍；综合的看，本发明在弧长参数化任务上的加速比为12.13倍。

本发明用于解决参数曲线弧长计算与弧长参数化这两个计算密集型任务。步骤包括：输入参数曲线数量、曲线参数；利用本发明提出的“双阶段的广度优先搜索”构建弧长表；从弧长表中取出曲线弧长；利用二分搜索方法在弧长表中进行搜索，确定待求参数的弧长的搜索区间，并利用“基于二进制位操作加速二叉树遍历回溯的深度优先搜索”将区间迭代地缩小至用户指定的误差；输出弧长对应的参数。本发明具有更高的可并行性，更为适合gpu的simt架构，相比传统的cpu算法和朴素的递归算法，计算耗时显著减少。

以上所述，仅为本发明较佳具体实施方式，但本发明保护范围并不局限于此，任何熟悉本技术领域的技术人员在本发明揭露的技术范围内，可轻易想到的变化或替换，都应涵盖在本发明的保护范围之内。