{

threshold = thresh;

octaves = octaves_in;

std::vector<float> rList;

std::vector<int> 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] = (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;

nList[2] = 14;

nList[3] = 15;

nList[4] = 20;

generateKernel(rList, nList, (float)(5.85 * patternScale), (float)(8.2 * patternScale));

std::vector<int> indexChange)

{

generateKernel(radiusList, numberList, dMax, dMin, indexChange);

float dMin, std::vector<int> indexChange)

{

dMax_ = dMax;

dMin_ = dMin;

// get the total number of points

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++)

{

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 = (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.3f;

for (unsigned int scale = 0; scale < scales_; ++scale)

{

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 * 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 * 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.5f;

}

else

{

patternIterator->sigma = (float)(sigma_scale * scaleList_[scale] * (double(radiusList[ring]))

* sin(CV_PI / numberList[ring]));

}

// adapt the sizeList if necessary

const unsigned int size = cvCeil(((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 = (unsigned int)indexChange.size();

if (indSize == 0)

{

indexChange.resize(points_ * (points_ - 1) / 2);

indSize = (unsigned int)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_<indSize);

// make sure the user passes something sensible

BriskShortPair& shortPair = shortPairs_[indexChange[noShortPairs_]];

shortPair.j = j;

shortPair.i = i;

++noShortPairs_;

}

}

}

// no bits:

strings_ = (int) ceil((float(noShortPairs_)) / 128.0) * 4 * 4;

const float key_y, const unsigned int scale, const unsigned int rot,

const unsigned int point) const

{

// get the float position

const BriskPatternPoint& briskPoint = patternPoints_[scale * n_rot_ * points_ + rot * points_ + point];

const float xf = briskPoint.x + key_x;

const float yf = briskPoint.y + key_y;

const int x = int(xf);

const int y = int(yf);

const int& imagecols = image.cols;

// get the sigma:

const float sigma_half = briskPoint.sigma;

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 = (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);

const uchar* ptr = &image.at<uchar>(y, x);

size_t step = image.step;

// just interpolate:

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 = (int)(4194304.0 / area);

const int scaling2 = int(float(scaling) * area / 1024.0);

// the integral image is larger:

const int integralcols = imagecols + 1;

// 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.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 = (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)

{

// 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;

{

const Point2f& pt = keyPt.pt;

return (pt.x < minX) || (pt.x >= maxX) || (pt.y < minY) || (pt.y >= maxY);

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)

{

doOrientation = true;

computeKeypointsNoOrientation(_image, _mask, keypoints);

}

//Remove keypoints very close to the border

size_t ksize = keypoints.size();

std::vector<int> kscales; // remember the scale per keypoint

kscales.resize(ksize);

static const float log2 = 0.693147180559945f;

static const float lb_scalerange = (float)(log(scalerange_) / (log2));

std::vector<cv::KeyPoint>::iterator beginning = keypoints.begin();

std::vector<int>::iterator beginningkscales = kscales.begin();

static const float basicSize06 = basicSize_ * 0.6f;

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((float)border, (float)border, (float)border_x, (float)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:

cv::Mat descriptors;

if (doDescriptors)

{

_descriptors.create((int)ksize, strings_, CV_8U);

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++)

{

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;

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 = (float)(atan2((float) direction1, (float) direction0) / CV_PI * 180.0);

if (kp.angle < 0)

kp.angle += 360.f;

}

if (!doDescriptors)

continue;

int theta;

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

int 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

delete[] _values;

{

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);

{

if (_octaves == 0)

layers_ = 1;

else

layers_ = 2 * _octaves;

{

// 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 (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));

}

{

// make sure keypoints is empty

keypoints.resize(0);

keypoints.reserve(2000);

// assign thresholds

int safeThreshold_ = (int)(threshold_ * safetyFactor_);

std::vector<std::vector<cv::KeyPoint> > agastPoints;

agastPoints.resize(layers_);

// go through the octaves and intra layers and calculate fast corner scores:

for (int 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 size_t num = agastPoints[0].size();

for (size_t n = 0; n < num; n++)

{

const cv::Point2f& point = agastPoints.at(0)[n].pt;

// first check if it is a maximum:

if (!isMax2D(0, (int)point.x, (int)point.y))

continue;

// let's do the subpixel and float scale refinement:

BriskLayer& l = pyramid_[0];

int s_0_0 = l.getAgastScore(point.x - 1, point.y - 1, 1);

int s_1_0 = l.getAgastScore(point.x, point.y - 1, 1);

int s_2_0 = l.getAgastScore(point.x + 1, point.y - 1, 1);

int s_2_1 = l.getAgastScore(point.x + 1, point.y, 1);

int s_1_1 = l.getAgastScore(point.x, point.y, 1);

int s_0_1 = l.getAgastScore(point.x - 1, point.y, 1);

int s_0_2 = l.getAgastScore(point.x - 1, point.y + 1, 1);

int s_1_2 = l.getAgastScore(point.x, point.y + 1, 1);

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 (int i = 0; i < layers_; i++)

{

BriskLayer& l = pyramid_[i];

const size_t num = agastPoints[i].size();

if (i == layers_ - 1)

{

for (size_t n = 0; n < num; n++)

{

const cv::Point2f& point = agastPoints.at(i)[n].pt;

// consider only 2D maxima...

if (!isMax2D(i, (int)point.x, (int)point.y))

continue;

bool ismax;

float dx, dy;

getScoreMaxBelow(i, (int)point.x, (int)point.y, l.getAgastScore(point.x, point.y, safeThreshold_), ismax, dx, dy);

if (!ismax)

continue;

// get the patch on this layer:

int s_0_0 = l.getAgastScore(point.x - 1, point.y - 1, 1);

int s_1_0 = l.getAgastScore(point.x, point.y - 1, 1);

int s_2_0 = l.getAgastScore(point.x + 1, point.y - 1, 1);

int s_2_1 = l.getAgastScore(point.x + 1, point.y, 1);

int s_1_1 = l.getAgastScore(point.x, point.y, 1);

int s_0_1 = l.getAgastScore(point.x - 1, point.y, 1);

int s_0_2 = l.getAgastScore(point.x - 1, point.y + 1, 1);

int s_1_2 = l.getAgastScore(point.x, point.y + 1, 1);

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 (size_t n = 0; n < num; n++)

{

const cv::Point2f& point = agastPoints.at(i)[n].pt;

// first check if it is a maximum:

if (!isMax2D(i, (int)point.x, (int)point.y))

continue;

// let's do the subpixel and float scale refinement:

bool ismax=false;

score = refine3D(i, (int)point.x, (int)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));

}

}

}

}

{

assert(layer<layers_-1);

const BriskLayer& l = pyramid_[layer + 1];

if (layer % 2 == 0)

{ // octave

const int sixths_x = 4 * x_layer - 1;

const int x_above = sixths_x / 6;

const int sixths_y = 4 * y_layer - 1;

const int y_above = sixths_y / 6;

const int r_x = (sixths_x % 6);

const int r_x_1 = 6 - r_x;

const int r_y = (sixths_y % 6);

const int r_y_1 = 6 - r_y;

uchar score = 0xFF

& ((r_x_1 * r_y_1 * l.getAgastScore(x_above, y_above, 1) + r_x * r_y_1

* l.getAgastScore(x_above + 1, y_above, 1)

+ r_x_1 * r_y * l.getAgastScore(x_above, y_above + 1, 1)

+ r_x * r_y * l.getAgastScore(x_above + 1, y_above + 1, 1) + 18)

/ 36);

return score;

}

else

{ // intra

const int eighths_x = 6 * x_layer - 1;

const int x_above = eighths_x / 8;

const int eighths_y = 6 * y_layer - 1;

const int y_above = eighths_y / 8;

const int r_x = (eighths_x % 8);

const int r_x_1 = 8 - r_x;

const int r_y = (eighths_y % 8);

const int r_y_1 = 8 - r_y;

uchar score = 0xFF

& ((r_x_1 * r_y_1 * l.getAgastScore(x_above, y_above, 1) + r_x * r_y_1

* l.getAgastScore(x_above + 1, y_above, 1)

+ r_x_1 * r_y * l.getAgastScore(x_above, y_above + 1, 1)

+ r_x * r_y * l.getAgastScore(x_above + 1, y_above + 1, 1) + 32)

/ 64);

return score;

}

{

assert(layer);

const BriskLayer& l = pyramid_[layer - 1];

int sixth_x;

int quarter_x;

float xf;

int sixth_y;

int quarter_y;

float yf;

// scaling:

float offs;

float area;

int scaling;

int scaling2;

if (layer % 2 == 0)

{ // octave

sixth_x = 8 * x_layer + 1;

xf = float(sixth_x) / 6.0f;

sixth_y = 8 * y_layer + 1;

yf = float(sixth_y) / 6.0f;

// scaling:

offs = 2.0f / 3.0f;

area = 4.0f * offs * offs;

scaling = (int)(4194304.0 / area);

scaling2 = (int)(float(scaling) * area);

}

else

{

quarter_x = 6 * x_layer + 1;

xf = float(quarter_x) / 4.0f;

quarter_y = 6 * y_layer + 1;

yf = float(quarter_y) / 4.0f;

// scaling:

offs = 3.0f / 4.0f;

area = 4.0f * offs * offs;

scaling = (int)(4194304.0 / area);