harris_response() [1/2]

MAYAFLUX_API void MayaFlux::Kinesis::Vision::harris_response	(	std::span< const float >	gray,
		std::span< float >	dx,
		std::span< float >	dy,
		std::span< float >	tmp,
		std::span< float >	ixx,
		std::span< float >	iyy,
		std::span< float >	ixy,
		std::span< float >	sxx,
		std::span< float >	syy,
		std::span< float >	sxy,
		std::span< float >	dst,
		uint32_t	w,
		uint32_t	h,
		float	k,
		std::span< const float >	kern
	)

Compute the Harris corner response map into caller-supplied buffers.

Writes the normalised response into dst. Intermediate buffers dx/dy/tmp/ixx/iyy/ixy/sxx/syy/sxy must each be size >= w * h and must not alias each other or gray.

Parameters

gray	Single-channel float span, size w * h.
dx	Scratch: horizontal gradient, size >= w * h.
dy	Scratch: vertical gradient, size >= w * h.
tmp	Scratch: filter horizontal pass, size >= w * h.
ixx	Scratch: structure tensor Ixx, size >= w * h.
iyy	Scratch: structure tensor Iyy, size >= w * h.
ixy	Scratch: structure tensor Ixy, size >= w * h.
sxx	Scratch: smoothed Sxx, size >= w * h.
syy	Scratch: smoothed Syy, size >= w * h.
sxy	Scratch: smoothed Sxy, size >= w * h.
dst	Output response map, size >= w * h, values in [0, 1].
w	Image width in pixels.
h	Image height in pixels.
k	Harris sensitivity parameter. Typical range [0.04, 0.06].
kern	Precomputed 1D Gaussian kernel for structure tensor smoothing.

Definition at line 66 of file Harris.cpp.

{
    const size_t n = static_cast<size_t>(w) * h;
 
    sobel(gray, dx, dy, tmp, w, h);
 
    {
        const float* dx_ptr = dx.data();
        const float* dy_ptr = dy.data();
        float* ixx_ptr = ixx.data();
        float* iyy_ptr = iyy.data();
        float* ixy_ptr = ixy.data();
 
        P::for_each(P::par_unseq,
            std::views::iota(size_t { 0 }, n).begin(),
            std::views::iota(size_t { 0 }, n).end(),
            [&](size_t i) {
                const float gx = dx_ptr[i];
                const float gy = dy_ptr[i];
                ixx_ptr[i] = gx * gx;
                iyy_ptr[i] = gy * gy;
                ixy_ptr[i] = gx * gy;
            });
    }
 
    const float* h_src[3] = { ixx.data(), iyy.data(), ixy.data() };
    float* h_dst[3] = { dx.data(), dy.data(), tmp.data() };
    filter_horizontal_planes({ h_src, 3 }, { h_dst, 3 }, w, h, kern);
 
    const float* v_src[3] = { dx.data(), dy.data(), tmp.data() };
    float* v_dst[3] = { sxx.data(), syy.data(), sxy.data() };
    filter_vertical_planes({ v_src, 3 }, { v_dst, 3 }, w, h, kern);
 
    {
        const float* sxx_ptr = sxx.data();
        const float* syy_ptr = syy.data();
        const float* sxy_ptr = sxy.data();
        float* dst_ptr = dst.data();
 
#ifdef MAYAFLUX_ARCH_X64
        const __m256 k_vec = _mm256_set1_ps(k);
        const __m256 zero = _mm256_setzero_ps();
 
        size_t i = 0;
        for (; i + 8 <= n; i += 8) {
            const __m256 s0 = _mm256_loadu_ps(sxx_ptr + i);
            const __m256 s1 = _mm256_loadu_ps(syy_ptr + i);
            const __m256 s01 = _mm256_loadu_ps(sxy_ptr + i);
 
            const __m256 det = _mm256_sub_ps(
                _mm256_mul_ps(s0, s1),
                _mm256_mul_ps(s01, s01));
            const __m256 trace = _mm256_add_ps(s0, s1);
            const __m256 resp = _mm256_max_ps(zero,
                _mm256_fnmadd_ps(k_vec,
                    _mm256_mul_ps(trace, trace),
                    det));
 
            _mm256_storeu_ps(dst_ptr + i, resp);
        }
        for (; i < n; ++i) {
            const float s0 = sxx_ptr[i];
            const float s1 = syy_ptr[i];
            const float s01 = sxy_ptr[i];
            const float det = s0 * s1 - s01 * s01;
            const float trace = s0 + s1;
            dst_ptr[i] = std::max(0.0F, det - k * trace * trace);
        }
 
#elif defined(MAYAFLUX_ARCH_ARM64)
        const float32x4_t k_vec = vdupq_n_f32(k);
        const float32x4_t zero = vdupq_n_f32(0.0F);
 
        size_t i = 0;
        for (; i + 4 <= n; i += 4) {
            const float32x4_t s0 = vld1q_f32(sxx_ptr + i);
            const float32x4_t s1 = vld1q_f32(syy_ptr + i);
            const float32x4_t s01 = vld1q_f32(sxy_ptr + i);
 
            const float32x4_t det = vsubq_f32(
                vmulq_f32(s0, s1),
                vmulq_f32(s01, s01));
            const float32x4_t trace = vaddq_f32(s0, s1);
            const float32x4_t resp = vmaxq_f32(zero,
                vmlsq_f32(det, k_vec, vmulq_f32(trace, trace)));
 
            vst1q_f32(dst_ptr + i, resp);
        }
        for (; i < n; ++i) {
            const float s0 = sxx_ptr[i];
            const float s1 = syy_ptr[i];
            const float s01 = sxy_ptr[i];
            const float det = s0 * s1 - s01 * s01;
            const float trace = s0 + s1;
            dst_ptr[i] = std::max(0.0F, det - k * trace * trace);
        }
 
#else
        for (size_t i = 0; i < n; ++i) {
            const float s0 = sxx_ptr[i];
            const float s1 = syy_ptr[i];
            const float s01 = sxy_ptr[i];
            const float det = s0 * s1 - s01 * s01;
            const float trace = s0 + s1;
            dst_ptr[i] = std::max(0.0F, det - k * trace * trace);
        }
#endif
    }
}

References filter_horizontal_planes(), filter_vertical_planes(), h, k, and sobel().

Here is the call graph for this function: