From 228070a74cf43f4f37e6c65110ea18587d76e4a6 Mon Sep 17 00:00:00 2001 From: Vincent Rabaud Date: Thu, 23 Aug 2012 14:33:11 +0200 Subject: [PATCH 1/8] split FAST in order to reuse it in BRISK --- modules/features2d/src/fast.cpp | 328 +--------------------- modules/features2d/src/fast_score.cpp | 374 ++++++++++++++++++++++++++ modules/features2d/src/fast_score.hpp | 64 +++++ 3 files changed, 439 insertions(+), 327 deletions(-) create mode 100644 modules/features2d/src/fast_score.cpp create mode 100644 modules/features2d/src/fast_score.hpp diff --git a/modules/features2d/src/fast.cpp b/modules/features2d/src/fast.cpp index b9834bbf0..9beba2145 100644 --- a/modules/features2d/src/fast.cpp +++ b/modules/features2d/src/fast.cpp @@ -42,335 +42,11 @@ The references are: */ #include "precomp.hpp" +#include "fast_score.hpp" namespace cv { -static void makeOffsets(int pixel[], int row_stride, int patternSize) -{ - switch(patternSize) { - case 16: - pixel[0] = 0 + row_stride * 3; - pixel[1] = 1 + row_stride * 3; - pixel[2] = 2 + row_stride * 2; - pixel[3] = 3 + row_stride * 1; - pixel[4] = 3 + row_stride * 0; - pixel[5] = 3 + row_stride * -1; - pixel[6] = 2 + row_stride * -2; - pixel[7] = 1 + row_stride * -3; - pixel[8] = 0 + row_stride * -3; - pixel[9] = -1 + row_stride * -3; - pixel[10] = -2 + row_stride * -2; - pixel[11] = -3 + row_stride * -1; - pixel[12] = -3 + row_stride * 0; - pixel[13] = -3 + row_stride * 1; - pixel[14] = -2 + row_stride * 2; - pixel[15] = -1 + row_stride * 3; - break; - case 12: - pixel[0] = 0 + row_stride * 2; - pixel[1] = 1 + row_stride * 2; - pixel[2] = 2 + row_stride * 1; - pixel[3] = 2 + row_stride * 0; - pixel[4] = 2 + row_stride * -1; - pixel[5] = 1 + row_stride * -2; - pixel[6] = 0 + row_stride * -2; - pixel[7] = -1 + row_stride * -2; - pixel[8] = -2 + row_stride * -1; - pixel[9] = -2 + row_stride * 0; - pixel[10] = -2 + row_stride * 1; - pixel[11] = -1 + row_stride * 2; - break; - case 8: - pixel[0] = 0 + row_stride * 1; - pixel[1] = 1 + row_stride * 1; - pixel[2] = 1 + row_stride * 0; - pixel[3] = 1 + row_stride * -1; - pixel[4] = 0 + row_stride * -1; - pixel[5] = -1 + row_stride * -1; - pixel[6] = 0 + row_stride * 0; - pixel[7] = 1 + row_stride * 1; - break; - } -} - -/*static void testCorner(const uchar* ptr, const int pixel[], int K, int N, int threshold) { - // check that with the computed "threshold" the pixel is still a corner - // and that with the increased-by-1 "threshold" the pixel is not a corner anymore - for( int delta = 0; delta <= 1; delta++ ) - { - int v0 = std::min(ptr[0] + threshold + delta, 255); - int v1 = std::max(ptr[0] - threshold - delta, 0); - int c0 = 0, c1 = 0; - - for( int k = 0; k < N; k++ ) - { - int x = ptr[pixel[k]]; - if(x > v0) - { - if( ++c0 > K ) - break; - c1 = 0; - } - else if( x < v1 ) - { - if( ++c1 > K ) - break; - c0 = 0; - } - else - { - c0 = c1 = 0; - } - } - CV_Assert( (delta == 0 && std::max(c0, c1) > K) || - (delta == 1 && std::max(c0, c1) <= K) ); - } -}*/ - -template -int cornerScore(const uchar* ptr, const int pixel[], int threshold); - -template<> -int cornerScore<16>(const uchar* ptr, const int pixel[], int threshold) -{ - const int K = 8, N = 16 + K + 1; - int k, v = ptr[0]; - short d[N]; - for( k = 0; k < N; k++ ) - d[k] = (short)(v - ptr[pixel[k]]); - -#if CV_SSE2 - __m128i q0 = _mm_set1_epi16(-1000), q1 = _mm_set1_epi16(1000); - for( k = 0; k < 16; k += 8 ) - { - __m128i v0 = _mm_loadu_si128((__m128i*)(d+k+1)); - __m128i v1 = _mm_loadu_si128((__m128i*)(d+k+2)); - __m128i a = _mm_min_epi16(v0, v1); - __m128i b = _mm_max_epi16(v0, v1); - v0 = _mm_loadu_si128((__m128i*)(d+k+3)); - a = _mm_min_epi16(a, v0); - b = _mm_max_epi16(b, v0); - v0 = _mm_loadu_si128((__m128i*)(d+k+4)); - a = _mm_min_epi16(a, v0); - b = _mm_max_epi16(b, v0); - v0 = _mm_loadu_si128((__m128i*)(d+k+5)); - a = _mm_min_epi16(a, v0); - b = _mm_max_epi16(b, v0); - v0 = _mm_loadu_si128((__m128i*)(d+k+6)); - a = _mm_min_epi16(a, v0); - b = _mm_max_epi16(b, v0); - v0 = _mm_loadu_si128((__m128i*)(d+k+7)); - a = _mm_min_epi16(a, v0); - b = _mm_max_epi16(b, v0); - v0 = _mm_loadu_si128((__m128i*)(d+k+8)); - a = _mm_min_epi16(a, v0); - b = _mm_max_epi16(b, v0); - v0 = _mm_loadu_si128((__m128i*)(d+k)); - q0 = _mm_max_epi16(q0, _mm_min_epi16(a, v0)); - q1 = _mm_min_epi16(q1, _mm_max_epi16(b, v0)); - v0 = _mm_loadu_si128((__m128i*)(d+k+9)); - q0 = _mm_max_epi16(q0, _mm_min_epi16(a, v0)); - q1 = _mm_min_epi16(q1, _mm_max_epi16(b, v0)); - } - q0 = _mm_max_epi16(q0, _mm_sub_epi16(_mm_setzero_si128(), q1)); - q0 = _mm_max_epi16(q0, _mm_unpackhi_epi64(q0, q0)); - q0 = _mm_max_epi16(q0, _mm_srli_si128(q0, 4)); - q0 = _mm_max_epi16(q0, _mm_srli_si128(q0, 2)); - threshold = (short)_mm_cvtsi128_si32(q0) - 1; -#else - int a0 = threshold; - for( k = 0; k < 16; k += 2 ) - { - int a = std::min((int)d[k+1], (int)d[k+2]); - a = std::min(a, (int)d[k+3]); - if( a <= a0 ) - continue; - a = std::min(a, (int)d[k+4]); - a = std::min(a, (int)d[k+5]); - a = std::min(a, (int)d[k+6]); - a = std::min(a, (int)d[k+7]); - a = std::min(a, (int)d[k+8]); - a0 = std::max(a0, std::min(a, (int)d[k])); - a0 = std::max(a0, std::min(a, (int)d[k+9])); - } - - int b0 = -a0; - for( k = 0; k < 16; k += 2 ) - { - int b = std::max((int)d[k+1], (int)d[k+2]); - b = std::max(b, (int)d[k+3]); - b = std::max(b, (int)d[k+4]); - b = std::max(b, (int)d[k+5]); - if( b >= b0 ) - continue; - b = std::max(b, (int)d[k+6]); - b = std::max(b, (int)d[k+7]); - b = std::max(b, (int)d[k+8]); - - b0 = std::min(b0, std::max(b, (int)d[k])); - b0 = std::min(b0, std::max(b, (int)d[k+9])); - } - - threshold = -b0-1; -#endif - -#if 0 - testCorner(ptr, pixel, K, N, threshold); -#endif - return threshold; -} - -template<> -int cornerScore<12>(const uchar* ptr, const int pixel[], int threshold) -{ - const int K = 6, N = 12 + K + 1; - int k, v = ptr[0]; - short d[N]; - for( k = 0; k < N; k++ ) - d[k] = (short)(v - ptr[pixel[k]]); - -#if CV_SSE2 - __m128i q0 = _mm_set1_epi16(-1000), q1 = _mm_set1_epi16(1000); - for( k = 0; k < 16; k += 8 ) - { - __m128i v0 = _mm_loadu_si128((__m128i*)(d+k+1)); - __m128i v1 = _mm_loadu_si128((__m128i*)(d+k+2)); - __m128i a = _mm_min_epi16(v0, v1); - __m128i b = _mm_max_epi16(v0, v1); - v0 = _mm_loadu_si128((__m128i*)(d+k+3)); - a = _mm_min_epi16(a, v0); - b = _mm_max_epi16(b, v0); - v0 = _mm_loadu_si128((__m128i*)(d+k+4)); - a = _mm_min_epi16(a, v0); - b = _mm_max_epi16(b, v0); - v0 = _mm_loadu_si128((__m128i*)(d+k+5)); - a = _mm_min_epi16(a, v0); - b = _mm_max_epi16(b, v0); - v0 = _mm_loadu_si128((__m128i*)(d+k+6)); - a = _mm_min_epi16(a, v0); - b = _mm_max_epi16(b, v0); - v0 = _mm_loadu_si128((__m128i*)(d+k)); - q0 = _mm_max_epi16(q0, _mm_min_epi16(a, v0)); - q1 = _mm_min_epi16(q1, _mm_max_epi16(b, v0)); - v0 = _mm_loadu_si128((__m128i*)(d+k+7)); - q0 = _mm_max_epi16(q0, _mm_min_epi16(a, v0)); - q1 = _mm_min_epi16(q1, _mm_max_epi16(b, v0)); - } - q0 = _mm_max_epi16(q0, _mm_sub_epi16(_mm_setzero_si128(), q1)); - q0 = _mm_max_epi16(q0, _mm_unpackhi_epi64(q0, q0)); - q0 = _mm_max_epi16(q0, _mm_srli_si128(q0, 4)); - q0 = _mm_max_epi16(q0, _mm_srli_si128(q0, 2)); - threshold = (short)_mm_cvtsi128_si32(q0) - 1; -#else - int a0 = threshold; - for( k = 0; k < 12; k += 2 ) - { - int a = std::min((int)d[k+1], (int)d[k+2]); - if( a <= a0 ) - continue; - a = std::min(a, (int)d[k+3]); - a = std::min(a, (int)d[k+4]); - a = std::min(a, (int)d[k+5]); - a = std::min(a, (int)d[k+6]); - a0 = std::max(a0, std::min(a, (int)d[k])); - a0 = std::max(a0, std::min(a, (int)d[k+7])); - } - - int b0 = -a0; - for( k = 0; k < 12; k += 2 ) - { - int b = std::max((int)d[k+1], (int)d[k+2]); - b = std::max(b, (int)d[k+3]); - b = std::max(b, (int)d[k+4]); - if( b >= b0 ) - continue; - b = std::max(b, (int)d[k+5]); - b = std::max(b, (int)d[k+6]); - - b0 = std::min(b0, std::max(b, (int)d[k])); - b0 = std::min(b0, std::max(b, (int)d[k+7])); - } - - threshold = -b0-1; -#endif - -#if 0 - testCorner(ptr, pixel, K, N, threshold); -#endif - return threshold; -} - -template<> -int cornerScore<8>(const uchar* ptr, const int pixel[], int threshold) -{ - const int K = 4, N = 8 + K + 1; - int k, v = ptr[0]; - short d[N]; - for( k = 0; k < N; k++ ) - d[k] = (short)(v - ptr[pixel[k]]); - -#if CV_SSE2 - __m128i q0 = _mm_set1_epi16(-1000), q1 = _mm_set1_epi16(1000); - for( k = 0; k < 16; k += 8 ) - { - __m128i v0 = _mm_loadu_si128((__m128i*)(d+k+1)); - __m128i v1 = _mm_loadu_si128((__m128i*)(d+k+2)); - __m128i a = _mm_min_epi16(v0, v1); - __m128i b = _mm_max_epi16(v0, v1); - v0 = _mm_loadu_si128((__m128i*)(d+k+3)); - a = _mm_min_epi16(a, v0); - b = _mm_max_epi16(b, v0); - v0 = _mm_loadu_si128((__m128i*)(d+k+4)); - a = _mm_min_epi16(a, v0); - b = _mm_max_epi16(b, v0); - v0 = _mm_loadu_si128((__m128i*)(d+k)); - q0 = _mm_max_epi16(q0, _mm_min_epi16(a, v0)); - q1 = _mm_min_epi16(q1, _mm_max_epi16(b, v0)); - v0 = _mm_loadu_si128((__m128i*)(d+k+5)); - q0 = _mm_max_epi16(q0, _mm_min_epi16(a, v0)); - q1 = _mm_min_epi16(q1, _mm_max_epi16(b, v0)); - } - q0 = _mm_max_epi16(q0, _mm_sub_epi16(_mm_setzero_si128(), q1)); - q0 = _mm_max_epi16(q0, _mm_unpackhi_epi64(q0, q0)); - q0 = _mm_max_epi16(q0, _mm_srli_si128(q0, 4)); - q0 = _mm_max_epi16(q0, _mm_srli_si128(q0, 2)); - threshold = (short)_mm_cvtsi128_si32(q0) - 1; -#else - int a0 = threshold; - for( k = 0; k < 8; k += 2 ) - { - int a = std::min((int)d[k+1], (int)d[k+2]); - if( a <= a0 ) - continue; - a = std::min(a, (int)d[k+3]); - a = std::min(a, (int)d[k+4]); - a0 = std::max(a0, std::min(a, (int)d[k])); - a0 = std::max(a0, std::min(a, (int)d[k+5])); - } - - int b0 = -a0; - for( k = 0; k < 8; k += 2 ) - { - int b = std::max((int)d[k+1], (int)d[k+2]); - b = std::max(b, (int)d[k+3]); - if( b >= b0 ) - continue; - b = std::max(b, (int)d[k+4]); - - b0 = std::min(b0, std::max(b, (int)d[k])); - b0 = std::min(b0, std::max(b, (int)d[k+5])); - } - - threshold = -b0-1; -#endif - -#if 0 - testCorner(ptr, pixel, K, N, threshold); -#endif - return threshold; -} - template void FAST_t(InputArray _img, std::vector& keypoints, int threshold, bool nonmax_suppression) { @@ -381,8 +57,6 @@ void FAST_t(InputArray _img, std::vector& keypoints, int threshold, bo #endif int i, j, k, pixel[25]; makeOffsets(pixel, (int)img.step, patternSize); - for(k = patternSize; k < 25; k++) - pixel[k] = pixel[k - patternSize]; keypoints.clear(); diff --git a/modules/features2d/src/fast_score.cpp b/modules/features2d/src/fast_score.cpp new file mode 100644 index 000000000..aa5881e2d --- /dev/null +++ b/modules/features2d/src/fast_score.cpp @@ -0,0 +1,374 @@ +/* This is FAST corner detector, contributed to OpenCV by the author, Edward Rosten. + Below is the original copyright and the references */ + +/* +Copyright (c) 2006, 2008 Edward Rosten +All rights reserved. + +Redistribution and use in source and binary forms, with or without +modification, are permitted provided that the following conditions +are met: + + *Redistributions of source code must retain the above copyright + notice, this list of conditions and the following disclaimer. + + *Redistributions in binary form must reproduce the above copyright + notice, this list of conditions and the following disclaimer in the + documentation and/or other materials provided with the distribution. + + *Neither the name of the University of Cambridge nor the names of + its contributors may be used to endorse or promote products derived + from this software without specific prior written permission. + +THIS SOFTWARE IS PROVIDED BY THE COPYRIGHT HOLDERS AND CONTRIBUTORS +"AS IS" AND ANY EXPRESS OR IMPLIED WARRANTIES, INCLUDING, BUT NOT +LIMITED TO, THE IMPLIED WARRANTIES OF MERCHANTABILITY AND FITNESS FOR +A PARTICULAR PURPOSE ARE DISCLAIMED. IN NO EVENT SHALL THE COPYRIGHT OWNER OR +CONTRIBUTORS BE LIABLE FOR ANY DIRECT, INDIRECT, INCIDENTAL, SPECIAL, +EXEMPLARY, OR CONSEQUENTIAL DAMAGES (INCLUDING, BUT NOT LIMITED TO, +PROCUREMENT OF SUBSTITUTE GOODS OR SERVICES; LOSS OF USE, DATA, OR +PROFITS; OR BUSINESS INTERRUPTION) HOWEVER CAUSED AND ON ANY THEORY OF +LIABILITY, WHETHER IN CONTRACT, STRICT LIABILITY, OR TORT (INCLUDING +NEGLIGENCE OR OTHERWISE) ARISING IN ANY WAY OUT OF THE USE OF THIS +SOFTWARE, EVEN IF ADVISED OF THE POSSIBILITY OF SUCH DAMAGE. +*/ + +/* +The references are: + * Machine learning for high-speed corner detection, + E. Rosten and T. Drummond, ECCV 2006 + * Faster and better: A machine learning approach to corner detection + E. Rosten, R. Porter and T. Drummond, PAMI, 2009 +*/ + +#include "fast_score.hpp" + +namespace cv +{ + +void makeOffsets(int pixel[25], int row_stride, int patternSize) +{ + CV_Assert(pixel != 0); + switch(patternSize) { + case 16: + pixel[0] = 0 + row_stride * 3; + pixel[1] = 1 + row_stride * 3; + pixel[2] = 2 + row_stride * 2; + pixel[3] = 3 + row_stride * 1; + pixel[4] = 3 + row_stride * 0; + pixel[5] = 3 + row_stride * -1; + pixel[6] = 2 + row_stride * -2; + pixel[7] = 1 + row_stride * -3; + pixel[8] = 0 + row_stride * -3; + pixel[9] = -1 + row_stride * -3; + pixel[10] = -2 + row_stride * -2; + pixel[11] = -3 + row_stride * -1; + pixel[12] = -3 + row_stride * 0; + pixel[13] = -3 + row_stride * 1; + pixel[14] = -2 + row_stride * 2; + pixel[15] = -1 + row_stride * 3; + break; + case 12: + pixel[0] = 0 + row_stride * 2; + pixel[1] = 1 + row_stride * 2; + pixel[2] = 2 + row_stride * 1; + pixel[3] = 2 + row_stride * 0; + pixel[4] = 2 + row_stride * -1; + pixel[5] = 1 + row_stride * -2; + pixel[6] = 0 + row_stride * -2; + pixel[7] = -1 + row_stride * -2; + pixel[8] = -2 + row_stride * -1; + pixel[9] = -2 + row_stride * 0; + pixel[10] = -2 + row_stride * 1; + pixel[11] = -1 + row_stride * 2; + break; + case 8: + pixel[0] = 0 + row_stride * 1; + pixel[1] = 1 + row_stride * 1; + pixel[2] = 1 + row_stride * 0; + pixel[3] = 1 + row_stride * -1; + pixel[4] = 0 + row_stride * -1; + pixel[5] = -1 + row_stride * -1; + pixel[6] = 0 + row_stride * 0; + pixel[7] = 1 + row_stride * 1; + break; + } + for(int k = patternSize; k < 25; k++) + pixel[k] = pixel[k - patternSize]; +} + +/*static void testCorner(const uchar* ptr, const int pixel[], int K, int N, int threshold) { + // check that with the computed "threshold" the pixel is still a corner + // and that with the increased-by-1 "threshold" the pixel is not a corner anymore + for( int delta = 0; delta <= 1; delta++ ) + { + int v0 = std::min(ptr[0] + threshold + delta, 255); + int v1 = std::max(ptr[0] - threshold - delta, 0); + int c0 = 0, c1 = 0; + + for( int k = 0; k < N; k++ ) + { + int x = ptr[pixel[k]]; + if(x > v0) + { + if( ++c0 > K ) + break; + c1 = 0; + } + else if( x < v1 ) + { + if( ++c1 > K ) + break; + c0 = 0; + } + else + { + c0 = c1 = 0; + } + } + CV_Assert( (delta == 0 && std::max(c0, c1) > K) || + (delta == 1 && std::max(c0, c1) <= K) ); + } +}*/ + +template<> +int cornerScore<16>(const uchar* ptr, const int pixel[], int threshold) +{ + const int K = 8, N = 16 + K + 1; + int k, v = ptr[0]; + short d[N]; + for( k = 0; k < N; k++ ) + d[k] = (short)(v - ptr[pixel[k]]); + +#if CV_SSE2 + __m128i q0 = _mm_set1_epi16(-1000), q1 = _mm_set1_epi16(1000); + for( k = 0; k < 16; k += 8 ) + { + __m128i v0 = _mm_loadu_si128((__m128i*)(d+k+1)); + __m128i v1 = _mm_loadu_si128((__m128i*)(d+k+2)); + __m128i a = _mm_min_epi16(v0, v1); + __m128i b = _mm_max_epi16(v0, v1); + v0 = _mm_loadu_si128((__m128i*)(d+k+3)); + a = _mm_min_epi16(a, v0); + b = _mm_max_epi16(b, v0); + v0 = _mm_loadu_si128((__m128i*)(d+k+4)); + a = _mm_min_epi16(a, v0); + b = _mm_max_epi16(b, v0); + v0 = _mm_loadu_si128((__m128i*)(d+k+5)); + a = _mm_min_epi16(a, v0); + b = _mm_max_epi16(b, v0); + v0 = _mm_loadu_si128((__m128i*)(d+k+6)); + a = _mm_min_epi16(a, v0); + b = _mm_max_epi16(b, v0); + v0 = _mm_loadu_si128((__m128i*)(d+k+7)); + a = _mm_min_epi16(a, v0); + b = _mm_max_epi16(b, v0); + v0 = _mm_loadu_si128((__m128i*)(d+k+8)); + a = _mm_min_epi16(a, v0); + b = _mm_max_epi16(b, v0); + v0 = _mm_loadu_si128((__m128i*)(d+k)); + q0 = _mm_max_epi16(q0, _mm_min_epi16(a, v0)); + q1 = _mm_min_epi16(q1, _mm_max_epi16(b, v0)); + v0 = _mm_loadu_si128((__m128i*)(d+k+9)); + q0 = _mm_max_epi16(q0, _mm_min_epi16(a, v0)); + q1 = _mm_min_epi16(q1, _mm_max_epi16(b, v0)); + } + q0 = _mm_max_epi16(q0, _mm_sub_epi16(_mm_setzero_si128(), q1)); + q0 = _mm_max_epi16(q0, _mm_unpackhi_epi64(q0, q0)); + q0 = _mm_max_epi16(q0, _mm_srli_si128(q0, 4)); + q0 = _mm_max_epi16(q0, _mm_srli_si128(q0, 2)); + threshold = (short)_mm_cvtsi128_si32(q0) - 1; +#else + int a0 = threshold; + for( k = 0; k < 16; k += 2 ) + { + int a = std::min((int)d[k+1], (int)d[k+2]); + a = std::min(a, (int)d[k+3]); + if( a <= a0 ) + continue; + a = std::min(a, (int)d[k+4]); + a = std::min(a, (int)d[k+5]); + a = std::min(a, (int)d[k+6]); + a = std::min(a, (int)d[k+7]); + a = std::min(a, (int)d[k+8]); + a0 = std::max(a0, std::min(a, (int)d[k])); + a0 = std::max(a0, std::min(a, (int)d[k+9])); + } + + int b0 = -a0; + for( k = 0; k < 16; k += 2 ) + { + int b = std::max((int)d[k+1], (int)d[k+2]); + b = std::max(b, (int)d[k+3]); + b = std::max(b, (int)d[k+4]); + b = std::max(b, (int)d[k+5]); + if( b >= b0 ) + continue; + b = std::max(b, (int)d[k+6]); + b = std::max(b, (int)d[k+7]); + b = std::max(b, (int)d[k+8]); + + b0 = std::min(b0, std::max(b, (int)d[k])); + b0 = std::min(b0, std::max(b, (int)d[k+9])); + } + + threshold = -b0-1; +#endif + +#if 0 + testCorner(ptr, pixel, K, N, threshold); +#endif + return threshold; +} + +template<> +int cornerScore<12>(const uchar* ptr, const int pixel[], int threshold) +{ + const int K = 6, N = 12 + K + 1; + int k, v = ptr[0]; + short d[N]; + for( k = 0; k < N; k++ ) + d[k] = (short)(v - ptr[pixel[k]]); + +#if CV_SSE2 + __m128i q0 = _mm_set1_epi16(-1000), q1 = _mm_set1_epi16(1000); + for( k = 0; k < 16; k += 8 ) + { + __m128i v0 = _mm_loadu_si128((__m128i*)(d+k+1)); + __m128i v1 = _mm_loadu_si128((__m128i*)(d+k+2)); + __m128i a = _mm_min_epi16(v0, v1); + __m128i b = _mm_max_epi16(v0, v1); + v0 = _mm_loadu_si128((__m128i*)(d+k+3)); + a = _mm_min_epi16(a, v0); + b = _mm_max_epi16(b, v0); + v0 = _mm_loadu_si128((__m128i*)(d+k+4)); + a = _mm_min_epi16(a, v0); + b = _mm_max_epi16(b, v0); + v0 = _mm_loadu_si128((__m128i*)(d+k+5)); + a = _mm_min_epi16(a, v0); + b = _mm_max_epi16(b, v0); + v0 = _mm_loadu_si128((__m128i*)(d+k+6)); + a = _mm_min_epi16(a, v0); + b = _mm_max_epi16(b, v0); + v0 = _mm_loadu_si128((__m128i*)(d+k)); + q0 = _mm_max_epi16(q0, _mm_min_epi16(a, v0)); + q1 = _mm_min_epi16(q1, _mm_max_epi16(b, v0)); + v0 = _mm_loadu_si128((__m128i*)(d+k+7)); + q0 = _mm_max_epi16(q0, _mm_min_epi16(a, v0)); + q1 = _mm_min_epi16(q1, _mm_max_epi16(b, v0)); + } + q0 = _mm_max_epi16(q0, _mm_sub_epi16(_mm_setzero_si128(), q1)); + q0 = _mm_max_epi16(q0, _mm_unpackhi_epi64(q0, q0)); + q0 = _mm_max_epi16(q0, _mm_srli_si128(q0, 4)); + q0 = _mm_max_epi16(q0, _mm_srli_si128(q0, 2)); + threshold = (short)_mm_cvtsi128_si32(q0) - 1; +#else + int a0 = threshold; + for( k = 0; k < 12; k += 2 ) + { + int a = std::min((int)d[k+1], (int)d[k+2]); + if( a <= a0 ) + continue; + a = std::min(a, (int)d[k+3]); + a = std::min(a, (int)d[k+4]); + a = std::min(a, (int)d[k+5]); + a = std::min(a, (int)d[k+6]); + a0 = std::max(a0, std::min(a, (int)d[k])); + a0 = std::max(a0, std::min(a, (int)d[k+7])); + } + + int b0 = -a0; + for( k = 0; k < 12; k += 2 ) + { + int b = std::max((int)d[k+1], (int)d[k+2]); + b = std::max(b, (int)d[k+3]); + b = std::max(b, (int)d[k+4]); + if( b >= b0 ) + continue; + b = std::max(b, (int)d[k+5]); + b = std::max(b, (int)d[k+6]); + + b0 = std::min(b0, std::max(b, (int)d[k])); + b0 = std::min(b0, std::max(b, (int)d[k+7])); + } + + threshold = -b0-1; +#endif + +#if 0 + testCorner(ptr, pixel, K, N, threshold); +#endif + return threshold; +} + +template<> +int cornerScore<8>(const uchar* ptr, const int pixel[], int threshold) +{ + const int K = 4, N = 8 + K + 1; + int k, v = ptr[0]; + short d[N]; + for( k = 0; k < N; k++ ) + d[k] = (short)(v - ptr[pixel[k]]); + +#if CV_SSE2 + __m128i q0 = _mm_set1_epi16(-1000), q1 = _mm_set1_epi16(1000); + for( k = 0; k < 16; k += 8 ) + { + __m128i v0 = _mm_loadu_si128((__m128i*)(d+k+1)); + __m128i v1 = _mm_loadu_si128((__m128i*)(d+k+2)); + __m128i a = _mm_min_epi16(v0, v1); + __m128i b = _mm_max_epi16(v0, v1); + v0 = _mm_loadu_si128((__m128i*)(d+k+3)); + a = _mm_min_epi16(a, v0); + b = _mm_max_epi16(b, v0); + v0 = _mm_loadu_si128((__m128i*)(d+k+4)); + a = _mm_min_epi16(a, v0); + b = _mm_max_epi16(b, v0); + v0 = _mm_loadu_si128((__m128i*)(d+k)); + q0 = _mm_max_epi16(q0, _mm_min_epi16(a, v0)); + q1 = _mm_min_epi16(q1, _mm_max_epi16(b, v0)); + v0 = _mm_loadu_si128((__m128i*)(d+k+5)); + q0 = _mm_max_epi16(q0, _mm_min_epi16(a, v0)); + q1 = _mm_min_epi16(q1, _mm_max_epi16(b, v0)); + } + q0 = _mm_max_epi16(q0, _mm_sub_epi16(_mm_setzero_si128(), q1)); + q0 = _mm_max_epi16(q0, _mm_unpackhi_epi64(q0, q0)); + q0 = _mm_max_epi16(q0, _mm_srli_si128(q0, 4)); + q0 = _mm_max_epi16(q0, _mm_srli_si128(q0, 2)); + threshold = (short)_mm_cvtsi128_si32(q0) - 1; +#else + int a0 = threshold; + for( k = 0; k < 8; k += 2 ) + { + int a = std::min((int)d[k+1], (int)d[k+2]); + if( a <= a0 ) + continue; + a = std::min(a, (int)d[k+3]); + a = std::min(a, (int)d[k+4]); + a0 = std::max(a0, std::min(a, (int)d[k])); + a0 = std::max(a0, std::min(a, (int)d[k+5])); + } + + int b0 = -a0; + for( k = 0; k < 8; k += 2 ) + { + int b = std::max((int)d[k+1], (int)d[k+2]); + b = std::max(b, (int)d[k+3]); + if( b >= b0 ) + continue; + b = std::max(b, (int)d[k+4]); + + b0 = std::min(b0, std::max(b, (int)d[k])); + b0 = std::min(b0, std::max(b, (int)d[k+5])); + } + + threshold = -b0-1; +#endif + +#if 0 + testCorner(ptr, pixel, K, N, threshold); +#endif + return threshold; +} + +} diff --git a/modules/features2d/src/fast_score.hpp b/modules/features2d/src/fast_score.hpp new file mode 100644 index 000000000..f6e9f6fc6 --- /dev/null +++ b/modules/features2d/src/fast_score.hpp @@ -0,0 +1,64 @@ +/* This is FAST corner detector, contributed to OpenCV by the author, Edward Rosten. + Below is the original copyright and the references */ + +/* +Copyright (c) 2006, 2008 Edward Rosten +All rights reserved. + +Redistribution and use in source and binary forms, with or without +modification, are permitted provided that the following conditions +are met: + + *Redistributions of source code must retain the above copyright + notice, this list of conditions and the following disclaimer. + + *Redistributions in binary form must reproduce the above copyright + notice, this list of conditions and the following disclaimer in the + documentation and/or other materials provided with the distribution. + + *Neither the name of the University of Cambridge nor the names of + its contributors may be used to endorse or promote products derived + from this software without specific prior written permission. + +THIS SOFTWARE IS PROVIDED BY THE COPYRIGHT HOLDERS AND CONTRIBUTORS +"AS IS" AND ANY EXPRESS OR IMPLIED WARRANTIES, INCLUDING, BUT NOT +LIMITED TO, THE IMPLIED WARRANTIES OF MERCHANTABILITY AND FITNESS FOR +A PARTICULAR PURPOSE ARE DISCLAIMED. IN NO EVENT SHALL THE COPYRIGHT OWNER OR +CONTRIBUTORS BE LIABLE FOR ANY DIRECT, INDIRECT, INCIDENTAL, SPECIAL, +EXEMPLARY, OR CONSEQUENTIAL DAMAGES (INCLUDING, BUT NOT LIMITED TO, +PROCUREMENT OF SUBSTITUTE GOODS OR SERVICES; LOSS OF USE, DATA, OR +PROFITS; OR BUSINESS INTERRUPTION) HOWEVER CAUSED AND ON ANY THEORY OF +LIABILITY, WHETHER IN CONTRACT, STRICT LIABILITY, OR TORT (INCLUDING +NEGLIGENCE OR OTHERWISE) ARISING IN ANY WAY OUT OF THE USE OF THIS +SOFTWARE, EVEN IF ADVISED OF THE POSSIBILITY OF SUCH DAMAGE. +*/ + +/* +The references are: + * Machine learning for high-speed corner detection, + E. Rosten and T. Drummond, ECCV 2006 + * Faster and better: A machine learning approach to corner detection + E. Rosten, R. Porter and T. Drummond, PAMI, 2009 +*/ + +#ifndef __OPENCV_FEATURES_2D_FAST_HPP__ +#define __OPENCV_FEATURES_2D_FAST_HPP__ + +#ifdef __cplusplus + +#include "precomp.hpp" + +namespace cv +{ + +void makeOffsets(int pixel[25], int row_stride, int patternSize); + +//static void testCorner(const uchar* ptr, const int pixel[], int K, int N, int threshold); + +template +int cornerScore(const uchar* ptr, const int pixel[], int threshold); + +} + +#endif +#endif From 13ded36ecb831d512a0f92a0baebdab8bedded6f Mon Sep 17 00:00:00 2001 From: Vincent Rabaud Date: Thu, 23 Aug 2012 14:52:01 +0200 Subject: [PATCH 2/8] initial addition of BRISK with some tests --- .../include/opencv2/features2d/features2d.hpp | 86 + modules/features2d/src/brisk.cpp | 2277 +++++++++++++++++ modules/features2d/src/features2d_init.cpp | 7 + modules/features2d/test/test_brisk.cpp | 95 + .../test/test_descriptors_regression.cpp | 7 + .../test/test_detectors_regression.cpp | 6 + 6 files changed, 2478 insertions(+) create mode 100755 modules/features2d/src/brisk.cpp create mode 100644 modules/features2d/test/test_brisk.cpp diff --git a/modules/features2d/include/opencv2/features2d/features2d.hpp b/modules/features2d/include/opencv2/features2d/features2d.hpp index ff52822ed..a081fd4ee 100644 --- a/modules/features2d/include/opencv2/features2d/features2d.hpp +++ b/modules/features2d/include/opencv2/features2d/features2d.hpp @@ -267,6 +267,92 @@ public: static Ptr create( const string& name ); }; +/*! + BRISK implementation +*/ +class CV_EXPORTS_W BRISK : public Feature2D +{ +public: + CV_WRAP explicit BRISK(int thresh=30, int octaves=3, float patternScale=1.0f); + + virtual ~BRISK(); + + // returns the descriptor size in bytes + int descriptorSize() const; + // returns the descriptor type + int descriptorType() const; + + // Compute the BRISK features on an image + void operator()(InputArray image, InputArray mask, vector& keypoints) const; + + // Compute the BRISK features and descriptors on an image + void operator()( InputArray image, InputArray mask, vector& keypoints, + OutputArray descriptors, bool useProvidedKeypoints=false ) const; + + AlgorithmInfo* info() const; + + // custom setup + CV_WRAP explicit BRISK(std::vector &radiusList, std::vector &numberList, + float dMax=5.85f, float dMin=8.2f, std::vector indexChange=std::vector()); + + // call this to generate the kernel: + // circle of radius r (pixels), with n points; + // short pairings with dMax, long pairings with dMin + CV_WRAP void generateKernel(std::vector &radiusList, + std::vector &numberList, float dMax=5.85f, float dMin=8.2f, + std::vector indexChange=std::vector()); + +protected: + + void computeImpl( const Mat& image, vector& keypoints, Mat& descriptors ) const; + void detectImpl( const Mat& image, vector& keypoints, const Mat& mask=Mat() ) const; + + // Feature parameters + CV_PROP_RW int threshold; + CV_PROP_RW int octaves; + + // some helper structures for the Brisk pattern representation + struct BriskPatternPoint{ + float x; // x coordinate relative to center + float y; // x coordinate relative to center + float sigma; // Gaussian smoothing sigma + }; + struct BriskShortPair{ + unsigned int i; // index of the first pattern point + unsigned int j; // index of other pattern point + }; + struct BriskLongPair{ + unsigned int i; // index of the first pattern point + unsigned int j; // index of other pattern point + int weighted_dx; // 1024.0/dx + int weighted_dy; // 1024.0/dy + }; + inline int smoothedIntensity(const cv::Mat& image, + const cv::Mat& integral,const float key_x, + const float key_y, const unsigned int scale, + const unsigned int rot, const unsigned int point) const; + // pattern properties + BriskPatternPoint* patternPoints_; //[i][rotation][scale] + unsigned int points_; // total number of collocation points + float* scaleList_; // lists the scaling per scale index [scale] + unsigned int* sizeList_; // lists the total pattern size per scale index [scale] + static const unsigned int scales_; // scales discretization + static const float scalerange_; // span of sizes 40->4 Octaves - else, this needs to be adjusted... + static const unsigned int n_rot_; // discretization of the rotation look-up + + // pairs + int strings_; // number of uchars the descriptor consists of + float dMax_; // short pair maximum distance + float dMin_; // long pair maximum distance + BriskShortPair* shortPairs_; // d<_dMax + BriskLongPair* longPairs_; // d>_dMin + unsigned int noShortPairs_; // number of shortParis + unsigned int noLongPairs_; // number of longParis + + // general + static const float basicSize_; +}; + /*! ORB implementation. diff --git a/modules/features2d/src/brisk.cpp b/modules/features2d/src/brisk.cpp new file mode 100755 index 000000000..415ed7821 --- /dev/null +++ b/modules/features2d/src/brisk.cpp @@ -0,0 +1,2277 @@ +/********************************************************************* + * Software License Agreement (BSD License) + * + * Copyright (C) 2011 The Autonomous Systems Lab (ASL), ETH Zurich, + * Stefan Leutenegger, Simon Lynen and Margarita Chli. + * Copyright (c) 2009, Willow Garage, Inc. + * All rights reserved. + * + * Redistribution and use in source and binary forms, with or without + * modification, are permitted provided that the following conditions + * are met: + * + * * Redistributions of source code must retain the above copyright + * notice, this list of conditions and the following disclaimer. + * * Redistributions in binary form must reproduce the above + * copyright notice, this list of conditions and the following + * disclaimer in the documentation and/or other materials provided + * with the distribution. + * * Neither the name of the Willow Garage nor the names of its + * contributors may be used to endorse or promote products derived + * from this software without specific prior written permission. + * + * THIS SOFTWARE IS PROVIDED BY THE COPYRIGHT HOLDERS AND CONTRIBUTORS + * "AS IS" AND ANY EXPRESS OR IMPLIED WARRANTIES, INCLUDING, BUT NOT + * LIMITED TO, THE IMPLIED WARRANTIES OF MERCHANTABILITY AND FITNESS + * FOR A PARTICULAR PURPOSE ARE DISCLAIMED. IN NO EVENT SHALL THE + * COPYRIGHT OWNER OR CONTRIBUTORS BE LIABLE FOR ANY DIRECT, INDIRECT, + * INCIDENTAL, SPECIAL, EXEMPLARY, OR CONSEQUENTIAL DAMAGES (INCLUDING, + * BUT NOT LIMITED TO, PROCUREMENT OF SUBSTITUTE GOODS OR SERVICES; + * LOSS OF USE, DATA, OR PROFITS; OR BUSINESS INTERRUPTION) HOWEVER + * CAUSED AND ON ANY THEORY OF LIABILITY, WHETHER IN CONTRACT, STRICT + * LIABILITY, OR TORT (INCLUDING NEGLIGENCE OR OTHERWISE) ARISING IN + * ANY WAY OUT OF THE USE OF THIS SOFTWARE, EVEN IF ADVISED OF THE + * POSSIBILITY OF SUCH DAMAGE. + *********************************************************************/ + +/* + BRISK - Binary Robust Invariant Scalable Keypoints + Reference implementation of + [1] Stefan Leutenegger,Margarita Chli and Roland Siegwart, BRISK: + Binary Robust Invariant Scalable Keypoints, in Proceedings of + the IEEE International Conference on Computer Vision (ICCV2011). + */ + +#include +#include +#include +#include +#include + +#include "fast_score.hpp" + +namespace cv +{ + +// a layer in the Brisk detector pyramid +class CV_EXPORTS BriskLayer +{ +public: + // constructor arguments + struct CV_EXPORTS CommonParams + { + static const int HALFSAMPLE = 0; + static const int TWOTHIRDSAMPLE = 1; + }; + // construct a base layer + BriskLayer(const cv::Mat& img, float scale = 1.0f, float offset = 0.0f); + // derive a layer + BriskLayer(const BriskLayer& layer, int mode); + + // Fast/Agast without non-max suppression + void + getAgastPoints(uint8_t threshold, std::vector& keypoints); + + // get scores - attention, this is in layer coordinates, not scale=1 coordinates! + inline uint8_t + getAgastScore(int x, int y, uint8_t threshold); + inline uint8_t + getAgastScore_5_8(int x, int y, uint8_t threshold); + inline uint8_t + getAgastScore(float xf, float yf, uint8_t threshold, float scale = 1.0f); + + // accessors + inline const cv::Mat& + img() const + { + return img_; + } + inline const cv::Mat& + scores() const + { + return scores_; + } + inline float + scale() const + { + return scale_; + } + inline float + offset() const + { + return offset_; + } + + // half sampling + static inline void + halfsample(const cv::Mat& srcimg, cv::Mat& dstimg); + // two third sampling + static inline void + twothirdsample(const cv::Mat& srcimg, cv::Mat& dstimg); + +private: + // access gray values (smoothed/interpolated) + __inline__ uint8_t + value(const cv::Mat& mat, float xf, float yf, float scale); + // the image + cv::Mat img_; + // its Fast scores + cv::Mat_ scores_; + // coordinate transformation + float scale_; + float offset_; + // agast + cv::Ptr fast_9_16_; + int pixel_5_8_[25]; + int pixel_9_16_[25]; +}; + +class CV_EXPORTS BriskScaleSpace +{ +public: + // construct telling the octaves number: + BriskScaleSpace(uint8_t _octaves = 3); + ~BriskScaleSpace(); + + // construct the image pyramids + void + constructPyramid(const cv::Mat& image); + + // get Keypoints + void + getKeypoints(const uint8_t _threshold, std::vector& keypoints); + +protected: + // nonmax suppression: + __inline__ bool + isMax2D(const uint8_t layer, const int x_layer, const int y_layer); + // 1D (scale axis) refinement: + __inline__ float + refine1D(const float s_05, const float s0, const float s05, float& max); // around octave + __inline__ float + refine1D_1(const float s_05, const float s0, const float s05, float& max); // around intra + __inline__ float + refine1D_2(const float s_05, const float s0, const float s05, float& max); // around octave 0 only + // 2D maximum refinement: + __inline__ float + subpixel2D(const int s_0_0, const int s_0_1, const int s_0_2, const int s_1_0, const int s_1_1, const int s_1_2, + const int s_2_0, const int s_2_1, const int s_2_2, float& delta_x, float& delta_y); + // 3D maximum refinement centered around (x_layer,y_layer) + __inline__ float + refine3D(const uint8_t layer, const int x_layer, const int y_layer, float& x, float& y, float& scale, bool& ismax); + + // interpolated score access with recalculation when needed: + __inline__ int + getScoreAbove(const uint8_t layer, const int x_layer, const int y_layer); + __inline__ int + getScoreBelow(const uint8_t layer, const int x_layer, const int y_layer); + + // return the maximum of score patches above or below + __inline__ float + getScoreMaxAbove(const uint8_t layer, const int x_layer, const int y_layer, const int threshold, bool& ismax, + float& dx, float& dy); + __inline__ float + getScoreMaxBelow(const uint8_t layer, const int x_layer, const int y_layer, const int threshold, bool& ismax, + float& dx, float& dy); + + // the image pyramids: + uint8_t layers_; + std::vector pyramid_; + + // some constant parameters: + static const float safetyFactor_; + static const float basicSize_; +}; + +using namespace cv; + +const float BRISK::basicSize_ = 12.0; +const unsigned int BRISK::scales_ = 64; +const float BRISK::scalerange_ = 30; // 40->4 Octaves - else, this needs to be adjusted... +const unsigned int BRISK::n_rot_ = 1024; // discretization of the rotation look-up + +const float BriskScaleSpace::safetyFactor_ = 1.0; +const float BriskScaleSpace::basicSize_ = 12.0; + +// constructors +BRISK::BRISK(int thresh, int octaves_in, float patternScale) +{ + threshold = thresh; + octaves = octaves_in; + + std::vector rList; + std::vector nList; + + // this is the standard pattern found to be suitable also + rList.resize(5); + nList.resize(5); + const double f = 0.85 * patternScale; + + rList[0] = f * 0; + rList[1] = f * 2.9; + rList[2] = f * 4.9; + rList[3] = f * 7.4; + rList[4] = f * 10.8; + + nList[0] = 1; + nList[1] = 10; + nList[2] = 14; + nList[3] = 15; + nList[4] = 20; + + generateKernel(rList, nList, 5.85 * patternScale, 8.2 * patternScale); + +} +BRISK::BRISK(std::vector &radiusList, std::vector &numberList, float dMax, float dMin, + std::vector indexChange) +{ + generateKernel(radiusList, numberList, dMax, dMin, indexChange); +} + +void +BRISK::generateKernel(std::vector &radiusList, std::vector &numberList, float dMax, + float dMin, std::vector indexChange) +{ + + dMax_ = dMax; + dMin_ = dMin; + + // get the total number of points + const int rings = radiusList.size(); + assert(radiusList.size()!=0&&radiusList.size()==numberList.size()); + points_ = 0; // remember the total number of points + for (int ring = 0; ring < rings; ring++) + { + points_ += numberList[ring]; + } + // set up the patterns + patternPoints_ = new BriskPatternPoint[points_ * scales_ * n_rot_]; + BriskPatternPoint* patternIterator = patternPoints_; + + // define the scale discretization: + static const float lb_scale = log(scalerange_) / log(2.0); + static const float lb_scale_step = lb_scale / (scales_); + + scaleList_ = new float[scales_]; + sizeList_ = new unsigned int[scales_]; + + const float sigma_scale = 1.3; + + for (unsigned int scale = 0; scale < scales_; ++scale) + { + scaleList_[scale] = pow((double) 2.0, (double) (scale * lb_scale_step)); + sizeList_[scale] = 0; + + // generate the pattern points look-up + double alpha, theta; + for (size_t rot = 0; rot < n_rot_; ++rot) + { + theta = double(rot) * 2 * M_PI / double(n_rot_); // this is the rotation of the feature + for (int ring = 0; ring < rings; ++ring) + { + for (int num = 0; num < numberList[ring]; ++num) + { + // the actual coordinates on the circle + alpha = (double(num)) * 2 * M_PI / double(numberList[ring]); + patternIterator->x = scaleList_[scale] * radiusList[ring] * cos(alpha + theta); // feature rotation plus angle of the point + patternIterator->y = scaleList_[scale] * radiusList[ring] * sin(alpha + theta); + // and the gaussian kernel sigma + if (ring == 0) + { + patternIterator->sigma = sigma_scale * scaleList_[scale] * 0.5; + } + else + { + patternIterator->sigma = sigma_scale * scaleList_[scale] * (double(radiusList[ring])) + * sin(M_PI / numberList[ring]); + } + // adapt the sizeList if necessary + const unsigned int size = ceil(((scaleList_[scale] * radiusList[ring]) + patternIterator->sigma)) + 1; + if (sizeList_[scale] < size) + { + sizeList_[scale] = size; + } + + // increment the iterator + ++patternIterator; + } + } + } + } + + // now also generate pairings + shortPairs_ = new BriskShortPair[points_ * (points_ - 1) / 2]; + longPairs_ = new BriskLongPair[points_ * (points_ - 1) / 2]; + noShortPairs_ = 0; + noLongPairs_ = 0; + + // fill indexChange with 0..n if empty + unsigned int indSize = indexChange.size(); + if (indSize == 0) + { + indexChange.resize(points_ * (points_ - 1) / 2); + indSize = indexChange.size(); + } + for (unsigned int i = 0; i < indSize; i++) + { + indexChange[i] = i; + } + const float dMin_sq = dMin_ * dMin_; + const float dMax_sq = dMax_ * dMax_; + for (unsigned int i = 1; i < points_; i++) + { + for (unsigned int j = 0; j < i; j++) + { //(find all the pairs) + // point pair distance: + const float dx = patternPoints_[j].x - patternPoints_[i].x; + const float dy = patternPoints_[j].y - patternPoints_[i].y; + const float norm_sq = (dx * dx + dy * dy); + if (norm_sq > dMin_sq) + { + // save to long pairs + BriskLongPair& longPair = longPairs_[noLongPairs_]; + longPair.weighted_dx = int((dx / (norm_sq)) * 2048.0 + 0.5); + longPair.weighted_dy = int((dy / (norm_sq)) * 2048.0 + 0.5); + longPair.i = i; + longPair.j = j; + ++noLongPairs_; + } + else if (norm_sq < dMax_sq) + { + // save to short pairs + assert(noShortPairs_ 2) + { + // now the calculation: + uchar* ptr = image.data + x_left + imagecols * y_top; + // first the corners: + ret_val = A * int(*ptr); + ptr += dx + 1; + ret_val += B * int(*ptr); + ptr += dy * imagecols + 1; + ret_val += C * int(*ptr); + ptr -= dx + 1; + ret_val += D * int(*ptr); + + // next the edges: + int* ptr_integral = (int*) integral.data + x_left + integralcols * y_top + 1; + // find a simple path through the different surface corners + const int tmp1 = (*ptr_integral); + ptr_integral += dx; + const int tmp2 = (*ptr_integral); + ptr_integral += integralcols; + const int tmp3 = (*ptr_integral); + ptr_integral++; + const int tmp4 = (*ptr_integral); + ptr_integral += dy * integralcols; + const int tmp5 = (*ptr_integral); + ptr_integral--; + const int tmp6 = (*ptr_integral); + ptr_integral += integralcols; + const int tmp7 = (*ptr_integral); + ptr_integral -= dx; + const int tmp8 = (*ptr_integral); + ptr_integral -= integralcols; + const int tmp9 = (*ptr_integral); + ptr_integral--; + const int tmp10 = (*ptr_integral); + ptr_integral -= dy * integralcols; + const int tmp11 = (*ptr_integral); + ptr_integral++; + const int tmp12 = (*ptr_integral); + + // assign the weighted surface integrals: + const int upper = (tmp3 - tmp2 + tmp1 - tmp12) * r_y_1_i; + const int middle = (tmp6 - tmp3 + tmp12 - tmp9) * scaling; + const int left = (tmp9 - tmp12 + tmp11 - tmp10) * r_x_1_i; + const int right = (tmp5 - tmp4 + tmp3 - tmp6) * r_x1_i; + const int bottom = (tmp7 - tmp6 + tmp9 - tmp8) * r_y1_i; + + return (ret_val + upper + middle + left + right + bottom + scaling2 / 2) / scaling2; + } + + // now the calculation: + uchar* ptr = image.data + x_left + imagecols * y_top; + // first row: + ret_val = A * int(*ptr); + ptr++; + const uchar* end1 = ptr + dx; + for (; ptr < end1; ptr++) + { + ret_val += r_y_1_i * int(*ptr); + } + ret_val += B * int(*ptr); + // middle ones: + ptr += imagecols - dx - 1; + uchar* end_j = ptr + dy * imagecols; + for (; ptr < end_j; ptr += imagecols - dx - 1) + { + ret_val += r_x_1_i * int(*ptr); + ptr++; + const uchar* end2 = ptr + dx; + for (; ptr < end2; ptr++) + { + ret_val += int(*ptr) * scaling; + } + ret_val += r_x1_i * int(*ptr); + } + // last row: + ret_val += D * int(*ptr); + ptr++; + const uchar* end3 = ptr + dx; + for (; ptr < end3; ptr++) + { + ret_val += r_y1_i * int(*ptr); + } + ret_val += C * int(*ptr); + + return (ret_val + scaling2 / 2) / scaling2; +} + +inline bool +RoiPredicate(const float minX, const float minY, const float maxX, const float maxY, const KeyPoint& keyPt) +{ + const Point2f& pt = keyPt.pt; + return (pt.x < minX) || (pt.x >= maxX) || (pt.y < minY) || (pt.y >= maxY); +} + +// computes the descriptor +void +BRISK::operator()( InputArray _image, InputArray _mask, vector& keypoints, + OutputArray _descriptors, bool useProvidedKeypoints) const +{ + Mat image = _image.getMat(), mask = _mask.getMat(); + if (!useProvidedKeypoints) + detectImpl(image, keypoints, mask); + + //Remove keypoints very close to the border + size_t ksize = keypoints.size(); + std::vector kscales; // remember the scale per keypoint + kscales.resize(ksize); + static const float log2 = 0.693147180559945; + static const float lb_scalerange = log(scalerange_) / (log2); + std::vector::iterator beginning = keypoints.begin(); + std::vector::iterator beginningkscales = kscales.begin(); + static const float basicSize06 = basicSize_ * 0.6; + for (size_t k = 0; k < ksize; k++) + { + unsigned int scale; + scale = std::max((int) (scales_ / lb_scalerange * (log(keypoints[k].size / (basicSize06)) / log2) + 0.5), 0); + // saturate + if (scale >= scales_) + scale = scales_ - 1; + kscales[k] = scale; + const int border = sizeList_[scale]; + const int border_x = image.cols - border; + const int border_y = image.rows - border; + if (RoiPredicate(border, border, border_x, border_y, keypoints[k])) + { + keypoints.erase(beginning + k); + kscales.erase(beginningkscales + k); + if (k == 0) + { + beginning = keypoints.begin(); + beginningkscales = kscales.begin(); + } + ksize--; + k--; + } + } + + // first, calculate the integral image over the whole image: + // current integral image + cv::Mat _integral; // the integral image + cv::integral(image, _integral); + + int* _values = new int[points_]; // for temporary use + + // resize the descriptors: + _descriptors.create(ksize, strings_, CV_8U); + cv::Mat descriptors = _descriptors.getMat(); + descriptors.setTo(0); + + // now do the extraction for all keypoints: + + // temporary variables containing gray values at sample points: + int t1; + int t2; + + // the feature orientation + uchar* ptr = descriptors.data; + for (size_t k = 0; k < ksize; k++) + { + int theta; + cv::KeyPoint& kp = keypoints[k]; + const int& scale = kscales[k]; + int shifter = 0; + int* pvalues = _values; + const float& x = kp.pt.x; + const float& y = kp.pt.y; + if (kp.angle==-1) + { + // don't compute the gradient direction, just assign a rotation of 0° + theta = 0; + } + else + { + theta = (int) (n_rot_ * (kp.angle / (360.0)) + 0.5); + if (theta < 0) + theta += n_rot_; + if (theta >= int(n_rot_)) + theta -= n_rot_; + } + + // now also extract the stuff for the actual direction: + // let us compute the smoothed values + shifter = 0; + + //unsigned int mean=0; + pvalues = _values; + // get the gray values in the rotated pattern + for (unsigned int i = 0; i < points_; i++) + { + *(pvalues++) = smoothedIntensity(image, _integral, x, y, scale, theta, i); + } + + // now iterate through all the pairings + unsigned int* ptr2 = (unsigned int*) ptr; + const BriskShortPair* max = shortPairs_ + noShortPairs_; + for (BriskShortPair* iter = shortPairs_; iter < max; ++iter) + { + t1 = *(_values + iter->i); + t2 = *(_values + iter->j); + if (t1 > t2) + { + *ptr2 |= ((1) << shifter); + + } // else already initialized with zero + // take care of the iterators: + ++shifter; + if (shifter == 32) + { + shifter = 0; + ++ptr2; + } + } + + ptr += strings_; + } + + // clean-up + _integral.release(); + delete[] _values; +} + +int +BRISK::descriptorSize() const +{ + return strings_; +} + +int +BRISK::descriptorType() const +{ + return CV_8U; +} + +BRISK::~BRISK() +{ + delete[] patternPoints_; + delete[] shortPairs_; + delete[] longPairs_; + delete[] scaleList_; + delete[] sizeList_; +} + +void +BRISK::operator()(InputArray _image, InputArray _mask, vector& keypoints) const +{ + Mat image = _image.getMat(), mask = _mask.getMat(); + if( image.type() != CV_8UC1 ) + cvtColor(_image, image, CV_BGR2GRAY); + + BriskScaleSpace briskScaleSpace(octaves); + briskScaleSpace.constructPyramid(image); + briskScaleSpace.getKeypoints(threshold, keypoints); + + // remove invalid points + removeInvalidPoints(mask, keypoints); + + // Compute the orientations of the keypoints + //Remove keypoints very close to the border + size_t ksize = keypoints.size(); + std::vector kscales; // remember the scale per keypoint + kscales.resize(ksize); + static const float log2 = 0.693147180559945; + static const float lb_scalerange = log(scalerange_) / (log2); + std::vector::iterator beginning = keypoints.begin(); + std::vector::iterator beginningkscales = kscales.begin(); + static const float basicSize06 = basicSize_ * 0.6; + for (size_t k = 0; k < ksize; k++) + { + unsigned int scale; + scale = std::max((int) (scales_ / lb_scalerange * (log(keypoints[k].size / (basicSize06)) / log2) + 0.5), 0); + // saturate + if (scale >= scales_) + scale = scales_ - 1; + kscales[k] = scale; + const int border = sizeList_[scale]; + const int border_x = image.cols - border; + const int border_y = image.rows - border; + if (RoiPredicate(border, border, border_x, border_y, keypoints[k])) + { + keypoints.erase(beginning + k); + kscales.erase(beginningkscales + k); + if (k == 0) + { + beginning = keypoints.begin(); + beginningkscales = kscales.begin(); + } + ksize--; + k--; + } + } + + // first, calculate the integral image over the whole image: + // current integral image + cv::Mat _integral; // the integral image + cv::integral(image, _integral); + + int* _values = new int[points_]; // for temporary use + + // now do the extraction for all keypoints: + + // temporary variables containing gray values at sample points: + int t1; + int t2; + + // the feature orientation + int direction0; + int direction1; + + for (size_t k = 0; k < ksize; k++) + { + cv::KeyPoint& kp = keypoints[k]; + const int& scale = kscales[k]; + int* pvalues = _values; + const float& x = kp.pt.x; + const float& y = kp.pt.y; + // get the gray values in the unrotated pattern + for (unsigned int i = 0; i < points_; i++) + { + *(pvalues++) = smoothedIntensity(image, _integral, x, y, scale, 0, i); + } + + direction0 = 0; + direction1 = 0; + // now iterate through the long pairings + const BriskLongPair* max = longPairs_ + noLongPairs_; + for (BriskLongPair* iter = longPairs_; iter < max; ++iter) + { + t1 = *(_values + iter->i); + t2 = *(_values + iter->j); + const int delta_t = (t1 - t2); + // update the direction: + const int tmp0 = delta_t * (iter->weighted_dx) / 1024; + const int tmp1 = delta_t * (iter->weighted_dy) / 1024; + direction0 += tmp0; + direction1 += tmp1; + } + kp.angle = atan2((float) direction1, (float) direction0) / M_PI * 180.0; + if (kp.angle < 0) + kp.angle += 360; + } +} + + +void +BRISK::detectImpl( const Mat& image, vector& keypoints, const Mat& mask) const +{ + (*this)(image, mask, keypoints); +} + +void +BRISK::computeImpl( const Mat& image, vector& keypoints, Mat& descriptors) const +{ + (*this)(image, Mat(), keypoints, descriptors, true); +} + +// construct telling the octaves number: +BriskScaleSpace::BriskScaleSpace(uint8_t _octaves) +{ + if (_octaves == 0) + layers_ = 1; + else + layers_ = 2 * _octaves; +} +BriskScaleSpace::~BriskScaleSpace() +{ + +} +// construct the image pyramids +void +BriskScaleSpace::constructPyramid(const cv::Mat& image) +{ + + // set correct size: + pyramid_.clear(); + + // fill the pyramid: + pyramid_.push_back(BriskLayer(image.clone())); + if (layers_ > 1) + { + pyramid_.push_back(BriskLayer(pyramid_.back(), BriskLayer::CommonParams::TWOTHIRDSAMPLE)); + } + const int octaves2 = layers_; + + for (uint8_t i = 2; i < octaves2; i += 2) + { + pyramid_.push_back(BriskLayer(pyramid_[i - 2], BriskLayer::CommonParams::HALFSAMPLE)); + pyramid_.push_back(BriskLayer(pyramid_[i - 1], BriskLayer::CommonParams::HALFSAMPLE)); + } +} + +void +BriskScaleSpace::getKeypoints(const uint8_t _threshold, std::vector& keypoints) +{ + // make sure keypoints is empty + keypoints.resize(0); + keypoints.reserve(2000); + + // assign thresholds + uint8_t threshold_ = _threshold; + uint8_t safeThreshold_ = threshold_ * safetyFactor_; + std::vector > agastPoints; + agastPoints.resize(layers_); + + // go through the octaves and intra layers and calculate fast corner scores: + for (uint8_t i = 0; i < layers_; i++) + { + // call OAST16_9 without nms + BriskLayer& l = pyramid_[i]; + l.getAgastPoints(safeThreshold_, agastPoints[i]); + } + + if (layers_ == 1) + { + // just do a simple 2d subpixel refinement... + const int num = agastPoints[0].size(); + for (int n = 0; n < num; n++) + { + const cv::Point2f& point = agastPoints.at(0)[n].pt; + // first check if it is a maximum: + if (!isMax2D(0, point.x, point.y)) + continue; + + // let's do the subpixel and float scale refinement: + BriskLayer& l = pyramid_[0]; + register int s_0_0 = l.getAgastScore(point.x - 1, point.y - 1, 1); + register int s_1_0 = l.getAgastScore(point.x, point.y - 1, 1); + register int s_2_0 = l.getAgastScore(point.x + 1, point.y - 1, 1); + register int s_2_1 = l.getAgastScore(point.x + 1, point.y, 1); + register int s_1_1 = l.getAgastScore(point.x, point.y, 1); + register int s_0_1 = l.getAgastScore(point.x - 1, point.y, 1); + register int s_0_2 = l.getAgastScore(point.x - 1, point.y + 1, 1); + register int s_1_2 = l.getAgastScore(point.x, point.y + 1, 1); + register int s_2_2 = l.getAgastScore(point.x + 1, point.y + 1, 1); + float delta_x, delta_y; + float max = subpixel2D(s_0_0, s_0_1, s_0_2, s_1_0, s_1_1, s_1_2, s_2_0, s_2_1, s_2_2, delta_x, delta_y); + + // store: + keypoints.push_back(cv::KeyPoint(float(point.x) + delta_x, float(point.y) + delta_y, basicSize_, -1, max, 0)); + + } + + return; + } + + float x, y, scale, score; + for (uint8_t i = 0; i < layers_; i++) + { + BriskLayer& l = pyramid_[i]; + const int num = agastPoints[i].size(); + if (i == layers_ - 1) + { + for (int n = 0; n < num; n++) + { + const cv::Point2f& point = agastPoints.at(i)[n].pt; + // consider only 2D maxima... + if (!isMax2D(i, point.x, point.y)) + continue; + + bool ismax; + float dx, dy; + getScoreMaxBelow(i, point.x, point.y, l.getAgastScore(point.x, point.y, safeThreshold_), ismax, dx, dy); + if (!ismax) + continue; + + // get the patch on this layer: + register int s_0_0 = l.getAgastScore(point.x - 1, point.y - 1, 1); + register int s_1_0 = l.getAgastScore(point.x, point.y - 1, 1); + register int s_2_0 = l.getAgastScore(point.x + 1, point.y - 1, 1); + register int s_2_1 = l.getAgastScore(point.x + 1, point.y, 1); + register int s_1_1 = l.getAgastScore(point.x, point.y, 1); + register int s_0_1 = l.getAgastScore(point.x - 1, point.y, 1); + register int s_0_2 = l.getAgastScore(point.x - 1, point.y + 1, 1); + register int s_1_2 = l.getAgastScore(point.x, point.y + 1, 1); + register int s_2_2 = l.getAgastScore(point.x + 1, point.y + 1, 1); + float delta_x, delta_y; + float max = subpixel2D(s_0_0, s_0_1, s_0_2, s_1_0, s_1_1, s_1_2, s_2_0, s_2_1, s_2_2, delta_x, delta_y); + + // store: + keypoints.push_back( + cv::KeyPoint((float(point.x) + delta_x) * l.scale() + l.offset(), + (float(point.y) + delta_y) * l.scale() + l.offset(), basicSize_ * l.scale(), -1, max, i)); + } + } + else + { + // not the last layer: + for (int n = 0; n < num; n++) + { + const cv::Point2f& point = agastPoints.at(i)[n].pt; + + // first check if it is a maximum: + if (!isMax2D(i, point.x, point.y)) + continue; + + // let's do the subpixel and float scale refinement: + bool ismax; + score = refine3D(i, point.x, point.y, x, y, scale, ismax); + if (!ismax) + { + continue; + } + + // finally store the detected keypoint: + if (score > float(threshold_)) + { + keypoints.push_back(cv::KeyPoint(x, y, basicSize_ * scale, -1, score, i)); + } + } + } + } +} + +// interpolated score access with recalculation when needed: +__inline__ int +BriskScaleSpace::getScoreAbove(const uint8_t layer, const int x_layer, const int y_layer) +{ + assert(layer delta; + // put together a list of 2d-offsets to where the maximum is also reached + if (center == s_1_1) + { + delta.push_back(-1); + delta.push_back(-1); + } + if (center == s0_1) + { + delta.push_back(0); + delta.push_back(-1); + } + if (center == s1_1) + { + delta.push_back(1); + delta.push_back(-1); + } + if (center == s_10) + { + delta.push_back(-1); + delta.push_back(0); + } + if (center == s10) + { + delta.push_back(1); + delta.push_back(0); + } + if (center == s_11) + { + delta.push_back(-1); + delta.push_back(1); + } + if (center == s01) + { + delta.push_back(0); + delta.push_back(1); + } + if (center == s11) + { + delta.push_back(1); + delta.push_back(1); + } + const unsigned int deltasize = delta.size(); + if (deltasize != 0) + { + // in this case, we have to analyze the situation more carefully: + // the values are gaussian blurred and then we really decide + data = scores.data + y_layer * scorescols + x_layer; + int smoothedcenter = 4 * center + 2 * (s_10 + s10 + s0_1 + s01) + s_1_1 + s1_1 + s_11 + s11; + for (unsigned int i = 0; i < deltasize; i += 2) + { + data = scores.data + (y_layer - 1 + delta[i + 1]) * scorescols + x_layer + delta[i] - 1; + int othercenter = *data; + data++; + othercenter += 2 * (*data); + data++; + othercenter += *data; + data += scorescols; + othercenter += 2 * (*data); + data--; + othercenter += 4 * (*data); + data--; + othercenter += 2 * (*data); + data += scorescols; + othercenter += *data; + data++; + othercenter += 2 * (*data); + data++; + othercenter += *data; + if (othercenter > smoothedcenter) + return false; + } + } + return true; +} + +// 3D maximum refinement centered around (x_layer,y_layer) +__inline__ float +BriskScaleSpace::refine3D(const uint8_t layer, const int x_layer, const int y_layer, float& x, float& y, float& scale, + bool& ismax) +{ + ismax = true; + BriskLayer& thisLayer = pyramid_[layer]; + const int center = thisLayer.getAgastScore(x_layer, y_layer, 1); + + // check and get above maximum: + float delta_x_above, delta_y_above; + float max_above = getScoreMaxAbove(layer, x_layer, y_layer, center, ismax, delta_x_above, delta_y_above); + + if (!ismax) + return 0.0; + + float max; // to be returned + + if (layer % 2 == 0) + { // on octave + // treat the patch below: + float delta_x_below, delta_y_below; + float max_below_float; + uchar max_below_uchar = 0; + if (layer == 0) + { + // guess the lower intra octave... + BriskLayer& l = pyramid_[0]; + register int s_0_0 = l.getAgastScore_5_8(x_layer - 1, y_layer - 1, 1); + max_below_uchar = s_0_0; + register int s_1_0 = l.getAgastScore_5_8(x_layer, y_layer - 1, 1); + if (s_1_0 > max_below_uchar) + max_below_uchar = s_1_0; + register int s_2_0 = l.getAgastScore_5_8(x_layer + 1, y_layer - 1, 1); + if (s_2_0 > max_below_uchar) + max_below_uchar = s_2_0; + register int s_2_1 = l.getAgastScore_5_8(x_layer + 1, y_layer, 1); + if (s_2_1 > max_below_uchar) + max_below_uchar = s_2_1; + register int s_1_1 = l.getAgastScore_5_8(x_layer, y_layer, 1); + if (s_1_1 > max_below_uchar) + max_below_uchar = s_1_1; + register int s_0_1 = l.getAgastScore_5_8(x_layer - 1, y_layer, 1); + if (s_0_1 > max_below_uchar) + max_below_uchar = s_0_1; + register int s_0_2 = l.getAgastScore_5_8(x_layer - 1, y_layer + 1, 1); + if (s_0_2 > max_below_uchar) + max_below_uchar = s_0_2; + register int s_1_2 = l.getAgastScore_5_8(x_layer, y_layer + 1, 1); + if (s_1_2 > max_below_uchar) + max_below_uchar = s_1_2; + register int s_2_2 = l.getAgastScore_5_8(x_layer + 1, y_layer + 1, 1); + if (s_2_2 > max_below_uchar) + max_below_uchar = s_2_2; + + max_below_float = subpixel2D(s_0_0, s_0_1, s_0_2, s_1_0, s_1_1, s_1_2, s_2_0, s_2_1, s_2_2, delta_x_below, + delta_y_below); + max_below_float = max_below_uchar; + } + else + { + max_below_float = getScoreMaxBelow(layer, x_layer, y_layer, center, ismax, delta_x_below, delta_y_below); + if (!ismax) + return 0; + } + + // get the patch on this layer: + register int s_0_0 = thisLayer.getAgastScore(x_layer - 1, y_layer - 1, 1); + register int s_1_0 = thisLayer.getAgastScore(x_layer, y_layer - 1, 1); + register int s_2_0 = thisLayer.getAgastScore(x_layer + 1, y_layer - 1, 1); + register int s_2_1 = thisLayer.getAgastScore(x_layer + 1, y_layer, 1); + register int s_1_1 = thisLayer.getAgastScore(x_layer, y_layer, 1); + register int s_0_1 = thisLayer.getAgastScore(x_layer - 1, y_layer, 1); + register int s_0_2 = thisLayer.getAgastScore(x_layer - 1, y_layer + 1, 1); + register int s_1_2 = thisLayer.getAgastScore(x_layer, y_layer + 1, 1); + register int s_2_2 = thisLayer.getAgastScore(x_layer + 1, y_layer + 1, 1); + float delta_x_layer, delta_y_layer; + float max_layer = subpixel2D(s_0_0, s_0_1, s_0_2, s_1_0, s_1_1, s_1_2, s_2_0, s_2_1, s_2_2, delta_x_layer, + delta_y_layer); + + // calculate the relative scale (1D maximum): + if (layer == 0) + { + scale = refine1D_2(max_below_float, std::max(float(center), max_layer), max_above, max); + } + else + scale = refine1D(max_below_float, std::max(float(center), max_layer), max_above, max); + + if (scale > 1.0) + { + // interpolate the position: + const float r0 = (1.5 - scale) / .5; + const float r1 = 1.0 - r0; + x = (r0 * delta_x_layer + r1 * delta_x_above + float(x_layer)) * thisLayer.scale() + thisLayer.offset(); + y = (r0 * delta_y_layer + r1 * delta_y_above + float(y_layer)) * thisLayer.scale() + thisLayer.offset(); + } + else + { + if (layer == 0) + { + // interpolate the position: + const float r0 = (scale - 0.5) / 0.5; + const float r_1 = 1.0 - r0; + x = r0 * delta_x_layer + r_1 * delta_x_below + float(x_layer); + y = r0 * delta_y_layer + r_1 * delta_y_below + float(y_layer); + } + else + { + // interpolate the position: + const float r0 = (scale - 0.75) / 0.25; + const float r_1 = 1.0 - r0; + x = (r0 * delta_x_layer + r_1 * delta_x_below + float(x_layer)) * thisLayer.scale() + thisLayer.offset(); + y = (r0 * delta_y_layer + r_1 * delta_y_below + float(y_layer)) * thisLayer.scale() + thisLayer.offset(); + } + } + } + else + { + // on intra + // check the patch below: + float delta_x_below, delta_y_below; + float max_below = getScoreMaxBelow(layer, x_layer, y_layer, center, ismax, delta_x_below, delta_y_below); + if (!ismax) + return 0.0; + + // get the patch on this layer: + register int s_0_0 = thisLayer.getAgastScore(x_layer - 1, y_layer - 1, 1); + register int s_1_0 = thisLayer.getAgastScore(x_layer, y_layer - 1, 1); + register int s_2_0 = thisLayer.getAgastScore(x_layer + 1, y_layer - 1, 1); + register int s_2_1 = thisLayer.getAgastScore(x_layer + 1, y_layer, 1); + register int s_1_1 = thisLayer.getAgastScore(x_layer, y_layer, 1); + register int s_0_1 = thisLayer.getAgastScore(x_layer - 1, y_layer, 1); + register int s_0_2 = thisLayer.getAgastScore(x_layer - 1, y_layer + 1, 1); + register int s_1_2 = thisLayer.getAgastScore(x_layer, y_layer + 1, 1); + register int s_2_2 = thisLayer.getAgastScore(x_layer + 1, y_layer + 1, 1); + float delta_x_layer, delta_y_layer; + float max_layer = subpixel2D(s_0_0, s_0_1, s_0_2, s_1_0, s_1_1, s_1_2, s_2_0, s_2_1, s_2_2, delta_x_layer, + delta_y_layer); + + // calculate the relative scale (1D maximum): + scale = refine1D_1(max_below, std::max(float(center), max_layer), max_above, max); + if (scale > 1.0) + { + // interpolate the position: + const float r0 = 4.0 - scale * 3.0; + const float r1 = 1.0 - r0; + x = (r0 * delta_x_layer + r1 * delta_x_above + float(x_layer)) * thisLayer.scale() + thisLayer.offset(); + y = (r0 * delta_y_layer + r1 * delta_y_above + float(y_layer)) * thisLayer.scale() + thisLayer.offset(); + } + else + { + // interpolate the position: + const float r0 = scale * 3.0 - 2.0; + const float r_1 = 1.0 - r0; + x = (r0 * delta_x_layer + r_1 * delta_x_below + float(x_layer)) * thisLayer.scale() + thisLayer.offset(); + y = (r0 * delta_y_layer + r_1 * delta_y_below + float(y_layer)) * thisLayer.scale() + thisLayer.offset(); + } + } + + // calculate the absolute scale: + scale *= thisLayer.scale(); + + // that's it, return the refined maximum: + return max; +} + +// return the maximum of score patches above or below +__inline__ float +BriskScaleSpace::getScoreMaxAbove(const uint8_t layer, const int x_layer, const int y_layer, const int threshold, + bool& ismax, float& dx, float& dy) +{ + + ismax = false; + // relevant floating point coordinates + float x_1; + float x1; + float y_1; + float y1; + + // the layer above + assert(layer+1 threshold) + return 0; + for (int x = x_1 + 1; x <= int(x1); x++) + { + tmp_max = layerAbove.getAgastScore(float(x), y_1, 1); + if (tmp_max > threshold) + return 0; + if (tmp_max > max) + { + max = tmp_max; + max_x = x; + } + } + tmp_max = layerAbove.getAgastScore(x1, y_1, 1); + if (tmp_max > threshold) + return 0; + if (tmp_max > max) + { + max = tmp_max; + max_x = int(x1); + } + + // middle rows + for (int y = y_1 + 1; y <= int(y1); y++) + { + tmp_max = layerAbove.getAgastScore(x_1, float(y), 1); + if (tmp_max > threshold) + return 0; + if (tmp_max > max) + { + max = tmp_max; + max_x = int(x_1 + 1); + max_y = y; + } + for (int x = x_1 + 1; x <= int(x1); x++) + { + tmp_max = layerAbove.getAgastScore(x, y, 1); + if (tmp_max > threshold) + return 0; + if (tmp_max > max) + { + max = tmp_max; + max_x = x; + max_y = y; + } + } + tmp_max = layerAbove.getAgastScore(x1, float(y), 1); + if (tmp_max > threshold) + return 0; + if (tmp_max > max) + { + max = tmp_max; + max_x = int(x1); + max_y = y; + } + } + + // bottom row + tmp_max = layerAbove.getAgastScore(x_1, y1, 1); + if (tmp_max > max) + { + max = tmp_max; + max_x = int(x_1 + 1); + max_y = int(y1); + } + for (int x = x_1 + 1; x <= int(x1); x++) + { + tmp_max = layerAbove.getAgastScore(float(x), y1, 1); + if (tmp_max > max) + { + max = tmp_max; + max_x = x; + max_y = int(y1); + } + } + tmp_max = layerAbove.getAgastScore(x1, y1, 1); + if (tmp_max > max) + { + max = tmp_max; + max_x = int(x1); + max_y = int(y1); + } + + //find dx/dy: + register int s_0_0 = layerAbove.getAgastScore(max_x - 1, max_y - 1, 1); + register int s_1_0 = layerAbove.getAgastScore(max_x, max_y - 1, 1); + register int s_2_0 = layerAbove.getAgastScore(max_x + 1, max_y - 1, 1); + register int s_2_1 = layerAbove.getAgastScore(max_x + 1, max_y, 1); + register int s_1_1 = layerAbove.getAgastScore(max_x, max_y, 1); + register int s_0_1 = layerAbove.getAgastScore(max_x - 1, max_y, 1); + register int s_0_2 = layerAbove.getAgastScore(max_x - 1, max_y + 1, 1); + register int s_1_2 = layerAbove.getAgastScore(max_x, max_y + 1, 1); + register int s_2_2 = layerAbove.getAgastScore(max_x + 1, max_y + 1, 1); + float dx_1, dy_1; + float refined_max = subpixel2D(s_0_0, s_0_1, s_0_2, s_1_0, s_1_1, s_1_2, s_2_0, s_2_1, s_2_2, dx_1, dy_1); + + // calculate dx/dy in above coordinates + float real_x = float(max_x) + dx_1; + float real_y = float(max_y) + dy_1; + bool returnrefined = true; + if (layer % 2 == 0) + { + dx = (real_x * 6.0f + 1.0f) / 4.0f - float(x_layer); + dy = (real_y * 6.0f + 1.0f) / 4.0f - float(y_layer); + } + else + { + dx = (real_x * 8.0 + 1.0) / 6.0 - float(x_layer); + dy = (real_y * 8.0 + 1.0) / 6.0 - float(y_layer); + } + + // saturate + if (dx > 1.0f) + { + dx = 1.0f; + returnrefined = false; + } + if (dx < -1.0f) + { + dx = -1.0f; + returnrefined = false; + } + if (dy > 1.0f) + { + dy = 1.0f; + returnrefined = false; + } + if (dy < -1.0f) + { + dy = -1.0f; + returnrefined = false; + } + + // done and ok. + ismax = true; + if (returnrefined) + { + return std::max(refined_max, max); + } + return max; +} + +__inline__ float +BriskScaleSpace::getScoreMaxBelow(const uint8_t layer, const int x_layer, const int y_layer, const int threshold, + bool& ismax, float& dx, float& dy) +{ + ismax = false; + + // relevant floating point coordinates + float x_1; + float x1; + float y_1; + float y1; + + if (layer % 2 == 0) + { + // octave + x_1 = float(8 * (x_layer) + 1 - 4) / 6.0; + x1 = float(8 * (x_layer) + 1 + 4) / 6.0; + y_1 = float(8 * (y_layer) + 1 - 4) / 6.0; + y1 = float(8 * (y_layer) + 1 + 4) / 6.0; + } + else + { + x_1 = float(6 * (x_layer) + 1 - 3) / 4.0; + x1 = float(6 * (x_layer) + 1 + 3) / 4.0; + y_1 = float(6 * (y_layer) + 1 - 3) / 4.0; + y1 = float(6 * (y_layer) + 1 + 3) / 4.0; + } + + // the layer below + assert(layer>0); + BriskLayer& layerBelow = pyramid_[layer - 1]; + + // check the first row + int max_x = x_1 + 1; + int max_y = y_1 + 1; + float tmp_max; + float max = layerBelow.getAgastScore(x_1, y_1, 1); + if (max > threshold) + return 0; + for (int x = x_1 + 1; x <= int(x1); x++) + { + tmp_max = layerBelow.getAgastScore(float(x), y_1, 1); + if (tmp_max > threshold) + return 0; + if (tmp_max > max) + { + max = tmp_max; + max_x = x; + } + } + tmp_max = layerBelow.getAgastScore(x1, y_1, 1); + if (tmp_max > threshold) + return 0; + if (tmp_max > max) + { + max = tmp_max; + max_x = int(x1); + } + + // middle rows + for (int y = y_1 + 1; y <= int(y1); y++) + { + tmp_max = layerBelow.getAgastScore(x_1, float(y), 1); + if (tmp_max > threshold) + return 0; + if (tmp_max > max) + { + max = tmp_max; + max_x = int(x_1 + 1); + max_y = y; + } + for (int x = x_1 + 1; x <= int(x1); x++) + { + tmp_max = layerBelow.getAgastScore(x, y, 1); + if (tmp_max > threshold) + return 0; + if (tmp_max == max) + { + const int t1 = 2 + * (layerBelow.getAgastScore(x - 1, y, 1) + layerBelow.getAgastScore(x + 1, y, 1) + + layerBelow.getAgastScore(x, y + 1, 1) + layerBelow.getAgastScore(x, y - 1, 1)) + + (layerBelow.getAgastScore(x + 1, y + 1, 1) + layerBelow.getAgastScore(x - 1, y + 1, 1) + + layerBelow.getAgastScore(x + 1, y - 1, 1) + layerBelow.getAgastScore(x - 1, y - 1, 1)); + const int t2 = 2 + * (layerBelow.getAgastScore(max_x - 1, max_y, 1) + layerBelow.getAgastScore(max_x + 1, max_y, 1) + + layerBelow.getAgastScore(max_x, max_y + 1, 1) + layerBelow.getAgastScore(max_x, max_y - 1, 1)) + + (layerBelow.getAgastScore(max_x + 1, max_y + 1, 1) + layerBelow.getAgastScore(max_x - 1, + max_y + 1, 1) + + layerBelow.getAgastScore(max_x + 1, max_y - 1, 1) + + layerBelow.getAgastScore(max_x - 1, max_y - 1, 1)); + if (t1 > t2) + { + max_x = x; + max_y = y; + } + } + if (tmp_max > max) + { + max = tmp_max; + max_x = x; + max_y = y; + } + } + tmp_max = layerBelow.getAgastScore(x1, float(y), 1); + if (tmp_max > threshold) + return 0; + if (tmp_max > max) + { + max = tmp_max; + max_x = int(x1); + max_y = y; + } + } + + // bottom row + tmp_max = layerBelow.getAgastScore(x_1, y1, 1); + if (tmp_max > max) + { + max = tmp_max; + max_x = int(x_1 + 1); + max_y = int(y1); + } + for (int x = x_1 + 1; x <= int(x1); x++) + { + tmp_max = layerBelow.getAgastScore(float(x), y1, 1); + if (tmp_max > max) + { + max = tmp_max; + max_x = x; + max_y = int(y1); + } + } + tmp_max = layerBelow.getAgastScore(x1, y1, 1); + if (tmp_max > max) + { + max = tmp_max; + max_x = int(x1); + max_y = int(y1); + } + + //find dx/dy: + register int s_0_0 = layerBelow.getAgastScore(max_x - 1, max_y - 1, 1); + register int s_1_0 = layerBelow.getAgastScore(max_x, max_y - 1, 1); + register int s_2_0 = layerBelow.getAgastScore(max_x + 1, max_y - 1, 1); + register int s_2_1 = layerBelow.getAgastScore(max_x + 1, max_y, 1); + register int s_1_1 = layerBelow.getAgastScore(max_x, max_y, 1); + register int s_0_1 = layerBelow.getAgastScore(max_x - 1, max_y, 1); + register int s_0_2 = layerBelow.getAgastScore(max_x - 1, max_y + 1, 1); + register int s_1_2 = layerBelow.getAgastScore(max_x, max_y + 1, 1); + register int s_2_2 = layerBelow.getAgastScore(max_x + 1, max_y + 1, 1); + float dx_1, dy_1; + float refined_max = subpixel2D(s_0_0, s_0_1, s_0_2, s_1_0, s_1_1, s_1_2, s_2_0, s_2_1, s_2_2, dx_1, dy_1); + + // calculate dx/dy in above coordinates + float real_x = float(max_x) + dx_1; + float real_y = float(max_y) + dy_1; + bool returnrefined = true; + if (layer % 2 == 0) + { + dx = (real_x * 6.0 + 1.0) / 8.0 - float(x_layer); + dy = (real_y * 6.0 + 1.0) / 8.0 - float(y_layer); + } + else + { + dx = (real_x * 4.0 - 1.0) / 6.0 - float(x_layer); + dy = (real_y * 4.0 - 1.0) / 6.0 - float(y_layer); + } + + // saturate + if (dx > 1.0) + { + dx = 1.0; + returnrefined = false; + } + if (dx < -1.0) + { + dx = -1.0; + returnrefined = false; + } + if (dy > 1.0) + { + dy = 1.0; + returnrefined = false; + } + if (dy < -1.0) + { + dy = -1.0; + returnrefined = false; + } + + // done and ok. + ismax = true; + if (returnrefined) + { + return std::max(refined_max, max); + } + return max; +} + +__inline__ float +BriskScaleSpace::refine1D(const float s_05, const float s0, const float s05, float& max) +{ + int i_05 = int(1024.0 * s_05 + 0.5); + int i0 = int(1024.0 * s0 + 0.5); + int i05 = int(1024.0 * s05 + 0.5); + + // 16.0000 -24.0000 8.0000 + // -40.0000 54.0000 -14.0000 + // 24.0000 -27.0000 6.0000 + + int three_a = 16 * i_05 - 24 * i0 + 8 * i05; + // second derivative must be negative: + if (three_a >= 0) + { + if (s0 >= s_05 && s0 >= s05) + { + max = s0; + return 1.0; + } + if (s_05 >= s0 && s_05 >= s05) + { + max = s_05; + return 0.75; + } + if (s05 >= s0 && s05 >= s_05) + { + max = s05; + return 1.5; + } + } + + int three_b = -40 * i_05 + 54 * i0 - 14 * i05; + // calculate max location: + float ret_val = -float(three_b) / float(2 * three_a); + // saturate and return + if (ret_val < 0.75) + ret_val = 0.75; + else if (ret_val > 1.5) + ret_val = 1.5; // allow to be slightly off bounds ...? + int three_c = +24 * i_05 - 27 * i0 + 6 * i05; + max = float(three_c) + float(three_a) * ret_val * ret_val + float(three_b) * ret_val; + max /= 3072.0; + return ret_val; +} + +__inline__ float +BriskScaleSpace::refine1D_1(const float s_05, const float s0, const float s05, float& max) +{ + int i_05 = int(1024.0 * s_05 + 0.5); + int i0 = int(1024.0 * s0 + 0.5); + int i05 = int(1024.0 * s05 + 0.5); + + // 4.5000 -9.0000 4.5000 + //-10.5000 18.0000 -7.5000 + // 6.0000 -8.0000 3.0000 + + int two_a = 9 * i_05 - 18 * i0 + 9 * i05; + // second derivative must be negative: + if (two_a >= 0) + { + if (s0 >= s_05 && s0 >= s05) + { + max = s0; + return 1.0; + } + if (s_05 >= s0 && s_05 >= s05) + { + max = s_05; + return 0.6666666666666666666666666667; + } + if (s05 >= s0 && s05 >= s_05) + { + max = s05; + return 1.3333333333333333333333333333; + } + } + + int two_b = -21 * i_05 + 36 * i0 - 15 * i05; + // calculate max location: + float ret_val = -float(two_b) / float(2 * two_a); + // saturate and return + if (ret_val < 0.6666666666666666666666666667) + ret_val = 0.666666666666666666666666667; + else if (ret_val > 1.33333333333333333333333333) + ret_val = 1.333333333333333333333333333; + int two_c = +12 * i_05 - 16 * i0 + 6 * i05; + max = float(two_c) + float(two_a) * ret_val * ret_val + float(two_b) * ret_val; + max /= 2048.0; + return ret_val; +} + +__inline__ float +BriskScaleSpace::refine1D_2(const float s_05, const float s0, const float s05, float& max) +{ + int i_05 = int(1024.0 * s_05 + 0.5); + int i0 = int(1024.0 * s0 + 0.5); + int i05 = int(1024.0 * s05 + 0.5); + + // 18.0000 -30.0000 12.0000 + // -45.0000 65.0000 -20.0000 + // 27.0000 -30.0000 8.0000 + + int a = 2 * i_05 - 4 * i0 + 2 * i05; + // second derivative must be negative: + if (a >= 0) + { + if (s0 >= s_05 && s0 >= s05) + { + max = s0; + return 1.0; + } + if (s_05 >= s0 && s_05 >= s05) + { + max = s_05; + return 0.7; + } + if (s05 >= s0 && s05 >= s_05) + { + max = s05; + return 1.5; + } + } + + int b = -5 * i_05 + 8 * i0 - 3 * i05; + // calculate max location: + float ret_val = -float(b) / float(2 * a); + // saturate and return + if (ret_val < 0.7) + ret_val = 0.7; + else if (ret_val > 1.5) + ret_val = 1.5; // allow to be slightly off bounds ...? + int c = +3 * i_05 - 3 * i0 + 1 * i05; + max = float(c) + float(a) * ret_val * ret_val + float(b) * ret_val; + max /= 1024; + return ret_val; +} + +__inline__ float +BriskScaleSpace::subpixel2D(const int s_0_0, const int s_0_1, const int s_0_2, const int s_1_0, const int s_1_1, + const int s_1_2, const int s_2_0, const int s_2_1, const int s_2_2, float& delta_x, + float& delta_y) +{ + + // the coefficients of the 2d quadratic function least-squares fit: + register int tmp1 = s_0_0 + s_0_2 - 2 * s_1_1 + s_2_0 + s_2_2; + register int coeff1 = 3 * (tmp1 + s_0_1 - ((s_1_0 + s_1_2) << 1) + s_2_1); + register int coeff2 = 3 * (tmp1 - ((s_0_1 + s_2_1) << 1) + s_1_0 + s_1_2); + register int tmp2 = s_0_2 - s_2_0; + register int tmp3 = (s_0_0 + tmp2 - s_2_2); + register int tmp4 = tmp3 - 2 * tmp2; + register int coeff3 = -3 * (tmp3 + s_0_1 - s_2_1); + register int coeff4 = -3 * (tmp4 + s_1_0 - s_1_2); + register int coeff5 = (s_0_0 - s_0_2 - s_2_0 + s_2_2) << 2; + register int coeff6 = -(s_0_0 + s_0_2 - ((s_1_0 + s_0_1 + s_1_2 + s_2_1) << 1) - 5 * s_1_1 + s_2_0 + s_2_2) << 1; + + // 2nd derivative test: + register int H_det = 4 * coeff1 * coeff2 - coeff5 * coeff5; + + if (H_det == 0) + { + delta_x = 0.0; + delta_y = 0.0; + return float(coeff6) / 18.0; + } + + if (!(H_det > 0 && coeff1 < 0)) + { + // The maximum must be at the one of the 4 patch corners. + int tmp_max = coeff3 + coeff4 + coeff5; + delta_x = 1.0; + delta_y = 1.0; + + int tmp = -coeff3 + coeff4 - coeff5; + if (tmp > tmp_max) + { + tmp_max = tmp; + delta_x = -1.0; + delta_y = 1.0; + } + tmp = coeff3 - coeff4 - coeff5; + if (tmp > tmp_max) + { + tmp_max = tmp; + delta_x = 1.0; + delta_y = -1.0; + } + tmp = -coeff3 - coeff4 + coeff5; + if (tmp > tmp_max) + { + tmp_max = tmp; + delta_x = -1.0; + delta_y = -1.0; + } + return float(tmp_max + coeff1 + coeff2 + coeff6) / 18.0; + } + + // this is hopefully the normal outcome of the Hessian test + delta_x = float(2 * coeff2 * coeff3 - coeff4 * coeff5) / float(-H_det); + delta_y = float(2 * coeff1 * coeff4 - coeff3 * coeff5) / float(-H_det); + // TODO: this is not correct, but easy, so perform a real boundary maximum search: + bool tx = false; + bool tx_ = false; + bool ty = false; + bool ty_ = false; + if (delta_x > 1.0) + tx = true; + else if (delta_x < -1.0) + tx_ = true; + if (delta_y > 1.0) + ty = true; + if (delta_y < -1.0) + ty_ = true; + + if (tx || tx_ || ty || ty_) + { + // get two candidates: + float delta_x1 = 0.0, delta_x2 = 0.0, delta_y1 = 0.0, delta_y2 = 0.0; + if (tx) + { + delta_x1 = 1.0; + delta_y1 = -float(coeff4 + coeff5) / float(2 * coeff2); + if (delta_y1 > 1.0) + delta_y1 = 1.0; + else if (delta_y1 < -1.0) + delta_y1 = -1.0; + } + else if (tx_) + { + delta_x1 = -1.0; + delta_y1 = -float(coeff4 - coeff5) / float(2 * coeff2); + if (delta_y1 > 1.0) + delta_y1 = 1.0; + else if (delta_y1 < -1.0) + delta_y1 = -1.0; + } + if (ty) + { + delta_y2 = 1.0; + delta_x2 = -float(coeff3 + coeff5) / float(2 * coeff1); + if (delta_x2 > 1.0) + delta_x2 = 1.0; + else if (delta_x2 < -1.0) + delta_x2 = -1.0; + } + else if (ty_) + { + delta_y2 = -1.0; + delta_x2 = -float(coeff3 - coeff5) / float(2 * coeff1); + if (delta_x2 > 1.0) + delta_x2 = 1.0; + else if (delta_x2 < -1.0) + delta_x2 = -1.0; + } + // insert both options for evaluation which to pick + float max1 = (coeff1 * delta_x1 * delta_x1 + coeff2 * delta_y1 * delta_y1 + coeff3 * delta_x1 + coeff4 * delta_y1 + + coeff5 * delta_x1 * delta_y1 + coeff6) + / 18.0; + float max2 = (coeff1 * delta_x2 * delta_x2 + coeff2 * delta_y2 * delta_y2 + coeff3 * delta_x2 + coeff4 * delta_y2 + + coeff5 * delta_x2 * delta_y2 + coeff6) + / 18.0; + if (max1 > max2) + { + delta_x = delta_x1; + delta_y = delta_x1; + return max1; + } + else + { + delta_x = delta_x2; + delta_y = delta_x2; + return max2; + } + } + + // this is the case of the maximum inside the boundaries: + return (coeff1 * delta_x * delta_x + coeff2 * delta_y * delta_y + coeff3 * delta_x + coeff4 * delta_y + + coeff5 * delta_x * delta_y + coeff6) + / 18.0; +} + +// construct a layer +BriskLayer::BriskLayer(const cv::Mat& img_in, float scale_in, float offset_in) +{ + img_ = img_in; + scores_ = cv::Mat_::zeros(img_in.rows, img_in.cols); + // attention: this means that the passed image reference must point to persistent memory + scale_ = scale_in; + offset_ = offset_in; + // create an agast detector + fast_9_16_ = new FastFeatureDetector(1, true, FastFeatureDetector::TYPE_9_16); + makeOffsets(pixel_5_8_, img_.step, 8); + makeOffsets(pixel_9_16_, img_.step, 16); +} +// derive a layer +BriskLayer::BriskLayer(const BriskLayer& layer, int mode) +{ + if (mode == CommonParams::HALFSAMPLE) + { + img_.create(layer.img().rows / 2, layer.img().cols / 2, CV_8U); + halfsample(layer.img(), img_); + scale_ = layer.scale() * 2; + offset_ = 0.5 * scale_ - 0.5; + } + else + { + img_.create(2 * (layer.img().rows / 3), 2 * (layer.img().cols / 3), CV_8U); + twothirdsample(layer.img(), img_); + scale_ = layer.scale() * 1.5; + offset_ = 0.5 * scale_ - 0.5; + } + scores_ = cv::Mat::zeros(img_.rows, img_.cols, CV_8U); + fast_9_16_ = new FastFeatureDetector(1, false, FastFeatureDetector::TYPE_9_16); + makeOffsets(pixel_5_8_, img_.step, 8); + makeOffsets(pixel_9_16_, img_.step, 16); +} + +// Fast/Agast +// wraps the agast class +void +BriskLayer::getAgastPoints(uint8_t threshold, std::vector& keypoints) +{ + fast_9_16_->set("threshold", threshold); + fast_9_16_->detect(img_, keypoints); + + // also write scores + const int num = keypoints.size(); + + for (int i = 0; i < num; i++) + scores_(keypoints[i].pt.y, keypoints[i].pt.x) = keypoints[i].response; +} +inline uint8_t +BriskLayer::getAgastScore(int x, int y, uint8_t threshold) +{ + if (x < 3 || y < 3) + return 0; + if (x >= img_.cols - 3 || y >= img_.rows - 3) + return 0; + uint8_t& score = *(scores_.data + x + y * scores_.cols); + if (score > 2) + { + return score; + } + score = cornerScore<16>(img_.data + x + y * img_.cols, pixel_9_16_, threshold - 1); + if (score < threshold) + score = 0; + return score; +} + +inline uint8_t +BriskLayer::getAgastScore_5_8(int x, int y, uint8_t threshold) +{ + if (x < 2 || y < 2) + return 0; + if (x >= img_.cols - 2 || y >= img_.rows - 2) + return 0; + uint8_t score = cornerScore<8>(img_.data + x + y * img_.cols, pixel_5_8_, threshold - 1); + if (score < threshold) + score = 0; + return score; +} + +inline uint8_t +BriskLayer::getAgastScore(float xf, float yf, uint8_t threshold_in, float scale_in) +{ + if (scale_in <= 1.0f) + { + // just do an interpolation inside the layer + const int x = int(xf); + const float rx1 = xf - float(x); + const float rx = 1.0f - rx1; + const int y = int(yf); + const float ry1 = yf - float(y); + const float ry = 1.0f - ry1; + + return rx * ry * getAgastScore(x, y, threshold_in) + rx1 * ry * getAgastScore(x + 1, y, threshold_in) + + rx * ry1 * getAgastScore(x, y + 1, threshold_in) + rx1 * ry1 * getAgastScore(x + 1, y + 1, threshold_in); + } + else + { + // this means we overlap area smoothing + const float halfscale = scale_in / 2.0f; + // get the scores first: + for (int x = int(xf - halfscale); x <= int(xf + halfscale + 1.0f); x++) + { + for (int y = int(yf - halfscale); y <= int(yf + halfscale + 1.0f); y++) + { + getAgastScore(x, y, threshold_in); + } + } + // get the smoothed value + return value(scores_, xf, yf, scale_in); + } +} + +// access gray values (smoothed/interpolated) +__inline__ uint8_t +BriskLayer::value(const cv::Mat& mat, float xf, float yf, float scale_in) +{ + assert(!mat.empty()); + // get the position + const int x = floor(xf); + const int y = floor(yf); + const cv::Mat& image = mat; + const int& imagecols = image.cols; + + // get the sigma_half: + const float sigma_half = scale_in / 2; + const float area = 4.0 * sigma_half * sigma_half; + // calculate output: + int ret_val; + if (sigma_half < 0.5) + { + //interpolation multipliers: + const int r_x = (xf - x) * 1024; + const int r_y = (yf - y) * 1024; + const int r_x_1 = (1024 - r_x); + const int r_y_1 = (1024 - r_y); + uchar* ptr = image.data + x + y * imagecols; + // just interpolate: + ret_val = (r_x_1 * r_y_1 * int(*ptr)); + ptr++; + ret_val += (r_x * r_y_1 * int(*ptr)); + ptr += imagecols; + ret_val += (r_x * r_y * int(*ptr)); + ptr--; + ret_val += (r_x_1 * r_y * int(*ptr)); + return 0xFF & ((ret_val + 512) / 1024 / 1024); + } + + // this is the standard case (simple, not speed optimized yet): + + // scaling: + const int scaling = 4194304.0 / area; + const int scaling2 = float(scaling) * area / 1024.0; + + // calculate borders + const float x_1 = xf - sigma_half; + const float x1 = xf + sigma_half; + const float y_1 = yf - sigma_half; + const float y1 = yf + sigma_half; + + const int x_left = int(x_1 + 0.5); + const int y_top = int(y_1 + 0.5); + const int x_right = int(x1 + 0.5); + const int y_bottom = int(y1 + 0.5); + + // overlap area - multiplication factors: + const float r_x_1 = float(x_left) - x_1 + 0.5; + const float r_y_1 = float(y_top) - y_1 + 0.5; + const float r_x1 = x1 - float(x_right) + 0.5; + const float r_y1 = y1 - float(y_bottom) + 0.5; + const int dx = x_right - x_left - 1; + const int dy = y_bottom - y_top - 1; + const int A = (r_x_1 * r_y_1) * scaling; + const int B = (r_x1 * r_y_1) * scaling; + const int C = (r_x1 * r_y1) * scaling; + const int D = (r_x_1 * r_y1) * scaling; + const int r_x_1_i = r_x_1 * scaling; + const int r_y_1_i = r_y_1 * scaling; + const int r_x1_i = r_x1 * scaling; + const int r_y1_i = r_y1 * scaling; + + // now the calculation: + uchar* ptr = image.data + x_left + imagecols * y_top; + // first row: + ret_val = A * int(*ptr); + ptr++; + const uchar* end1 = ptr + dx; + for (; ptr < end1; ptr++) + { + ret_val += r_y_1_i * int(*ptr); + } + ret_val += B * int(*ptr); + // middle ones: + ptr += imagecols - dx - 1; + uchar* end_j = ptr + dy * imagecols; + for (; ptr < end_j; ptr += imagecols - dx - 1) + { + ret_val += r_x_1_i * int(*ptr); + ptr++; + const uchar* end2 = ptr + dx; + for (; ptr < end2; ptr++) + { + ret_val += int(*ptr) * scaling; + } + ret_val += r_x1_i * int(*ptr); + } + // last row: + ret_val += D * int(*ptr); + ptr++; + const uchar* end3 = ptr + dx; + for (; ptr < end3; ptr++) + { + ret_val += r_y1_i * int(*ptr); + } + ret_val += C * int(*ptr); + + return 0xFF & ((ret_val + scaling2 / 2) / scaling2 / 1024); +} + +// half sampling +inline void +BriskLayer::halfsample(const cv::Mat& srcimg, cv::Mat& dstimg) +{ + // make sure the destination image is of the right size: + assert(srcimg.cols/2==dstimg.cols); + assert(srcimg.rows/2==dstimg.rows); + + // handle non-SSE case + resize(srcimg, dstimg, dstimg.size(), 0, 0, INTER_AREA); +} + +inline void +BriskLayer::twothirdsample(const cv::Mat& srcimg, cv::Mat& dstimg) +{ + // make sure the destination image is of the right size: + assert((srcimg.cols/3)*2==dstimg.cols); + assert((srcimg.rows/3)*2==dstimg.rows); + + resize(srcimg, dstimg, dstimg.size(), 0, 0, INTER_AREA); +} + +} diff --git a/modules/features2d/src/features2d_init.cpp b/modules/features2d/src/features2d_init.cpp index c9abfefa9..fa11b809c 100644 --- a/modules/features2d/src/features2d_init.cpp +++ b/modules/features2d/src/features2d_init.cpp @@ -51,6 +51,12 @@ using namespace cv; Otherwise, linker may throw away some seemingly unused stuff. */ +CV_INIT_ALGORITHM(BRISK, "Feature2D.BRISK", + obj.info()->addParam(obj, "thres", obj.threshold); + obj.info()->addParam(obj, "octaves", obj.octaves)); + +/////////////////////////////////////////////////////////////////////////////////////////////////////////// + CV_INIT_ALGORITHM(BriefDescriptorExtractor, "Feature2D.BRIEF", obj.info()->addParam(obj, "bytes", obj.bytes_)); @@ -154,6 +160,7 @@ bool cv::initModule_features2d(void) { bool all = true; all &= !BriefDescriptorExtractor_info_auto.name().empty(); + all &= !BRISK_info_auto.name().empty(); all &= !FastFeatureDetector_info_auto.name().empty(); all &= !StarDetector_info_auto.name().empty(); all &= !MSER_info_auto.name().empty(); diff --git a/modules/features2d/test/test_brisk.cpp b/modules/features2d/test/test_brisk.cpp new file mode 100644 index 000000000..10e71802a --- /dev/null +++ b/modules/features2d/test/test_brisk.cpp @@ -0,0 +1,95 @@ +/*M/////////////////////////////////////////////////////////////////////////////////////// +// +// IMPORTANT: READ BEFORE DOWNLOADING, COPYING, INSTALLING OR USING. +// +// By downloading, copying, installing or using the software you agree to this license. +// If you do not agree to this license, do not download, install, +// copy or use the software. +// +// +// License Agreement +// For Open Source Computer Vision Library +// +// Copyright (C) 2000-2008, Intel Corporation, all rights reserved. +// Copyright (C) 2009, Willow Garage Inc., all rights reserved. +// Third party copyrights are property of their respective owners. +// +// Redistribution and use in source and binary forms, with or without modification, +// are permitted provided that the following conditions are met: +// +// * Redistribution's of source code must retain the above copyright notice, +// this list of conditions and the following disclaimer. +// +// * Redistribution's in binary form must reproduce the above copyright notice, +// this list of conditions and the following disclaimer in the documentation +// and/or other materials provided with the distribution. +// +// * The name of the copyright holders may not be used to endorse or promote products +// derived from this software without specific prior written permission. +// +// This software is provided by the copyright holders and contributors "as is" and +// any express or implied warranties, including, but not limited to, the implied +// warranties of merchantability and fitness for a particular purpose are disclaimed. +// In no event shall the Intel Corporation or contributors be liable for any direct, +// indirect, incidental, special, exemplary, or consequential damages +// (including, but not limited to, procurement of substitute goods or services; +// loss of use, data, or profits; or business interruption) however caused +// and on any theory of liability, whether in contract, strict liability, +// or tort (including negligence or otherwise) arising in any way out of +// the use of this software, even if advised of the possibility of such damage. +// +//M*/ + +#include "test_precomp.hpp" + +using namespace cv; + +class CV_BRISKTest : public cvtest::BaseTest +{ +public: + CV_BRISKTest(); + ~CV_BRISKTest(); +protected: + void run(int); +}; + +CV_BRISKTest::CV_BRISKTest() {} +CV_BRISKTest::~CV_BRISKTest() {} + +void CV_BRISKTest::run( int ) +{ + Mat image1 = imread(string(ts->get_data_path()) + "inpaint/orig.jpg"); + Mat image2 = imread(string(ts->get_data_path()) + "cameracalibration/chess9.jpg"); + + if (image1.empty() || image2.empty()) + { + ts->set_failed_test_info( cvtest::TS::FAIL_INVALID_TEST_DATA ); + return; + } + + Mat gray1, gray2; + cvtColor(image1, gray1, CV_BGR2GRAY); + cvtColor(image2, gray2, CV_BGR2GRAY); + + Ptr detector = Algorithm::create("Feature2D.BRISK"); + + vector keypoints1; + vector keypoints2; + detector->detect(image1, keypoints1); + detector->detect(image2, keypoints2); + + for(size_t i = 0; i < keypoints1.size(); ++i) + { + const KeyPoint& kp = keypoints1[i]; + ASSERT_NE(kp.angle, -1); + } + + for(size_t i = 0; i < keypoints2.size(); ++i) + { + const KeyPoint& kp = keypoints2[i]; + ASSERT_NE(kp.angle, -1); + } +} + +TEST(Features2d_BRISK, regression) { CV_BRISKTest test; test.safe_run(); } + diff --git a/modules/features2d/test/test_descriptors_regression.cpp b/modules/features2d/test/test_descriptors_regression.cpp index b53ef57f7..2185625ae 100644 --- a/modules/features2d/test/test_descriptors_regression.cpp +++ b/modules/features2d/test/test_descriptors_regression.cpp @@ -301,6 +301,13 @@ private: * Tests registrations * \****************************************************************************************/ +TEST( Features2d_DescriptorExtractor_BRISK, regression ) +{ + CV_DescriptorExtractorTest test( "descriptor-brisk", (CV_DescriptorExtractorTest::DistanceType)2.f, + DescriptorExtractor::create("BRISK") ); + test.safe_run(); +} + TEST( Features2d_DescriptorExtractor_ORB, regression ) { // TODO adjust the parameters below diff --git a/modules/features2d/test/test_detectors_regression.cpp b/modules/features2d/test/test_detectors_regression.cpp index 0077f84cc..e5e8712ce 100644 --- a/modules/features2d/test/test_detectors_regression.cpp +++ b/modules/features2d/test/test_detectors_regression.cpp @@ -247,6 +247,12 @@ void CV_FeatureDetectorTest::run( int /*start_from*/ ) * Tests registrations * \****************************************************************************************/ +TEST( Features2d_Detector_BRISK, regression ) +{ + CV_FeatureDetectorTest test( "detector-brisk", FeatureDetector::create("BRISK") ); + test.safe_run(); +} + TEST( Features2d_Detector_FAST, regression ) { CV_FeatureDetectorTest test( "detector-fast", FeatureDetector::create("FAST") ); From b325b2f9a42568c1db6fd586ca5c2075103f0036 Mon Sep 17 00:00:00 2001 From: Vincent Rabaud Date: Fri, 24 Aug 2012 02:31:37 +0200 Subject: [PATCH 3/8] add more test to BRISK --- modules/features2d/test/test_keypoints.cpp | 6 ++++++ 1 file changed, 6 insertions(+) diff --git a/modules/features2d/test/test_keypoints.cpp b/modules/features2d/test/test_keypoints.cpp index c3f8c6fb5..b5d01a737 100644 --- a/modules/features2d/test/test_keypoints.cpp +++ b/modules/features2d/test/test_keypoints.cpp @@ -119,6 +119,12 @@ protected: // Registration of tests +TEST(Features2d_Detector_Keypoints_BRISK, validation) +{ + CV_FeatureDetectorKeypointsTest test(Algorithm::create("Feature2D.BRISK")); + test.safe_run(); +} + TEST(Features2d_Detector_Keypoints_FAST, validation) { CV_FeatureDetectorKeypointsTest test(Algorithm::create("Feature2D.FAST")); From 92da6d381bd069751c452cff7cc1890eae4670e9 Mon Sep 17 00:00:00 2001 From: Vincent Rabaud Date: Fri, 24 Aug 2012 03:42:58 +0200 Subject: [PATCH 4/8] add rotation tests for the descriptors --- modules/features2d/src/brisk.cpp | 3 ++ .../test_rotation_and_scale_invariance.cpp | 36 +++++++++++++++++++ 2 files changed, 39 insertions(+) diff --git a/modules/features2d/src/brisk.cpp b/modules/features2d/src/brisk.cpp index 415ed7821..8d5fe0d2e 100755 --- a/modules/features2d/src/brisk.cpp +++ b/modules/features2d/src/brisk.cpp @@ -530,6 +530,9 @@ BRISK::operator()( InputArray _image, InputArray _mask, vector& keypoi OutputArray _descriptors, bool useProvidedKeypoints) const { Mat image = _image.getMat(), mask = _mask.getMat(); + if( image.type() != CV_8UC1 ) + cvtColor(image, image, CV_BGR2GRAY); + if (!useProvidedKeypoints) detectImpl(image, keypoints, mask); diff --git a/modules/features2d/test/test_rotation_and_scale_invariance.cpp b/modules/features2d/test/test_rotation_and_scale_invariance.cpp index 48dddd3ee..ac3094663 100644 --- a/modules/features2d/test/test_rotation_and_scale_invariance.cpp +++ b/modules/features2d/test/test_rotation_and_scale_invariance.cpp @@ -592,6 +592,15 @@ protected: /* * Detector's rotation invariance check */ + +TEST(Features2d_RotationInvariance_Detector_BRISK, regression) +{ + DetectorRotationInvarianceTest test(Algorithm::create("Feature2D.BRISK"), + 0.32f, + 0.81f); + test.safe_run(); +} + TEST(Features2d_RotationInvariance_Detector_ORB, regression) { DetectorRotationInvarianceTest test(Algorithm::create("Feature2D.ORB"), @@ -603,6 +612,16 @@ TEST(Features2d_RotationInvariance_Detector_ORB, regression) /* * Descriptors's rotation invariance check */ + +TEST(Features2d_RotationInvariance_Descriptor_BRISK, regression) +{ + DescriptorRotationInvarianceTest test(Algorithm::create("Feature2D.BRISK"), + Algorithm::create("Feature2D.BRISK"), + NORM_HAMMING, + 0.99f); + test.safe_run(); +} + TEST(Features2d_RotationInvariance_Descriptor_ORB, regression) { DescriptorRotationInvarianceTest test(Algorithm::create("Feature2D.ORB"), @@ -625,6 +644,14 @@ TEST(Features2d_RotationInvariance_Descriptor_ORB, regression) * Detector's scale invariance check */ +//TEST(Features2d_ScaleInvariance_Detector_BRISK, regression) +//{ +// DetectorScaleInvarianceTest test(Algorithm::create("Feature2D.BRISK"), +// 0.09f, +// 0.52f); +// test.safe_run(); +//} + //TEST(Features2d_ScaleInvariance_Detector_ORB, regression) //{ // DetectorScaleInvarianceTest test(Algorithm::create("Feature2D.ORB"), @@ -637,6 +664,15 @@ TEST(Features2d_RotationInvariance_Descriptor_ORB, regression) * Descriptor's scale invariance check */ +//TEST(Features2d_ScaleInvariance_Descriptor_BRISK, regression) +//{ +// DescriptorScaleInvarianceTest test(Algorithm::create("Feature2D.BRISK"), +// Algorithm::create("Feature2D.BRISK"), +// NORM_HAMMING, +// 0.99f); +// test.safe_run(); +//} + //TEST(Features2d_ScaleInvariance_Descriptor_ORB, regression) //{ // DescriptorScaleInvarianceTest test(Algorithm::create("Feature2D.ORB"), From da1921b2fc727062049974a61537976e44c7037a Mon Sep 17 00:00:00 2001 From: Vincent Rabaud Date: Fri, 24 Aug 2012 04:03:06 +0200 Subject: [PATCH 5/8] add const correctness, replace __inline__ and remote a useless release() --- modules/features2d/src/brisk.cpp | 113 +++++++++++++++---------------- 1 file changed, 56 insertions(+), 57 deletions(-) diff --git a/modules/features2d/src/brisk.cpp b/modules/features2d/src/brisk.cpp index 8d5fe0d2e..537b8f797 100755 --- a/modules/features2d/src/brisk.cpp +++ b/modules/features2d/src/brisk.cpp @@ -74,11 +74,11 @@ public: // get scores - attention, this is in layer coordinates, not scale=1 coordinates! inline uint8_t - getAgastScore(int x, int y, uint8_t threshold); + getAgastScore(int x, int y, uint8_t threshold) const; inline uint8_t - getAgastScore_5_8(int x, int y, uint8_t threshold); + getAgastScore_5_8(int x, int y, uint8_t threshold) const; inline uint8_t - getAgastScore(float xf, float yf, uint8_t threshold, float scale = 1.0f); + getAgastScore(float xf, float yf, uint8_t threshold, float scale = 1.0f) const; // accessors inline const cv::Mat& @@ -111,8 +111,8 @@ public: private: // access gray values (smoothed/interpolated) - __inline__ uint8_t - value(const cv::Mat& mat, float xf, float yf, float scale); + inline uint8_t + value(const cv::Mat& mat, float xf, float yf, float scale) const; // the image cv::Mat img_; // its Fast scores @@ -139,40 +139,40 @@ public: // get Keypoints void - getKeypoints(const uint8_t _threshold, std::vector& keypoints); + getKeypoints(const uint8_t _threshold, std::vector& keypoints); protected: // nonmax suppression: - __inline__ bool + inline bool isMax2D(const uint8_t layer, const int x_layer, const int y_layer); // 1D (scale axis) refinement: - __inline__ float - refine1D(const float s_05, const float s0, const float s05, float& max); // around octave - __inline__ float - refine1D_1(const float s_05, const float s0, const float s05, float& max); // around intra - __inline__ float - refine1D_2(const float s_05, const float s0, const float s05, float& max); // around octave 0 only + inline float + refine1D(const float s_05, const float s0, const float s05, float& max) const; // around octave + inline float + refine1D_1(const float s_05, const float s0, const float s05, float& max) const; // around intra + inline float + refine1D_2(const float s_05, const float s0, const float s05, float& max) const; // around octave 0 only // 2D maximum refinement: - __inline__ float + inline float subpixel2D(const int s_0_0, const int s_0_1, const int s_0_2, const int s_1_0, const int s_1_1, const int s_1_2, - const int s_2_0, const int s_2_1, const int s_2_2, float& delta_x, float& delta_y); + const int s_2_0, const int s_2_1, const int s_2_2, float& delta_x, float& delta_y) const; // 3D maximum refinement centered around (x_layer,y_layer) - __inline__ float - refine3D(const uint8_t layer, const int x_layer, const int y_layer, float& x, float& y, float& scale, bool& ismax); + inline float + refine3D(const uint8_t layer, const int x_layer, const int y_layer, float& x, float& y, float& scale, bool& ismax) const; // interpolated score access with recalculation when needed: - __inline__ int - getScoreAbove(const uint8_t layer, const int x_layer, const int y_layer); - __inline__ int - getScoreBelow(const uint8_t layer, const int x_layer, const int y_layer); + inline int + getScoreAbove(const uint8_t layer, const int x_layer, const int y_layer) const; + inline int + getScoreBelow(const uint8_t layer, const int x_layer, const int y_layer) const; // return the maximum of score patches above or below - __inline__ float + inline float getScoreMaxAbove(const uint8_t layer, const int x_layer, const int y_layer, const int threshold, bool& ismax, - float& dx, float& dy); - __inline__ float + float& dx, float& dy) const; + inline float getScoreMaxBelow(const uint8_t layer, const int x_layer, const int y_layer, const int threshold, bool& ismax, - float& dx, float& dy); + float& dx, float& dy) const; // the image pyramids: uint8_t layers_; @@ -354,7 +354,7 @@ BRISK::generateKernel(std::vector &radiusList, std::vector &numberLi } // simple alternative: -__inline__ int +inline int BRISK::smoothedIntensity(const cv::Mat& image, const cv::Mat& integral, const float key_x, const float key_y, const unsigned int scale, const unsigned int rot, const unsigned int point) const @@ -650,7 +650,6 @@ BRISK::operator()( InputArray _image, InputArray _mask, vector& keypoi } // clean-up - _integral.release(); delete[] _values; } @@ -946,11 +945,11 @@ BriskScaleSpace::getKeypoints(const uint8_t _threshold, std::vector0); - BriskLayer& layerBelow = pyramid_[layer - 1]; + const BriskLayer& layerBelow = pyramid_[layer - 1]; // check the first row int max_x = x_1 + 1; @@ -1754,8 +1753,8 @@ BriskScaleSpace::getScoreMaxBelow(const uint8_t layer, const int x_layer, const return max; } -__inline__ float -BriskScaleSpace::refine1D(const float s_05, const float s0, const float s05, float& max) +inline float +BriskScaleSpace::refine1D(const float s_05, const float s0, const float s05, float& max) const { int i_05 = int(1024.0 * s_05 + 0.5); int i0 = int(1024.0 * s0 + 0.5); @@ -1800,8 +1799,8 @@ BriskScaleSpace::refine1D(const float s_05, const float s0, const float s05, flo return ret_val; } -__inline__ float -BriskScaleSpace::refine1D_1(const float s_05, const float s0, const float s05, float& max) +inline float +BriskScaleSpace::refine1D_1(const float s_05, const float s0, const float s05, float& max) const { int i_05 = int(1024.0 * s_05 + 0.5); int i0 = int(1024.0 * s0 + 0.5); @@ -1846,8 +1845,8 @@ BriskScaleSpace::refine1D_1(const float s_05, const float s0, const float s05, f return ret_val; } -__inline__ float -BriskScaleSpace::refine1D_2(const float s_05, const float s0, const float s05, float& max) +inline float +BriskScaleSpace::refine1D_2(const float s_05, const float s0, const float s05, float& max) const { int i_05 = int(1024.0 * s_05 + 0.5); int i0 = int(1024.0 * s0 + 0.5); @@ -1892,10 +1891,10 @@ BriskScaleSpace::refine1D_2(const float s_05, const float s0, const float s05, f return ret_val; } -__inline__ float +inline float BriskScaleSpace::subpixel2D(const int s_0_0, const int s_0_1, const int s_0_2, const int s_1_0, const int s_1_1, const int s_1_2, const int s_2_0, const int s_2_1, const int s_2_2, float& delta_x, - float& delta_y) + float& delta_y) const { // the coefficients of the 2d quadratic function least-squares fit: @@ -2086,7 +2085,7 @@ BriskLayer::getAgastPoints(uint8_t threshold, std::vector& keypoints) scores_(keypoints[i].pt.y, keypoints[i].pt.x) = keypoints[i].response; } inline uint8_t -BriskLayer::getAgastScore(int x, int y, uint8_t threshold) +BriskLayer::getAgastScore(int x, int y, uint8_t threshold) const { if (x < 3 || y < 3) return 0; @@ -2104,7 +2103,7 @@ BriskLayer::getAgastScore(int x, int y, uint8_t threshold) } inline uint8_t -BriskLayer::getAgastScore_5_8(int x, int y, uint8_t threshold) +BriskLayer::getAgastScore_5_8(int x, int y, uint8_t threshold) const { if (x < 2 || y < 2) return 0; @@ -2117,7 +2116,7 @@ BriskLayer::getAgastScore_5_8(int x, int y, uint8_t threshold) } inline uint8_t -BriskLayer::getAgastScore(float xf, float yf, uint8_t threshold_in, float scale_in) +BriskLayer::getAgastScore(float xf, float yf, uint8_t threshold_in, float scale_in) const { if (scale_in <= 1.0f) { @@ -2150,8 +2149,8 @@ BriskLayer::getAgastScore(float xf, float yf, uint8_t threshold_in, float scale_ } // access gray values (smoothed/interpolated) -__inline__ uint8_t -BriskLayer::value(const cv::Mat& mat, float xf, float yf, float scale_in) +inline uint8_t +BriskLayer::value(const cv::Mat& mat, float xf, float yf, float scale_in) const { assert(!mat.empty()); // get the position From 84c4797030e981f46257b34a1260892cceabb0a1 Mon Sep 17 00:00:00 2001 From: Vincent Rabaud Date: Fri, 24 Aug 2012 16:23:07 +0200 Subject: [PATCH 6/8] revert orientation computation in jeypoint detection for efficiency (like done originally) --- .../include/opencv2/features2d/features2d.hpp | 5 + modules/features2d/src/brisk.cpp | 158 +++++++----------- 2 files changed, 70 insertions(+), 93 deletions(-) diff --git a/modules/features2d/include/opencv2/features2d/features2d.hpp b/modules/features2d/include/opencv2/features2d/features2d.hpp index a081fd4ee..774b01c18 100644 --- a/modules/features2d/include/opencv2/features2d/features2d.hpp +++ b/modules/features2d/include/opencv2/features2d/features2d.hpp @@ -307,6 +307,11 @@ protected: void computeImpl( const Mat& image, vector& keypoints, Mat& descriptors ) const; void detectImpl( const Mat& image, vector& keypoints, const Mat& mask=Mat() ) const; + void computeKeypointsNoOrientation(InputArray image, InputArray mask, vector& keypoints) const; + void computeDescriptorsAndOrOrientation(InputArray image, InputArray mask, vector& keypoints, + OutputArray descriptors, bool doDescriptors, bool doOrientation, + bool useProvidedKeypoints) const; + // Feature parameters CV_PROP_RW int threshold; CV_PROP_RW int octaves; diff --git a/modules/features2d/src/brisk.cpp b/modules/features2d/src/brisk.cpp index 537b8f797..72fc9a8ab 100755 --- a/modules/features2d/src/brisk.cpp +++ b/modules/features2d/src/brisk.cpp @@ -528,13 +528,28 @@ RoiPredicate(const float minX, const float minY, const float maxX, const float m void BRISK::operator()( InputArray _image, InputArray _mask, vector& keypoints, OutputArray _descriptors, bool useProvidedKeypoints) const +{ + bool doOrientation; + if (useProvidedKeypoints) + doOrientation = false; + computeDescriptorsAndOrOrientation(_image, _mask, keypoints, _descriptors, true, doOrientation, + useProvidedKeypoints); +} + +void +BRISK::computeDescriptorsAndOrOrientation(InputArray _image, InputArray _mask, vector& keypoints, + OutputArray _descriptors, bool doDescriptors, bool doOrientation, + bool useProvidedKeypoints) const { Mat image = _image.getMat(), mask = _mask.getMat(); if( image.type() != CV_8UC1 ) cvtColor(image, image, CV_BGR2GRAY); if (!useProvidedKeypoints) - detectImpl(image, keypoints, mask); + { + doOrientation = true; + computeKeypointsNoOrientation(_image, _mask, keypoints); + } //Remove keypoints very close to the border size_t ksize = keypoints.size(); @@ -578,9 +593,13 @@ BRISK::operator()( InputArray _image, InputArray _mask, vector& keypoi int* _values = new int[points_]; // for temporary use // resize the descriptors: - _descriptors.create(ksize, strings_, CV_8U); - cv::Mat descriptors = _descriptors.getMat(); - descriptors.setTo(0); + cv::Mat descriptors; + if (doDescriptors) + { + _descriptors.create(ksize, strings_, CV_8U); + descriptors = _descriptors.getMat(); + descriptors.setTo(0); + } // now do the extraction for all keypoints: @@ -592,13 +611,44 @@ BRISK::operator()( InputArray _image, InputArray _mask, vector& keypoi uchar* ptr = descriptors.data; for (size_t k = 0; k < ksize; k++) { - int theta; cv::KeyPoint& kp = keypoints[k]; const int& scale = kscales[k]; - int shifter = 0; int* pvalues = _values; const float& x = kp.pt.x; const float& y = kp.pt.y; + + if (doOrientation) + { + // get the gray values in the unrotated pattern + for (unsigned int i = 0; i < points_; i++) + { + *(pvalues++) = smoothedIntensity(image, _integral, x, y, scale, 0, i); + } + + int direction0 = 0; + int direction1 = 0; + // now iterate through the long pairings + const BriskLongPair* max = longPairs_ + noLongPairs_; + for (BriskLongPair* iter = longPairs_; iter < max; ++iter) + { + t1 = *(_values + iter->i); + t2 = *(_values + iter->j); + const int delta_t = (t1 - t2); + // update the direction: + const int tmp0 = delta_t * (iter->weighted_dx) / 1024; + const int tmp1 = delta_t * (iter->weighted_dy) / 1024; + direction0 += tmp0; + direction1 += tmp1; + } + kp.angle = atan2((float) direction1, (float) direction0) / M_PI * 180.0; + if (kp.angle < 0) + kp.angle += 360; + } + + if (!doDescriptors) + continue; + + int theta; if (kp.angle==-1) { // don't compute the gradient direction, just assign a rotation of 0° @@ -615,7 +665,7 @@ BRISK::operator()( InputArray _image, InputArray _mask, vector& keypoi // now also extract the stuff for the actual direction: // let us compute the smoothed values - shifter = 0; + int shifter = 0; //unsigned int mean=0; pvalues = _values; @@ -675,7 +725,14 @@ BRISK::~BRISK() } void -BRISK::operator()(InputArray _image, InputArray _mask, vector& keypoints) const +BRISK::operator()(InputArray image, InputArray mask, vector& keypoints) const +{ + computeKeypointsNoOrientation(image, mask, keypoints); + computeDescriptorsAndOrOrientation(image, mask, keypoints, cv::noArray(), false, true, true); +} + +void +BRISK::computeKeypointsNoOrientation(InputArray _image, InputArray _mask, vector& keypoints) const { Mat image = _image.getMat(), mask = _mask.getMat(); if( image.type() != CV_8UC1 ) @@ -687,91 +744,6 @@ BRISK::operator()(InputArray _image, InputArray _mask, vector& keypoin // remove invalid points removeInvalidPoints(mask, keypoints); - - // Compute the orientations of the keypoints - //Remove keypoints very close to the border - size_t ksize = keypoints.size(); - std::vector kscales; // remember the scale per keypoint - kscales.resize(ksize); - static const float log2 = 0.693147180559945; - static const float lb_scalerange = log(scalerange_) / (log2); - std::vector::iterator beginning = keypoints.begin(); - std::vector::iterator beginningkscales = kscales.begin(); - static const float basicSize06 = basicSize_ * 0.6; - for (size_t k = 0; k < ksize; k++) - { - unsigned int scale; - scale = std::max((int) (scales_ / lb_scalerange * (log(keypoints[k].size / (basicSize06)) / log2) + 0.5), 0); - // saturate - if (scale >= scales_) - scale = scales_ - 1; - kscales[k] = scale; - const int border = sizeList_[scale]; - const int border_x = image.cols - border; - const int border_y = image.rows - border; - if (RoiPredicate(border, border, border_x, border_y, keypoints[k])) - { - keypoints.erase(beginning + k); - kscales.erase(beginningkscales + k); - if (k == 0) - { - beginning = keypoints.begin(); - beginningkscales = kscales.begin(); - } - ksize--; - k--; - } - } - - // first, calculate the integral image over the whole image: - // current integral image - cv::Mat _integral; // the integral image - cv::integral(image, _integral); - - int* _values = new int[points_]; // for temporary use - - // now do the extraction for all keypoints: - - // temporary variables containing gray values at sample points: - int t1; - int t2; - - // the feature orientation - int direction0; - int direction1; - - for (size_t k = 0; k < ksize; k++) - { - cv::KeyPoint& kp = keypoints[k]; - const int& scale = kscales[k]; - int* pvalues = _values; - const float& x = kp.pt.x; - const float& y = kp.pt.y; - // get the gray values in the unrotated pattern - for (unsigned int i = 0; i < points_; i++) - { - *(pvalues++) = smoothedIntensity(image, _integral, x, y, scale, 0, i); - } - - direction0 = 0; - direction1 = 0; - // now iterate through the long pairings - const BriskLongPair* max = longPairs_ + noLongPairs_; - for (BriskLongPair* iter = longPairs_; iter < max; ++iter) - { - t1 = *(_values + iter->i); - t2 = *(_values + iter->j); - const int delta_t = (t1 - t2); - // update the direction: - const int tmp0 = delta_t * (iter->weighted_dx) / 1024; - const int tmp1 = delta_t * (iter->weighted_dy) / 1024; - direction0 += tmp0; - direction1 += tmp1; - } - kp.angle = atan2((float) direction1, (float) direction0) / M_PI * 180.0; - if (kp.angle < 0) - kp.angle += 360; - } } From 3ca0cc22531473c2aa0adb3acf8a06d9334a2b07 Mon Sep 17 00:00:00 2001 From: Vincent Rabaud Date: Fri, 24 Aug 2012 16:35:35 +0200 Subject: [PATCH 7/8] add a scale invariance test for the detector --- .../test/test_rotation_and_scale_invariance.cpp | 14 +++++++------- 1 file changed, 7 insertions(+), 7 deletions(-) diff --git a/modules/features2d/test/test_rotation_and_scale_invariance.cpp b/modules/features2d/test/test_rotation_and_scale_invariance.cpp index ac3094663..1441c87fd 100644 --- a/modules/features2d/test/test_rotation_and_scale_invariance.cpp +++ b/modules/features2d/test/test_rotation_and_scale_invariance.cpp @@ -644,13 +644,13 @@ TEST(Features2d_RotationInvariance_Descriptor_ORB, regression) * Detector's scale invariance check */ -//TEST(Features2d_ScaleInvariance_Detector_BRISK, regression) -//{ -// DetectorScaleInvarianceTest test(Algorithm::create("Feature2D.BRISK"), -// 0.09f, -// 0.52f); -// test.safe_run(); -//} +TEST(Features2d_ScaleInvariance_Detector_BRISK, regression) +{ + DetectorScaleInvarianceTest test(Algorithm::create("Feature2D.BRISK"), + 0.08f, + 0.54f); + test.safe_run(); +} //TEST(Features2d_ScaleInvariance_Detector_ORB, regression) //{ From 014106783d854cb4d89015171b26e79dda3c2632 Mon Sep 17 00:00:00 2001 From: Vadim Pisarevsky Date: Thu, 30 Aug 2012 17:32:47 +0400 Subject: [PATCH 8/8] fixed building BRISK on Windows --- modules/features2d/src/brisk.cpp | 1096 +++++++++++++++--------------- modules/imgproc/src/imgwarp.cpp | 2 +- 2 files changed, 542 insertions(+), 556 deletions(-) diff --git a/modules/features2d/src/brisk.cpp b/modules/features2d/src/brisk.cpp index 72fc9a8ab..2a54d7c1c 100755 --- a/modules/features2d/src/brisk.cpp +++ b/modules/features2d/src/brisk.cpp @@ -70,15 +70,15 @@ public: // Fast/Agast without non-max suppression void - getAgastPoints(uint8_t threshold, std::vector& keypoints); + getAgastPoints(int threshold, std::vector& keypoints); // get scores - attention, this is in layer coordinates, not scale=1 coordinates! - inline uint8_t - getAgastScore(int x, int y, uint8_t threshold) const; - inline uint8_t - getAgastScore_5_8(int x, int y, uint8_t threshold) const; - inline uint8_t - getAgastScore(float xf, float yf, uint8_t threshold, float scale = 1.0f) const; + inline int + getAgastScore(int x, int y, int threshold) const; + inline int + getAgastScore_5_8(int x, int y, int threshold) const; + inline int + getAgastScore(float xf, float yf, int threshold, float scale = 1.0f) const; // accessors inline const cv::Mat& @@ -111,7 +111,7 @@ public: private: // access gray values (smoothed/interpolated) - inline uint8_t + inline int value(const cv::Mat& mat, float xf, float yf, float scale) const; // the image cv::Mat img_; @@ -130,7 +130,7 @@ class CV_EXPORTS BriskScaleSpace { public: // construct telling the octaves number: - BriskScaleSpace(uint8_t _octaves = 3); + BriskScaleSpace(int _octaves = 3); ~BriskScaleSpace(); // construct the image pyramids @@ -139,12 +139,12 @@ public: // get Keypoints void - getKeypoints(const uint8_t _threshold, std::vector& keypoints); + getKeypoints(const int _threshold, std::vector& keypoints); protected: // nonmax suppression: inline bool - isMax2D(const uint8_t layer, const int x_layer, const int y_layer); + isMax2D(const int layer, const int x_layer, const int y_layer); // 1D (scale axis) refinement: inline float refine1D(const float s_05, const float s0, const float s05, float& max) const; // around octave @@ -158,24 +158,24 @@ protected: const int s_2_0, const int s_2_1, const int s_2_2, float& delta_x, float& delta_y) const; // 3D maximum refinement centered around (x_layer,y_layer) inline float - refine3D(const uint8_t layer, const int x_layer, const int y_layer, float& x, float& y, float& scale, bool& ismax) const; + refine3D(const int layer, const int x_layer, const int y_layer, float& x, float& y, float& scale, bool& ismax) const; // interpolated score access with recalculation when needed: inline int - getScoreAbove(const uint8_t layer, const int x_layer, const int y_layer) const; + getScoreAbove(const int layer, const int x_layer, const int y_layer) const; inline int - getScoreBelow(const uint8_t layer, const int x_layer, const int y_layer) const; + getScoreBelow(const int layer, const int x_layer, const int y_layer) const; // return the maximum of score patches above or below inline float - getScoreMaxAbove(const uint8_t layer, const int x_layer, const int y_layer, const int threshold, bool& ismax, + getScoreMaxAbove(const int layer, const int x_layer, const int y_layer, const int threshold, bool& ismax, float& dx, float& dy) const; inline float - getScoreMaxBelow(const uint8_t layer, const int x_layer, const int y_layer, const int threshold, bool& ismax, + getScoreMaxBelow(const int layer, const int x_layer, const int y_layer, const int threshold, bool& ismax, float& dx, float& dy) const; // the image pyramids: - uint8_t layers_; + int layers_; std::vector pyramid_; // some constant parameters: @@ -183,15 +183,13 @@ protected: static const float basicSize_; }; -using namespace cv; - -const float BRISK::basicSize_ = 12.0; +const float BRISK::basicSize_ = 12.0f; const unsigned int BRISK::scales_ = 64; -const float BRISK::scalerange_ = 30; // 40->4 Octaves - else, this needs to be adjusted... +const float BRISK::scalerange_ = 30.f; // 40->4 Octaves - else, this needs to be adjusted... const unsigned int BRISK::n_rot_ = 1024; // discretization of the rotation look-up -const float BriskScaleSpace::safetyFactor_ = 1.0; -const float BriskScaleSpace::basicSize_ = 12.0; +const float BriskScaleSpace::safetyFactor_ = 1.0f; +const float BriskScaleSpace::basicSize_ = 12.0f; // constructors BRISK::BRISK(int thresh, int octaves_in, float patternScale) @@ -207,11 +205,11 @@ BRISK::BRISK(int thresh, int octaves_in, float patternScale) nList.resize(5); const double f = 0.85 * patternScale; - rList[0] = f * 0; - rList[1] = f * 2.9; - rList[2] = f * 4.9; - rList[3] = f * 7.4; - rList[4] = f * 10.8; + rList[0] = (float)(f * 0.); + rList[1] = (float)(f * 2.9); + rList[2] = (float)(f * 4.9); + rList[3] = (float)(f * 7.4); + rList[4] = (float)(f * 10.8); nList[0] = 1; nList[1] = 10; @@ -219,7 +217,7 @@ BRISK::BRISK(int thresh, int octaves_in, float patternScale) nList[3] = 15; nList[4] = 20; - generateKernel(rList, nList, 5.85 * patternScale, 8.2 * patternScale); + generateKernel(rList, nList, (float)(5.85 * patternScale), (float)(8.2 * patternScale)); } BRISK::BRISK(std::vector &radiusList, std::vector &numberList, float dMax, float dMin, @@ -237,7 +235,7 @@ BRISK::generateKernel(std::vector &radiusList, std::vector &numberLi dMin_ = dMin; // get the total number of points - const int rings = radiusList.size(); + const int rings = (int)radiusList.size(); assert(radiusList.size()!=0&&radiusList.size()==numberList.size()); points_ = 0; // remember the total number of points for (int ring = 0; ring < rings; ring++) @@ -249,44 +247,44 @@ BRISK::generateKernel(std::vector &radiusList, std::vector &numberLi BriskPatternPoint* patternIterator = patternPoints_; // define the scale discretization: - static const float lb_scale = log(scalerange_) / log(2.0); + static const float lb_scale = (float)(log(scalerange_) / log(2.0)); static const float lb_scale_step = lb_scale / (scales_); scaleList_ = new float[scales_]; sizeList_ = new unsigned int[scales_]; - const float sigma_scale = 1.3; + const float sigma_scale = 1.3f; for (unsigned int scale = 0; scale < scales_; ++scale) { - scaleList_[scale] = pow((double) 2.0, (double) (scale * lb_scale_step)); + scaleList_[scale] = (float)pow((double) 2.0, (double) (scale * lb_scale_step)); sizeList_[scale] = 0; // generate the pattern points look-up double alpha, theta; for (size_t rot = 0; rot < n_rot_; ++rot) { - theta = double(rot) * 2 * M_PI / double(n_rot_); // this is the rotation of the feature + theta = double(rot) * 2 * CV_PI / double(n_rot_); // this is the rotation of the feature for (int ring = 0; ring < rings; ++ring) { for (int num = 0; num < numberList[ring]; ++num) { // the actual coordinates on the circle - alpha = (double(num)) * 2 * M_PI / double(numberList[ring]); - patternIterator->x = scaleList_[scale] * radiusList[ring] * cos(alpha + theta); // feature rotation plus angle of the point - patternIterator->y = scaleList_[scale] * radiusList[ring] * sin(alpha + theta); + alpha = (double(num)) * 2 * CV_PI / double(numberList[ring]); + patternIterator->x = (float)(scaleList_[scale] * radiusList[ring] * cos(alpha + theta)); // feature rotation plus angle of the point + patternIterator->y = (float)(scaleList_[scale] * radiusList[ring] * sin(alpha + theta)); // and the gaussian kernel sigma if (ring == 0) { - patternIterator->sigma = sigma_scale * scaleList_[scale] * 0.5; + patternIterator->sigma = sigma_scale * scaleList_[scale] * 0.5f; } else { - patternIterator->sigma = sigma_scale * scaleList_[scale] * (double(radiusList[ring])) - * sin(M_PI / numberList[ring]); + patternIterator->sigma = (float)(sigma_scale * scaleList_[scale] * (double(radiusList[ring])) + * sin(CV_PI / numberList[ring])); } // adapt the sizeList if necessary - const unsigned int size = ceil(((scaleList_[scale] * radiusList[ring]) + patternIterator->sigma)) + 1; + const unsigned int size = cvCeil(((scaleList_[scale] * radiusList[ring]) + patternIterator->sigma)) + 1; if (sizeList_[scale] < size) { sizeList_[scale] = size; @@ -306,11 +304,11 @@ BRISK::generateKernel(std::vector &radiusList, std::vector &numberLi noLongPairs_ = 0; // fill indexChange with 0..n if empty - unsigned int indSize = indexChange.size(); + unsigned int indSize = (unsigned int)indexChange.size(); if (indSize == 0) { indexChange.resize(points_ * (points_ - 1) / 2); - indSize = indexChange.size(); + indSize = (unsigned int)indexChange.size(); } for (unsigned int i = 0; i < indSize; i++) { @@ -370,34 +368,30 @@ BRISK::smoothedIntensity(const cv::Mat& image, const cv::Mat& integral, const fl // get the sigma: const float sigma_half = briskPoint.sigma; - const float area = 4.0 * sigma_half * sigma_half; + const float area = 4.0f * sigma_half * sigma_half; // calculate output: int ret_val; if (sigma_half < 0.5) { //interpolation multipliers: - const int r_x = (xf - x) * 1024; - const int r_y = (yf - y) * 1024; + const int r_x = (int)((xf - x) * 1024); + const int r_y = (int)((yf - y) * 1024); const int r_x_1 = (1024 - r_x); const int r_y_1 = (1024 - r_y); - uchar* ptr = image.data + x + y * imagecols; + const uchar* ptr = &image.at(y, x); + size_t step = image.step; // just interpolate: - ret_val = (r_x_1 * r_y_1 * int(*ptr)); - ptr++; - ret_val += (r_x * r_y_1 * int(*ptr)); - ptr += imagecols; - ret_val += (r_x * r_y * int(*ptr)); - ptr--; - ret_val += (r_x_1 * r_y * int(*ptr)); + ret_val = r_x_1 * r_y_1 * ptr[0] + r_x * r_y_1 * ptr[1] + + r_x * r_y * ptr[step] + r_x_1 * r_y * ptr[step+1]; return (ret_val + 512) / 1024; } // this is the standard case (simple, not speed optimized yet): // scaling: - const int scaling = 4194304.0 / area; - const int scaling2 = float(scaling) * area / 1024.0; + const int scaling = (int)(4194304.0 / area); + const int scaling2 = int(float(scaling) * area / 1024.0); // the integral image is larger: const int integralcols = imagecols + 1; @@ -414,20 +408,20 @@ BRISK::smoothedIntensity(const cv::Mat& image, const cv::Mat& integral, const fl const int y_bottom = int(y1 + 0.5); // overlap area - multiplication factors: - const float r_x_1 = float(x_left) - x_1 + 0.5; - const float r_y_1 = float(y_top) - y_1 + 0.5; - const float r_x1 = x1 - float(x_right) + 0.5; - const float r_y1 = y1 - float(y_bottom) + 0.5; + const float r_x_1 = float(x_left) - x_1 + 0.5f; + const float r_y_1 = float(y_top) - y_1 + 0.5f; + const float r_x1 = x1 - float(x_right) + 0.5f; + const float r_y1 = y1 - float(y_bottom) + 0.5f; const int dx = x_right - x_left - 1; const int dy = y_bottom - y_top - 1; - const int A = (r_x_1 * r_y_1) * scaling; - const int B = (r_x1 * r_y_1) * scaling; - const int C = (r_x1 * r_y1) * scaling; - const int D = (r_x_1 * r_y1) * scaling; - const int r_x_1_i = r_x_1 * scaling; - const int r_y_1_i = r_y_1 * scaling; - const int r_x1_i = r_x1 * scaling; - const int r_y1_i = r_y1 * scaling; + const int A = (int)((r_x_1 * r_y_1) * scaling); + const int B = (int)((r_x1 * r_y_1) * scaling); + const int C = (int)((r_x1 * r_y1) * scaling); + const int D = (int)((r_x_1 * r_y1) * scaling); + const int r_x_1_i = (int)(r_x_1 * scaling); + const int r_y_1_i = (int)(r_y_1 * scaling); + const int r_x1_i = (int)(r_x1 * scaling); + const int r_y1_i = (int)(r_y1 * scaling); if (dx + dy > 2) { @@ -529,7 +523,7 @@ void BRISK::operator()( InputArray _image, InputArray _mask, vector& keypoints, OutputArray _descriptors, bool useProvidedKeypoints) const { - bool doOrientation; + bool doOrientation=true; if (useProvidedKeypoints) doOrientation = false; computeDescriptorsAndOrOrientation(_image, _mask, keypoints, _descriptors, true, doOrientation, @@ -555,11 +549,11 @@ BRISK::computeDescriptorsAndOrOrientation(InputArray _image, InputArray _mask, v size_t ksize = keypoints.size(); std::vector kscales; // remember the scale per keypoint kscales.resize(ksize); - static const float log2 = 0.693147180559945; - static const float lb_scalerange = log(scalerange_) / (log2); + static const float log2 = 0.693147180559945f; + static const float lb_scalerange = (float)(log(scalerange_) / (log2)); std::vector::iterator beginning = keypoints.begin(); std::vector::iterator beginningkscales = kscales.begin(); - static const float basicSize06 = basicSize_ * 0.6; + static const float basicSize06 = basicSize_ * 0.6f; for (size_t k = 0; k < ksize; k++) { unsigned int scale; @@ -571,7 +565,7 @@ BRISK::computeDescriptorsAndOrOrientation(InputArray _image, InputArray _mask, v const int border = sizeList_[scale]; const int border_x = image.cols - border; const int border_y = image.rows - border; - if (RoiPredicate(border, border, border_x, border_y, keypoints[k])) + if (RoiPredicate((float)border, (float)border, (float)border_x, (float)border_y, keypoints[k])) { keypoints.erase(beginning + k); kscales.erase(beginningkscales + k); @@ -596,7 +590,7 @@ BRISK::computeDescriptorsAndOrOrientation(InputArray _image, InputArray _mask, v cv::Mat descriptors; if (doDescriptors) { - _descriptors.create(ksize, strings_, CV_8U); + _descriptors.create((int)ksize, strings_, CV_8U); descriptors = _descriptors.getMat(); descriptors.setTo(0); } @@ -640,9 +634,9 @@ BRISK::computeDescriptorsAndOrOrientation(InputArray _image, InputArray _mask, v direction0 += tmp0; direction1 += tmp1; } - kp.angle = atan2((float) direction1, (float) direction0) / M_PI * 180.0; + kp.angle = (float)(atan2((float) direction1, (float) direction0) / CV_PI * 180.0); if (kp.angle < 0) - kp.angle += 360; + kp.angle += 360.f; } if (!doDescriptors) @@ -760,7 +754,7 @@ BRISK::computeImpl( const Mat& image, vector& keypoints, Mat& descript } // construct telling the octaves number: -BriskScaleSpace::BriskScaleSpace(uint8_t _octaves) +BriskScaleSpace::BriskScaleSpace(int _octaves) { if (_octaves == 0) layers_ = 1; @@ -787,7 +781,7 @@ BriskScaleSpace::constructPyramid(const cv::Mat& image) } const int octaves2 = layers_; - for (uint8_t i = 2; i < octaves2; i += 2) + for (uchar i = 2; i < octaves2; i += 2) { pyramid_.push_back(BriskLayer(pyramid_[i - 2], BriskLayer::CommonParams::HALFSAMPLE)); pyramid_.push_back(BriskLayer(pyramid_[i - 1], BriskLayer::CommonParams::HALFSAMPLE)); @@ -795,20 +789,19 @@ BriskScaleSpace::constructPyramid(const cv::Mat& image) } void -BriskScaleSpace::getKeypoints(const uint8_t _threshold, std::vector& keypoints) +BriskScaleSpace::getKeypoints(const int threshold_, std::vector& keypoints) { // make sure keypoints is empty keypoints.resize(0); keypoints.reserve(2000); // assign thresholds - uint8_t threshold_ = _threshold; - uint8_t safeThreshold_ = threshold_ * safetyFactor_; + int safeThreshold_ = (int)(threshold_ * safetyFactor_); std::vector > agastPoints; agastPoints.resize(layers_); // go through the octaves and intra layers and calculate fast corner scores: - for (uint8_t i = 0; i < layers_; i++) + for (int i = 0; i < layers_; i++) { // call OAST16_9 without nms BriskLayer& l = pyramid_[i]; @@ -818,25 +811,25 @@ BriskScaleSpace::getKeypoints(const uint8_t _threshold, std::vector max_below_uchar) - max_below_uchar = s_1_0; - register int s_2_0 = l.getAgastScore_5_8(x_layer + 1, y_layer - 1, 1); - if (s_2_0 > max_below_uchar) - max_below_uchar = s_2_0; - register int s_2_1 = l.getAgastScore_5_8(x_layer + 1, y_layer, 1); - if (s_2_1 > max_below_uchar) - max_below_uchar = s_2_1; - register int s_1_1 = l.getAgastScore_5_8(x_layer, y_layer, 1); - if (s_1_1 > max_below_uchar) - max_below_uchar = s_1_1; - register int s_0_1 = l.getAgastScore_5_8(x_layer - 1, y_layer, 1); - if (s_0_1 > max_below_uchar) - max_below_uchar = s_0_1; - register int s_0_2 = l.getAgastScore_5_8(x_layer - 1, y_layer + 1, 1); - if (s_0_2 > max_below_uchar) - max_below_uchar = s_0_2; - register int s_1_2 = l.getAgastScore_5_8(x_layer, y_layer + 1, 1); - if (s_1_2 > max_below_uchar) - max_below_uchar = s_1_2; - register int s_2_2 = l.getAgastScore_5_8(x_layer + 1, y_layer + 1, 1); - if (s_2_2 > max_below_uchar) - max_below_uchar = s_2_2; + int s_0_0 = l.getAgastScore_5_8(x_layer - 1, y_layer - 1, 1); + max_below = s_0_0; + int s_1_0 = l.getAgastScore_5_8(x_layer, y_layer - 1, 1); + max_below = std::max(s_1_0, max_below); + int s_2_0 = l.getAgastScore_5_8(x_layer + 1, y_layer - 1, 1); + max_below = std::max(s_2_0, max_below); + int s_2_1 = l.getAgastScore_5_8(x_layer + 1, y_layer, 1); + max_below = std::max(s_2_1, max_below); + int s_1_1 = l.getAgastScore_5_8(x_layer, y_layer, 1); + max_below = std::max(s_1_1, max_below); + int s_0_1 = l.getAgastScore_5_8(x_layer - 1, y_layer, 1); + max_below = std::max(s_0_1, max_below); + int s_0_2 = l.getAgastScore_5_8(x_layer - 1, y_layer + 1, 1); + max_below = std::max(s_0_2, max_below); + int s_1_2 = l.getAgastScore_5_8(x_layer, y_layer + 1, 1); + max_below = std::max(s_1_2, max_below); + int s_2_2 = l.getAgastScore_5_8(x_layer + 1, y_layer + 1, 1); + max_below = std::max(s_2_2, max_below); max_below_float = subpixel2D(s_0_0, s_0_1, s_0_2, s_1_0, s_1_1, s_1_2, s_2_0, s_2_1, s_2_2, delta_x_below, delta_y_below); - max_below_float = max_below_uchar; + max_below_float = (float)max_below; } else { @@ -1246,15 +1231,15 @@ BriskScaleSpace::refine3D(const uint8_t layer, const int x_layer, const int y_la } // get the patch on this layer: - register int s_0_0 = thisLayer.getAgastScore(x_layer - 1, y_layer - 1, 1); - register int s_1_0 = thisLayer.getAgastScore(x_layer, y_layer - 1, 1); - register int s_2_0 = thisLayer.getAgastScore(x_layer + 1, y_layer - 1, 1); - register int s_2_1 = thisLayer.getAgastScore(x_layer + 1, y_layer, 1); - register int s_1_1 = thisLayer.getAgastScore(x_layer, y_layer, 1); - register int s_0_1 = thisLayer.getAgastScore(x_layer - 1, y_layer, 1); - register int s_0_2 = thisLayer.getAgastScore(x_layer - 1, y_layer + 1, 1); - register int s_1_2 = thisLayer.getAgastScore(x_layer, y_layer + 1, 1); - register int s_2_2 = thisLayer.getAgastScore(x_layer + 1, y_layer + 1, 1); + int s_0_0 = thisLayer.getAgastScore(x_layer - 1, y_layer - 1, 1); + int s_1_0 = thisLayer.getAgastScore(x_layer, y_layer - 1, 1); + int s_2_0 = thisLayer.getAgastScore(x_layer + 1, y_layer - 1, 1); + int s_2_1 = thisLayer.getAgastScore(x_layer + 1, y_layer, 1); + int s_1_1 = thisLayer.getAgastScore(x_layer, y_layer, 1); + int s_0_1 = thisLayer.getAgastScore(x_layer - 1, y_layer, 1); + int s_0_2 = thisLayer.getAgastScore(x_layer - 1, y_layer + 1, 1); + int s_1_2 = thisLayer.getAgastScore(x_layer, y_layer + 1, 1); + int s_2_2 = thisLayer.getAgastScore(x_layer + 1, y_layer + 1, 1); float delta_x_layer, delta_y_layer; float max_layer = subpixel2D(s_0_0, s_0_1, s_0_2, s_1_0, s_1_1, s_1_2, s_2_0, s_2_1, s_2_2, delta_x_layer, delta_y_layer); @@ -1270,8 +1255,8 @@ BriskScaleSpace::refine3D(const uint8_t layer, const int x_layer, const int y_la if (scale > 1.0) { // interpolate the position: - const float r0 = (1.5 - scale) / .5; - const float r1 = 1.0 - r0; + const float r0 = (1.5f - scale) / .5f; + const float r1 = 1.0f - r0; x = (r0 * delta_x_layer + r1 * delta_x_above + float(x_layer)) * thisLayer.scale() + thisLayer.offset(); y = (r0 * delta_y_layer + r1 * delta_y_above + float(y_layer)) * thisLayer.scale() + thisLayer.offset(); } @@ -1280,16 +1265,16 @@ BriskScaleSpace::refine3D(const uint8_t layer, const int x_layer, const int y_la if (layer == 0) { // interpolate the position: - const float r0 = (scale - 0.5) / 0.5; - const float r_1 = 1.0 - r0; + const float r0 = (scale - 0.5f) / 0.5f; + const float r_1 = 1.0f - r0; x = r0 * delta_x_layer + r_1 * delta_x_below + float(x_layer); y = r0 * delta_y_layer + r_1 * delta_y_below + float(y_layer); } else { // interpolate the position: - const float r0 = (scale - 0.75) / 0.25; - const float r_1 = 1.0 - r0; + const float r0 = (scale - 0.75f) / 0.25f; + const float r_1 = 1.0f - r0; x = (r0 * delta_x_layer + r_1 * delta_x_below + float(x_layer)) * thisLayer.scale() + thisLayer.offset(); y = (r0 * delta_y_layer + r_1 * delta_y_below + float(y_layer)) * thisLayer.scale() + thisLayer.offset(); } @@ -1302,18 +1287,18 @@ BriskScaleSpace::refine3D(const uint8_t layer, const int x_layer, const int y_la float delta_x_below, delta_y_below; float max_below = getScoreMaxBelow(layer, x_layer, y_layer, center, ismax, delta_x_below, delta_y_below); if (!ismax) - return 0.0; + return 0.0f; // get the patch on this layer: - register int s_0_0 = thisLayer.getAgastScore(x_layer - 1, y_layer - 1, 1); - register int s_1_0 = thisLayer.getAgastScore(x_layer, y_layer - 1, 1); - register int s_2_0 = thisLayer.getAgastScore(x_layer + 1, y_layer - 1, 1); - register int s_2_1 = thisLayer.getAgastScore(x_layer + 1, y_layer, 1); - register int s_1_1 = thisLayer.getAgastScore(x_layer, y_layer, 1); - register int s_0_1 = thisLayer.getAgastScore(x_layer - 1, y_layer, 1); - register int s_0_2 = thisLayer.getAgastScore(x_layer - 1, y_layer + 1, 1); - register int s_1_2 = thisLayer.getAgastScore(x_layer, y_layer + 1, 1); - register int s_2_2 = thisLayer.getAgastScore(x_layer + 1, y_layer + 1, 1); + int s_0_0 = thisLayer.getAgastScore(x_layer - 1, y_layer - 1, 1); + int s_1_0 = thisLayer.getAgastScore(x_layer, y_layer - 1, 1); + int s_2_0 = thisLayer.getAgastScore(x_layer + 1, y_layer - 1, 1); + int s_2_1 = thisLayer.getAgastScore(x_layer + 1, y_layer, 1); + int s_1_1 = thisLayer.getAgastScore(x_layer, y_layer, 1); + int s_0_1 = thisLayer.getAgastScore(x_layer - 1, y_layer, 1); + int s_0_2 = thisLayer.getAgastScore(x_layer - 1, y_layer + 1, 1); + int s_1_2 = thisLayer.getAgastScore(x_layer, y_layer + 1, 1); + int s_2_2 = thisLayer.getAgastScore(x_layer + 1, y_layer + 1, 1); float delta_x_layer, delta_y_layer; float max_layer = subpixel2D(s_0_0, s_0_1, s_0_2, s_1_0, s_1_1, s_1_2, s_2_0, s_2_1, s_2_2, delta_x_layer, delta_y_layer); @@ -1323,16 +1308,16 @@ BriskScaleSpace::refine3D(const uint8_t layer, const int x_layer, const int y_la if (scale > 1.0) { // interpolate the position: - const float r0 = 4.0 - scale * 3.0; - const float r1 = 1.0 - r0; + const float r0 = 4.0f - scale * 3.0f; + const float r1 = 1.0f - r0; x = (r0 * delta_x_layer + r1 * delta_x_above + float(x_layer)) * thisLayer.scale() + thisLayer.offset(); y = (r0 * delta_y_layer + r1 * delta_y_above + float(y_layer)) * thisLayer.scale() + thisLayer.offset(); } else { // interpolate the position: - const float r0 = scale * 3.0 - 2.0; - const float r_1 = 1.0 - r0; + const float r0 = scale * 3.0f - 2.0f; + const float r_1 = 1.0f - r0; x = (r0 * delta_x_layer + r_1 * delta_x_below + float(x_layer)) * thisLayer.scale() + thisLayer.offset(); y = (r0 * delta_y_layer + r_1 * delta_y_below + float(y_layer)) * thisLayer.scale() + thisLayer.offset(); } @@ -1347,7 +1332,7 @@ BriskScaleSpace::refine3D(const uint8_t layer, const int x_layer, const int y_la // return the maximum of score patches above or below inline float -BriskScaleSpace::getScoreMaxAbove(const uint8_t layer, const int x_layer, const int y_layer, const int threshold, +BriskScaleSpace::getScoreMaxAbove(const int layer, const int x_layer, const int y_layer, const int threshold, bool& ismax, float& dx, float& dy) const { @@ -1365,10 +1350,10 @@ BriskScaleSpace::getScoreMaxAbove(const uint8_t layer, const int x_layer, const if (layer % 2 == 0) { // octave - x_1 = float(4 * (x_layer) - 1 - 2) / 6.0; - x1 = float(4 * (x_layer) - 1 + 2) / 6.0; - y_1 = float(4 * (y_layer) - 1 - 2) / 6.0; - y1 = float(4 * (y_layer) - 1 + 2) / 6.0; + x_1 = float(4 * (x_layer) - 1 - 2) / 6.0f; + x1 = float(4 * (x_layer) - 1 + 2) / 6.0f; + y_1 = float(4 * (y_layer) - 1 - 2) / 6.0f; + y1 = float(4 * (y_layer) - 1 + 2) / 6.0f; } else { @@ -1380,103 +1365,103 @@ BriskScaleSpace::getScoreMaxAbove(const uint8_t layer, const int x_layer, const } // check the first row - int max_x = x_1 + 1; - int max_y = y_1 + 1; + int max_x = (int)x_1 + 1; + int max_y = (int)y_1 + 1; float tmp_max; - float max = layerAbove.getAgastScore(x_1, y_1, 1); - if (max > threshold) + float maxval = (float)layerAbove.getAgastScore(x_1, y_1, 1); + if (maxval > threshold) return 0; - for (int x = x_1 + 1; x <= int(x1); x++) + for (int x = (int)x_1 + 1; x <= int(x1); x++) { - tmp_max = layerAbove.getAgastScore(float(x), y_1, 1); + tmp_max = (float)layerAbove.getAgastScore(float(x), y_1, 1); if (tmp_max > threshold) return 0; - if (tmp_max > max) + if (tmp_max > maxval) { - max = tmp_max; + maxval = tmp_max; max_x = x; } } - tmp_max = layerAbove.getAgastScore(x1, y_1, 1); + tmp_max = (float)layerAbove.getAgastScore(x1, y_1, 1); if (tmp_max > threshold) return 0; - if (tmp_max > max) + if (tmp_max > maxval) { - max = tmp_max; + maxval = tmp_max; max_x = int(x1); } // middle rows - for (int y = y_1 + 1; y <= int(y1); y++) + for (int y = (int)y_1 + 1; y <= int(y1); y++) { - tmp_max = layerAbove.getAgastScore(x_1, float(y), 1); + tmp_max = (float)layerAbove.getAgastScore(x_1, float(y), 1); if (tmp_max > threshold) return 0; - if (tmp_max > max) + if (tmp_max > maxval) { - max = tmp_max; + maxval = tmp_max; max_x = int(x_1 + 1); max_y = y; } - for (int x = x_1 + 1; x <= int(x1); x++) + for (int x = (int)x_1 + 1; x <= int(x1); x++) { - tmp_max = layerAbove.getAgastScore(x, y, 1); + tmp_max = (float)layerAbove.getAgastScore(x, y, 1); if (tmp_max > threshold) return 0; - if (tmp_max > max) + if (tmp_max > maxval) { - max = tmp_max; + maxval = tmp_max; max_x = x; max_y = y; } } - tmp_max = layerAbove.getAgastScore(x1, float(y), 1); + tmp_max = (float)layerAbove.getAgastScore(x1, float(y), 1); if (tmp_max > threshold) return 0; - if (tmp_max > max) + if (tmp_max > maxval) { - max = tmp_max; + maxval = tmp_max; max_x = int(x1); max_y = y; } } // bottom row - tmp_max = layerAbove.getAgastScore(x_1, y1, 1); - if (tmp_max > max) + tmp_max = (float)layerAbove.getAgastScore(x_1, y1, 1); + if (tmp_max > maxval) { - max = tmp_max; + maxval = tmp_max; max_x = int(x_1 + 1); max_y = int(y1); } - for (int x = x_1 + 1; x <= int(x1); x++) + for (int x = (int)x_1 + 1; x <= int(x1); x++) { - tmp_max = layerAbove.getAgastScore(float(x), y1, 1); - if (tmp_max > max) + tmp_max = (float)layerAbove.getAgastScore(float(x), y1, 1); + if (tmp_max > maxval) { - max = tmp_max; + maxval = tmp_max; max_x = x; max_y = int(y1); } } - tmp_max = layerAbove.getAgastScore(x1, y1, 1); - if (tmp_max > max) + tmp_max = (float)layerAbove.getAgastScore(x1, y1, 1); + if (tmp_max > maxval) { - max = tmp_max; + maxval = tmp_max; max_x = int(x1); max_y = int(y1); } //find dx/dy: - register int s_0_0 = layerAbove.getAgastScore(max_x - 1, max_y - 1, 1); - register int s_1_0 = layerAbove.getAgastScore(max_x, max_y - 1, 1); - register int s_2_0 = layerAbove.getAgastScore(max_x + 1, max_y - 1, 1); - register int s_2_1 = layerAbove.getAgastScore(max_x + 1, max_y, 1); - register int s_1_1 = layerAbove.getAgastScore(max_x, max_y, 1); - register int s_0_1 = layerAbove.getAgastScore(max_x - 1, max_y, 1); - register int s_0_2 = layerAbove.getAgastScore(max_x - 1, max_y + 1, 1); - register int s_1_2 = layerAbove.getAgastScore(max_x, max_y + 1, 1); - register int s_2_2 = layerAbove.getAgastScore(max_x + 1, max_y + 1, 1); + int s_0_0 = layerAbove.getAgastScore(max_x - 1, max_y - 1, 1); + int s_1_0 = layerAbove.getAgastScore(max_x, max_y - 1, 1); + int s_2_0 = layerAbove.getAgastScore(max_x + 1, max_y - 1, 1); + int s_2_1 = layerAbove.getAgastScore(max_x + 1, max_y, 1); + int s_1_1 = layerAbove.getAgastScore(max_x, max_y, 1); + int s_0_1 = layerAbove.getAgastScore(max_x - 1, max_y, 1); + int s_0_2 = layerAbove.getAgastScore(max_x - 1, max_y + 1, 1); + int s_1_2 = layerAbove.getAgastScore(max_x, max_y + 1, 1); + int s_2_2 = layerAbove.getAgastScore(max_x + 1, max_y + 1, 1); float dx_1, dy_1; float refined_max = subpixel2D(s_0_0, s_0_1, s_0_2, s_1_0, s_1_1, s_1_2, s_2_0, s_2_1, s_2_2, dx_1, dy_1); @@ -1491,8 +1476,8 @@ BriskScaleSpace::getScoreMaxAbove(const uint8_t layer, const int x_layer, const } else { - dx = (real_x * 8.0 + 1.0) / 6.0 - float(x_layer); - dy = (real_y * 8.0 + 1.0) / 6.0 - float(y_layer); + dx = (real_x * 8.0f + 1.0f) / 6.0f - float(x_layer); + dy = (real_y * 8.0f + 1.0f) / 6.0f - float(y_layer); } // saturate @@ -1517,6 +1502,205 @@ BriskScaleSpace::getScoreMaxAbove(const uint8_t layer, const int x_layer, const returnrefined = false; } + // done and ok. + ismax = true; + if (returnrefined) + { + return std::max(refined_max, maxval); + } + return maxval; +} + +inline float +BriskScaleSpace::getScoreMaxBelow(const int layer, const int x_layer, const int y_layer, const int threshold, + bool& ismax, float& dx, float& dy) const +{ + ismax = false; + + // relevant floating point coordinates + float x_1; + float x1; + float y_1; + float y1; + + if (layer % 2 == 0) + { + // octave + x_1 = float(8 * (x_layer) + 1 - 4) / 6.0f; + x1 = float(8 * (x_layer) + 1 + 4) / 6.0f; + y_1 = float(8 * (y_layer) + 1 - 4) / 6.0f; + y1 = float(8 * (y_layer) + 1 + 4) / 6.0f; + } + else + { + x_1 = float(6 * (x_layer) + 1 - 3) / 4.0f; + x1 = float(6 * (x_layer) + 1 + 3) / 4.0f; + y_1 = float(6 * (y_layer) + 1 - 3) / 4.0f; + y1 = float(6 * (y_layer) + 1 + 3) / 4.0f; + } + + // the layer below + assert(layer>0); + const BriskLayer& layerBelow = pyramid_[layer - 1]; + + // check the first row + int max_x = (int)x_1 + 1; + int max_y = (int)y_1 + 1; + float tmp_max; + float max = (float)layerBelow.getAgastScore(x_1, y_1, 1); + if (max > threshold) + return 0; + for (int x = (int)x_1 + 1; x <= int(x1); x++) + { + tmp_max = (float)layerBelow.getAgastScore(float(x), y_1, 1); + if (tmp_max > threshold) + return 0; + if (tmp_max > max) + { + max = tmp_max; + max_x = x; + } + } + tmp_max = (float)layerBelow.getAgastScore(x1, y_1, 1); + if (tmp_max > threshold) + return 0; + if (tmp_max > max) + { + max = tmp_max; + max_x = int(x1); + } + + // middle rows + for (int y = (int)y_1 + 1; y <= int(y1); y++) + { + tmp_max = (float)layerBelow.getAgastScore(x_1, float(y), 1); + if (tmp_max > threshold) + return 0; + if (tmp_max > max) + { + max = tmp_max; + max_x = int(x_1 + 1); + max_y = y; + } + for (int x = (int)x_1 + 1; x <= int(x1); x++) + { + tmp_max = (float)layerBelow.getAgastScore(x, y, 1); + if (tmp_max > threshold) + return 0; + if (tmp_max == max) + { + const int t1 = 2 + * (layerBelow.getAgastScore(x - 1, y, 1) + layerBelow.getAgastScore(x + 1, y, 1) + + layerBelow.getAgastScore(x, y + 1, 1) + layerBelow.getAgastScore(x, y - 1, 1)) + + (layerBelow.getAgastScore(x + 1, y + 1, 1) + layerBelow.getAgastScore(x - 1, y + 1, 1) + + layerBelow.getAgastScore(x + 1, y - 1, 1) + layerBelow.getAgastScore(x - 1, y - 1, 1)); + const int t2 = 2 + * (layerBelow.getAgastScore(max_x - 1, max_y, 1) + layerBelow.getAgastScore(max_x + 1, max_y, 1) + + layerBelow.getAgastScore(max_x, max_y + 1, 1) + layerBelow.getAgastScore(max_x, max_y - 1, 1)) + + (layerBelow.getAgastScore(max_x + 1, max_y + 1, 1) + layerBelow.getAgastScore(max_x - 1, + max_y + 1, 1) + + layerBelow.getAgastScore(max_x + 1, max_y - 1, 1) + + layerBelow.getAgastScore(max_x - 1, max_y - 1, 1)); + if (t1 > t2) + { + max_x = x; + max_y = y; + } + } + if (tmp_max > max) + { + max = tmp_max; + max_x = x; + max_y = y; + } + } + tmp_max = (float)layerBelow.getAgastScore(x1, float(y), 1); + if (tmp_max > threshold) + return 0; + if (tmp_max > max) + { + max = tmp_max; + max_x = int(x1); + max_y = y; + } + } + + // bottom row + tmp_max = (float)layerBelow.getAgastScore(x_1, y1, 1); + if (tmp_max > max) + { + max = tmp_max; + max_x = int(x_1 + 1); + max_y = int(y1); + } + for (int x = (int)x_1 + 1; x <= int(x1); x++) + { + tmp_max = (float)layerBelow.getAgastScore(float(x), y1, 1); + if (tmp_max > max) + { + max = tmp_max; + max_x = x; + max_y = int(y1); + } + } + tmp_max = (float)layerBelow.getAgastScore(x1, y1, 1); + if (tmp_max > max) + { + max = tmp_max; + max_x = int(x1); + max_y = int(y1); + } + + //find dx/dy: + int s_0_0 = layerBelow.getAgastScore(max_x - 1, max_y - 1, 1); + int s_1_0 = layerBelow.getAgastScore(max_x, max_y - 1, 1); + int s_2_0 = layerBelow.getAgastScore(max_x + 1, max_y - 1, 1); + int s_2_1 = layerBelow.getAgastScore(max_x + 1, max_y, 1); + int s_1_1 = layerBelow.getAgastScore(max_x, max_y, 1); + int s_0_1 = layerBelow.getAgastScore(max_x - 1, max_y, 1); + int s_0_2 = layerBelow.getAgastScore(max_x - 1, max_y + 1, 1); + int s_1_2 = layerBelow.getAgastScore(max_x, max_y + 1, 1); + int s_2_2 = layerBelow.getAgastScore(max_x + 1, max_y + 1, 1); + float dx_1, dy_1; + float refined_max = subpixel2D(s_0_0, s_0_1, s_0_2, s_1_0, s_1_1, s_1_2, s_2_0, s_2_1, s_2_2, dx_1, dy_1); + + // calculate dx/dy in above coordinates + float real_x = float(max_x) + dx_1; + float real_y = float(max_y) + dy_1; + bool returnrefined = true; + if (layer % 2 == 0) + { + dx = (float)((real_x * 6.0 + 1.0) / 8.0) - float(x_layer); + dy = (float)((real_y * 6.0 + 1.0) / 8.0) - float(y_layer); + } + else + { + dx = (float)((real_x * 4.0 - 1.0) / 6.0) - float(x_layer); + dy = (float)((real_y * 4.0 - 1.0) / 6.0) - float(y_layer); + } + + // saturate + if (dx > 1.0) + { + dx = 1.0f; + returnrefined = false; + } + if (dx < -1.0f) + { + dx = -1.0f; + returnrefined = false; + } + if (dy > 1.0f) + { + dy = 1.0f; + returnrefined = false; + } + if (dy < -1.0f) + { + dy = -1.0f; + returnrefined = false; + } + // done and ok. ismax = true; if (returnrefined) @@ -1526,205 +1710,6 @@ BriskScaleSpace::getScoreMaxAbove(const uint8_t layer, const int x_layer, const return max; } -inline float -BriskScaleSpace::getScoreMaxBelow(const uint8_t layer, const int x_layer, const int y_layer, const int threshold, - bool& ismax, float& dx, float& dy) const -{ - ismax = false; - - // relevant floating point coordinates - float x_1; - float x1; - float y_1; - float y1; - - if (layer % 2 == 0) - { - // octave - x_1 = float(8 * (x_layer) + 1 - 4) / 6.0; - x1 = float(8 * (x_layer) + 1 + 4) / 6.0; - y_1 = float(8 * (y_layer) + 1 - 4) / 6.0; - y1 = float(8 * (y_layer) + 1 + 4) / 6.0; - } - else - { - x_1 = float(6 * (x_layer) + 1 - 3) / 4.0; - x1 = float(6 * (x_layer) + 1 + 3) / 4.0; - y_1 = float(6 * (y_layer) + 1 - 3) / 4.0; - y1 = float(6 * (y_layer) + 1 + 3) / 4.0; - } - - // the layer below - assert(layer>0); - const BriskLayer& layerBelow = pyramid_[layer - 1]; - - // check the first row - int max_x = x_1 + 1; - int max_y = y_1 + 1; - float tmp_max; - float max = layerBelow.getAgastScore(x_1, y_1, 1); - if (max > threshold) - return 0; - for (int x = x_1 + 1; x <= int(x1); x++) - { - tmp_max = layerBelow.getAgastScore(float(x), y_1, 1); - if (tmp_max > threshold) - return 0; - if (tmp_max > max) - { - max = tmp_max; - max_x = x; - } - } - tmp_max = layerBelow.getAgastScore(x1, y_1, 1); - if (tmp_max > threshold) - return 0; - if (tmp_max > max) - { - max = tmp_max; - max_x = int(x1); - } - - // middle rows - for (int y = y_1 + 1; y <= int(y1); y++) - { - tmp_max = layerBelow.getAgastScore(x_1, float(y), 1); - if (tmp_max > threshold) - return 0; - if (tmp_max > max) - { - max = tmp_max; - max_x = int(x_1 + 1); - max_y = y; - } - for (int x = x_1 + 1; x <= int(x1); x++) - { - tmp_max = layerBelow.getAgastScore(x, y, 1); - if (tmp_max > threshold) - return 0; - if (tmp_max == max) - { - const int t1 = 2 - * (layerBelow.getAgastScore(x - 1, y, 1) + layerBelow.getAgastScore(x + 1, y, 1) - + layerBelow.getAgastScore(x, y + 1, 1) + layerBelow.getAgastScore(x, y - 1, 1)) - + (layerBelow.getAgastScore(x + 1, y + 1, 1) + layerBelow.getAgastScore(x - 1, y + 1, 1) - + layerBelow.getAgastScore(x + 1, y - 1, 1) + layerBelow.getAgastScore(x - 1, y - 1, 1)); - const int t2 = 2 - * (layerBelow.getAgastScore(max_x - 1, max_y, 1) + layerBelow.getAgastScore(max_x + 1, max_y, 1) - + layerBelow.getAgastScore(max_x, max_y + 1, 1) + layerBelow.getAgastScore(max_x, max_y - 1, 1)) - + (layerBelow.getAgastScore(max_x + 1, max_y + 1, 1) + layerBelow.getAgastScore(max_x - 1, - max_y + 1, 1) - + layerBelow.getAgastScore(max_x + 1, max_y - 1, 1) - + layerBelow.getAgastScore(max_x - 1, max_y - 1, 1)); - if (t1 > t2) - { - max_x = x; - max_y = y; - } - } - if (tmp_max > max) - { - max = tmp_max; - max_x = x; - max_y = y; - } - } - tmp_max = layerBelow.getAgastScore(x1, float(y), 1); - if (tmp_max > threshold) - return 0; - if (tmp_max > max) - { - max = tmp_max; - max_x = int(x1); - max_y = y; - } - } - - // bottom row - tmp_max = layerBelow.getAgastScore(x_1, y1, 1); - if (tmp_max > max) - { - max = tmp_max; - max_x = int(x_1 + 1); - max_y = int(y1); - } - for (int x = x_1 + 1; x <= int(x1); x++) - { - tmp_max = layerBelow.getAgastScore(float(x), y1, 1); - if (tmp_max > max) - { - max = tmp_max; - max_x = x; - max_y = int(y1); - } - } - tmp_max = layerBelow.getAgastScore(x1, y1, 1); - if (tmp_max > max) - { - max = tmp_max; - max_x = int(x1); - max_y = int(y1); - } - - //find dx/dy: - register int s_0_0 = layerBelow.getAgastScore(max_x - 1, max_y - 1, 1); - register int s_1_0 = layerBelow.getAgastScore(max_x, max_y - 1, 1); - register int s_2_0 = layerBelow.getAgastScore(max_x + 1, max_y - 1, 1); - register int s_2_1 = layerBelow.getAgastScore(max_x + 1, max_y, 1); - register int s_1_1 = layerBelow.getAgastScore(max_x, max_y, 1); - register int s_0_1 = layerBelow.getAgastScore(max_x - 1, max_y, 1); - register int s_0_2 = layerBelow.getAgastScore(max_x - 1, max_y + 1, 1); - register int s_1_2 = layerBelow.getAgastScore(max_x, max_y + 1, 1); - register int s_2_2 = layerBelow.getAgastScore(max_x + 1, max_y + 1, 1); - float dx_1, dy_1; - float refined_max = subpixel2D(s_0_0, s_0_1, s_0_2, s_1_0, s_1_1, s_1_2, s_2_0, s_2_1, s_2_2, dx_1, dy_1); - - // calculate dx/dy in above coordinates - float real_x = float(max_x) + dx_1; - float real_y = float(max_y) + dy_1; - bool returnrefined = true; - if (layer % 2 == 0) - { - dx = (real_x * 6.0 + 1.0) / 8.0 - float(x_layer); - dy = (real_y * 6.0 + 1.0) / 8.0 - float(y_layer); - } - else - { - dx = (real_x * 4.0 - 1.0) / 6.0 - float(x_layer); - dy = (real_y * 4.0 - 1.0) / 6.0 - float(y_layer); - } - - // saturate - if (dx > 1.0) - { - dx = 1.0; - returnrefined = false; - } - if (dx < -1.0) - { - dx = -1.0; - returnrefined = false; - } - if (dy > 1.0) - { - dy = 1.0; - returnrefined = false; - } - if (dy < -1.0) - { - dy = -1.0; - returnrefined = false; - } - - // done and ok. - ismax = true; - if (returnrefined) - { - return std::max(refined_max, max); - } - return max; -} - inline float BriskScaleSpace::refine1D(const float s_05, const float s0, const float s05, float& max) const { @@ -1743,17 +1728,17 @@ BriskScaleSpace::refine1D(const float s_05, const float s0, const float s05, flo if (s0 >= s_05 && s0 >= s05) { max = s0; - return 1.0; + return 1.0f; } if (s_05 >= s0 && s_05 >= s05) { max = s_05; - return 0.75; + return 0.75f; } if (s05 >= s0 && s05 >= s_05) { max = s05; - return 1.5; + return 1.5f; } } @@ -1767,7 +1752,7 @@ BriskScaleSpace::refine1D(const float s_05, const float s0, const float s05, flo ret_val = 1.5; // allow to be slightly off bounds ...? int three_c = +24 * i_05 - 27 * i0 + 6 * i05; max = float(three_c) + float(three_a) * ret_val * ret_val + float(three_b) * ret_val; - max /= 3072.0; + max /= 3072.0f; return ret_val; } @@ -1789,17 +1774,17 @@ BriskScaleSpace::refine1D_1(const float s_05, const float s0, const float s05, f if (s0 >= s_05 && s0 >= s05) { max = s0; - return 1.0; + return 1.0f; } if (s_05 >= s0 && s_05 >= s05) { max = s_05; - return 0.6666666666666666666666666667; + return 0.6666666666666666666666666667f; } if (s05 >= s0 && s05 >= s_05) { max = s05; - return 1.3333333333333333333333333333; + return 1.3333333333333333333333333333f; } } @@ -1807,13 +1792,13 @@ BriskScaleSpace::refine1D_1(const float s_05, const float s0, const float s05, f // calculate max location: float ret_val = -float(two_b) / float(2 * two_a); // saturate and return - if (ret_val < 0.6666666666666666666666666667) - ret_val = 0.666666666666666666666666667; - else if (ret_val > 1.33333333333333333333333333) - ret_val = 1.333333333333333333333333333; + if (ret_val < 0.6666666666666666666666666667f) + ret_val = 0.666666666666666666666666667f; + else if (ret_val > 1.33333333333333333333333333f) + ret_val = 1.333333333333333333333333333f; int two_c = +12 * i_05 - 16 * i0 + 6 * i05; max = float(two_c) + float(two_a) * ret_val * ret_val + float(two_b) * ret_val; - max /= 2048.0; + max /= 2048.0f; return ret_val; } @@ -1835,17 +1820,17 @@ BriskScaleSpace::refine1D_2(const float s_05, const float s0, const float s05, f if (s0 >= s_05 && s0 >= s05) { max = s0; - return 1.0; + return 1.0f; } if (s_05 >= s0 && s_05 >= s05) { max = s_05; - return 0.7; + return 0.7f; } if (s05 >= s0 && s05 >= s_05) { max = s05; - return 1.5; + return 1.5f; } } @@ -1853,10 +1838,10 @@ BriskScaleSpace::refine1D_2(const float s_05, const float s0, const float s05, f // calculate max location: float ret_val = -float(b) / float(2 * a); // saturate and return - if (ret_val < 0.7) - ret_val = 0.7; - else if (ret_val > 1.5) - ret_val = 1.5; // allow to be slightly off bounds ...? + if (ret_val < 0.7f) + ret_val = 0.7f; + else if (ret_val > 1.5f) + ret_val = 1.5f; // allow to be slightly off bounds ...? int c = +3 * i_05 - 3 * i0 + 1 * i05; max = float(c) + float(a) * ret_val * ret_val + float(b) * ret_val; max /= 1024; @@ -1870,56 +1855,56 @@ BriskScaleSpace::subpixel2D(const int s_0_0, const int s_0_1, const int s_0_2, c { // the coefficients of the 2d quadratic function least-squares fit: - register int tmp1 = s_0_0 + s_0_2 - 2 * s_1_1 + s_2_0 + s_2_2; - register int coeff1 = 3 * (tmp1 + s_0_1 - ((s_1_0 + s_1_2) << 1) + s_2_1); - register int coeff2 = 3 * (tmp1 - ((s_0_1 + s_2_1) << 1) + s_1_0 + s_1_2); - register int tmp2 = s_0_2 - s_2_0; - register int tmp3 = (s_0_0 + tmp2 - s_2_2); - register int tmp4 = tmp3 - 2 * tmp2; - register int coeff3 = -3 * (tmp3 + s_0_1 - s_2_1); - register int coeff4 = -3 * (tmp4 + s_1_0 - s_1_2); - register int coeff5 = (s_0_0 - s_0_2 - s_2_0 + s_2_2) << 2; - register int coeff6 = -(s_0_0 + s_0_2 - ((s_1_0 + s_0_1 + s_1_2 + s_2_1) << 1) - 5 * s_1_1 + s_2_0 + s_2_2) << 1; + int tmp1 = s_0_0 + s_0_2 - 2 * s_1_1 + s_2_0 + s_2_2; + int coeff1 = 3 * (tmp1 + s_0_1 - ((s_1_0 + s_1_2) << 1) + s_2_1); + int coeff2 = 3 * (tmp1 - ((s_0_1 + s_2_1) << 1) + s_1_0 + s_1_2); + int tmp2 = s_0_2 - s_2_0; + int tmp3 = (s_0_0 + tmp2 - s_2_2); + int tmp4 = tmp3 - 2 * tmp2; + int coeff3 = -3 * (tmp3 + s_0_1 - s_2_1); + int coeff4 = -3 * (tmp4 + s_1_0 - s_1_2); + int coeff5 = (s_0_0 - s_0_2 - s_2_0 + s_2_2) << 2; + int coeff6 = -(s_0_0 + s_0_2 - ((s_1_0 + s_0_1 + s_1_2 + s_2_1) << 1) - 5 * s_1_1 + s_2_0 + s_2_2) << 1; // 2nd derivative test: - register int H_det = 4 * coeff1 * coeff2 - coeff5 * coeff5; + int H_det = 4 * coeff1 * coeff2 - coeff5 * coeff5; if (H_det == 0) { - delta_x = 0.0; - delta_y = 0.0; - return float(coeff6) / 18.0; + delta_x = 0.0f; + delta_y = 0.0f; + return float(coeff6) / 18.0f; } if (!(H_det > 0 && coeff1 < 0)) { // The maximum must be at the one of the 4 patch corners. int tmp_max = coeff3 + coeff4 + coeff5; - delta_x = 1.0; - delta_y = 1.0; + delta_x = 1.0f; + delta_y = 1.0f; int tmp = -coeff3 + coeff4 - coeff5; if (tmp > tmp_max) { tmp_max = tmp; - delta_x = -1.0; - delta_y = 1.0; + delta_x = -1.0f; + delta_y = 1.0f; } tmp = coeff3 - coeff4 - coeff5; if (tmp > tmp_max) { tmp_max = tmp; - delta_x = 1.0; - delta_y = -1.0; + delta_x = 1.0f; + delta_y = -1.0f; } tmp = -coeff3 - coeff4 + coeff5; if (tmp > tmp_max) { tmp_max = tmp; - delta_x = -1.0; - delta_y = -1.0; + delta_x = -1.0f; + delta_y = -1.0f; } - return float(tmp_max + coeff1 + coeff2 + coeff6) / 18.0; + return float(tmp_max + coeff1 + coeff2 + coeff6) / 18.0f; } // this is hopefully the normal outcome of the Hessian test @@ -1942,50 +1927,50 @@ BriskScaleSpace::subpixel2D(const int s_0_0, const int s_0_1, const int s_0_2, c if (tx || tx_ || ty || ty_) { // get two candidates: - float delta_x1 = 0.0, delta_x2 = 0.0, delta_y1 = 0.0, delta_y2 = 0.0; + float delta_x1 = 0.0f, delta_x2 = 0.0f, delta_y1 = 0.0f, delta_y2 = 0.0f; if (tx) { - delta_x1 = 1.0; + delta_x1 = 1.0f; delta_y1 = -float(coeff4 + coeff5) / float(2 * coeff2); - if (delta_y1 > 1.0) - delta_y1 = 1.0; - else if (delta_y1 < -1.0) - delta_y1 = -1.0; + if (delta_y1 > 1.0f) + delta_y1 = 1.0f; + else if (delta_y1 < -1.0f) + delta_y1 = -1.0f; } else if (tx_) { - delta_x1 = -1.0; + delta_x1 = -1.0f; delta_y1 = -float(coeff4 - coeff5) / float(2 * coeff2); - if (delta_y1 > 1.0) - delta_y1 = 1.0; + if (delta_y1 > 1.0f) + delta_y1 = 1.0f; else if (delta_y1 < -1.0) - delta_y1 = -1.0; + delta_y1 = -1.0f; } if (ty) { - delta_y2 = 1.0; + delta_y2 = 1.0f; delta_x2 = -float(coeff3 + coeff5) / float(2 * coeff1); - if (delta_x2 > 1.0) - delta_x2 = 1.0; - else if (delta_x2 < -1.0) - delta_x2 = -1.0; + if (delta_x2 > 1.0f) + delta_x2 = 1.0f; + else if (delta_x2 < -1.0f) + delta_x2 = -1.0f; } else if (ty_) { - delta_y2 = -1.0; + delta_y2 = -1.0f; delta_x2 = -float(coeff3 - coeff5) / float(2 * coeff1); - if (delta_x2 > 1.0) - delta_x2 = 1.0; - else if (delta_x2 < -1.0) - delta_x2 = -1.0; + if (delta_x2 > 1.0f) + delta_x2 = 1.0f; + else if (delta_x2 < -1.0f) + delta_x2 = -1.0f; } // insert both options for evaluation which to pick float max1 = (coeff1 * delta_x1 * delta_x1 + coeff2 * delta_y1 * delta_y1 + coeff3 * delta_x1 + coeff4 * delta_y1 + coeff5 * delta_x1 * delta_y1 + coeff6) - / 18.0; + / 18.0f; float max2 = (coeff1 * delta_x2 * delta_x2 + coeff2 * delta_y2 * delta_y2 + coeff3 * delta_x2 + coeff4 * delta_y2 + coeff5 * delta_x2 * delta_y2 + coeff6) - / 18.0; + / 18.0f; if (max1 > max2) { delta_x = delta_x1; @@ -2003,7 +1988,7 @@ BriskScaleSpace::subpixel2D(const int s_0_0, const int s_0_1, const int s_0_2, c // this is the case of the maximum inside the boundaries: return (coeff1 * delta_x * delta_x + coeff2 * delta_y * delta_y + coeff3 * delta_x + coeff4 * delta_y + coeff5 * delta_x * delta_y + coeff6) - / 18.0; + / 18.0f; } // construct a layer @@ -2016,8 +2001,8 @@ BriskLayer::BriskLayer(const cv::Mat& img_in, float scale_in, float offset_in) offset_ = offset_in; // create an agast detector fast_9_16_ = new FastFeatureDetector(1, true, FastFeatureDetector::TYPE_9_16); - makeOffsets(pixel_5_8_, img_.step, 8); - makeOffsets(pixel_9_16_, img_.step, 16); + makeOffsets(pixel_5_8_, (int)img_.step, 8); + makeOffsets(pixel_9_16_, (int)img_.step, 16); } // derive a layer BriskLayer::BriskLayer(const BriskLayer& layer, int mode) @@ -2027,68 +2012,69 @@ BriskLayer::BriskLayer(const BriskLayer& layer, int mode) img_.create(layer.img().rows / 2, layer.img().cols / 2, CV_8U); halfsample(layer.img(), img_); scale_ = layer.scale() * 2; - offset_ = 0.5 * scale_ - 0.5; + offset_ = 0.5f * scale_ - 0.5f; } else { img_.create(2 * (layer.img().rows / 3), 2 * (layer.img().cols / 3), CV_8U); twothirdsample(layer.img(), img_); - scale_ = layer.scale() * 1.5; - offset_ = 0.5 * scale_ - 0.5; + scale_ = layer.scale() * 1.5f; + offset_ = 0.5f * scale_ - 0.5f; } scores_ = cv::Mat::zeros(img_.rows, img_.cols, CV_8U); fast_9_16_ = new FastFeatureDetector(1, false, FastFeatureDetector::TYPE_9_16); - makeOffsets(pixel_5_8_, img_.step, 8); - makeOffsets(pixel_9_16_, img_.step, 16); + makeOffsets(pixel_5_8_, (int)img_.step, 8); + makeOffsets(pixel_9_16_, (int)img_.step, 16); } // Fast/Agast // wraps the agast class void -BriskLayer::getAgastPoints(uint8_t threshold, std::vector& keypoints) +BriskLayer::getAgastPoints(int threshold, std::vector& keypoints) { fast_9_16_->set("threshold", threshold); fast_9_16_->detect(img_, keypoints); // also write scores - const int num = keypoints.size(); + const size_t num = keypoints.size(); - for (int i = 0; i < num; i++) - scores_(keypoints[i].pt.y, keypoints[i].pt.x) = keypoints[i].response; + for (size_t i = 0; i < num; i++) + scores_((int)keypoints[i].pt.y, (int)keypoints[i].pt.x) = saturate_cast(keypoints[i].response); } -inline uint8_t -BriskLayer::getAgastScore(int x, int y, uint8_t threshold) const + +inline int +BriskLayer::getAgastScore(int x, int y, int threshold) const { if (x < 3 || y < 3) return 0; if (x >= img_.cols - 3 || y >= img_.rows - 3) return 0; - uint8_t& score = *(scores_.data + x + y * scores_.cols); + uchar& score = (uchar&)scores_(y, x); if (score > 2) { return score; } - score = cornerScore<16>(img_.data + x + y * img_.cols, pixel_9_16_, threshold - 1); + score = (uchar)cornerScore<16>(&img_.at(y, x), pixel_9_16_, threshold - 1); if (score < threshold) score = 0; return score; } -inline uint8_t -BriskLayer::getAgastScore_5_8(int x, int y, uint8_t threshold) const +inline int +BriskLayer::getAgastScore_5_8(int x, int y, int threshold) const { if (x < 2 || y < 2) return 0; if (x >= img_.cols - 2 || y >= img_.rows - 2) return 0; - uint8_t score = cornerScore<8>(img_.data + x + y * img_.cols, pixel_5_8_, threshold - 1); + int score = cornerScore<8>(&img_.at(y, x), pixel_5_8_, threshold - 1); if (score < threshold) score = 0; return score; } -inline uint8_t -BriskLayer::getAgastScore(float xf, float yf, uint8_t threshold_in, float scale_in) const +inline int +BriskLayer::getAgastScore(float xf, float yf, int threshold_in, float scale_in) const { if (scale_in <= 1.0f) { @@ -2100,8 +2086,8 @@ BriskLayer::getAgastScore(float xf, float yf, uint8_t threshold_in, float scale_ const float ry1 = yf - float(y); const float ry = 1.0f - ry1; - return rx * ry * getAgastScore(x, y, threshold_in) + rx1 * ry * getAgastScore(x + 1, y, threshold_in) - + rx * ry1 * getAgastScore(x, y + 1, threshold_in) + rx1 * ry1 * getAgastScore(x + 1, y + 1, threshold_in); + return (uchar)(rx * ry * getAgastScore(x, y, threshold_in) + rx1 * ry * getAgastScore(x + 1, y, threshold_in) + + rx * ry1 * getAgastScore(x, y + 1, threshold_in) + rx1 * ry1 * getAgastScore(x + 1, y + 1, threshold_in)); } else { @@ -2121,26 +2107,26 @@ BriskLayer::getAgastScore(float xf, float yf, uint8_t threshold_in, float scale_ } // access gray values (smoothed/interpolated) -inline uint8_t +inline int BriskLayer::value(const cv::Mat& mat, float xf, float yf, float scale_in) const { assert(!mat.empty()); // get the position - const int x = floor(xf); - const int y = floor(yf); + const int x = cvFloor(xf); + const int y = cvFloor(yf); const cv::Mat& image = mat; const int& imagecols = image.cols; // get the sigma_half: const float sigma_half = scale_in / 2; - const float area = 4.0 * sigma_half * sigma_half; + const float area = 4.0f * sigma_half * sigma_half; // calculate output: int ret_val; if (sigma_half < 0.5) { //interpolation multipliers: - const int r_x = (xf - x) * 1024; - const int r_y = (yf - y) * 1024; + const int r_x = (int)((xf - x) * 1024); + const int r_y = (int)((yf - y) * 1024); const int r_x_1 = (1024 - r_x); const int r_y_1 = (1024 - r_y); uchar* ptr = image.data + x + y * imagecols; @@ -2158,8 +2144,8 @@ BriskLayer::value(const cv::Mat& mat, float xf, float yf, float scale_in) const // this is the standard case (simple, not speed optimized yet): // scaling: - const int scaling = 4194304.0 / area; - const int scaling2 = float(scaling) * area / 1024.0; + const int scaling = (int)(4194304.0f / area); + const int scaling2 = (int)(float(scaling) * area / 1024.0f); // calculate borders const float x_1 = xf - sigma_half; @@ -2173,20 +2159,20 @@ BriskLayer::value(const cv::Mat& mat, float xf, float yf, float scale_in) const const int y_bottom = int(y1 + 0.5); // overlap area - multiplication factors: - const float r_x_1 = float(x_left) - x_1 + 0.5; - const float r_y_1 = float(y_top) - y_1 + 0.5; - const float r_x1 = x1 - float(x_right) + 0.5; - const float r_y1 = y1 - float(y_bottom) + 0.5; + const float r_x_1 = float(x_left) - x_1 + 0.5f; + const float r_y_1 = float(y_top) - y_1 + 0.5f; + const float r_x1 = x1 - float(x_right) + 0.5f; + const float r_y1 = y1 - float(y_bottom) + 0.5f; const int dx = x_right - x_left - 1; const int dy = y_bottom - y_top - 1; - const int A = (r_x_1 * r_y_1) * scaling; - const int B = (r_x1 * r_y_1) * scaling; - const int C = (r_x1 * r_y1) * scaling; - const int D = (r_x_1 * r_y1) * scaling; - const int r_x_1_i = r_x_1 * scaling; - const int r_y_1_i = r_y_1 * scaling; - const int r_x1_i = r_x1 * scaling; - const int r_y1_i = r_y1 * scaling; + const int A = (int)((r_x_1 * r_y_1) * scaling); + const int B = (int)((r_x1 * r_y_1) * scaling); + const int C = (int)((r_x1 * r_y1) * scaling); + const int D = (int)((r_x_1 * r_y1) * scaling); + const int r_x_1_i = (int)(r_x_1 * scaling); + const int r_y_1_i = (int)(r_y_1 * scaling); + const int r_x1_i = (int)(r_x1 * scaling); + const int r_y1_i = (int)(r_y1 * scaling); // now the calculation: uchar* ptr = image.data + x_left + imagecols * y_top; diff --git a/modules/imgproc/src/imgwarp.cpp b/modules/imgproc/src/imgwarp.cpp index ac0bb61f1..cbe9e29b1 100644 --- a/modules/imgproc/src/imgwarp.cpp +++ b/modules/imgproc/src/imgwarp.cpp @@ -1314,7 +1314,7 @@ public: ssize.width *= cn; int dy, dx, k = 0; - VecOp vop(scale_x, scale_y, src.channels(), src.step/*, area_ofs*/); + VecOp vop(scale_x, scale_y, src.channels(), (int)src.step/*, area_ofs*/); for( dy = range.start; dy < range.end; dy++ ) {