Function GoldComputeHistogram

tests/cvcuda/system/TestOpSIFT.cpp:467–523

}

inline void GoldComputeHistogram(float (&histogram)[kHistogramBins], float featRadius,
                                 const RawPyramidType &srcGaussianPyramid, const long3 &currStrides,
                                 const long3 &currShape, int octave, int layer, int currBatch, int r, int c)
{
    std::vector<float> tempHistogram(kHistogramBins + 4, 0.f);

    int radius = std::round(featRadius * kOrientationRadius);

    float weightScale = -1.f / (2.f * (featRadius * kOrientationSigma) * (featRadius * kOrientationSigma));

    auto gaussVal = [&srcGaussianPyramid, &currStrides, &octave, &layer, &currBatch](int row, int col)
    {
        return util::ValueAt<WT>(srcGaussianPyramid[octave][layer], currStrides, long3{currBatch, row, col});
    };

    for (int i = -radius; i <= radius; i++)
    {
        if (r + i <= 0 || r + i >= currShape.y - 1)
        {
            continue;
        }

        for (int j = -radius; j <= radius; j++)
        {
            if (c + j <= 0 || c + j >= currShape.z - 1)
            {
                continue;
            }

            float dx = gaussVal(r + i + 0, c + j + 1) - gaussVal(r + i + 0, c + j - 1);
            float dy = gaussVal(r + i - 1, c + j + 0) - gaussVal(r + i + 1, c + j + 0);

            float angle     = std::atan2(dy, dx) * 180.f / M_PI;
            float weight    = std::exp2f((i * i + j * j) * weightScale);
            float magnitude = std::sqrt(dx * dx + dy * dy);

            int bin = std::round(angle * kHistogramBins / 360.f);

            bin = (bin >= kHistogramBins ? bin - kHistogramBins : (bin < 0 ? bin + kHistogramBins : bin));

            tempHistogram[2 + bin] += weight * magnitude;
        }
    }

    tempHistogram[0] = tempHistogram[2 + kHistogramBins - 2];
    tempHistogram[1] = tempHistogram[2 + kHistogramBins - 1];

    tempHistogram[2 + kHistogramBins + 0] = tempHistogram[2 + 0];
    tempHistogram[2 + kHistogramBins + 1] = tempHistogram[2 + 1];

    for (int i = 0; i < kHistogramBins; i++)
    {
        histogram[i] = (tempHistogram[2 + i - 2] + tempHistogram[2 + i + 2]) * 1.f / 16
                     + (tempHistogram[2 + i - 1] + tempHistogram[2 + i + 1]) * 4.f / 16
                     + (tempHistogram[2 + i + 0]) * 6.f / 16;
    }
}