Skip to content
New issue

Have a question about this project? Sign up for a free GitHub account to open an issue and contact its maintainers and the community.

Sign up for GitHub

By clicking “Sign up for GitHub”, you agree to our terms of service and privacy statement. We’ll occasionally send you account related emails.

Already on GitHub? Sign in to your account

Jump to bottom

Updating STAR detector and FREAK descriptor to work with large and/or 16-bit images #1932

Merged
merged 4 commits into from Apr 10, 2014
Merged

Updating STAR detector and FREAK descriptor to work with large and/or 16-bit images #1932

Show file tree
Hide file tree
Changes from all commits
Commits
Show all changes
4 commits
Select commit Hold shift + click to select a range
2cc11e2
Updating STAR detector and FREAK descriptor to work with large and/or…
Dec 5, 2013
73b0827
Don't automatically convert to grayscale (keep previous behavior.)
Dec 5, 2013
1c15776
Merge branch 'master' of https://github.com/Itseez/opencv
Dec 11, 2013
75ed2f5
Merge branch 'master' of https://github.com/Itseez/opencv
Mar 6, 2014
File filter

Filter by extension

Filter by extension
Conversations
Failed to load comments.
Jump to
Jump to file
Failed to load files.
Diff view
Diff view
12 changes: 10 additions & 2 deletions modules/features2d/include/opencv2/features2d.hpp
View file
Open in desktop
Expand Up @@ -405,8 +405,16 @@ class CV_EXPORTS FREAK : public DescriptorExtractor
protected:
virtual void computeImpl( InputArray image, std::vector<KeyPoint>& keypoints, OutputArray descriptors ) const;
void buildPattern();
uchar meanIntensity( InputArray image, InputArray integral, const float kp_x, const float kp_y,
const unsigned int scale, const unsigned int rot, const unsigned int point ) const;

template <typename imgType, typename iiType>
imgType meanIntensity( InputArray image, InputArray integral, const float kp_x, const float kp_y,
const unsigned int scale, const unsigned int rot, const unsigned int point ) const;

template <typename srcMatType, typename iiMatType>
void computeDescriptors( InputArray image, std::vector<KeyPoint>& keypoints, OutputArray descriptors ) const;

template <typename srcMatType>
void extractDescriptor(srcMatType *pointsValue, void ** ptr) const;

bool orientationNormalized; //true if the orientation is normalized, false otherwise
bool scaleNormalized; //true if the scale is normalized, false otherwise
Expand Down
271 changes: 158 additions & 113 deletions modules/features2d/src/freak.cpp
View file
Open in desktop
Expand Up @@ -239,13 +239,129 @@ void FREAK::computeImpl( InputArray _image, std::vector<KeyPoint>& keypoints, Ou

((FREAK*)this)->buildPattern();

// Convert to gray if not already
Mat grayImage = image;
// if( image.channels() > 1 )
Copy link

Choose a reason for hiding this comment

The reason will be displayed to describe this comment to others. Learn more.

Why these lines are commented?

Copy link
Author

Choose a reason for hiding this comment

The reason will be displayed to describe this comment to others. Learn more.

Because they were causing one of the tests to fail. It's mentioned in here:
#1654

Copy link

Choose a reason for hiding this comment

The reason will be displayed to describe this comment to others. Learn more.

So, may be they should be deleted or there should be some notes, why they are still be in request?

Copy link
Author

Choose a reason for hiding this comment

The reason will be displayed to describe this comment to others. Learn more.

There really wasn't any conclusion (or response) in #1654, so I just commented in order to pass the tests. However, keeping the lines in strikes me as correct behavior. So, the code is currently in logic purgatory.

// cvtColor( image, grayImage, COLOR_BGR2GRAY );

// Use 32-bit integers if we won't overflow in the integral image
if ((image.depth() == CV_8U || image.depth() == CV_8S) &&
(image.rows * image.cols) < 8388608 ) // 8388608 = 2 ^ (32 - 8(bit depth) - 1(sign bit))
{
// Create the integral image appropriate for our type & usage
if (image.depth() == CV_8U)
computeDescriptors<uchar, int>(grayImage, keypoints, _descriptors);
else if (image.depth() == CV_8S)
computeDescriptors<char, int>(grayImage, keypoints, _descriptors);
else
CV_Error( Error::StsUnsupportedFormat, "" );
} else {
// Create the integral image appropriate for our type & usage
if ( image.depth() == CV_8U )
computeDescriptors<uchar, double>(grayImage, keypoints, _descriptors);
else if ( image.depth() == CV_8S )
computeDescriptors<char, double>(grayImage, keypoints, _descriptors);
else if ( image.depth() == CV_16U )
computeDescriptors<ushort, double>(grayImage, keypoints, _descriptors);
else if ( image.depth() == CV_16S )
computeDescriptors<short, double>(grayImage, keypoints, _descriptors);
else
CV_Error( Error::StsUnsupportedFormat, "" );
}
}

