| 214 |
214 |
descriptionPairs[i] = allPairs[FREAK_DEF_PAIRS[i]];
|
| 215 |
215 |
}
|
| 216 |
216 |
}
|
| 217 |
|
|
|
217 |
|
| 218 |
218 |
void FREAK::computeImpl( const Mat& image, std::vector<KeyPoint>& keypoints, Mat& descriptors ) const {
|
| 219 |
|
|
|
219 |
|
| 220 |
220 |
if( image.empty() )
|
| 221 |
221 |
return;
|
| 222 |
222 |
if( keypoints.empty() )
|
| 223 |
223 |
return;
|
| 224 |
|
|
|
224 |
|
| 225 |
225 |
((FREAK*)this)->buildPattern();
|
|
226 |
|
|
227 |
// Use 32-bit integers if we won't overflow in the integral image
|
|
228 |
if ((image.type() == CV_8U || image.type() == CV_8S) &&
|
|
229 |
(image.rows * image.cols) < 8388608 ) // 8388608 = 2 ^ (32 - 8(bit depth) - 1(sign bit))
|
|
230 |
{
|
|
231 |
// Create the integral image appropriate for our type & usage
|
|
232 |
if (image.type() == CV_8U)
|
|
233 |
computeDescriptors<uchar, int>(image, keypoints, descriptors);
|
|
234 |
else if (image.type() == CV_8S)
|
|
235 |
computeDescriptors<char, int>(image, keypoints, descriptors);
|
|
236 |
else
|
|
237 |
CV_Error( CV_StsUnsupportedFormat, "" );
|
|
238 |
} else {
|
|
239 |
// Create the integral image appropriate for our type & usage
|
|
240 |
if ( image.type() == CV_8U )
|
|
241 |
computeDescriptors<uchar, double>(image, keypoints, descriptors);
|
|
242 |
else if ( image.type() == CV_8S )
|
|
243 |
computeDescriptors<char, double>(image, keypoints, descriptors);
|
|
244 |
else if ( image.type() == CV_16U )
|
|
245 |
computeDescriptors<ushort, double>(image, keypoints, descriptors);
|
|
246 |
else if ( image.type() == CV_16S )
|
|
247 |
computeDescriptors<short, double>(image, keypoints, descriptors);
|
|
248 |
else
|
|
249 |
CV_Error( CV_StsUnsupportedFormat, "" );
|
|
250 |
}
|
|
251 |
}
|
| 226 |
252 |
|
|
253 |
template <typename srcMatType, typename iiMatType>
|
|
254 |
void FREAK::computeDescriptors( const Mat& image, std::vector<KeyPoint>& keypoints, Mat& descriptors ) const {
|
|
255 |
|
| 227 |
256 |
Mat imgIntegral;
|
| 228 |
|
integral(image, imgIntegral);
|
|
257 |
if (sizeof(iiMatType) == sizeof(double))
|
|
258 |
integral(image, imgIntegral, CV_64F);
|
|
259 |
else
|
|
260 |
integral(image, imgIntegral, CV_32S);
|
|
261 |
|
| 229 |
262 |
std::vector<int> kpScaleIdx(keypoints.size()); // used to save pattern scale index corresponding to each keypoints
|
| 230 |
263 |
const std::vector<int>::iterator ScaleIdxBegin = kpScaleIdx.begin(); // used in std::vector erase function
|
| 231 |
264 |
const std::vector<cv::KeyPoint>::iterator kpBegin = keypoints.begin(); // used in std::vector erase function
|
| 232 |
265 |
const float sizeCst = static_cast<float>(FREAK_NB_SCALES/(FREAK_LOG2* nOctaves));
|
| 233 |
|
uchar pointsValue[FREAK_NB_POINTS];
|
|
266 |
srcMatType pointsValue[FREAK_NB_POINTS];
|
| 234 |
267 |
int thetaIdx = 0;
|
| 235 |
268 |
int direction0;
|
| 236 |
269 |
int direction1;
|
| ... | ... | |
| 276 |
309 |
// extract the best comparisons only
|
| 277 |
310 |
descriptors = cv::Mat::zeros((int)keypoints.size(), FREAK_NB_PAIRS/8, CV_8U);
|
| 278 |
311 |
#if CV_SSE2
|
| 279 |
|
__m128i* ptr= (__m128i*) (descriptors.data+(keypoints.size()-1)*descriptors.step[0]);
|
| 280 |
|
#else
|
| 281 |
|
std::bitset<FREAK_NB_PAIRS>* ptr = (std::bitset<FREAK_NB_PAIRS>*) (descriptors.data+(keypoints.size()-1)*descriptors.step[0]);
|
|
312 |
__m128i* ptrSSE = (__m128i*) (descriptors.data+(keypoints.size()-1)*descriptors.step[0]);
|
| 282 |
313 |
#endif
|
|
314 |
std::bitset<FREAK_NB_PAIRS>* ptrScalar = (std::bitset<FREAK_NB_PAIRS>*) (descriptors.data+(keypoints.size()-1)*descriptors.step[0]);
|
|
315 |
|
| 283 |
316 |
for( size_t k = keypoints.size(); k--; ) {
|
| 284 |
317 |
// estimate orientation (gradient)
|
| 285 |
318 |
if( !orientationNormalized ) {
|
| ... | ... | |
| 289 |
322 |
else {
|
| 290 |
323 |
// get the points intensity value in the un-rotated pattern
|
| 291 |
324 |
for( int i = FREAK_NB_POINTS; i--; ) {
|
| 292 |
|
pointsValue[i] = meanIntensity(image, imgIntegral, keypoints[k].pt.x,keypoints[k].pt.y, kpScaleIdx[k], 0, i);
|
|
325 |
pointsValue[i] = meanIntensity<srcMatType, iiMatType>(image, imgIntegral,
|
|
326 |
keypoints[k].pt.x, keypoints[k].pt.y,
|
|
327 |
kpScaleIdx[k], 0, i);
|
| 293 |
328 |
}
|
| 294 |
329 |
direction0 = 0;
|
| 295 |
330 |
direction1 = 0;
|
| ... | ... | |
| 310 |
345 |
}
|
| 311 |
346 |
// extract descriptor at the computed orientation
|
| 312 |
347 |
for( int i = FREAK_NB_POINTS; i--; ) {
|
| 313 |
|
pointsValue[i] = meanIntensity(image, imgIntegral, keypoints[k].pt.x,keypoints[k].pt.y, kpScaleIdx[k], thetaIdx, i);
|
|
348 |
pointsValue[i] = meanIntensity<srcMatType, iiMatType>(image, imgIntegral,
|
|
349 |
keypoints[k].pt.x, keypoints[k].pt.y,
|
|
350 |
kpScaleIdx[k], thetaIdx, i);
|
| 314 |
351 |
}
|
| 315 |
352 |
#if CV_SSE2
|
| 316 |
|
// note that comparisons order is modified in each block (but first 128 comparisons remain globally the same-->does not affect the 128,384 bits segmanted matching strategy)
|
| 317 |
|
int cnt = 0;
|
| 318 |
|
for( int n = FREAK_NB_PAIRS/128; n-- ; )
|
| 319 |
|
{
|
| 320 |
|
__m128i result128 = _mm_setzero_si128();
|
| 321 |
|
for( int m = 128/16; m--; cnt += 16 )
|
|
353 |
if (sizeof(srcMatType) == sizeof(char)) {
|
|
354 |
// note that comparisons order is modified in each block (but first 128 comparisons remain globally the same-->does not affect the 128,384 bits segmanted matching strategy)
|
|
355 |
int cnt = 0;
|
|
356 |
for( int n = FREAK_NB_PAIRS/128; n-- ; )
|
| 322 |
357 |
{
|
| 323 |
|
__m128i operand1 = _mm_set_epi8(
|
| 324 |
|
pointsValue[descriptionPairs[cnt+0].i],
|
| 325 |
|
pointsValue[descriptionPairs[cnt+1].i],
|
| 326 |
|
pointsValue[descriptionPairs[cnt+2].i],
|
| 327 |
|
pointsValue[descriptionPairs[cnt+3].i],
|
| 328 |
|
pointsValue[descriptionPairs[cnt+4].i],
|
| 329 |
|
pointsValue[descriptionPairs[cnt+5].i],
|
| 330 |
|
pointsValue[descriptionPairs[cnt+6].i],
|
| 331 |
|
pointsValue[descriptionPairs[cnt+7].i],
|
| 332 |
|
pointsValue[descriptionPairs[cnt+8].i],
|
| 333 |
|
pointsValue[descriptionPairs[cnt+9].i],
|
| 334 |
|
pointsValue[descriptionPairs[cnt+10].i],
|
| 335 |
|
pointsValue[descriptionPairs[cnt+11].i],
|
| 336 |
|
pointsValue[descriptionPairs[cnt+12].i],
|
| 337 |
|
pointsValue[descriptionPairs[cnt+13].i],
|
| 338 |
|
pointsValue[descriptionPairs[cnt+14].i],
|
| 339 |
|
pointsValue[descriptionPairs[cnt+15].i]);
|
| 340 |
|
|
| 341 |
|
__m128i operand2 = _mm_set_epi8(
|
| 342 |
|
pointsValue[descriptionPairs[cnt+0].j],
|
| 343 |
|
pointsValue[descriptionPairs[cnt+1].j],
|
| 344 |
|
pointsValue[descriptionPairs[cnt+2].j],
|
| 345 |
|
pointsValue[descriptionPairs[cnt+3].j],
|
| 346 |
|
pointsValue[descriptionPairs[cnt+4].j],
|
| 347 |
|
pointsValue[descriptionPairs[cnt+5].j],
|
| 348 |
|
pointsValue[descriptionPairs[cnt+6].j],
|
| 349 |
|
pointsValue[descriptionPairs[cnt+7].j],
|
| 350 |
|
pointsValue[descriptionPairs[cnt+8].j],
|
| 351 |
|
pointsValue[descriptionPairs[cnt+9].j],
|
| 352 |
|
pointsValue[descriptionPairs[cnt+10].j],
|
| 353 |
|
pointsValue[descriptionPairs[cnt+11].j],
|
| 354 |
|
pointsValue[descriptionPairs[cnt+12].j],
|
| 355 |
|
pointsValue[descriptionPairs[cnt+13].j],
|
| 356 |
|
pointsValue[descriptionPairs[cnt+14].j],
|
| 357 |
|
pointsValue[descriptionPairs[cnt+15].j]);
|
|
358 |
__m128i result128 = _mm_setzero_si128();
|
|
359 |
for( int m = 128/16; m--; cnt += 16 )
|
|
360 |
{
|
|
361 |
__m128i operand1 = _mm_set_epi8(
|
|
362 |
pointsValue[descriptionPairs[cnt+0].i],
|
|
363 |
pointsValue[descriptionPairs[cnt+1].i],
|
|
364 |
pointsValue[descriptionPairs[cnt+2].i],
|
|
365 |
pointsValue[descriptionPairs[cnt+3].i],
|
|
366 |
pointsValue[descriptionPairs[cnt+4].i],
|
|
367 |
pointsValue[descriptionPairs[cnt+5].i],
|
|
368 |
pointsValue[descriptionPairs[cnt+6].i],
|
|
369 |
pointsValue[descriptionPairs[cnt+7].i],
|
|
370 |
pointsValue[descriptionPairs[cnt+8].i],
|
|
371 |
pointsValue[descriptionPairs[cnt+9].i],
|
|
372 |
pointsValue[descriptionPairs[cnt+10].i],
|
|
373 |
pointsValue[descriptionPairs[cnt+11].i],
|
|
374 |
pointsValue[descriptionPairs[cnt+12].i],
|
|
375 |
pointsValue[descriptionPairs[cnt+13].i],
|
|
376 |
pointsValue[descriptionPairs[cnt+14].i],
|
|
377 |
pointsValue[descriptionPairs[cnt+15].i]);
|
|
378 |
|
|
379 |
__m128i operand2 = _mm_set_epi8(
|
|
380 |
pointsValue[descriptionPairs[cnt+0].j],
|
|
381 |
pointsValue[descriptionPairs[cnt+1].j],
|
|
382 |
pointsValue[descriptionPairs[cnt+2].j],
|
|
383 |
pointsValue[descriptionPairs[cnt+3].j],
|
|
384 |
pointsValue[descriptionPairs[cnt+4].j],
|
|
385 |
pointsValue[descriptionPairs[cnt+5].j],
|
|
386 |
pointsValue[descriptionPairs[cnt+6].j],
|
|
387 |
pointsValue[descriptionPairs[cnt+7].j],
|
|
388 |
pointsValue[descriptionPairs[cnt+8].j],
|
|
389 |
pointsValue[descriptionPairs[cnt+9].j],
|
|
390 |
pointsValue[descriptionPairs[cnt+10].j],
|
|
391 |
pointsValue[descriptionPairs[cnt+11].j],
|
|
392 |
pointsValue[descriptionPairs[cnt+12].j],
|
|
393 |
pointsValue[descriptionPairs[cnt+13].j],
|
|
394 |
pointsValue[descriptionPairs[cnt+14].j],
|
|
395 |
pointsValue[descriptionPairs[cnt+15].j]);
|
| 358 |
396 |
|
| 359 |
|
__m128i workReg = _mm_min_epu8(operand1, operand2); // emulated "not less than" for 8-bit UNSIGNED integers
|
| 360 |
|
workReg = _mm_cmpeq_epi8(workReg, operand2); // emulated "not less than" for 8-bit UNSIGNED integers
|
|
397 |
__m128i workReg = _mm_min_epu8(operand1, operand2); // emulated "not less than" for 8-bit UNSIGNED integers
|
|
398 |
workReg = _mm_cmpeq_epi8(workReg, operand2); // emulated "not less than" for 8-bit UNSIGNED integers
|
| 361 |
399 |
|
| 362 |
|
workReg = _mm_and_si128(_mm_set1_epi16(short(0x8080 >> m)), workReg); // merge the last 16 bits with the 128bits std::vector until full
|
| 363 |
|
result128 = _mm_or_si128(result128, workReg);
|
|
400 |
workReg = _mm_and_si128(_mm_set1_epi16(short(0x8080 >> m)), workReg); // merge the last 16 bits with the 128bits std::vector until full
|
|
401 |
result128 = _mm_or_si128(result128, workReg);
|
|
402 |
}
|
|
403 |
(*ptrSSE) = result128;
|
|
404 |
++ptrSSE;
|
| 364 |
405 |
}
|
| 365 |
|
(*ptr) = result128;
|
| 366 |
|
++ptr;
|
| 367 |
|
}
|
| 368 |
|
ptr -= 8;
|
| 369 |
|
#else
|
| 370 |
|
// extracting descriptor preserving the order of SSE version
|
| 371 |
|
int cnt = 0;
|
| 372 |
|
for( int n = 7; n < FREAK_NB_PAIRS; n += 128)
|
|
406 |
ptrSSE -= 8;
|
|
407 |
} else
|
|
408 |
#endif
|
| 373 |
409 |
{
|
| 374 |
|
for( int m = 8; m--; )
|
|
410 |
// extracting descriptor preserving the order of SSE version
|
|
411 |
int cnt = 0;
|
|
412 |
for( int n = 7; n < FREAK_NB_PAIRS; n += 128)
|
| 375 |
413 |
{
|
| 376 |
|
int nm = n-m;
|
| 377 |
|
for(int kk = nm+15*8; kk >= nm; kk-=8, ++cnt)
|
|
414 |
for( int m = 8; m--; )
|
| 378 |
415 |
{
|
| 379 |
|
ptr->set(kk, pointsValue[descriptionPairs[cnt].i] >= pointsValue[descriptionPairs[cnt].j]);
|
|
416 |
int nm = n-m;
|
|
417 |
for(int kk = nm+15*8; kk >= nm; kk-=8, ++cnt)
|
|
418 |
{
|
|
419 |
ptrScalar->set(kk, pointsValue[descriptionPairs[cnt].i] >= pointsValue[descriptionPairs[cnt].j]);
|
|
420 |
}
|
| 380 |
421 |
}
|
| 381 |
422 |
}
|
|
423 |
--ptrScalar;
|
| 382 |
424 |
}
|
| 383 |
|
--ptr;
|
| 384 |
|
#endif
|
| 385 |
425 |
}
|
| 386 |
426 |
}
|
| 387 |
427 |
else { // extract all possible comparisons for selection
|
| ... | ... | |
| 397 |
437 |
else {
|
| 398 |
438 |
//get the points intensity value in the un-rotated pattern
|
| 399 |
439 |
for( int i = FREAK_NB_POINTS;i--; )
|
| 400 |
|
pointsValue[i] = meanIntensity(image, imgIntegral, keypoints[k].pt.x,keypoints[k].pt.y, kpScaleIdx[k], 0, i);
|
|
440 |
pointsValue[i] = meanIntensity<srcMatType, iiMatType>(image, imgIntegral,
|
|
441 |
keypoints[k].pt.x,keypoints[k].pt.y,
|
|
442 |
kpScaleIdx[k], 0, i);
|
| 401 |
443 |
|
| 402 |
444 |
direction0 = 0;
|
| 403 |
445 |
direction1 = 0;
|
| ... | ... | |
| 419 |
461 |
}
|
| 420 |
462 |
// get the points intensity value in the rotated pattern
|
| 421 |
463 |
for( int i = FREAK_NB_POINTS; i--; ) {
|
| 422 |
|
pointsValue[i] = meanIntensity(image, imgIntegral, keypoints[k].pt.x,
|
| 423 |
|
keypoints[k].pt.y, kpScaleIdx[k], thetaIdx, i);
|
|
464 |
pointsValue[i] = meanIntensity<srcMatType, iiMatType>(image, imgIntegral,
|
|
465 |
keypoints[k].pt.x, keypoints[k].pt.y,
|
|
466 |
kpScaleIdx[k], thetaIdx, i);
|
| 424 |
467 |
}
|
| 425 |
468 |
|
| 426 |
469 |
int cnt(0);
|
| ... | ... | |
| 437 |
480 |
}
|
| 438 |
481 |
|
| 439 |
482 |
// simply take average on a square patch, not even gaussian approx
|
| 440 |
|
uchar FREAK::meanIntensity( const cv::Mat& image, const cv::Mat& integral,
|
|
483 |
template <typename imgType, typename iiType>
|
|
484 |
int FREAK::meanIntensity( const cv::Mat& image, const cv::Mat& integral,
|
| 441 |
485 |
const float kp_x,
|
| 442 |
486 |
const float kp_y,
|
| 443 |
487 |
const unsigned int scale,
|
| ... | ... | |
| 461 |
505 |
const int r_y = static_cast<int>((yf-y)*1024);
|
| 462 |
506 |
const int r_x_1 = (1024-r_x);
|
| 463 |
507 |
const int r_y_1 = (1024-r_y);
|
| 464 |
|
uchar* ptr = image.data+x+y*imagecols;
|
| 465 |
508 |
unsigned int ret_val;
|
| 466 |
509 |
// linear interpolation:
|
| 467 |
|
ret_val = (r_x_1*r_y_1*int(*ptr));
|
| 468 |
|
ptr++;
|
| 469 |
|
ret_val += (r_x*r_y_1*int(*ptr));
|
| 470 |
|
ptr += imagecols;
|
| 471 |
|
ret_val += (r_x*r_y*int(*ptr));
|
| 472 |
|
ptr--;
|
| 473 |
|
ret_val += (r_x_1*r_y*int(*ptr));
|
|
510 |
ret_val = r_x_1*r_y_1*int(image.at<imgType>(y , x ))
|
|
511 |
+ r_x *r_y_1*int(image.at<imgType>(y , x+1))
|
|
512 |
+ r_x_1*r_y *int(image.at<imgType>(y+1, x ))
|
|
513 |
+ r_x *r_y *int(image.at<imgType>(y+1, x+1));
|
| 474 |
514 |
//return the rounded mean
|
| 475 |
515 |
ret_val += 2 * 1024 * 1024;
|
| 476 |
|
return static_cast<uchar>(ret_val / (4 * 1024 * 1024));
|
|
516 |
return static_cast<int>(ret_val / (4 * 1024 * 1024));
|
| 477 |
517 |
}
|
| 478 |
518 |
|
| 479 |
519 |
// expected case:
|
| ... | ... | |
| 483 |
523 |
const int y_top = int(yf-radius+0.5);
|
| 484 |
524 |
const int x_right = int(xf+radius+1.5);//integral image is 1px wider
|
| 485 |
525 |
const int y_bottom = int(yf+radius+1.5);//integral image is 1px higher
|
| 486 |
|
int ret_val;
|
|
526 |
iiType ret_val;
|
| 487 |
527 |
|
| 488 |
|
ret_val = integral.at<int>(y_bottom,x_right);//bottom right corner
|
| 489 |
|
ret_val -= integral.at<int>(y_bottom,x_left);
|
| 490 |
|
ret_val += integral.at<int>(y_top,x_left);
|
| 491 |
|
ret_val -= integral.at<int>(y_top,x_right);
|
|
528 |
ret_val = integral.at<iiType>(y_bottom,x_right);//bottom right corner
|
|
529 |
ret_val -= integral.at<iiType>(y_bottom,x_left);
|
|
530 |
ret_val += integral.at<iiType>(y_top,x_left);
|
|
531 |
ret_val -= integral.at<iiType>(y_top,x_right);
|
| 492 |
532 |
ret_val = ret_val/( (x_right-x_left)* (y_bottom-y_top) );
|
| 493 |
533 |
//~ std::cout<<integral.step[1]<<std::endl;
|
| 494 |
|
return static_cast<uchar>(ret_val);
|
|
534 |
return static_cast<int>(ret_val);
|
| 495 |
535 |
}
|
| 496 |
536 |
|
| 497 |
537 |
// pair selection algorithm from a set of training images and corresponding keypoints
|