diff --git a/modules/features/perf/opencl/perf_sift.cpp b/modules/features/perf/opencl/perf_sift.cpp new file mode 100644 index 0000000000..43fd733ffa --- /dev/null +++ b/modules/features/perf/opencl/perf_sift.cpp @@ -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 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::create(); + vector 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::create(); + vector 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 diff --git a/modules/features/src/opencl/sift.cl b/modules/features/src/opencl/sift.cl new file mode 100644 index 0000000000..5ab42c0992 --- /dev/null +++ b/modules/features/src/opencl/sift.cl @@ -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); + } + } +} diff --git a/modules/features/src/sift.dispatch.cpp b/modules/features/src/sift.dispatch.cpp index f7c6116093..6a45d0ec84 100644 --- a/modules/features/src/sift.dispatch.cpp +++ b/modules/features/src/sift.dispatch.cpp @@ -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& 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& 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 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(0, 0) = 0.5f; + H.at(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& gauss_packs, + int nOctaves, int nOctaveLayers, double sigma) +{ + const int levelsPerOctave = nOctaveLayers + 3; + std::vector 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 > gk_coeffs(n_kernels); + std::vector gk_radius(n_kernels); + for (int i = 0; i < n_kernels; i++) + computeGaussianKernel1D(sig[i + 1], gk_coeffs[i], gk_radius[i]); + + std::vector 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& gauss_packs, + std::vector& 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(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& gauss_packs, + const std::vector& 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 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 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(pos); + dst[0] = ptf.x; dst[1] = ptf.y; dst[2] = angle; dst[3] = scl; + dst[4] = (float)layer; + all_rows.ptr(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& 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 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& 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(); diff --git a/modules/features/test/ocl/test_feature2d.cpp b/modules/features/test/ocl/test_feature2d.cpp new file mode 100644 index 0000000000..ccfddd712b --- /dev/null +++ b/modules/features/test/ocl/test_feature2d.cpp @@ -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 + +#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()>, std::string, double) +{ + std::string filename; + double desc_eps; + Mat image, descriptors; + vector keypoints; + UMat uimage, udescriptors; + vector ukeypoints; + Ptr 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 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(r, c) - (double)udesc.at(r, c)); + else + d = std::abs((double)descriptors.at(r, c) - (double)udesc.at(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::create(0, 4); + vector 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 keypoints; + UMat descriptorsFloat, descriptorsUchar; + + Ptr 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 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 kps { + KeyPoint(154.076813f, 136.160904f, 111.078636f, 216.195618f, 0.00000899323549f, 7) + }; + Ptr 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 imgs; + vector descriptors; + vector > keypoints; + int n = 6; + Ptr 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