CUDA 从非均匀采样中重新采样
CUDA resample from non-even sampling
我想从非均匀采样中重新采样（插值）一个序列。我认为 tex 行不通，因为假设您的样本是均匀的，它基本上会进行插值？搜索会不会太慢？
我应该做推力吗？任何指针表示赞赏。任何例子都会很有帮助。
更新：
说圆圈标记的线是我的样本。我知道每个圆点的价值。显然，样本在水平轴上是均匀分布的。
现在，我想知道采样线下方直线上每个 x 标记处的值。 x标记沿直线均匀分布。
---o--------o----o------o------o------o------ （抽样）
--X-----X-----X-----X-----X-----X-----X---（已知插值）
所以我想知道如何使用 CUDA 获取每个 x 标记位置的值？很明显，使用C/C++的最基本的算法就是针对每个x标记位置，搜索最近的两个圆圈位置，然后进行线性插值。但是在这种情况下，您需要先对两个序列进行排序，然后循环遍历 x 标记，并针对每个 x 标记进行搜索。这听起来很膨胀。
我想知道我们应该如何在 CUDA 中做到这一点？谢谢。
可能有多种方法。例如，您可以使用基本的 cuda binary search in a thread-parallel fashion. I'll demonstrate a thrust 实现。
出于讨论的目的，我假设两个数据集（已知样本点和所需样本点）都是任意放置的（即我不假设任何一个样本集都是均匀分布的）。但是我会规定或要求所需的样本点完全包含在已知样本点内。我相信这是明智的，因为通常的线性插值需要在所需样本点的任一侧都有一个已知的样本点。
因此我们将使用这样的数据集：
o: 1,3,7
f(o): 3,7,15
x: 1.5, 2.5, 4.5, 5.0, 6.0, 6.5
f(x): ?, ?, ?, ?, ?, ?
我们看到f是已知的函数值，对应于f(o) = 2o+1，在这种情况下是一条直线（虽然这种方法不需要已知样本点来拟合任何特定的函数). x 表示我们希望根据已知值 (f(o)) 插入函数值的索引。我们的愿望是通过从最近的 f(o) 点进行插值来计算 f(x)。请注意，我们的数据集是 x 的所有值都位于最小值 (1) 和最大值 (7) o 之间。这就是我之前说的stipulation/requirement
我们的推力方法是使用矢量化二进制搜索，使用 thrust::upper_bound，定位 "insertion point"，其中每个所需的 x 值都适合 o 序列.这为我们提供了用于插值的右邻居和左邻居 (right-1)。一旦我们知道插入点，就可以对该算法进行微不足道的扩展来选择例如two left 和 two right neighbors（或更多）如果我们想使用线性插值以外的东西。
插入点然后给我们我们的左右邻居，我们使用此信息将 x 向量（所需的插值点）和 thrust::transform 向量（所需的插值点）以及 thrust::tuple（通过 thrust::zip_iterator）提供：
- 右邻索引
- 右邻函数值
- 左邻居索引
- 左邻域函数值
有了这些数量，再加上所需的索引 (x)，插值就很简单了。
编辑： 受另一个答案的启发，我决定包含一种避免并行二进制搜索的方法，而是使用前缀和方法来识别插入索引o 数据中的 x 数据。此方法假定 x 和 o 序列都已排序。
我们将从 merge_by_key 操作开始。我们将 x 与 o 合并，以建立排序（这似乎比二进制搜索更有效）。 x 和 o 数量将是 "keys"，并且 o 的值将全部为 1，而 x 的值将全部为 0。然后使用我们的样本数据，merge_by_key 将产生这个：
o keys: 1,3,7
o vals: 1,1,1
x keys: 1.5,2.5,4.5,5.0,6.0,6.5
x vals: 0, 0, 0, 0, 0, 0
merged keys: 1, 1.5, 2.5, 3, 4.5, 5.0, 6.0, 6.5, 7
merged vals: 1, 0, 0, 1, 0, 0, 0, 0, 1
当我们对合并的值进行前缀和（包含扫描）时，我们得到：
ins. ind.: 1, 1, 1, 2, 2, 2, 2, 2, 3
然后我们可以执行 copy_if 操作以仅提取与 x vals（其合并的 vals 为零）关联的插入索引，以生成与步骤 1 中相同的插入索引序列：
d_i: 1, 1, 2, 2, 2, 2
方法 2 的其余部分可以使用与方法 1 中使用的完全相同的剩余插值代码 (thrust::transform)。
这是一个展示这两种方法的完整示例：
$ cat t1224.cu
#include <thrust/device_vector.h>
#include <thrust/binary_search.h>
#include <thrust/transform.h>
#include <thrust/copy.h>
#include <thrust/iterator/zip_iterator.h>
#include <thrust/iterator/permutation_iterator.h>
#include <thrust/iterator/transform_iterator.h>
#include <iostream>
#include <thrust/merge.h>
#include <thrust/iterator/constant_iterator.h>
#include <thrust/scan.h>
struct interp_func
{
template <typename T>
__host__ __device__
float operator()(float t1, T t2){ // m = (y1-y0)/(x1-x0) y = m(x-x0) + y0
return ((thrust::get<1>(t2) - thrust::get<3>(t2))/(thrust::get<0>(t2) - thrust::get<2>(t2)))*(t1 - thrust::get<2>(t2)) + thrust::get<3>(t2);
}
};
using namespace thrust::placeholders;
int main(){
// sample data
float o[] = {1.0f, 3.0f, 7.0f}; // unevenly spaced sample points for function f
float f[] = {3.0f, 7.0f, 15.0f}; // f(o) = 2o+1
float x[] = {1.5f, 2.5f, 4.5f, 5.0f, 6.0f, 6.5f}; // additional desired sample points for f
int so = sizeof(o)/sizeof(o[0]);
int sx = sizeof(x)/sizeof(x[0]);
// setup data on device
thrust::device_vector<float> d_o(o, o+so);
thrust::device_vector<float> d_f(f, f+so);
thrust::device_vector<float> d_x(x, x+sx);
thrust::device_vector<int> d_i(sx); // insertion indices
thrust::device_vector<float> d_r(sx); // results
// method 1: binary search
// perform search for insertion indices
thrust::upper_bound(d_o.begin(), d_o.end(), d_x.begin(), d_x.end(), d_i.begin());
// then perform linear interpolation based on left and right neighbors
std::cout << "Method 1 insertion indices:" << std::endl;
thrust::copy(d_i.begin(), d_i.end(), std::ostream_iterator<int>(std::cout, ","));
std::cout << std::endl;
thrust::transform(d_x.begin(), d_x.end(), thrust::make_zip_iterator(thrust::make_tuple(thrust::make_permutation_iterator(d_o.begin(), d_i.begin()), thrust::make_permutation_iterator(d_f.begin(), d_i.begin()), thrust::make_permutation_iterator(d_o.begin(), thrust::make_transform_iterator(d_i.begin(), _1-1)), thrust::make_permutation_iterator(d_f.begin(), thrust::make_transform_iterator(d_i.begin(), _1-1)))), d_r.begin(), interp_func());
// output results
std::cout << "Interpolation points:" << std::endl;
thrust::copy(d_x.begin(), d_x.end(), std::ostream_iterator<float>(std::cout, ","));
std::cout << std::endl << "Interpolated values:" << std::endl;
thrust::copy(d_r.begin(), d_r.end(), std::ostream_iterator<float>(std::cout, ","));
std::cout << std::endl << "Expected values:" << std::endl;
for (int i = 0; i < sx; i++) std::cout << 2*x[i]+1 << ",";
std::cout << std::endl;
//method 2: merge + prefix sum
thrust::device_vector<float> d_kr(sx+so);
thrust::device_vector<int> d_vr(sx+so);
thrust::device_vector<int> d_s(sx+so);
thrust::merge_by_key(d_o.begin(), d_o.end(), d_x.begin(), d_x.end(), thrust::constant_iterator<int>(1), thrust::constant_iterator<int>(0), d_kr.begin(), d_vr.begin());
thrust::inclusive_scan(d_vr.begin(), d_vr.end(), d_s.begin());
thrust::copy_if(d_s.begin(), d_s.end(), d_vr.begin(), d_i.begin(), _1 == 0);
std::cout << "Method 2 insertion indices:" << std::endl;
thrust::copy(d_i.begin(), d_i.end(), std::ostream_iterator<int>(std::cout, ","));
std::cout << std::endl;
// remainder of solution method would be identical to end of method 1 starting with the thrust::transform
return 0;
}
$ nvcc -o t1224 t1224.cu
$ ./t1224
Method 1 insertion indices:
1,1,2,2,2,2,
Interpolation points:
1.5,2.5,4.5,5,6,6.5,
Interpolated values:
4,6,10,11,13,14,
Expected values:
4,6,10,11,13,14,
Method 2 insertion indices:
1,1,2,2,2,2,
$
同样，一旦我们知道了插入点，选择 2 个右邻点和 2 个左邻点来进行更复杂的插值将是一个微不足道的扩展。我们只需修改传递给变换（插值）仿函数的 zip 迭代器，并修改仿函数本身以实现所需的算法。
另请注意，此方法假定输入 o 序列已经排序。如果不是，则需要添加 o（键）和 f（值）的按键排序。 x序列对于方法1不需要排序，但是对于方法2必须排序（合并需要两个序列都排序）。
最佳方法的细节取决于所涉及的大小（即它是一大批短序列还是单个巨大序列等），但在高层次上你可以只用一个（并行可能是 O(N)) 类型的输入序列和并行前缀和。特别是您可以避免任何二进制搜索。查看 modernGPU "intervalExpand" 背后的思想：https://nvlabs.github.io/moderngpu/intervalmove.html
伪代码简要说明：
1: sort the input sequence
2: for each input point seq[i]:
let count[i] = number of output points in the interval [seq[i], seq[i+1])
3: let indices = exclusive prefix-sum of count
4: use intervalExpand() to go from seq, count, indices to the desired output.
你可以在第 4 步中使用任何你想要的插值公式，包括线性、三次等。重要的是 intervalExpand 会告诉你每个 output 索引，它们是将输出夹在中间的正确 输入 索引。
同样，如果您正在处理一大批较小的序列，二分查找实际上可能比 运行 更快并且更容易编写。否则，您应该能够使用 modernGPU 库中的模板化代码相对轻松地完成此操作。
希望对您有所帮助。
我想从非均匀采样中重新采样（插值）一个序列。我认为 tex 行不通，因为假设您的样本是均匀的，它基本上会进行插值？搜索会不会太慢？
我应该做推力吗？任何指针表示赞赏。任何例子都会很有帮助。
更新：
说圆圈标记的线是我的样本。我知道每个圆点的价值。显然，样本在水平轴上是均匀分布的。
现在，我想知道采样线下方直线上每个 x 标记处的值。 x标记沿直线均匀分布。
---o--------o----o------o------o------o------ （抽样）
--X-----X-----X-----X-----X-----X-----X---（已知插值）
所以我想知道如何使用 CUDA 获取每个 x 标记位置的值？很明显，使用C/C++的最基本的算法就是针对每个x标记位置，搜索最近的两个圆圈位置，然后进行线性插值。但是在这种情况下，您需要先对两个序列进行排序，然后循环遍历 x 标记，并针对每个 x 标记进行搜索。这听起来很膨胀。
我想知道我们应该如何在 CUDA 中做到这一点？谢谢。
可能有多种方法。例如，您可以使用基本的 cuda binary search in a thread-parallel fashion. I'll demonstrate a thrust 实现。
出于讨论的目的，我假设两个数据集（已知样本点和所需样本点）都是任意放置的（即我不假设任何一个样本集都是均匀分布的）。但是我会规定或要求所需的样本点完全包含在已知样本点内。我相信这是明智的，因为通常的线性插值需要在所需样本点的任一侧都有一个已知的样本点。
因此我们将使用这样的数据集：
o: 1,3,7
f(o): 3,7,15
x: 1.5, 2.5, 4.5, 5.0, 6.0, 6.5
f(x): ?, ?, ?, ?, ?, ?
我们看到f是已知的函数值，对应于f(o) = 2o+1，在这种情况下是一条直线（虽然这种方法不需要已知样本点来拟合任何特定的函数). x 表示我们希望根据已知值 (f(o)) 插入函数值的索引。我们的愿望是通过从最近的 f(o) 点进行插值来计算 f(x)。请注意，我们的数据集是 x 的所有值都位于最小值 (1) 和最大值 (7) o 之间。这就是我之前说的stipulation/requirement
我们的推力方法是使用矢量化二进制搜索，使用 thrust::upper_bound，定位 "insertion point"，其中每个所需的 x 值都适合 o 序列.这为我们提供了用于插值的右邻居和左邻居 (right-1)。一旦我们知道插入点，就可以对该算法进行微不足道的扩展来选择例如two left 和 two right neighbors（或更多）如果我们想使用线性插值以外的东西。
插入点然后给我们我们的左右邻居，我们使用此信息将 x 向量（所需的插值点）和 thrust::transform 向量（所需的插值点）以及 thrust::tuple（通过 thrust::zip_iterator）提供：
- 右邻索引
- 右邻函数值
- 左邻居索引
- 左邻域函数值
有了这些数量，再加上所需的索引 (x)，插值就很简单了。
编辑： 受另一个答案的启发，我决定包含一种避免并行二进制搜索的方法，而是使用前缀和方法来识别插入索引o 数据中的 x 数据。此方法假定 x 和 o 序列都已排序。
我们将从 merge_by_key 操作开始。我们将 x 与 o 合并，以建立排序（这似乎比二进制搜索更有效）。 x 和 o 数量将是 "keys"，并且 o 的值将全部为 1，而 x 的值将全部为 0。然后使用我们的样本数据，merge_by_key 将产生这个：
o keys: 1,3,7
o vals: 1,1,1
x keys: 1.5,2.5,4.5,5.0,6.0,6.5
x vals: 0, 0, 0, 0, 0, 0
merged keys: 1, 1.5, 2.5, 3, 4.5, 5.0, 6.0, 6.5, 7
merged vals: 1, 0, 0, 1, 0, 0, 0, 0, 1
当我们对合并的值进行前缀和（包含扫描）时，我们得到：
ins. ind.: 1, 1, 1, 2, 2, 2, 2, 2, 3
然后我们可以执行 copy_if 操作以仅提取与 x vals（其合并的 vals 为零）关联的插入索引，以生成与步骤 1 中相同的插入索引序列：
d_i: 1, 1, 2, 2, 2, 2
方法 2 的其余部分可以使用与方法 1 中使用的完全相同的剩余插值代码 (thrust::transform)。
这是一个展示这两种方法的完整示例：
$ cat t1224.cu
#include <thrust/device_vector.h>
#include <thrust/binary_search.h>
#include <thrust/transform.h>
#include <thrust/copy.h>
#include <thrust/iterator/zip_iterator.h>
#include <thrust/iterator/permutation_iterator.h>
#include <thrust/iterator/transform_iterator.h>
#include <iostream>
#include <thrust/merge.h>
#include <thrust/iterator/constant_iterator.h>
#include <thrust/scan.h>
struct interp_func
{
template <typename T>
__host__ __device__
float operator()(float t1, T t2){ // m = (y1-y0)/(x1-x0) y = m(x-x0) + y0
return ((thrust::get<1>(t2) - thrust::get<3>(t2))/(thrust::get<0>(t2) - thrust::get<2>(t2)))*(t1 - thrust::get<2>(t2)) + thrust::get<3>(t2);
}
};
using namespace thrust::placeholders;
int main(){
// sample data
float o[] = {1.0f, 3.0f, 7.0f}; // unevenly spaced sample points for function f
float f[] = {3.0f, 7.0f, 15.0f}; // f(o) = 2o+1
float x[] = {1.5f, 2.5f, 4.5f, 5.0f, 6.0f, 6.5f}; // additional desired sample points for f
int so = sizeof(o)/sizeof(o[0]);
int sx = sizeof(x)/sizeof(x[0]);
// setup data on device
thrust::device_vector<float> d_o(o, o+so);
thrust::device_vector<float> d_f(f, f+so);
thrust::device_vector<float> d_x(x, x+sx);
thrust::device_vector<int> d_i(sx); // insertion indices
thrust::device_vector<float> d_r(sx); // results
// method 1: binary search
// perform search for insertion indices
thrust::upper_bound(d_o.begin(), d_o.end(), d_x.begin(), d_x.end(), d_i.begin());
// then perform linear interpolation based on left and right neighbors
std::cout << "Method 1 insertion indices:" << std::endl;
thrust::copy(d_i.begin(), d_i.end(), std::ostream_iterator<int>(std::cout, ","));
std::cout << std::endl;
thrust::transform(d_x.begin(), d_x.end(), thrust::make_zip_iterator(thrust::make_tuple(thrust::make_permutation_iterator(d_o.begin(), d_i.begin()), thrust::make_permutation_iterator(d_f.begin(), d_i.begin()), thrust::make_permutation_iterator(d_o.begin(), thrust::make_transform_iterator(d_i.begin(), _1-1)), thrust::make_permutation_iterator(d_f.begin(), thrust::make_transform_iterator(d_i.begin(), _1-1)))), d_r.begin(), interp_func());
// output results
std::cout << "Interpolation points:" << std::endl;
thrust::copy(d_x.begin(), d_x.end(), std::ostream_iterator<float>(std::cout, ","));
std::cout << std::endl << "Interpolated values:" << std::endl;
thrust::copy(d_r.begin(), d_r.end(), std::ostream_iterator<float>(std::cout, ","));
std::cout << std::endl << "Expected values:" << std::endl;
for (int i = 0; i < sx; i++) std::cout << 2*x[i]+1 << ",";
std::cout << std::endl;
//method 2: merge + prefix sum
thrust::device_vector<float> d_kr(sx+so);
thrust::device_vector<int> d_vr(sx+so);
thrust::device_vector<int> d_s(sx+so);
thrust::merge_by_key(d_o.begin(), d_o.end(), d_x.begin(), d_x.end(), thrust::constant_iterator<int>(1), thrust::constant_iterator<int>(0), d_kr.begin(), d_vr.begin());
thrust::inclusive_scan(d_vr.begin(), d_vr.end(), d_s.begin());
thrust::copy_if(d_s.begin(), d_s.end(), d_vr.begin(), d_i.begin(), _1 == 0);
std::cout << "Method 2 insertion indices:" << std::endl;
thrust::copy(d_i.begin(), d_i.end(), std::ostream_iterator<int>(std::cout, ","));
std::cout << std::endl;
// remainder of solution method would be identical to end of method 1 starting with the thrust::transform
return 0;
}
$ nvcc -o t1224 t1224.cu
$ ./t1224
Method 1 insertion indices:
1,1,2,2,2,2,
Interpolation points:
1.5,2.5,4.5,5,6,6.5,
Interpolated values:
4,6,10,11,13,14,
Expected values:
4,6,10,11,13,14,
Method 2 insertion indices:
1,1,2,2,2,2,
$
同样，一旦我们知道了插入点，选择 2 个右邻点和 2 个左邻点来进行更复杂的插值将是一个微不足道的扩展。我们只需修改传递给变换（插值）仿函数的 zip 迭代器，并修改仿函数本身以实现所需的算法。
另请注意，此方法假定输入 o 序列已经排序。如果不是，则需要添加 o（键）和 f（值）的按键排序。 x序列对于方法1不需要排序，但是对于方法2必须排序（合并需要两个序列都排序）。
最佳方法的细节取决于所涉及的大小（即它是一大批短序列还是单个巨大序列等），但在高层次上你可以只用一个（并行可能是 O(N)) 类型的输入序列和并行前缀和。特别是您可以避免任何二进制搜索。查看 modernGPU "intervalExpand" 背后的思想：https://nvlabs.github.io/moderngpu/intervalmove.html
伪代码简要说明：
1: sort the input sequence
2: for each input point seq[i]:
let count[i] = number of output points in the interval [seq[i], seq[i+1])
3: let indices = exclusive prefix-sum of count
4: use intervalExpand() to go from seq, count, indices to the desired output.
你可以在第 4 步中使用任何你想要的插值公式，包括线性、三次等。重要的是 intervalExpand 会告诉你每个 output 索引，它们是将输出夹在中间的正确 输入 索引。
同样，如果您正在处理一大批较小的序列，二分查找实际上可能比 运行 更快并且更容易编写。否则，您应该能够使用 modernGPU 库中的模板化代码相对轻松地完成此操作。
希望对您有所帮助。