| 300 | // http://research.microsoft.com/en-us/um/people/zhang/INRIA/Publis/Tutorial-Estim/node25.html |
| 301 | template<typename T> |
| 302 | __global__ void computeEvalHomography( |
| 303 | Param<unsigned> inliers, Param<unsigned> idx, Param<T> H, Param<float> err, |
| 304 | CParam<float> x_src, CParam<float> y_src, CParam<float> x_dst, |
| 305 | CParam<float> y_dst, CParam<float> rnd, const unsigned iterations, |
| 306 | const unsigned nsamples, const float inlier_thr, |
| 307 | const af_homography_type htype) { |
| 308 | unsigned bid_x = blockIdx.x; |
| 309 | unsigned tid_x = threadIdx.x; |
| 310 | unsigned i = bid_x * blockDim.x + tid_x; |
| 311 | |
| 312 | __shared__ unsigned s_inliers[256]; |
| 313 | __shared__ unsigned s_idx[256]; |
| 314 | |
| 315 | s_inliers[tid_x] = 0; |
| 316 | s_idx[tid_x] = 0; |
| 317 | __syncthreads(); |
| 318 | |
| 319 | if (i < iterations) { |
| 320 | const unsigned Hidx = H.dims[0] * i; |
| 321 | T* H_ptr = H.ptr + Hidx; |
| 322 | T H_tmp[9]; |
| 323 | for (int h = 0; h < 9; h++) H_tmp[h] = H_ptr[h]; |
| 324 | |
| 325 | if (htype == AF_HOMOGRAPHY_RANSAC) { |
| 326 | // Compute inliers |
| 327 | unsigned inliers_count = 0; |
| 328 | for (unsigned j = 0; j < nsamples; j++) { |
| 329 | float z = H_tmp[6] * x_src.ptr[j] + H_tmp[7] * y_src.ptr[j] + |
| 330 | H_tmp[8]; |
| 331 | float x = (H_tmp[0] * x_src.ptr[j] + H_tmp[1] * y_src.ptr[j] + |
| 332 | H_tmp[2]) / |
| 333 | z; |
| 334 | float y = (H_tmp[3] * x_src.ptr[j] + H_tmp[4] * y_src.ptr[j] + |
| 335 | H_tmp[5]) / |
| 336 | z; |
| 337 | |
| 338 | float dist = sq(x_dst.ptr[j] - x) + sq(y_dst.ptr[j] - y); |
| 339 | if (dist < inlier_thr * inlier_thr) inliers_count++; |
| 340 | } |
| 341 | |
| 342 | s_inliers[tid_x] = inliers_count; |
| 343 | s_idx[tid_x] = i; |
| 344 | } else if (htype == AF_HOMOGRAPHY_LMEDS) { |
| 345 | // Compute error |
| 346 | for (unsigned j = 0; j < nsamples; j++) { |
| 347 | float z = H_tmp[6] * x_src.ptr[j] + H_tmp[7] * y_src.ptr[j] + |
| 348 | H_tmp[8]; |
| 349 | float x = (H_tmp[0] * x_src.ptr[j] + H_tmp[1] * y_src.ptr[j] + |
| 350 | H_tmp[2]) / |
| 351 | z; |
| 352 | float y = (H_tmp[3] * x_src.ptr[j] + H_tmp[4] * y_src.ptr[j] + |
| 353 | H_tmp[5]) / |
| 354 | z; |
| 355 | |
| 356 | float dist = sq(x_dst.ptr[j] - x) + sq(y_dst.ptr[j] - y); |
| 357 | err.ptr[i * err.dims[0] + j] = sqrt(dist); |
| 358 | } |
| 359 | } |