用cuda实现ransac圆拟合

上一篇实现了ransac的圆拟合，一般来说基本够用，但如果采样点较多时，考虑实时性要求做加速也是必要的。c++移植到cuda还是比较顺利的。重点理解原子操作。对比下来，速度提升5-8倍。不过精度浮动略大。

在CUDA中，可以通过原子操作轻松避免竞态条件。原子操作能够读取、修改并写回一个值到内存，而不会有其他线程的干扰。这些操作通常适用于共享内存和全局内存。共享内存中的原子操作通常用于防止同一线程块内不同线程之间的竞态条件，而全局内存中的原子操作用于防止不同线程块中的线程之间的竞态条件。需要注意的是，共享内存通常比全局内存快得多。

CUDA支持的原子操作

CUDA提供了多种原子操作函数，包括但不限于：

• atomicAdd()：加法操作，用于将一个值添加到指定的内存地址。
• atomicSub()：减法操作，用于从指定的内存地址减去一个值。
• atomicExch()：交换操作，用于交换内存地址中的值和给定的值。
• atomicMin() 和 atomicMax()：用于计算并存储最小值或最大值。
• atomicInc() 和 atomicDec()：增量和减量操作，用于增加或减少内存地址中的值。
• atomicCAS()：比较并交换操作，用于比较内存地址中的值，如果等于给定值，则替换为新值。
• atomicAnd()、atomicOr() 和 atomicXor()：用于执行位与、位或和位异或操作。

#include <cuda_runtime.h>
#include "device_launch_parameters.h"

#include <curand_kernel.h>
#include <stdio.h>
#include  <device_atomic_functions.h>


__global__ void circleRANSACKernel(const Point2f* points, int cloudSize, float maxErrorThreshold, float* center_x, float* center_y, float* radius, int* bestConsensusCnt, int* modelMeanErrors) 
{

    // 随机选择三个点
    curandState randState;
    curand_init(blockIdx.x * blockDim.x + threadIdx.x, 0, 0, &randState);
    int indices[3];
    for (int i = 0; i < 3; ++i) {
        int index = curand(&randState) % cloudSize;
        indices[i] = index;
    }

    // 计算圆心和半径
    float x1 = points[indices[0]].x;
    float y1 = points[indices[0]].y;
    float x2 = points[indices[1]].x;
    float y2 = points[indices[1]].y;
    float x3 = points[indices[2]].x;
    float y3 = points[indices[2]].y;

    double xa = (x1 + x2) / 2.0, ya = (y1 + y2) / 2.0;
    double xb = (x1 + x3) / 2.0, yb = (y1 + y3) / 2.0;

    double ka = (x1 - x2) / (y2 - y1);
    double kb = (x1 - x3) / (y3 - y1);

    float centerX = (yb - ya + ka * xa - kb * xb) / (ka - kb);
    float centerY = ka * centerX - ka * xa + ya;
    float R = sqrt((centerX - xa) * (centerX - xa) + (centerY - ya) * (centerY - ya));

    // 计算内点数量和平均误差
    int sums = 0;
    double meanError = 0;
    for (int j = 0; j < cloudSize; ++j) {
        Point2f point = points[j];
        double distance = abs(R - sqrt((point.x - centerX) * (point.x - centerX) + (point.y - centerY) * (point.y - centerY)));
        if (distance < maxErrorThreshold) {
            sums++;
            meanError += distance;
        }
    }
    double currentMeanError0 = sums > 0 ? meanError / sums : meanError;
    int currentMeanError = (int)10000*currentMeanError0;
    // 更新最佳模型
    if (sums > cloudSize * 0.5 && currentMeanError < *modelMeanErrors) {
        atomicMax(&bestConsensusCnt[0], sums);
        atomicMin(modelMeanErrors, currentMeanError);
        atomicExch(center_x, centerX);
        atomicExch(center_y, centerY);
        atomicExch(radius, R);
    }
}