| 357 | // Based on Section 4 of Lowe's paper. |
| 358 | template<typename T> |
| 359 | __global__ void interpolateExtrema( |
| 360 | float* x_out, float* y_out, unsigned* layer_out, float* response_out, |
| 361 | float* size_out, unsigned* counter, const float* x_in, const float* y_in, |
| 362 | const unsigned* layer_in, const unsigned extrema_feat, |
| 363 | const CParam<T> dog_octave, const unsigned max_feat, const unsigned octave, |
| 364 | const unsigned n_layers, const float contrast_thr, const float edge_thr, |
| 365 | const float sigma, const float img_scale) { |
| 366 | const unsigned f = blockIdx.x * blockDim.x + threadIdx.x; |
| 367 | |
| 368 | if (f < extrema_feat) { |
| 369 | const float first_deriv_scale = img_scale * 0.5f; |
| 370 | const float second_deriv_scale = img_scale; |
| 371 | const float cross_deriv_scale = img_scale * 0.25f; |
| 372 | |
| 373 | float xl = 0, xy = 0, xx = 0, contr = 0; |
| 374 | int i = 0; |
| 375 | |
| 376 | unsigned x = x_in[f]; |
| 377 | unsigned y = y_in[f]; |
| 378 | unsigned layer = layer_in[f]; |
| 379 | |
| 380 | const int dim0 = dog_octave.dims[0]; |
| 381 | const int dim1 = dog_octave.dims[1]; |
| 382 | const int imel = dim0 * dim1; |
| 383 | |
| 384 | const T* prev_ptr = dog_octave.ptr + (layer - 1) * imel; |
| 385 | const T* center_ptr = dog_octave.ptr + (layer)*imel; |
| 386 | const T* next_ptr = dog_octave.ptr + (layer + 1) * imel; |
| 387 | |
| 388 | for (i = 0; i < MAX_INTERP_STEPS; i++) { |
| 389 | float dD[3] = { |
| 390 | (float)(CPTR(x + 1, y) - CPTR(x - 1, y)) * first_deriv_scale, |
| 391 | (float)(CPTR(x, y + 1) - CPTR(x, y - 1)) * first_deriv_scale, |
| 392 | (float)(NPTR(x, y) - PPTR(x, y)) * first_deriv_scale}; |
| 393 | |
| 394 | float d2 = CPTR(x, y) * 2.f; |
| 395 | float dxx = |
| 396 | (CPTR(x + 1, y) + CPTR(x - 1, y) - d2) * second_deriv_scale; |
| 397 | float dyy = |
| 398 | (CPTR(x, y + 1) + CPTR(x, y - 1) - d2) * second_deriv_scale; |
| 399 | float dss = (NPTR(x, y) + PPTR(x, y) - d2) * second_deriv_scale; |
| 400 | float dxy = (CPTR(x + 1, y + 1) - CPTR(x - 1, y + 1) - |
| 401 | CPTR(x + 1, y - 1) + CPTR(x - 1, y - 1)) * |
| 402 | cross_deriv_scale; |
| 403 | float dxs = (NPTR(x + 1, y) - NPTR(x - 1, y) - PPTR(x + 1, y) + |
| 404 | PPTR(x - 1, y)) * |
| 405 | cross_deriv_scale; |
| 406 | float dys = (NPTR(x, y + 1) - NPTR(x - 1, y - 1) - PPTR(x, y - 1) + |
| 407 | PPTR(x - 1, y - 1)) * |
| 408 | cross_deriv_scale; |
| 409 | |
| 410 | float H[9] = {dxx, dxy, dxs, dxy, dyy, dys, dxs, dys, dss}; |
| 411 | |
| 412 | float X[3]; |
| 413 | gaussianElimination<3>(H, dD, X); |
| 414 | |
| 415 | xl = -X[2]; |
| 416 | xy = -X[1]; |