potok/opencv
1
0
Fork 0
mirror of https://github.com/opencv/opencv.git synced 2026-07-21 19:33:03 +04:00
Code Packages Releases Wiki Activity

Merge pull request #29450 from manand881:feat/ocl-sift-descriptor

Browse Source 
Feat: add OpenCL SIFT detector and descriptor #29450

Add OpenCL implementation of SIFT detector and
descriptor extractor, triggered via T-API when
the caller passes a UMat. Targets the features
module (5.x directory layout).

Implementation:
- 4 OpenCL kernels in sift.cl: gaussian_blur_h, gaussian_blur_v, detect_and_orient, and compute_descriptor.
- Host-side dispatch in sift.dispatch.cpp using CV_OCL_RUN_ macro, pure additions to the existing CPU path.
- T-API entry: detectAndCompute, detect, compute.

Gating and fallback:
- 64-bit only. 32-bit falls back to CPU via sizeof(void*) > 4 guard.
- OpenCL 1.1 or later required. Devices below 1.1 fall back to CPU.
- No-OpenCL builds fall back to CPU. Verified with WITH_OPENCL=OFF: 79 tests pass, same 3 pre-existing failures (DISK/AFFINE_FEATURE missing .npy data files).

Correctness:
- nOctaveLayers > 3 OOB fix: gk_coeffs and gk_radius resized to nOctaveLayers+2 via std::vector (were fixed size 5).
- Regression test for nOctaveLayers > 3.
- 10 OCL SIFT tests pass, 60 OCL tests pass, 139 features tests pass (3 pre-existing failures unrelated to SIFT).

Performance optimizations:
- native_exp, native_sqrt, native_recip for hardware approximations (OpenCL 1.1 builtins).
- Eliminate intermediate rawDst[128] buffer; normalize in-place from hist[].
- Non-blocking kernel launches with single ocl::finish() before host reads keypoint count.
- Custom separable Gaussian blur for init image bypassing T-API GaussianBlur dispatch overhead.
- exp2() replaces pow(2.0f, x); scl_octv reused.
- Hoist UMat tmp allocation out of pyramid loop.
- Interleaved keypoint output buffer (6 arrays to 1) for coalesced writes and fewer copies.
- Consolidated descriptor keypoint copies (2N to 2 host-to-device transfers).
- std::map replaced with flat vector indexed by pyramid level (O(1) vs O(log n)).
- DoG computed on-the-fly in detect kernel via READ_DOG macro, eliminating 45 subtract() calls and dog_pack allocation.
- Shared gauss_packs between detect and descriptor, eliminating duplicate copyTo.
- Cross-level packing per octave: one descriptor kernel launch per octave instead of per level. Keypoint buffer expanded to 5 floats (added layer index). Levels packed into single UMat. Reduces kernel launch overhead and improves GPU utilization for octaves with few keypoints.
- __local memory tiling for Gaussian blur kernels with 16x16 workgroups and cooperative halo loading. Reduces global memory traffic.
- Direct uchar descriptor output for CV_8U: kernel writes convert_uchar_sat_rte directly instead of float buffer + host convertTo pass. Eliminates temp allocation and extra kernel launch.

Performance result (stitching/s2.jpg resized):
- DetectAndCompute: OCL wins at >=960x540 (1.32x at 960x540, 1.80x at 1280x720, 1.87x at 1920x1080).
- Compute-only: OCL wins at >=1280x720 (1.60x at 1280x720, 1.78x at 1920x1080).

Tests:
- modules/features/test/ocl/test_feature2d.cpp with OCL SIFT correctness tests on real images (leuven img1.png, a3.png, s2.jpg).
- DescriptorType, Regression_26139, and Batch tests mirror CPU coverage (CV_8U descriptor type, single-keypoint edge case, 6-image detect+compute).
- modules/features/perf/opencl/perf_sift.cpp with real image and 5-size sweep.
- CPU SIFT_Scaled fixture added to perf_feature2d.cpp for aligned OCL vs CPU comparison.

### Pull Request Readiness Checklist

See details at https://github.com/opencv/opencv/wiki/How_to_contribute#making-a-good-pull-request

- [x] I agree to contribute to the project under Apache 2 License.
- [x] To the best of my knowledge, the proposed patch is not based on a code under GPL or another license that is incompatible with OpenCV
- [x] The PR is proposed to the proper branch
- [ ] There is a reference to the original bug report and related work
- [x] There is accuracy test, performance test and test data in opencv_extra repository, if applicable
      Patch to opencv_extra has the same branch name.