template <typename srcMatType>
void FREAK::extractDescriptor(srcMatType *pointsValue, void ** ptr) const
{
std::bitset<FREAK_NB_PAIRS>** ptrScalar = (std::bitset<FREAK_NB_PAIRS>**) ptr;

// extracting descriptor preserving the order of SSE version
int cnt = 0;
for( int n = 7; n < FREAK_NB_PAIRS; n += 128)
{
for( int m = 8; m--; )
{
int nm = n-m;
for(int kk = nm+15*8; kk >= nm; kk-=8, ++cnt)
{
(*ptrScalar)->set(kk, pointsValue[descriptionPairs[cnt].i] >= pointsValue[descriptionPairs[cnt].j]);
}
}
}
--(*ptrScalar);
}

#if CV_SSE2
template <>
void FREAK::extractDescriptor(uchar *pointsValue, void ** ptr) const
{
__m128i** ptrSSE = (__m128i**) ptr;

// note that comparisons order is modified in each block (but first 128 comparisons remain globally the same-->does not affect the 128,384 bits segmanted matching strategy)
int cnt = 0;
for( int n = FREAK_NB_PAIRS/128; n-- ; )
{
__m128i result128 = _mm_setzero_si128();
for( int m = 128/16; m--; cnt += 16 )
{
__m128i operand1 = _mm_set_epi8(pointsValue[descriptionPairs[cnt+0].i],
pointsValue[descriptionPairs[cnt+1].i],
pointsValue[descriptionPairs[cnt+2].i],
pointsValue[descriptionPairs[cnt+3].i],
pointsValue[descriptionPairs[cnt+4].i],
pointsValue[descriptionPairs[cnt+5].i],
pointsValue[descriptionPairs[cnt+6].i],
pointsValue[descriptionPairs[cnt+7].i],
pointsValue[descriptionPairs[cnt+8].i],
pointsValue[descriptionPairs[cnt+9].i],
pointsValue[descriptionPairs[cnt+10].i],
pointsValue[descriptionPairs[cnt+11].i],
pointsValue[descriptionPairs[cnt+12].i],
pointsValue[descriptionPairs[cnt+13].i],
pointsValue[descriptionPairs[cnt+14].i],
pointsValue[descriptionPairs[cnt+15].i]);

__m128i operand2 = _mm_set_epi8(pointsValue[descriptionPairs[cnt+0].j],
pointsValue[descriptionPairs[cnt+1].j],
pointsValue[descriptionPairs[cnt+2].j],
pointsValue[descriptionPairs[cnt+3].j],
pointsValue[descriptionPairs[cnt+4].j],
pointsValue[descriptionPairs[cnt+5].j],
pointsValue[descriptionPairs[cnt+6].j],
pointsValue[descriptionPairs[cnt+7].j],
pointsValue[descriptionPairs[cnt+8].j],
pointsValue[descriptionPairs[cnt+9].j],
pointsValue[descriptionPairs[cnt+10].j],
pointsValue[descriptionPairs[cnt+11].j],
pointsValue[descriptionPairs[cnt+12].j],
pointsValue[descriptionPairs[cnt+13].j],
pointsValue[descriptionPairs[cnt+14].j],
pointsValue[descriptionPairs[cnt+15].j]);

__m128i workReg = _mm_min_epu8(operand1, operand2); // emulated "not less than" for 8-bit UNSIGNED integers
workReg = _mm_cmpeq_epi8(workReg, operand2); // emulated "not less than" for 8-bit UNSIGNED integers

workReg = _mm_and_si128(_mm_set1_epi16(short(0x8080 >> m)), workReg); // merge the last 16 bits with the 128bits std::vector until full
result128 = _mm_or_si128(result128, workReg);
}
(**ptrSSE) = result128;
++(*ptrSSE);
}
(*ptrSSE) -= 8;
}
#endif

