Function computeDescriptor

src/backend/cuda/kernel/sift.hpp:686–819

// of Lowe's paper.
template<typename T>
__global__ void computeDescriptor(
    float* desc_out, const unsigned desc_len, const unsigned histsz,
    const float* x_in, const float* y_in, const unsigned* layer_in,
    const float* response_in, const float* size_in, const float* ori_in,
    const unsigned total_feat, const CParam<T> gauss_octave, const int d,
    const int n, const float scale, const int n_layers) {
    const int tid_x = threadIdx.x;
    const int tid_y = threadIdx.y;
    const int bsz_x = blockDim.x;
    const int bsz_y = blockDim.y;

    const int f = blockIdx.y * bsz_y + tid_y;

    SharedMemory<float> shared;
    float* shrdMem = shared.getPointer();
    float* desc    = shrdMem;
    float* accum   = shrdMem + desc_len * histsz;

    for (int i = tid_x; i < desc_len * histsz; i += bsz_x)
        desc[tid_y * desc_len + i] = 0.f;
    __syncthreads();

    if (f < total_feat) {
        const unsigned layer = layer_in[f];
        float ori            = (360.f - ori_in[f]) * PI_VAL / 180.f;
        ori                  = (ori > PI_VAL) ? ori - PI_VAL * 2 : ori;
        const float size     = size_in[f];
        const int fx         = round(x_in[f] * scale);
        const int fy         = round(y_in[f] * scale);

        const int dim0 = gauss_octave.dims[0];
        const int dim1 = gauss_octave.dims[1];
        const int imel = dim0 * dim1;

        // Points img to correct Gaussian pyramid layer
        const T* img_ptr = gauss_octave.ptr + layer * imel;

        float cos_t        = cosf(ori);
        float sin_t        = sinf(ori);
        float bins_per_rad = n / (PI_VAL * 2.f);
        float exp_denom    = d * d * 0.5f;
        float hist_width   = DESCR_SCL_FCTR * size * scale * 0.5f;
        int radius         = hist_width * sqrtf(2.f) * (d + 1.f) * 0.5f + 0.5f;

        int len            = radius * 2 + 1;
        const int hist_off = (tid_x % histsz) * desc_len;

        // Calculate orientation histogram
        for (int l = tid_x; l < len * len; l += bsz_x) {
            int i = l / len - radius;
            int j = l % len - radius;

            int y = fy + i;
            int x = fx + j;

            float x_rot = (j * cos_t - i * sin_t) / hist_width;
            float y_rot = (j * sin_t + i * cos_t) / hist_width;
            float xbin  = x_rot + d / 2 - 0.5f;