- [x] The feature is well documented and sample code can be built with the project CMake
This commit is contained in:
Anand Mahesh
2026-07-15 14:59:43 +05:30
committed by GitHub
parent 5c7b438295
commit 60c705e93b
4 changed files with 1379 additions and 0 deletions
Download Patch File Download Diff File Expand all files Collapse all files
modules/features/perf/opencl/perf_sift.cpp
+63
View File
@@ -0,0 +1,63 @@
// This file is part of OpenCV project.
// It is subject to the license terms in the LICENSE file found in the top-level directory
// of this distribution and at http://opencv.org/license.html
#include "../perf_precomp.hpp"
#include "opencv2/ts/ocl_perf.hpp"
#ifdef HAVE_OPENCL
namespace opencv_test {
namespace ocl {
typedef TestBaseWithParam<Size> SIFTOclFixture;
OCL_PERF_TEST_P(SIFTOclFixture, DetectAndCompute, ::testing::Values(
Size(320, 240), Size(640, 480), Size(960, 540),
Size(1280, 720), Size(1920, 1080)))
{
Size sz = GetParam();
Mat msrc = imread(getDataPath("stitching/s2.jpg"), IMREAD_GRAYSCALE);
ASSERT_FALSE(msrc.empty()) << "Failed to load stitching/s2.jpg";
Mat mimg;
resize(msrc, mimg, sz, 0, 0, INTER_LINEAR);
UMat img = mimg.getUMat(ACCESS_READ), mask;
declare.in(img, WARMUP_READ);
Ptr<SIFT> sift = SIFT::create();
vector<KeyPoint> points;
UMat descriptors;
OCL_TEST_CYCLE() sift->detectAndCompute(img, mask, points, descriptors, false);
EXPECT_GT(points.size(), 20u);
EXPECT_EQ((size_t)descriptors.rows, points.size());
SANITY_CHECK_NOTHING();
}
OCL_PERF_TEST_P(SIFTOclFixture, Compute, ::testing::Values(
Size(640, 480), Size(1280, 720), Size(1920, 1080)))
{
Size sz = GetParam();
Mat msrc = imread(getDataPath("stitching/s2.jpg"), IMREAD_GRAYSCALE);
ASSERT_FALSE(msrc.empty()) << "Failed to load stitching/s2.jpg";
Mat mimg;
resize(msrc, mimg, sz, 0, 0, INTER_LINEAR);
UMat img = mimg.getUMat(ACCESS_READ);
declare.in(img, WARMUP_READ);
Ptr<SIFT> sift = SIFT::create();
vector<KeyPoint> points;
sift->detect(img, points);
UMat descriptors;
OCL_TEST_CYCLE() sift->compute(img, points, descriptors);
EXPECT_EQ((size_t)descriptors.rows, points.size());
SANITY_CHECK_NOTHING();
}
} // ocl
} // opencv_test
#endif // HAVE_OPENCL
modules/features/src/opencl/sift.cl
+659
View File
@@ -0,0 +1,659 @@
// OpenCL port of the SIFT detector and descriptor extractor
// See sift.simd.hpp for the CPU reference.
// Patent US6711293 expired in March 2020.
#define SIFT_DESCR_WIDTH 4
#define SIFT_DESCR_HIST_BINS 8
#define SIFT_DESCR_SCL_FCTR 3.0f
#define SIFT_DESCR_MAG_THR 0.2f
#define SIFT_INT_DESCR_FCTR 512.0f
#define SIFT_D SIFT_DESCR_WIDTH
#define SIFT_N SIFT_DESCR_HIST_BINS
#define SIFT_DH (SIFT_D + 2)
#define SIFT_NH (SIFT_N + 2)
#define SIFT_HIST_LEN (SIFT_DH * SIFT_DH * SIFT_NH)
#define SIFT_DESCR_LEN (SIFT_D * SIFT_D * SIFT_N)
#define SIFT_PI 3.14159265358979323846f
#define SIFT_RAD2DEG (180.0f / SIFT_PI)
#define SIFT_FLT_EPS 1.1920928955078125e-7f
#define SIFT_IMG_BORDER 5
#define SIFT_MAX_INTERP_STEPS 5
#define SIFT_ORI_HIST_BINS 36
#define SIFT_ORI_SIG_FCTR 1.5f
#define SIFT_ORI_RADIUS 4.5f
#define SIFT_ORI_PEAK_RATIO 0.8f
#define SIFT_BLUR_WG_X 16
#define SIFT_BLUR_WG_Y 16
#define SIFT_BLUR_MAX_RADIUS 32
#define SIFT_DESC_KP_PER_WG 4
#define SIFT_DESC_THREADS_PER_KP 16
#define SIFT_DESC_WG_SIZE (SIFT_DESC_KP_PER_WG * SIFT_DESC_THREADS_PER_KP)
#define SIFT_DESC_BINS (SIFT_N + 2)
#define READ_F32(buf, step, ofs, x, y) \
(*(__global const float*)((buf) + (ofs) + (y) * (step) + (x) * 4))
#define READ_DOG(gp, gs, go, x, y, lc, rows) \
(READ_F32(gp, gs, go, x, (y) + (lc+1)*(rows)) - \
READ_F32(gp, gs, go, x, (y) + (lc)*(rows)))
__kernel void
SIFT_gaussian_blur_h(
__global const uchar* src, int src_step, int src_ofs,
__global uchar* dst, int dst_step, int dst_ofs,
int cols, int rows,
__global const float* coeffs, int radius)
{
int x = get_global_id(0);
int y = get_global_id(1);
int lx = get_local_id(0);
int ly = get_local_id(1);
int gx = get_group_id(0) * SIFT_BLUR_WG_X;
int gy = get_group_id(1) * SIFT_BLUR_WG_Y;
__local float lbuf[SIFT_BLUR_WG_Y][SIFT_BLUR_WG_X + 2 * SIFT_BLUR_MAX_RADIUS];
int tile_w = SIFT_BLUR_WG_X + 2 * radius;
int src_y = gy + ly;
if (src_y < rows)
{
for (int col = lx; col < tile_w; col += SIFT_BLUR_WG_X)
{
int xx = gx - radius + col;
if (xx < 0) xx = -xx;
if (xx >= cols) xx = 2 * cols - xx - 2;
lbuf[ly][col] = READ_F32(src, src_step, src_ofs, xx, src_y);
}
}
barrier(CLK_LOCAL_MEM_FENCE);
if (x >= cols || y >= rows)
return;
float sum = 0.0f;
for (int k = -radius; k <= radius; k++)
{
sum += lbuf[ly][lx + radius + k] * coeffs[k + radius];
}
*((__global float*)(dst + dst_ofs + y * dst_step + x * 4)) = sum;
}
__kernel void
SIFT_gaussian_blur_v(
__global const uchar* src, int src_step, int src_ofs,
__global uchar* dst, int dst_step, int dst_ofs,
int cols, int rows,
__global const float* coeffs, int radius)
{
int x = get_global_id(0);
int y = get_global_id(1);
int lx = get_local_id(0);
int ly = get_local_id(1);
int gx = get_group_id(0) * SIFT_BLUR_WG_X;
int gy = get_group_id(1) * SIFT_BLUR_WG_Y;
__local float lbuf[SIFT_BLUR_WG_Y + 2 * SIFT_BLUR_MAX_RADIUS][SIFT_BLUR_WG_X];
int tile_h = SIFT_BLUR_WG_Y + 2 * radius;
for (int row = ly; row < tile_h; row += SIFT_BLUR_WG_Y)
{
int xx = gx + lx;
int yy = gy - radius + row;
if (xx >= cols) xx = cols - 1;
if (yy < 0) yy = -yy;
if (yy >= rows) yy = 2 * rows - yy - 2;
lbuf[row][lx] = READ_F32(src, src_step, src_ofs, xx, yy);
}
barrier(CLK_LOCAL_MEM_FENCE);
if (x >= cols || y >= rows)
return;
float sum = 0.0f;
for (int k = -radius; k <= radius; k++)
{
sum += lbuf[ly + radius + k][lx] * coeffs[k + radius];
}
*((__global float*)(dst + dst_ofs + y * dst_step + x * 4)) = sum;
}
#define SIFT_SCATTER(cr, cc, w0, w1) \
if ((cr) >= 0 && (cr) < SIFT_D && (cc) >= 0 && (cc) < SIFT_D) { \
int dri = (cr) - ti, dci = (cc) - tj; \
if (dri == 0 && dci == 0) { \
my_bins[o0] += (w0); \
my_bins[o0 + 1] += (w1); \
} else if (dri == 1 && dci == 0) { \
border_down[kid_in_wg][tid][o0] += (w0); \
border_down[kid_in_wg][tid][o0 + 1] += (w1); \
} else if (dri == 0 && dci == 1) { \
border_right[kid_in_wg][tid][o0] += (w0); \
border_right[kid_in_wg][tid][o0 + 1] += (w1); \
} else { \
border_down_right[kid_in_wg][tid][o0] += (w0); \
border_down_right[kid_in_wg][tid][o0 + 1] += (w1); \
} \
}
__kernel void
SIFT_compute_descriptor(
__global const uchar* img, int img_step, int img_offset, int img_cols, int img_rows,
float diag,
__global const uchar* kpts, int kpts_step, int kpts_offset,
__global const int* out_rows, int row_start,
int nkeypoints,
__global uchar* desc, int desc_step, int desc_offset,
int descriptor_type)
{
int lid = get_local_id(0);
int kid_in_wg = lid >> 4;
int tid = lid & 15;
int ti = tid >> 2;
int tj = tid & 3;
int kid = get_group_id(0) * SIFT_DESC_KP_PER_WG + kid_in_wg;
bool valid = (kid < nkeypoints);
__local float border_down[SIFT_DESC_KP_PER_WG][16][SIFT_DESC_BINS];
__local float border_right[SIFT_DESC_KP_PER_WG][16][SIFT_DESC_BINS];
__local float border_down_right[SIFT_DESC_KP_PER_WG][16][SIFT_DESC_BINS];
__local float nrm2_buf[SIFT_DESC_KP_PER_WG][16];
for (int k = 0; k < SIFT_DESC_BINS; k++)
{
border_down[kid_in_wg][tid][k] = 0.0f;
border_right[kid_in_wg][tid][k] = 0.0f;
border_down_right[kid_in_wg][tid][k] = 0.0f;
}
float my_bins[SIFT_DESC_BINS];
for (int k = 0; k < SIFT_DESC_BINS; k++)
my_bins[k] = 0.0f;
float cos_t = 0.0f, sin_t = 0.0f;
float bins_per_rad = 0.0f, exp_scale = 0.0f;
int radius = 0, pt_x = 0, pt_y = 0;
float ori = 0.0f;
__global const uchar* img_base = 0;
int row = 0;
if (valid)
{
__global const float* kpt = (__global const float*)(kpts + kpts_offset + kid * kpts_step);
float ptx = kpt[0];
float pty = kpt[1];
ori = kpt[2];
float scl = kpt[3];
int level_idx = (int)kpt[4];
pt_x = convert_int_rte(ptx);
pt_y = convert_int_rte(pty);
cos_t = cos(ori * SIFT_PI / 180.0f);
sin_t = sin(ori * SIFT_PI / 180.0f);
bins_per_rad = (float)SIFT_N / 360.0f;
exp_scale = -1.0f / ((float)SIFT_D * (float)SIFT_D * 0.5f);
float hist_width = SIFT_DESCR_SCL_FCTR * scl;
radius = convert_int_rte(hist_width * 1.4142135623730951f * (float)(SIFT_D + 1) * 0.5f);
if ((float)radius > diag)
radius = convert_int_rte(diag);
float inv_hist_width = native_recip(hist_width);
cos_t *= inv_hist_width;
sin_t *= inv_hist_width;
img_base = img + img_offset + (long)level_idx * img_rows * img_step;
row = out_rows[kid + row_start];
}
barrier(CLK_LOCAL_MEM_FENCE);
if (valid)
{
float r_min = (ti == 0) ? -1.0f : (float)ti;
float r_max = (ti == 0) ? 1.0f : (float)(ti + 1);
float c_min = (tj == 0) ? -1.0f : (float)tj;
float c_max = (tj == 0) ? 1.0f : (float)(tj + 1);
float rr0 = r_min - 1.5f, rr1 = r_max - 1.5f;
float cc0 = c_min - 1.5f, cc1 = c_max - 1.5f;
float det = cos_t * cos_t + sin_t * sin_t;
float inv_det = native_recip(det);
float i1 = (-cc0 * sin_t + rr0 * cos_t) * inv_det;
float i2 = (-cc0 * sin_t + rr1 * cos_t) * inv_det;
float i3 = (-cc1 * sin_t + rr0 * cos_t) * inv_det;
float i4 = (-cc1 * sin_t + rr1 * cos_t) * inv_det;
float j1 = (cc0 * cos_t + rr0 * sin_t) * inv_det;
float j2 = (cc0 * cos_t + rr1 * sin_t) * inv_det;
float j3 = (cc1 * cos_t + rr0 * sin_t) * inv_det;
float j4 = (cc1 * cos_t + rr1 * sin_t) * inv_det;
int i_lo = max(-radius, (int)floor(fmin(fmin(i1, i2), fmin(i3, i4))) - 1);
int i_hi = min(radius, (int)ceil(fmax(fmax(i1, i2), fmax(i3, i4))) + 1);
int j_lo = max(-radius, (int)floor(fmin(fmin(j1, j2), fmin(j3, j4))) - 1);
int j_hi = min(radius, (int)ceil(fmax(fmax(j1, j2), fmax(j3, j4))) + 1);
for (int i = i_lo; i <= i_hi; i++)
{
for (int j = j_lo; j <= j_hi; j++)
{
float c_rot = (float)j * cos_t - (float)i * sin_t;
float r_rot = (float)j * sin_t + (float)i * cos_t;
float rbin = r_rot + (float)(SIFT_D / 2) - 0.5f;
float cbin = c_rot + (float)(SIFT_D / 2) - 0.5f;
int r = pt_y + i;
int c = pt_x + j;
if (rbin > -1.0f && rbin < (float)SIFT_D &&
cbin > -1.0f && cbin < (float)SIFT_D &&
r > 0 && r < img_rows - 1 && c > 0 && c < img_cols - 1)
{
int r0 = (int)floor(rbin);
int c0 = (int)floor(cbin);
int my_r0 = (r0 > 0) ? r0 : 0;
int my_c0 = (c0 > 0) ? c0 : 0;
if (my_r0 != ti || my_c0 != tj)
continue;
__global const uchar* row_base = img_base + r * img_step;
int c4 = c * 4;
float dx = *(__global const float*)(row_base + c4 + 4) -
*(__global const float*)(row_base + c4 - 4);
float dy = *(__global const float*)(row_base - img_step + c4) -
*(__global const float*)(row_base + img_step + c4);
float w = (c_rot * c_rot + r_rot * r_rot) * exp_scale;
float ang = atan2(dy, dx) * SIFT_RAD2DEG;
if (ang < 0.0f)
ang += 360.0f;
float mag = native_sqrt(dx * dx + dy * dy);
float W = native_exp(w);
float obin = (ang - ori) * bins_per_rad;
float magW = mag * W;
rbin -= (float)r0;
cbin -= (float)c0;
int o0 = (int)floor(obin);
obin -= (float)o0;
if (o0 < 0)
o0 += SIFT_N;
if (o0 >= SIFT_N)
o0 -= SIFT_N;
float v_r1 = magW * rbin, v_r0 = magW - v_r1;
float v_rc11 = v_r1 * cbin, v_rc10 = v_r1 - v_rc11;
float v_rc01 = v_r0 * cbin, v_rc00 = v_r0 - v_rc01;
float v_rco111 = v_rc11 * obin, v_rco110 = v_rc11 - v_rco111;
float v_rco101 = v_rc10 * obin, v_rco100 = v_rc10 - v_rco101;
float v_rco011 = v_rc01 * obin, v_rco010 = v_rc01 - v_rco011;
float v_rco001 = v_rc00 * obin, v_rco000 = v_rc00 - v_rco001;
SIFT_SCATTER(r0, c0, v_rco000, v_rco001);
SIFT_SCATTER(r0, c0 + 1, v_rco010, v_rco011);
SIFT_SCATTER(r0 + 1, c0, v_rco100, v_rco101);
SIFT_SCATTER(r0 + 1, c0 + 1, v_rco110, v_rco111);
}
}
}
}
barrier(CLK_LOCAL_MEM_FENCE);
if (valid)
{
if (ti > 0)
{
int above_tid = ((ti - 1) << 2) | tj;
for (int k = 0; k < SIFT_DESC_BINS; k++)
my_bins[k] += border_down[kid_in_wg][above_tid][k];
}
if (tj > 0)
{
int left_tid = (ti << 2) | (tj - 1);
for (int k = 0; k < SIFT_DESC_BINS; k++)
my_bins[k] += border_right[kid_in_wg][left_tid][k];
}
if (ti > 0 && tj > 0)
{
int ul_tid = ((ti - 1) << 2) | (tj - 1);
for (int k = 0; k < SIFT_DESC_BINS; k++)
my_bins[k] += border_down_right[kid_in_wg][ul_tid][k];
}
my_bins[0] += my_bins[SIFT_N];
my_bins[1] += my_bins[SIFT_N + 1];
float my_nrm2 = 0.0f;
for (int k = 0; k < SIFT_N; k++)
my_nrm2 += my_bins[k] * my_bins[k];
nrm2_buf[kid_in_wg][tid] = my_nrm2;
}
else
{
nrm2_buf[kid_in_wg][tid] = 0.0f;
}
barrier(CLK_LOCAL_MEM_FENCE);
if (valid)
{
float total_nrm2 = 0.0f;
for (int k = 0; k < 16; k++)
total_nrm2 += nrm2_buf[kid_in_wg][k];
float thr = native_sqrt(total_nrm2) * SIFT_DESCR_MAG_THR;
float my_nrm2 = 0.0f;
for (int k = 0; k < SIFT_N; k++)
{
float val = my_bins[k];
if (val > thr)
val = thr;
my_bins[k] = val;
my_nrm2 += val * val;
}
nrm2_buf[kid_in_wg][tid] = my_nrm2;
}
else
{
nrm2_buf[kid_in_wg][tid] = 0.0f;
}
barrier(CLK_LOCAL_MEM_FENCE);
if (valid)
{
float total_nrm2 = 0.0f;
for (int k = 0; k < 16; k++)
total_nrm2 += nrm2_buf[kid_in_wg][k];
float inv = SIFT_INT_DESCR_FCTR / fmax(native_sqrt(total_nrm2), SIFT_FLT_EPS);
int didx = (ti * SIFT_D + tj) * SIFT_N;
if (descriptor_type == 0)
{
__global uchar* dst = desc + desc_offset + row * desc_step;
for (int k = 0; k < SIFT_N; k++)
dst[didx + k] = convert_uchar_sat_rte(my_bins[k] * inv);
}
else
{
__global float* dst = (__global float*)(desc + desc_offset + row * desc_step);
for (int k = 0; k < SIFT_N; k++)
{
float v = round(my_bins[k] * inv);
dst[didx + k] = clamp(v, 0.0f, 255.0f);
}
}
}
}
#undef SIFT_SCATTER
__kernel void
SIFT_detect_and_orient(
__global const uchar* gauss_pack, int gauss_step, int gauss_offset,
int dog_cols, int dog_rows_per_layer,
int threshold,
float contrastThreshold, float edgeThreshold, float sigma,
int nOctaveLayers, int octave, int layer,
__global int* output_count,
__global float* output_kpts,
int max_output)
{
int c = get_global_id(0);
int r = get_global_id(1);
if (c < SIFT_IMG_BORDER || c >= dog_cols - SIFT_IMG_BORDER ||
r < SIFT_IMG_BORDER || r >= dog_rows_per_layer - SIFT_IMG_BORDER)
return;
int lc = layer;
float val = READ_DOG(gauss_pack, gauss_step, gauss_offset, c, r, lc, dog_rows_per_layer);
if (fabs(val) <= (float)threshold)
return;
float cv00 = READ_DOG(gauss_pack, gauss_step, gauss_offset, c-1, r-1, lc, dog_rows_per_layer);
float cv01 = READ_DOG(gauss_pack, gauss_step, gauss_offset, c, r-1, lc, dog_rows_per_layer);
float cv02 = READ_DOG(gauss_pack, gauss_step, gauss_offset, c+1, r-1, lc, dog_rows_per_layer);
float cv10 = READ_DOG(gauss_pack, gauss_step, gauss_offset, c-1, r, lc, dog_rows_per_layer);
float cv12 = READ_DOG(gauss_pack, gauss_step, gauss_offset, c+1, r, lc, dog_rows_per_layer);
float cv20 = READ_DOG(gauss_pack, gauss_step, gauss_offset, c-1, r+1, lc, dog_rows_per_layer);
float cv21 = READ_DOG(gauss_pack, gauss_step, gauss_offset, c, r+1, lc, dog_rows_per_layer);
float cv22 = READ_DOG(gauss_pack, gauss_step, gauss_offset, c+1, r+1, lc, dog_rows_per_layer);
float pv00 = READ_DOG(gauss_pack, gauss_step, gauss_offset, c-1, r-1, lc-1, dog_rows_per_layer);
float pv01 = READ_DOG(gauss_pack, gauss_step, gauss_offset, c, r-1, lc-1, dog_rows_per_layer);
float pv02 = READ_DOG(gauss_pack, gauss_step, gauss_offset, c+1, r-1, lc-1, dog_rows_per_layer);
float pv10 = READ_DOG(gauss_pack, gauss_step, gauss_offset, c-1, r, lc-1, dog_rows_per_layer);
float pv11 = READ_DOG(gauss_pack, gauss_step, gauss_offset, c, r, lc-1, dog_rows_per_layer);
float pv12 = READ_DOG(gauss_pack, gauss_step, gauss_offset, c+1, r, lc-1, dog_rows_per_layer);
float pv20 = READ_DOG(gauss_pack, gauss_step, gauss_offset, c-1, r+1, lc-1, dog_rows_per_layer);
float pv21 = READ_DOG(gauss_pack, gauss_step, gauss_offset, c, r+1, lc-1, dog_rows_per_layer);
float pv22 = READ_DOG(gauss_pack, gauss_step, gauss_offset, c+1, r+1, lc-1, dog_rows_per_layer);
float nv00 = READ_DOG(gauss_pack, gauss_step, gauss_offset, c-1, r-1, lc+1, dog_rows_per_layer);
float nv01 = READ_DOG(gauss_pack, gauss_step, gauss_offset, c, r-1, lc+1, dog_rows_per_layer);
float nv02 = READ_DOG(gauss_pack, gauss_step, gauss_offset, c+1, r-1, lc+1, dog_rows_per_layer);
float nv10 = READ_DOG(gauss_pack, gauss_step, gauss_offset, c-1, r, lc+1, dog_rows_per_layer);
float nv11 = READ_DOG(gauss_pack, gauss_step, gauss_offset, c, r, lc+1, dog_rows_per_layer);
float nv12 = READ_DOG(gauss_pack, gauss_step, gauss_offset, c+1, r, lc+1, dog_rows_per_layer);
float nv20 = READ_DOG(gauss_pack, gauss_step, gauss_offset, c-1, r+1, lc+1, dog_rows_per_layer);
float nv21 = READ_DOG(gauss_pack, gauss_step, gauss_offset, c, r+1, lc+1, dog_rows_per_layer);
float nv22 = READ_DOG(gauss_pack, gauss_step, gauss_offset, c+1, r+1, lc+1, dog_rows_per_layer);
bool is_extremum = false;
if (val > 0.0f)
{
float vmax = fmax(fmax(fmax(cv00, cv01), fmax(cv02, cv10)), fmax(fmax(cv12, cv20), fmax(cv21, cv22)));
vmax = fmax(vmax, fmax(fmax(fmax(pv00, pv01), fmax(pv02, pv10)), fmax(fmax(pv12, pv20), fmax(pv21, pv22))));
vmax = fmax(vmax, fmax(fmax(fmax(nv00, nv01), fmax(nv02, nv10)), fmax(fmax(nv12, nv20), fmax(nv21, nv22))));
vmax = fmax(vmax, fmax(pv11, nv11));
is_extremum = (val >= vmax);
}
else
{
float vmin = fmin(fmin(fmin(cv00, cv01), fmin(cv02, cv10)), fmin(fmin(cv12, cv20), fmin(cv21, cv22)));
vmin = fmin(vmin, fmin(fmin(fmin(pv00, pv01), fmin(pv02, pv10)), fmin(fmin(pv12, pv20), fmin(pv21, pv22))));
vmin = fmin(vmin, fmin(fmin(fmin(nv00, nv01), fmin(nv02, nv10)), fmin(fmin(nv12, nv20), fmin(nv21, nv22))));
vmin = fmin(vmin, fmin(pv11, nv11));
is_extremum = (val <= vmin);
}
if (!is_extremum)
return;
const float img_scale = 1.0f / 255.0f;
const float deriv_scale = img_scale * 0.5f;
const float second_deriv_scale = img_scale;
const float cross_deriv_scale = img_scale * 0.25f;
int rc = r, cc = c;
float xc = 0.0f, xr = 0.0f, xi = 0.0f, contr = 0.0f;
int iter = 0;
bool converged = false;
for (; iter < SIFT_MAX_INTERP_STEPS; iter++)
{
float cur_v = READ_DOG(gauss_pack, gauss_step, gauss_offset, cc, rc, lc, dog_rows_per_layer);
float prv_c = READ_DOG(gauss_pack, gauss_step, gauss_offset, cc, rc, lc-1, dog_rows_per_layer);
float nxt_c = READ_DOG(gauss_pack, gauss_step, gauss_offset, cc, rc, lc+1, dog_rows_per_layer);
float dD0 = (READ_DOG(gauss_pack, gauss_step, gauss_offset, cc+1, rc, lc, dog_rows_per_layer) -
READ_DOG(gauss_pack, gauss_step, gauss_offset, cc-1, rc, lc, dog_rows_per_layer)) * deriv_scale;
float dD1 = (READ_DOG(gauss_pack, gauss_step, gauss_offset, cc, rc+1, lc, dog_rows_per_layer) -
READ_DOG(gauss_pack, gauss_step, gauss_offset, cc, rc-1, lc, dog_rows_per_layer)) * deriv_scale;
float dD2 = (nxt_c - prv_c) * deriv_scale;
float v2 = cur_v * 2.0f;
float dxx = (READ_DOG(gauss_pack, gauss_step, gauss_offset, cc+1, rc, lc, dog_rows_per_layer) +
READ_DOG(gauss_pack, gauss_step, gauss_offset, cc-1, rc, lc, dog_rows_per_layer) - v2) * second_deriv_scale;
float dyy = (READ_DOG(gauss_pack, gauss_step, gauss_offset, cc, rc+1, lc, dog_rows_per_layer) +
READ_DOG(gauss_pack, gauss_step, gauss_offset, cc, rc-1, lc, dog_rows_per_layer) - v2) * second_deriv_scale;
float dss = (nxt_c + prv_c - v2) * second_deriv_scale;
float dxy = (READ_DOG(gauss_pack, gauss_step, gauss_offset, cc+1, rc+1, lc, dog_rows_per_layer) -
READ_DOG(gauss_pack, gauss_step, gauss_offset, cc-1, rc+1, lc, dog_rows_per_layer) -
READ_DOG(gauss_pack, gauss_step, gauss_offset, cc+1, rc-1, lc, dog_rows_per_layer) +
READ_DOG(gauss_pack, gauss_step, gauss_offset, cc-1, rc-1, lc, dog_rows_per_layer)) * cross_deriv_scale;
float dxs = (READ_DOG(gauss_pack, gauss_step, gauss_offset, cc+1, rc, lc+1, dog_rows_per_layer) -
READ_DOG(gauss_pack, gauss_step, gauss_offset, cc-1, rc, lc+1, dog_rows_per_layer) -
READ_DOG(gauss_pack, gauss_step, gauss_offset, cc+1, rc, lc-1, dog_rows_per_layer) +
READ_DOG(gauss_pack, gauss_step, gauss_offset, cc-1, rc, lc-1, dog_rows_per_layer)) * cross_deriv_scale;
float dys = (READ_DOG(gauss_pack, gauss_step, gauss_offset, cc, rc+1, lc+1, dog_rows_per_layer) -
READ_DOG(gauss_pack, gauss_step, gauss_offset, cc, rc-1, lc+1, dog_rows_per_layer) -
READ_DOG(gauss_pack, gauss_step, gauss_offset, cc, rc+1, lc-1, dog_rows_per_layer) +
READ_DOG(gauss_pack, gauss_step, gauss_offset, cc, rc-1, lc-1, dog_rows_per_layer)) * cross_deriv_scale;
float detH = dxx*(dyy*dss - dys*dys) - dxy*(dxy*dss - dys*dxs) + dxs*(dxy*dys - dyy*dxs);
if (fabs(detH) < SIFT_FLT_EPS)
{
xc = 0.0f; xr = 0.0f; xi = 0.0f;
converged = true;
contr = cur_v * img_scale;
break;
}
float inv_det = 1.0f / detH;
xc = -(dD0*(dyy*dss - dys*dys) - dxy*(dD1*dss - dys*dD2) + dxs*(dD1*dys - dyy*dD2)) * inv_det;
xr = -(dxx*(dD1*dss - dys*dD2) - dD0*(dxy*dss - dys*dxs) + dxs*(dxy*dD2 - dD1*dxs)) * inv_det;
xi = -(dxx*(dyy*dD2 - dD1*dys) - dxy*(dxy*dD2 - dD1*dxs) + dD0*(dxy*dys - dyy*dxs)) * inv_det;
if (fabs(xi) < 0.5f && fabs(xr) < 0.5f && fabs(xc) < 0.5f)
{
converged = true;
contr = cur_v * img_scale + (dD0*xc + dD1*xr + dD2*xi) * 0.5f;
break;
}
if (fabs(xc) > (float)(INT_MAX/3) || fabs(xr) > (float)(INT_MAX/3) || fabs(xi) > (float)(INT_MAX/3))
return;
cc += convert_int_rte(xc);
rc += convert_int_rte(xr);
lc += convert_int_rte(xi);
if (lc < 1 || lc > nOctaveLayers ||
cc < SIFT_IMG_BORDER || cc >= dog_cols - SIFT_IMG_BORDER ||
rc < SIFT_IMG_BORDER || rc >= dog_rows_per_layer - SIFT_IMG_BORDER)
return;
}
if (!converged)
return;
if (fabs(contr) * (float)nOctaveLayers < contrastThreshold)
return;
float cur_v = READ_DOG(gauss_pack, gauss_step, gauss_offset, cc, rc, lc, dog_rows_per_layer);
float v2 = cur_v * 2.0f;
float dxx = (READ_DOG(gauss_pack, gauss_step, gauss_offset, cc+1, rc, lc, dog_rows_per_layer) +
READ_DOG(gauss_pack, gauss_step, gauss_offset, cc-1, rc, lc, dog_rows_per_layer) - v2) * second_deriv_scale;
float dyy = (READ_DOG(gauss_pack, gauss_step, gauss_offset, cc, rc+1, lc, dog_rows_per_layer) +
READ_DOG(gauss_pack, gauss_step, gauss_offset, cc, rc-1, lc, dog_rows_per_layer) - v2) * second_deriv_scale;
float dxy = (READ_DOG(gauss_pack, gauss_step, gauss_offset, cc+1, rc+1, lc, dog_rows_per_layer) -
READ_DOG(gauss_pack, gauss_step, gauss_offset, cc-1, rc+1, lc, dog_rows_per_layer) -
READ_DOG(gauss_pack, gauss_step, gauss_offset, cc+1, rc-1, lc, dog_rows_per_layer) +
READ_DOG(gauss_pack, gauss_step, gauss_offset, cc-1, rc-1, lc, dog_rows_per_layer)) * cross_deriv_scale;
float tr = dxx + dyy;
float det = dxx * dyy - dxy * dxy;
if (det <= 0.0f || tr*tr*edgeThreshold >= (edgeThreshold + 1.0f)*(edgeThreshold + 1.0f)*det)
return;
float scl_octv = sigma * exp2(((float)lc + xi) / (float)nOctaveLayers);
int radius = convert_int_rte(SIFT_ORI_RADIUS * scl_octv);
float ori_sig = SIFT_ORI_SIG_FCTR * scl_octv;
float expf_scale = -1.0f / (2.0f * ori_sig * ori_sig);
float hist[SIFT_ORI_HIST_BINS];
for (int b = 0; b < SIFT_ORI_HIST_BINS; b++)
hist[b] = 0.0f;
int gauss_row_base = lc * dog_rows_per_layer;
for (int i = -radius; i <= radius; i++)
{
int gy = rc + i;
if (gy <= 0 || gy >= dog_rows_per_layer - 1)
continue;
int grow = gauss_row_base + gy;
for (int j = -radius; j <= radius; j++)
{
int gx = cc + j;
if (gx <= 0 || gx >= dog_cols - 1)
continue;
float dx = READ_F32(gauss_pack, gauss_step, gauss_offset, gx+1, grow) -
READ_F32(gauss_pack, gauss_step, gauss_offset, gx-1, grow);
float dy = READ_F32(gauss_pack, gauss_step, gauss_offset, gx, grow-1) -
READ_F32(gauss_pack, gauss_step, gauss_offset, gx, grow+1);
float ang = atan2(dy, dx) * SIFT_RAD2DEG;
if (ang < 0.0f) ang += 360.0f;
float mag = native_sqrt(dx*dx + dy*dy);
float w = native_exp((float)(i*i + j*j) * expf_scale);
int bin = convert_int_rte((float)SIFT_ORI_HIST_BINS / 360.0f * ang);
if (bin >= SIFT_ORI_HIST_BINS) bin -= SIFT_ORI_HIST_BINS;
if (bin < 0) bin += SIFT_ORI_HIST_BINS;
hist[bin] += w * mag;
}
}
float sm[SIFT_ORI_HIST_BINS];
for (int b = 0; b < SIFT_ORI_HIST_BINS; b++)
{
int bm2 = (b + SIFT_ORI_HIST_BINS - 2) % SIFT_ORI_HIST_BINS;
int bm1 = (b + SIFT_ORI_HIST_BINS - 1) % SIFT_ORI_HIST_BINS;
int bp1 = (b + 1) % SIFT_ORI_HIST_BINS;
int bp2 = (b + 2) % SIFT_ORI_HIST_BINS;
sm[b] = (hist[bm2] + hist[bp2]) * (1.0f/16.0f) +
(hist[bm1] + hist[bp1]) * (4.0f/16.0f) +
hist[b] * (6.0f/16.0f);
}
float maxval = sm[0];
for (int b = 1; b < SIFT_ORI_HIST_BINS; b++)
maxval = fmax(maxval, sm[b]);
float mag_thr = maxval * SIFT_ORI_PEAK_RATIO;
float kpt_x = (cc + xc) * (float)(1 << octave);
float kpt_y = (rc + xr) * (float)(1 << octave);
float kpt_size = scl_octv * (float)(1 << octave) * 2.0f;
float kpt_resp = fabs(contr);
int kpt_oct = octave + (lc << 8) + (convert_int_rte((xi + 0.5f) * 255.0f) << 16);
for (int b = 0; b < SIFT_ORI_HIST_BINS; b++)
{
int bm1 = (b + SIFT_ORI_HIST_BINS - 1) % SIFT_ORI_HIST_BINS;
int bp1 = (b + 1) % SIFT_ORI_HIST_BINS;
if (sm[b] > sm[bm1] && sm[b] > sm[bp1] && sm[b] >= mag_thr)
{
float denom = sm[bm1] - 2.0f*sm[b] + sm[bp1];
if (fabs(denom) < SIFT_FLT_EPS)
continue;
float bin = (float)b + 0.5f * (sm[bm1] - sm[bp1]) / denom;
if (bin < 0.0f) bin += (float)SIFT_ORI_HIST_BINS;
if (bin >= (float)SIFT_ORI_HIST_BINS) bin -= (float)SIFT_ORI_HIST_BINS;
float angle = 360.0f - (360.0f / (float)SIFT_ORI_HIST_BINS) * bin;
if (fabs(angle - 360.0f) < SIFT_FLT_EPS)
angle = 0.0f;
int out_idx = atomic_inc(output_count);
if (out_idx >= max_output)
return;
int base = out_idx * 6;
output_kpts[base + 0] = kpt_x;
output_kpts[base + 1] = kpt_y;
output_kpts[base + 2] = angle;
output_kpts[base + 3] = kpt_size;
output_kpts[base + 4] = kpt_resp;
output_kpts[base + 5] = as_float(kpt_oct);
}
}
}
modules/features/src/sift.dispatch.cpp
+466
View File
@@ -76,6 +76,7 @@
#include "sift.simd.hpp"
#include "sift.simd_declarations.hpp" // defines CV_CPU_DISPATCH_MODES_ALL=AVX2,...,BASELINE based on CMakeLists.txt content
#include "opencl_kernels_features.hpp"
namespace cv {
@@ -131,6 +132,13 @@ public:
void setSigma(double sigma_) CV_OVERRIDE { sigma = sigma_; }
double getSigma() const CV_OVERRIDE { return sigma; }
private:
#ifdef HAVE_OPENCL
bool ocl_detectAndCompute(InputArray _image, InputArray _mask,
std::vector<KeyPoint>& keypoints,
OutputArray _descriptors, bool useProvidedKeypoints);
#endif
protected:
CV_PROP_RW int nfeatures;
CV_PROP_RW int nOctaveLayers;
@@ -497,6 +505,450 @@ int SIFT_Impl::defaultNorm() const
return NORM_L2;
}
#ifdef HAVE_OPENCL
namespace {
static void computeGaussianKernel1D(double sigma, std::vector<float>& coeffs, int& radius)
{
radius = (int)std::ceil(sigma * 3.0);
coeffs.resize(2 * radius + 1);
double sum = 0.0;
for (int i = -radius; i <= radius; i++)
{
double v = std::exp(-0.5 * (double)(i*i) / (sigma * sigma));
coeffs[i + radius] = (float)v;
sum += v;
}
for (int i = 0; i < (int)coeffs.size(); i++)
coeffs[i] = (float)(coeffs[i] / sum);
}
static bool ocl_gaussianBlurSep(const UMat& src, UMat& dst, double sigma)
{
std::vector<float> coeffs;
int radius;
computeGaussianKernel1D(sigma, coeffs, radius);
UMat uc = Mat(coeffs, true).getUMat(ACCESS_READ);
dst.create(src.size(), src.type());
UMat tmp(src.size(), src.type());
size_t localSize[2] = {16, 16};
size_t globalSize[2] = {
(size_t)((src.cols + 15) / 16) * 16,
(size_t)((src.rows + 15) / 16) * 16
};
ocl::Kernel kerH("SIFT_gaussian_blur_h", ocl::features::sift_oclsrc);
if (kerH.empty())
return false;
if (!kerH.args(
ocl::KernelArg::ReadOnlyNoSize(src),
ocl::KernelArg::WriteOnlyNoSize(tmp),
src.cols, src.rows,
ocl::KernelArg::PtrReadOnly(uc), radius
).run(2, globalSize, localSize, false))
return false;
ocl::Kernel kerV("SIFT_gaussian_blur_v", ocl::features::sift_oclsrc);
if (kerV.empty())
return false;
return kerV.args(
ocl::KernelArg::ReadOnlyNoSize(tmp),
ocl::KernelArg::WriteOnlyNoSize(dst),
src.cols, src.rows,
ocl::KernelArg::PtrReadOnly(uc), radius
).run(2, globalSize, localSize, false);
}
static UMat ocl_createInitialImage(const UMat& img, bool doubleImageSize, float sigma, bool enable_precise_upscale)
{
UMat gray, gray_fpt;
if (img.channels() == 3 || img.channels() == 4)
{
cvtColor(img, gray, COLOR_BGR2GRAY);
gray.convertTo(gray_fpt, CV_32F, SIFT_FIXPT_SCALE, 0);
}
else
img.convertTo(gray_fpt, CV_32F, SIFT_FIXPT_SCALE, 0);
float sig_diff;
if (doubleImageSize)
{
sig_diff = sqrtf(std::max(sigma * sigma - SIFT_INIT_SIGMA * SIFT_INIT_SIGMA * 4, 0.01f));
UMat dbl;
if (enable_precise_upscale) {
dbl.create(Size(gray_fpt.cols*2, gray_fpt.rows*2), gray_fpt.type());
Mat H = Mat::zeros(2, 3, CV_32F);
H.at<float>(0, 0) = 0.5f;
H.at<float>(1, 1) = 0.5f;
warpAffine(gray_fpt, dbl, H, dbl.size(), INTER_LINEAR | WARP_INVERSE_MAP, BORDER_REFLECT);
} else {
resize(gray_fpt, dbl, Size(gray_fpt.cols*2, gray_fpt.rows*2), 0, 0, INTER_LINEAR);
}
UMat result;
if (!ocl_gaussianBlurSep(dbl, result, (double)sig_diff))
return UMat();
return result;
}
else
{
sig_diff = sqrtf(std::max(sigma * sigma - SIFT_INIT_SIGMA * SIFT_INIT_SIGMA, 0.01f));
UMat result;
if (!ocl_gaussianBlurSep(gray_fpt, result, (double)sig_diff))
return UMat();
return result;
}
}
static bool ocl_buildGaussianPyramid(const UMat& base, std::vector<UMat>& gauss_packs,
int nOctaves, int nOctaveLayers, double sigma)
{
const int levelsPerOctave = nOctaveLayers + 3;
std::vector<double> sig(nOctaveLayers + 3);
gauss_packs.resize(nOctaves);
sig[0] = sigma;
double k = std::pow(2., 1. / nOctaveLayers);
for (int i = 1; i < nOctaveLayers + 3; i++)
{
double sig_prev = std::pow(k, i-1)*sigma;
double sig_total = sig_prev*k;
sig[i] = std::sqrt(sig_total*sig_total - sig_prev*sig_prev);
}
int n_kernels = nOctaveLayers + 2;
std::vector<std::vector<float> > gk_coeffs(n_kernels);
std::vector<int> gk_radius(n_kernels);
for (int i = 0; i < n_kernels; i++)
computeGaussianKernel1D(sig[i + 1], gk_coeffs[i], gk_radius[i]);
std::vector<UMat> uc(n_kernels);
for (int i = 0; i < n_kernels; i++)
uc[i] = Mat(gk_coeffs[i], true).getUMat(ACCESS_READ);
for (int o = 0; o < nOctaves; o++)
{
int levelW = base.cols >> o;
int levelH = base.rows >> o;
gauss_packs[o].create(levelsPerOctave * levelH, levelW, base.type());
}
for (int o = 0; o < nOctaves; o++)
{
int img_cols = gauss_packs[o].cols;
int img_rows = gauss_packs[o].rows / levelsPerOctave;
UMat level0 = gauss_packs[o].rowRange(0, img_rows);
if (o == 0)
{
base.copyTo(level0);
}
else
{
int prev_rows = gauss_packs[o-1].rows / levelsPerOctave;
UMat src = gauss_packs[o-1].rowRange(nOctaveLayers * prev_rows,
(nOctaveLayers + 1) * prev_rows);
resize(src, level0, Size(img_cols, img_rows), 0, 0, INTER_NEAREST);
}
UMat tmp(Size(img_cols, img_rows), base.type());
UMat prev = level0;
size_t localSize[2] = {16, 16};
for (int i = 1; i < levelsPerOctave; i++)
{
int lev_idx = i - 1;
UMat dst = gauss_packs[o].rowRange(i * img_rows, (i + 1) * img_rows);
size_t globalSize[2] = {
(size_t)((img_cols + 15) / 16) * 16,
(size_t)((img_rows + 15) / 16) * 16
};
ocl::Kernel kerH("SIFT_gaussian_blur_h", ocl::features::sift_oclsrc);
if (kerH.empty())
return false;
if (!kerH.args(
ocl::KernelArg::ReadOnlyNoSize(prev),
ocl::KernelArg::WriteOnlyNoSize(tmp),
img_cols, img_rows,
ocl::KernelArg::PtrReadOnly(uc[lev_idx]), gk_radius[lev_idx]
).run(2, globalSize, localSize, false))
return false;
ocl::Kernel kerV("SIFT_gaussian_blur_v", ocl::features::sift_oclsrc);
if (kerV.empty())
return false;
if (!kerV.args(
ocl::KernelArg::ReadOnlyNoSize(tmp),
ocl::KernelArg::WriteOnlyNoSize(dst),
img_cols, img_rows,
ocl::KernelArg::PtrReadOnly(uc[lev_idx]), gk_radius[lev_idx]
).run(2, globalSize, localSize, false))
return false;
prev = dst;
}
}
return true;
}
static bool ocl_detectAndOrient(
const std::vector<UMat>& gauss_packs,
std::vector<KeyPoint>& keypoints,
int nOctaves, int nOctaveLayers,
double contrastThreshold, double edgeThreshold, double sigma)
{
const int threshold = cvFloor(0.5 * contrastThreshold / nOctaveLayers * 255 * SIFT_FIXPT_SCALE);
const int max_kpts_per_call = 100000;
UMat count_umat(1, 1, CV_32S, Scalar::all(0));
UMat kp_buf(1, max_kpts_per_call * 6, CV_32F);
for (int o = 0; o < nOctaves; o++)
{
int dog_cols = gauss_packs[o].cols;
int dog_rows = gauss_packs[o].rows / (nOctaveLayers + 3);
for (int i = 1; i <= nOctaveLayers; i++)
{
size_t globalSize[2] = {(size_t)dog_cols, (size_t)dog_rows};
ocl::Kernel ker("SIFT_detect_and_orient", ocl::features::sift_oclsrc);
if (ker.empty())
return false;
bool ok = ker.args(
ocl::KernelArg::ReadOnlyNoSize(gauss_packs[o]),
dog_cols, dog_rows,
threshold,
(float)contrastThreshold, (float)edgeThreshold, (float)sigma,
nOctaveLayers, o, i,
ocl::KernelArg::PtrWriteOnly(count_umat),
ocl::KernelArg::PtrWriteOnly(kp_buf),
max_kpts_per_call
).run(2, globalSize, 0, false);
if (!ok)
return false;
}
}
ocl::finish();
Mat count_mat = count_umat.getMat(ACCESS_READ);
int total = count_mat.at<int>(0, 0);
if (total > max_kpts_per_call)
total = max_kpts_per_call;
Mat kp = kp_buf.getMat(ACCESS_READ);
const float* kdata = (const float*)kp.data;
keypoints.resize(total);
for (int i = 0; i < total; i++)
{
KeyPoint& kpt = keypoints[i];
int base = i * 6;
kpt.pt.x = kdata[base + 0];
kpt.pt.y = kdata[base + 1];
kpt.angle = kdata[base + 2];
kpt.size = kdata[base + 3];
kpt.response = kdata[base + 4];
memcpy(&kpt.octave, &kdata[base + 5], sizeof(int));
kpt.class_id = -1;
}
return true;
}
static bool ocl_computeSIFTDescriptors(
const std::vector<UMat>& gauss_packs,
const std::vector<KeyPoint>& keypoints,
OutputArray _descriptors,
int nOctaveLayers, int firstOctave, int descriptor_type)
{
const int dsize = SIFT_DESCR_WIDTH*SIFT_DESCR_WIDTH*SIFT_DESCR_HIST_BINS;
const int nkp = (int)keypoints.size();
if (nkp == 0)
{
_descriptors.create(0, dsize, descriptor_type);
return true;
}
const int levelsPerOctave = nOctaveLayers + 3;
const int nOctaves = (int)gauss_packs.size();
struct KptGroup { int start; int count; };
std::vector<KptGroup> groups(nOctaves, {0, 0});
for (int i = 0; i < nkp; i++)
{
int octave, layer; float scale;
unpackOctave(keypoints[i], octave, layer, scale);
int oct_idx = octave - firstOctave;
if (oct_idx < 0 || oct_idx >= nOctaves)
return false;
groups[oct_idx].count++;
}
int offset = 0;
for (int o = 0; o < nOctaves; o++)
{
groups[o].start = offset;
offset += groups[o].count;
}
Mat all_kpts(nkp, 5, CV_32F);
Mat all_rows(nkp, 1, CV_32S);
std::vector<int> cursors(nOctaves, 0);
for (int i = 0; i < nkp; i++)
{
int octave, layer; float scale;
unpackOctave(keypoints[i], octave, layer, scale);
int oct_idx = octave - firstOctave;
float size = keypoints[i].size * scale;
Point2f ptf(keypoints[i].pt.x * scale, keypoints[i].pt.y * scale);
float angle = 360.f - keypoints[i].angle;
if (std::abs(angle - 360.f) < FLT_EPSILON)
angle = 0.f;
float scl = size * 0.5f;
int pos = groups[oct_idx].start + cursors[oct_idx]++;
float* dst = all_kpts.ptr<float>(pos);
dst[0] = ptf.x; dst[1] = ptf.y; dst[2] = angle; dst[3] = scl;
dst[4] = (float)layer;
all_rows.ptr<int>(pos)[0] = i;
}
UMat kpts_buf; all_kpts.copyTo(kpts_buf);
UMat rows_buf; all_rows.copyTo(rows_buf);
_descriptors.create(nkp, dsize, descriptor_type);
UMat desc_out = _descriptors.getUMat();
for (int o = 0; o < nOctaves; o++)
{
int nGroup = groups[o].count;
if (nGroup == 0)
continue;
const UMat& packed = gauss_packs[o];
int img_cols = packed.cols;
int img_rows = packed.rows / levelsPerOctave;
float diag = sqrt((float)img_cols * img_cols + (float)img_rows * img_rows);
UMat kpts_roi = kpts_buf.rowRange(groups[o].start,
groups[o].start + nGroup);
ocl::Kernel ker("SIFT_compute_descriptor", ocl::features::sift_oclsrc);
if (ker.empty())
return false;
const int KP_PER_WG = 4;
const int WG_SIZE = KP_PER_WG * 16;
size_t localsz = WG_SIZE;
size_t globalsz = (size_t)(((nGroup + KP_PER_WG - 1) / KP_PER_WG) * WG_SIZE);
bool ok = ker.args(
ocl::KernelArg::ReadOnlyNoSize(packed), img_cols, img_rows, diag,
ocl::KernelArg::ReadOnlyNoSize(kpts_roi),
ocl::KernelArg::PtrReadOnly(rows_buf),
groups[o].start,
nGroup,
ocl::KernelArg::WriteOnlyNoSize(desc_out),
descriptor_type
).run(1, &globalsz, &localsz, false);
if (!ok)
return false;
}
return true;
}
} // namespace
bool SIFT_Impl::ocl_detectAndCompute(InputArray _image, InputArray _mask,
std::vector<KeyPoint>& keypoints,
OutputArray _descriptors, bool useProvidedKeypoints)
{
UMat image = _image.getUMat();
if (image.empty() || image.depth() != CV_8U)
return false;
int firstOctave = -1;
int nOctaves = 0;
if (useProvidedKeypoints)
{
firstOctave = 0;
int maxOctave = INT_MIN;
int actualNLayers = 0;
for (size_t i = 0; i < keypoints.size(); i++)
{
int octave, layer; float scale;
unpackOctave(keypoints[i], octave, layer, scale);
firstOctave = std::min(firstOctave, octave);
maxOctave = std::max(maxOctave, octave);
actualNLayers = std::max(actualNLayers, layer-2);
}
firstOctave = std::min(firstOctave, 0);
if (firstOctave < -1 || actualNLayers > nOctaveLayers)
return false;
nOctaves = maxOctave - firstOctave + 1;
}
UMat base = ocl_createInitialImage(image, firstOctave < 0, (float)sigma, enable_precise_upscale);
if (base.empty())
return false;
if (!useProvidedKeypoints)
nOctaves = cvRound(std::log((double)std::min(base.cols, base.rows)) / std::log(2.) - 2) - firstOctave;
if (nOctaves <= 0)
return false;
std::vector<UMat> gauss_packs;
if (!ocl_buildGaussianPyramid(base, gauss_packs, nOctaves, nOctaveLayers, sigma))
return false;
if (!useProvidedKeypoints)
{
if (!ocl_detectAndOrient(gauss_packs, keypoints, nOctaves, nOctaveLayers,
contrastThreshold, edgeThreshold, sigma))
return false;
KeyPointsFilter::removeDuplicatedSorted(keypoints);
if (nfeatures > 0)
KeyPointsFilter::retainBest(keypoints, nfeatures);
const float scale = 1.f/(float)(1 << -firstOctave);
for (size_t i = 0; i < keypoints.size(); i++)
{
KeyPoint& kpt = keypoints[i];
kpt.octave = (kpt.octave & ~255) | ((kpt.octave + firstOctave) & 255);
kpt.pt *= scale;
kpt.size *= scale;
}
if (!_mask.empty())
{
Mat mask = _mask.getMat();
KeyPointsFilter::runByPixelsMask(keypoints, mask);
}
if (_descriptors.needed())
{
if (!ocl_computeSIFTDescriptors(gauss_packs, keypoints, _descriptors, nOctaveLayers, firstOctave, descriptor_type))
return false;
}
}
else
{
if (_descriptors.needed())
{
if (!ocl_computeSIFTDescriptors(gauss_packs, keypoints, _descriptors, nOctaveLayers, firstOctave, descriptor_type))
return false;
}
}
return true;
}
#endif // HAVE_OPENCL
void SIFT_Impl::detectAndCompute(InputArray _image, InputArray _mask,
std::vector<KeyPoint>& keypoints,
@@ -505,6 +957,20 @@ void SIFT_Impl::detectAndCompute(InputArray _image, InputArray _mask,
{
CV_TRACE_FUNCTION();
#ifdef HAVE_OPENCL
if (ocl::isOpenCLActivated() && _image.isUMat() && sizeof(void*) > 4)
{
ocl::Device d = ocl::Device::getDefault();
bool isDiscreteGPU = (d.type() == ocl::Device::TYPE_GPU) && !d.hostUnifiedMemory();
bool ocl11OrLater = (d.deviceVersionMajor() > 1) || (d.deviceVersionMajor() == 1 && d.deviceVersionMinor() >= 1);
if (isDiscreteGPU && ocl11OrLater && ocl_detectAndCompute(_image, _mask, keypoints, _descriptors, useProvidedKeypoints))
{
CV_IMPL_ADD(CV_IMPL_OCL);
return;
}
}
#endif
int firstOctave = -1, actualNOctaves = 0, actualNLayers = 0;
Mat image = _image.getMat(), mask = _mask.getMat();
modules/features/test/ocl/test_feature2d.cpp
+191
View File
@@ -0,0 +1,191 @@
// This file is part of OpenCV project.
// It is subject to the license terms in the LICENSE file found in the top-level directory
// of this distribution and at http://opencv.org/license.html
#include "../test_precomp.hpp"
#include "opencv2/ts/ocl_test.hpp"
#include <functional>
#ifdef HAVE_OPENCL
namespace opencv_test {
namespace ocl {
#define TEST_IMAGES testing::Values(\
"detectors_descriptors_evaluation/images_datasets/leuven/img1.png",\
"../stitching/a3.png", \
"../stitching/s2.jpg")
PARAM_TEST_CASE(Feature2DFixture, std::function<Ptr<Feature2D>()>, std::string, double)
{
std::string filename;
double desc_eps;
Mat image, descriptors;
vector<KeyPoint> keypoints;
UMat uimage, udescriptors;
vector<KeyPoint> ukeypoints;
Ptr<Feature2D> feature;
virtual void SetUp()
{
feature = GET_PARAM(0)();
filename = GET_PARAM(1);
desc_eps = GET_PARAM(2);
image = readImage(filename);
ASSERT_FALSE(image.empty());
image.copyTo(uimage);
OCL_OFF(feature->detect(image, keypoints));
OCL_ON(feature->detect(uimage, ukeypoints));
OCL_OFF(feature->compute(image, keypoints, descriptors));
OCL_ON(feature->compute(uimage, keypoints, udescriptors));
}
};
OCL_TEST_P(Feature2DFixture, KeypointsSame)
{
size_t count_diff = (keypoints.size() > ukeypoints.size()) ?
keypoints.size() - ukeypoints.size() : ukeypoints.size() - keypoints.size();
EXPECT_LE(count_diff, (size_t)(std::min(keypoints.size(), ukeypoints.size()) * 0.20 + 1));
std::vector<KeyPoint> cpu_sorted = keypoints, ocl_sorted = ukeypoints;
std::sort(cpu_sorted.begin(), cpu_sorted.end(), [](const KeyPoint& a, const KeyPoint& b) {
return a.pt.x < b.pt.x || (a.pt.x == b.pt.x && a.pt.y < b.pt.y);
});
std::sort(ocl_sorted.begin(), ocl_sorted.end(), [](const KeyPoint& a, const KeyPoint& b) {
return a.pt.x < b.pt.x || (a.pt.x == b.pt.x && a.pt.y < b.pt.y);
});
int matched = 0;
size_t j = 0;
for (size_t i = 0; i < cpu_sorted.size(); i++)
{
while (j < ocl_sorted.size() && ocl_sorted[j].pt.x < cpu_sorted[i].pt.x - 2.0f)
j++;
for (size_t k = j; k < ocl_sorted.size() && ocl_sorted[k].pt.x <= cpu_sorted[i].pt.x + 2.0f; k++)
{
if (std::abs(ocl_sorted[k].pt.y - cpu_sorted[i].pt.y) < 2.0f)
{
matched++;
break;
}
}
}
size_t n = std::min(cpu_sorted.size(), ocl_sorted.size());
EXPECT_GE(matched, (int)(n * 0.70));
}
OCL_TEST_P(Feature2DFixture, DescriptorsSame)
{
EXPECT_EQ(descriptors.size(), udescriptors.size());
ASSERT_EQ(descriptors.type(), udescriptors.type());
Mat udesc = udescriptors.getMat(ACCESS_READ);
double max_diff = 0.0;
for (int r = 0; r < descriptors.rows; r++)
{
for (int c = 0; c < descriptors.cols; c++)
{
double d;
if (descriptors.type() == CV_8U)
d = std::abs((double)descriptors.at<uchar>(r, c) - (double)udesc.at<uchar>(r, c));
else
d = std::abs((double)descriptors.at<float>(r, c) - (double)udesc.at<float>(r, c));
max_diff = std::max(max_diff, d);
}
}
EXPECT_LE(max_diff, desc_eps);
}
OCL_INSTANTIATE_TEST_CASE_P(SIFT, Feature2DFixture,
testing::Combine(testing::Values([]() { return SIFT::create(); }), TEST_IMAGES, testing::Values(2.0)));
OCL_TEST(SIFT, NonDefaultNOctaveLayers)
{
Mat image = TestUtils::readImage("../stitching/s2.jpg", IMREAD_GRAYSCALE);
ASSERT_FALSE(image.empty());
UMat uimage;
image.copyTo(uimage);
Ptr<SIFT> sift = SIFT::create(0, 4);
vector<KeyPoint> keypoints;
UMat descriptors;
ASSERT_NO_THROW(sift->detectAndCompute(uimage, noArray(), keypoints, descriptors, false));
EXPECT_GT(keypoints.size(), 20u);
EXPECT_EQ((size_t)descriptors.rows, keypoints.size());
}
OCL_TEST(SIFT, DescriptorType)
{
Mat image = imread(cvtest::findDataFile("features2d/tsukuba.png"), IMREAD_GRAYSCALE);
ASSERT_FALSE(image.empty());
UMat uimage;
image.copyTo(uimage);
vector<KeyPoint> keypoints;
UMat descriptorsFloat, descriptorsUchar;
Ptr<SIFT> siftFloat = SIFT::create(0, 3, 0.04, 10, 1.6, CV_32F);
siftFloat->detectAndCompute(uimage, noArray(), keypoints, descriptorsFloat, false);
ASSERT_EQ(descriptorsFloat.type(), CV_32F) << "type mismatch";
Ptr<SIFT> siftUchar = SIFT::create(0, 3, 0.04, 10, 1.6, CV_8U);
siftUchar->detectAndCompute(uimage, noArray(), keypoints, descriptorsUchar, false);
ASSERT_EQ(descriptorsUchar.type(), CV_8U) << "type mismatch";
Mat df = descriptorsFloat.getMat(ACCESS_READ);
Mat du = descriptorsUchar.getMat(ACCESS_READ);
Mat descriptorsFloat2;
du.assignTo(descriptorsFloat2, CV_32F);
Mat diff = df != descriptorsFloat2;
EXPECT_EQ(countNonZero(diff), 0) << "descriptors are not identical";
}
OCL_TEST(SIFT, Regression_26139)
{
UMat uimage(Size(300, 300), CV_8UC1, Scalar::all(0));
std::vector<KeyPoint> kps {
KeyPoint(154.076813f, 136.160904f, 111.078636f, 216.195618f, 0.00000899323549f, 7)
};
Ptr<SIFT> extractor = SIFT::create();
UMat descriptors;
extractor->compute(uimage, kps, descriptors);
ASSERT_EQ(descriptors.size(), Size(128, 1));
}
OCL_TEST(SIFT, Batch)
{
string path = cvtest::TS::ptr()->get_data_path() + "detectors_descriptors_evaluation/images_datasets/graf";
vector<UMat> imgs;
vector<UMat> descriptors;
vector<vector<KeyPoint> > keypoints;
int n = 6;
Ptr<SIFT> sift = SIFT::create();
for (int i = 0; i < n; i++)
{
string imgname = format("%s/img%d.png", path.c_str(), i + 1);
Mat img = imread(imgname, IMREAD_GRAYSCALE);
ASSERT_FALSE(img.empty()) << "Failed to load " << imgname;
UMat uimg;
img.copyTo(uimg);
imgs.push_back(uimg);
}
sift->detect(imgs, keypoints);
sift->compute(imgs, keypoints, descriptors);
ASSERT_EQ((int)keypoints.size(), n);
ASSERT_EQ((int)descriptors.size(), n);
for (int i = 0; i < n; i++)
{
EXPECT_GT((int)keypoints[i].size(), 100);
EXPECT_GT(descriptors[i].rows, 100);
}
}
}//ocl
}//opencv_test
#endif //HAVE_OPENCL