template <typename srcMatType, typename iiMatType>
void FREAK::computeDescriptors( InputArray _image, std::vector<KeyPoint>& keypoints, OutputArray _descriptors ) const {

Mat image = _image.getMat();
Mat imgIntegral;
integral(image, imgIntegral);
integral(image, imgIntegral, DataType<iiMatType>::type);
std::vector<int> kpScaleIdx(keypoints.size()); // used to save pattern scale index corresponding to each keypoints
const std::vector<int>::iterator ScaleIdxBegin = kpScaleIdx.begin(); // used in std::vector erase function
const std::vector<cv::KeyPoint>::iterator kpBegin = keypoints.begin(); // used in std::vector erase function
const float sizeCst = static_cast<float>(FREAK_NB_SCALES/(FREAK_LOG2* nOctaves));
uchar pointsValue[FREAK_NB_POINTS];
srcMatType pointsValue[FREAK_NB_POINTS];
int thetaIdx = 0;
int direction0;
int direction1;
Expand Down Expand Up @@ -300,13 +416,10 @@ void FREAK::computeImpl( InputArray _image, std::vector<KeyPoint>& keypoints, Ou
_descriptors.create((int)keypoints.size(), FREAK_NB_PAIRS/8, CV_8U);
_descriptors.setTo(Scalar::all(0));
Mat descriptors = _descriptors.getMat();
#if CV_SSE2
__m128i* ptr= (__m128i*) (descriptors.data+(keypoints.size()-1)*descriptors.step[0]);
#else
std::bitset<FREAK_NB_PAIRS>* ptr = (std::bitset<FREAK_NB_PAIRS>*) (descriptors.data+(keypoints.size()-1)*descriptors.step[0]);
#endif
for( size_t k = keypoints.size(); k--; )
{

void *ptr = descriptors.data+(keypoints.size()-1)*descriptors.step[0];

for( size_t k = keypoints.size(); k--; ) {
// estimate orientation (gradient)
if( !orientationNormalized )
{
Expand All @@ -316,9 +429,10 @@ void FREAK::computeImpl( InputArray _image, std::vector<KeyPoint>& keypoints, Ou
else
{
// get the points intensity value in the un-rotated pattern
for( int i = FREAK_NB_POINTS; i--; )
{
pointsValue[i] = meanIntensity(image, imgIntegral, keypoints[k].pt.x,keypoints[k].pt.y, kpScaleIdx[k], 0, i);
for( int i = FREAK_NB_POINTS; i--; ) {
pointsValue[i] = meanIntensity<srcMatType, iiMatType>(image, imgIntegral,
keypoints[k].pt.x, keypoints[k].pt.y,
kpScaleIdx[k], 0, i);
}
direction0 = 0;
direction1 = 0;
Expand All @@ -339,80 +453,14 @@ void FREAK::computeImpl( InputArray _image, std::vector<KeyPoint>& keypoints, Ou
thetaIdx -= FREAK_NB_ORIENTATION;
}
// extract descriptor at the computed orientation
for( int i = FREAK_NB_POINTS; i--; )
{
pointsValue[i] = meanIntensity(image, imgIntegral, keypoints[k].pt.x,keypoints[k].pt.y, kpScaleIdx[k], thetaIdx, i);
for( int i = FREAK_NB_POINTS; i--; ) {
pointsValue[i] = meanIntensity<srcMatType, iiMatType>(image, imgIntegral,
keypoints[k].pt.x, keypoints[k].pt.y,
kpScaleIdx[k], thetaIdx, i);
}
#if CV_SSE2
// note that comparisons order is modified in each block (but first 128 comparisons remain globally the same-->does not affect the 128,384 bits segmanted matching strategy)
int cnt = 0;
for( int n = FREAK_NB_PAIRS/128; n-- ; )
{
__m128i result128 = _mm_setzero_si128();
for( int m = 128/16; m--; cnt += 16 )
{
__m128i operand1 = _mm_set_epi8(
pointsValue[descriptionPairs[cnt+0].i],
pointsValue[descriptionPairs[cnt+1].i],
pointsValue[descriptionPairs[cnt+2].i],
pointsValue[descriptionPairs[cnt+3].i],
pointsValue[descriptionPairs[cnt+4].i],
pointsValue[descriptionPairs[cnt+5].i],
pointsValue[descriptionPairs[cnt+6].i],
pointsValue[descriptionPairs[cnt+7].i],
pointsValue[descriptionPairs[cnt+8].i],
pointsValue[descriptionPairs[cnt+9].i],
pointsValue[descriptionPairs[cnt+10].i],
pointsValue[descriptionPairs[cnt+11].i],
pointsValue[descriptionPairs[cnt+12].i],
pointsValue[descriptionPairs[cnt+13].i],
pointsValue[descriptionPairs[cnt+14].i],
pointsValue[descriptionPairs[cnt+15].i]);

__m128i operand2 = _mm_set_epi8(
pointsValue[descriptionPairs[cnt+0].j],
pointsValue[descriptionPairs[cnt+1].j],
pointsValue[descriptionPairs[cnt+2].j],
pointsValue[descriptionPairs[cnt+3].j],
pointsValue[descriptionPairs[cnt+4].j],
pointsValue[descriptionPairs[cnt+5].j],
pointsValue[descriptionPairs[cnt+6].j],
pointsValue[descriptionPairs[cnt+7].j],
pointsValue[descriptionPairs[cnt+8].j],
pointsValue[descriptionPairs[cnt+9].j],
pointsValue[descriptionPairs[cnt+10].j],
pointsValue[descriptionPairs[cnt+11].j],
pointsValue[descriptionPairs[cnt+12].j],
pointsValue[descriptionPairs[cnt+13].j],
pointsValue[descriptionPairs[cnt+14].j],
pointsValue[descriptionPairs[cnt+15].j]);

__m128i workReg = _mm_min_epu8(operand1, operand2); // emulated "not less than" for 8-bit UNSIGNED integers
workReg = _mm_cmpeq_epi8(workReg, operand2); // emulated "not less than" for 8-bit UNSIGNED integers

workReg = _mm_and_si128(_mm_set1_epi16(short(0x8080 >> m)), workReg); // merge the last 16 bits with the 128bits std::vector until full
result128 = _mm_or_si128(result128, workReg);
}
(*ptr) = result128;
++ptr;
}
ptr -= 8;
#else
// extracting descriptor preserving the order of SSE version
int cnt = 0;
for( int n = 7; n < FREAK_NB_PAIRS; n += 128)
{
for( int m = 8; m--; )
{
int nm = n-m;
for(int kk = nm+15*8; kk >= nm; kk-=8, ++cnt)
{
ptr->set(kk, pointsValue[descriptionPairs[cnt].i] >= pointsValue[descriptionPairs[cnt].j]);
}
}
}
--ptr;
#endif

// Extract descriptor
extractDescriptor<srcMatType>(pointsValue, &ptr);
}
}
else // extract all possible comparisons for selection
Expand All @@ -434,7 +482,9 @@ void FREAK::computeImpl( InputArray _image, std::vector<KeyPoint>& keypoints, Ou
{
//get the points intensity value in the un-rotated pattern
for( int i = FREAK_NB_POINTS;i--; )
pointsValue[i] = meanIntensity(image, imgIntegral, keypoints[k].pt.x,keypoints[k].pt.y, kpScaleIdx[k], 0, i);
pointsValue[i] = meanIntensity<srcMatType, iiMatType>(image, imgIntegral,
keypoints[k].pt.x,keypoints[k].pt.y,
kpScaleIdx[k], 0, i);

direction0 = 0;
direction1 = 0;
Expand All @@ -456,10 +506,10 @@ void FREAK::computeImpl( InputArray _image, std::vector<KeyPoint>& keypoints, Ou
thetaIdx -= FREAK_NB_ORIENTATION;
}
// get the points intensity value in the rotated pattern
for( int i = FREAK_NB_POINTS; i--; )
{
pointsValue[i] = meanIntensity(image, imgIntegral, keypoints[k].pt.x,
keypoints[k].pt.y, kpScaleIdx[k], thetaIdx, i);
for( int i = FREAK_NB_POINTS; i--; ) {
pointsValue[i] = meanIntensity<srcMatType, iiMatType>(image, imgIntegral,
keypoints[k].pt.x, keypoints[k].pt.y,
kpScaleIdx[k], thetaIdx, i);
}

int cnt(0);
Expand All @@ -478,21 +528,20 @@ void FREAK::computeImpl( InputArray _image, std::vector<KeyPoint>& keypoints, Ou
}

// simply take average on a square patch, not even gaussian approx
uchar FREAK::meanIntensity( InputArray _image, InputArray _integral,
const float kp_x,
const float kp_y,
const unsigned int scale,
const unsigned int rot,
const unsigned int point) const
{
template <typename imgType, typename iiType>
imgType FREAK::meanIntensity( InputArray _image, InputArray _integral,
const float kp_x,
const float kp_y,
const unsigned int scale,
const unsigned int rot,
const unsigned int point) const {
Mat image = _image.getMat(), integral = _integral.getMat();
// get point position in image
const PatternPoint& FreakPoint = patternLookup[scale*FREAK_NB_ORIENTATION*FREAK_NB_POINTS + rot*FREAK_NB_POINTS + point];
const float xf = FreakPoint.x+kp_x;
const float yf = FreakPoint.y+kp_y;
const int x = int(xf);
const int y = int(yf);
const int& imagecols = image.cols;

// get the sigma:
const float radius = FreakPoint.sigma;
Expand All @@ -505,19 +554,15 @@ uchar FREAK::meanIntensity( InputArray _image, InputArray _integral,
const int r_y = static_cast<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;
unsigned int ret_val;
// linear interpolation:
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*int(image.at<imgType>(y , x ))
+ r_x *r_y_1*int(image.at<imgType>(y , x+1))
+ r_x_1*r_y *int(image.at<imgType>(y+1, x ))
+ r_x *r_y *int(image.at<imgType>(y+1, x+1));
//return the rounded mean
ret_val += 2 * 1024 * 1024;
return static_cast<uchar>(ret_val / (4 * 1024 * 1024));
return static_cast<imgType>(ret_val / (4 * 1024 * 1024));
}

// expected case:
Expand All @@ -527,15 +572,15 @@ uchar FREAK::meanIntensity( InputArray _image, InputArray _integral,
const int y_top = int(yf-radius+0.5);
const int x_right = int(xf+radius+1.5);//integral image is 1px wider
const int y_bottom = int(yf+radius+1.5);//integral image is 1px higher
int ret_val;
iiType ret_val;

ret_val = integral.at<int>(y_bottom,x_right);//bottom right corner
ret_val -= integral.at<int>(y_bottom,x_left);
ret_val += integral.at<int>(y_top,x_left);
ret_val -= integral.at<int>(y_top,x_right);
ret_val = integral.at<iiType>(y_bottom,x_right);//bottom right corner
ret_val -= integral.at<iiType>(y_bottom,x_left);
ret_val += integral.at<iiType>(y_top,x_left);
ret_val -= integral.at<iiType>(y_top,x_right);
ret_val = ret_val/( (x_right-x_left)* (y_bottom-y_top) );
//~ std::cout<<integral.step[1]<<std::endl;
return static_cast<uchar>(ret_val);
return static_cast<imgType>(ret_val);
}

// pair selection algorithm from a set of training images and corresponding keypoints
Expand Down