ncycles_ = 0;

reordering_ = true;

if (options_.descriptor_size > 0 && options_.descriptor >= AKAZE::DESCRIPTOR_MLDB_UPRIGHT) {

generateDescriptorSubsample(descriptorSamples_, descriptorBits_, options_.descriptor_size,

options_.descriptor_pattern_size, options_.descriptor_channels);

}

Allocate_Memory_Evolution();

CV_INSTRUMENT_REGION()

float rfactor = 0.0f;

int level_height = 0, level_width = 0;

// maximum size of the area for the descriptor computation

float smax = 0.0;

if (options_.descriptor == AKAZE::DESCRIPTOR_MLDB_UPRIGHT || options_.descriptor == AKAZE::DESCRIPTOR_MLDB) {

smax = 10.0f*sqrtf(2.0f);

}

else if (options_.descriptor == AKAZE::DESCRIPTOR_KAZE_UPRIGHT || options_.descriptor == AKAZE::DESCRIPTOR_KAZE) {

smax = 12.0f*sqrtf(2.0f);

}

// Allocate the dimension of the matrices for the evolution

for (int i = 0, power = 1; i <= options_.omax - 1; i++, power *= 2) {

rfactor = 1.0f / power;

level_height = (int)(options_.img_height*rfactor);

level_width = (int)(options_.img_width*rfactor);

// Smallest possible octave and allow one scale if the image is small

if ((level_width < 80 || level_height < 40) && i != 0) {

options_.omax = i;

break;

}

for (int j = 0; j < options_.nsublevels; j++) {

MEvolution step;

step.size = Size(level_width, level_height);

step.esigma = options_.soffset*pow(2.f, (float)(j) / (float)(options_.nsublevels) + i);

step.sigma_size = cvRound(step.esigma * options_.derivative_factor / power); // In fact sigma_size only depends on j

step.etime = 0.5f * (step.esigma * step.esigma);

step.octave = i;

step.sublevel = j;

step.octave_ratio = (float)power;

step.border = cvRound(smax * step.sigma_size) + 1;

evolution_.push_back(step);

}

}

// Allocate memory for the number of cycles and time steps

for (size_t i = 1; i < evolution_.size(); i++) {

int naux = 0;

vector<float> tau;

float ttime = 0.0f;

ttime = evolution_[i].etime - evolution_[i - 1].etime;

naux = fed_tau_by_process_time(ttime, 1, 0.25f, reordering_, tau);

nsteps_.push_back(naux);

tsteps_.push_back(tau);

ncycles_++;

}

// Compute an appropriate kernel size according to the specified sigma

int ksize = (int)cvCeil(2.0f*(1.0f + (sigma - 0.8f) / (0.3f)));

ksize |= 1; // kernel should be odd

return ksize;

{

CV_INSTRUMENT_REGION()

/* The labeling scheme for this five star stencil:

Lstep.create(Lt.size(), Lt.type());

const int cols = Lt.cols - 2;

int row = row_begin;

const float *lt_a, *lt_c, *lt_b;

const float *lf_a, *lf_c, *lf_b;

float *dst;

float step_r = 0.f;

// Process the top row

if (row == 0) {

lt_c = Lt.ptr<float>(0) + 1; /* Skip the left-most column by +1 */

lf_c = Lf.ptr<float>(0) + 1;

lt_b = Lt.ptr<float>(1) + 1;

lf_b = Lf.ptr<float>(1) + 1;

// fill the corner to prevent uninitialized values

dst = Lstep.ptr<float>(0);

dst[0] = 0.0f;

++dst;

for (int j = 0; j < cols; j++) {

step_r = (lf_c[j] + lf_c[j + 1])*(lt_c[j + 1] - lt_c[j]) +

(lf_c[j] + lf_c[j - 1])*(lt_c[j - 1] - lt_c[j]) +

(lf_c[j] + lf_b[j ])*(lt_b[j ] - lt_c[j]);

dst[j] = step_r * step_size;

}

// fill the corner to prevent uninitialized values

dst[cols] = 0.0f;

++row;

}

// Process the middle rows

int middle_end = std::min(Lt.rows - 1, row_end);

for (; row < middle_end; ++row)

{

lt_a = Lt.ptr<float>(row - 1);

lf_a = Lf.ptr<float>(row - 1);

lt_c = Lt.ptr<float>(row );

lf_c = Lf.ptr<float>(row );

lt_b = Lt.ptr<float>(row + 1);

lf_b = Lf.ptr<float>(row + 1);

dst = Lstep.ptr<float>(row);

// The left-most column

step_r = (lf_c[0] + lf_c[1])*(lt_c[1] - lt_c[0]) +

(lf_c[0] + lf_b[0])*(lt_b[0] - lt_c[0]) +

(lf_c[0] + lf_a[0])*(lt_a[0] - lt_c[0]);

dst[0] = step_r * step_size;

lt_a++; lt_c++; lt_b++;

lf_a++; lf_c++; lf_b++;

dst++;

// The middle columns

for (int j = 0; j < cols; j++)

{

step_r = (lf_c[j] + lf_c[j + 1])*(lt_c[j + 1] - lt_c[j]) +

(lf_c[j] + lf_c[j - 1])*(lt_c[j - 1] - lt_c[j]) +

(lf_c[j] + lf_b[j ])*(lt_b[j ] - lt_c[j]) +

(lf_c[j] + lf_a[j ])*(lt_a[j ] - lt_c[j]);

dst[j] = step_r * step_size;

}

// The right-most column

step_r = (lf_c[cols] + lf_c[cols - 1])*(lt_c[cols - 1] - lt_c[cols]) +

(lf_c[cols] + lf_b[cols ])*(lt_b[cols ] - lt_c[cols]) +

(lf_c[cols] + lf_a[cols ])*(lt_a[cols ] - lt_c[cols]);

dst[cols] = step_r * step_size;

}

// Process the bottom row (row == Lt.rows - 1)

if (row_end == Lt.rows) {

lt_a = Lt.ptr<float>(row - 1) + 1; /* Skip the left-most column by +1 */

lf_a = Lf.ptr<float>(row - 1) + 1;

lt_c = Lt.ptr<float>(row ) + 1;

lf_c = Lf.ptr<float>(row ) + 1;

// fill the corner to prevent uninitialized values

dst = Lstep.ptr<float>(row);

dst[0] = 0.0f;

++dst;

for (int j = 0; j < cols; j++) {

step_r = (lf_c[j] + lf_c[j + 1])*(lt_c[j + 1] - lt_c[j]) +

(lf_c[j] + lf_c[j - 1])*(lt_c[j - 1] - lt_c[j]) +

(lf_c[j] + lf_a[j ])*(lt_a[j ] - lt_c[j]);

dst[j] = step_r * step_size;

}

// fill the corner to prevent uninitialized values

dst[cols] = 0.0f;

}

{

NonLinearScalarDiffusionStep(const Mat& Lt, const Mat& Lf, Mat& Lstep, float step_size)

: Lt_(&Lt), Lf_(&Lf), Lstep_(&Lstep), step_size_(step_size)

{}

void operator()(const Range& range) const

{

nld_step_scalar_one_lane(*Lt_, *Lf_, *Lstep_, step_size_, range.start, range.end);

}

const Mat* Lt_;

const Mat* Lf_;

Mat* Lstep_;

float step_size_;

{

if(!Lt_.isContinuous())

return false;

UMat Lt = Lt_.getUMat();

UMat Lf = Lf_.getUMat();

UMat Lstep = Lstep_.getUMat();

size_t globalSize[] = {(size_t)Lt.cols, (size_t)Lt.rows};

ocl::Kernel ker("AKAZE_nld_step_scalar", ocl::features2d::akaze_oclsrc);

if( ker.empty() )

return false;

return ker.args(

ocl::KernelArg::ReadOnly(Lt),

ocl::KernelArg::PtrReadOnly(Lf),

ocl::KernelArg::PtrWriteOnly(Lstep),

step_size).run(2, globalSize, 0, true);

{

CV_INSTRUMENT_REGION()

Lstep_.create(Lt_.size(), Lt_.type());

CV_OCL_RUN(Lt_.isUMat() && Lf_.isUMat() && Lstep_.isUMat(),

ocl_non_linear_diffusion_step(Lt_, Lf_, Lstep_, step_size));

Mat Lt = Lt_.getMat();

Mat Lf = Lf_.getMat();

Mat Lstep = Lstep_.getMat();

parallel_for_(Range(0, Lt.rows), NonLinearScalarDiffusionStep(Lt, Lf, Lstep, step_size));

{

CV_INSTRUMENT_REGION()

CV_Assert(nbins > 2);

CV_Assert(!Lx_.empty());

Mat Lx = Lx_.getMat();

Mat Ly = Ly_.getMat();

// temporary square roots of dot product

Mat modgs (Lx.rows - 2, Lx.cols - 2, CV_32F);

const int total = modgs.cols * modgs.rows;

float *modg = modgs.ptr<float>();

float hmax = 0.0f;

for (int i = 1; i < Lx.rows - 1; i++) {

const float *lx = Lx.ptr<float>(i) + 1;

const float *ly = Ly.ptr<float>(i) + 1;

const int cols = Lx.cols - 2;

for (int j = 0; j < cols; j++) {

float dist = sqrtf(lx[j] * lx[j] + ly[j] * ly[j]);

*modg++ = dist;

hmax = std::max(hmax, dist);

}

}

modg = modgs.ptr<float>();

if (hmax == 0.0f)

return 0.03f; // e.g. a blank image

// Compute the bin numbers: the value range [0, hmax] -> [0, nbins-1]

modgs *= (nbins - 1) / hmax;

// Count up histogram

std::vector<int> hist(nbins, 0);

for (int i = 0; i < total; i++)

hist[(int)modg[i]]++;

// Now find the perc of the histogram percentile

const int nthreshold = (int)((total - hist[0]) * perc); // Exclude hist[0] as background

int nelements = 0;

for (int k = 1; k < nbins; k++) {

if (nelements >= nthreshold)

return (float)hmax * k / nbins;

nelements += hist[k];

}

return 0.03f;

{

UMat Lx = Lx_.getUMat();

UMat Ly = Ly_.getUMat();

UMat Lflow = Lflow_.getUMat();

int total = Lx.rows * Lx.cols;

size_t globalSize[] = {(size_t)total};

ocl::Kernel ker("AKAZE_pm_g2", ocl::features2d::akaze_oclsrc);

if( ker.empty() )

return false;

return ker.args(

ocl::KernelArg::PtrReadOnly(Lx),

ocl::KernelArg::PtrReadOnly(Ly),

ocl::KernelArg::PtrWriteOnly(Lflow),

kcontrast, total).run(1, globalSize, 0, true);

{

CV_INSTRUMENT_REGION()

Lflow.create(Lx.size(), Lx.type());

switch (diffusivity) {

case KAZE::DIFF_PM_G1:

pm_g1(Lx, Ly, Lflow, kcontrast);

break;

case KAZE::DIFF_PM_G2:

CV_OCL_RUN(Lx.isUMat() && Ly.isUMat() && Lflow.isUMat(), ocl_pm_g2(Lx, Ly, Lflow, kcontrast));

pm_g2(Lx, Ly, Lflow, kcontrast);

break;

case KAZE::DIFF_WEICKERT:

weickert_diffusivity(Lx, Ly, Lflow, kcontrast);

break;

case KAZE::DIFF_CHARBONNIER:

charbonnier_diffusivity(Lx, Ly, Lflow, kcontrast);

break;

default:

CV_Error(diffusivity, "Diffusivity is not supported");

break;

}

{

Mat img = image.getMat();

if (img.channels() > 1)

cvtColor(image, img, COLOR_BGR2GRAY);

if ( img.depth() == CV_32F )

dst.assign(img);

else if ( img.depth() == CV_8U )

img.convertTo(dst, CV_32F, 1.0 / 255.0, 0);

else if ( img.depth() == CV_16U )

img.convertTo(dst, CV_32F, 1.0 / 65535.0, 0);

const std::vector<std::vector<float > > &tsteps_evolution, std::vector<Evolution<MatType> > &evolution)

{

CV_INSTRUMENT_REGION()

CV_Assert(evolution.size() > 0);

// convert input to grayscale float image if needed

MatType img;

prepareInputImage(image, img);

// create first level of the evolution

int ksize = getGaussianKernelSize(options.soffset);

GaussianBlur(img, evolution[0].Lsmooth, Size(ksize, ksize), options.soffset, options.soffset, BORDER_REPLICATE);

evolution[0].Lsmooth.copyTo(evolution[0].Lt);

if (evolution.size() == 1) {

// we don't need to compute kcontrast factor

Compute_Determinant_Hessian_Response(evolution);

return;

}

// derivatives, flow and diffusion step

MatType Lx, Ly, Lsmooth, Lflow, Lstep;

// compute derivatives for computing k contrast

GaussianBlur(img, Lsmooth, Size(5, 5), 1.0f, 1.0f, BORDER_REPLICATE);

Scharr(Lsmooth, Lx, CV_32F, 1, 0, 1, 0, BORDER_DEFAULT);

Scharr(Lsmooth, Ly, CV_32F, 0, 1, 1, 0, BORDER_DEFAULT);

Lsmooth.release();

// compute the kcontrast factor

float kcontrast = compute_kcontrast(Lx, Ly, options.kcontrast_percentile, options.kcontrast_nbins);

// Now generate the rest of evolution levels

for (size_t i = 1; i < evolution.size(); i++) {

Evolution<MatType> &e = evolution[i];

if (e.octave > evolution[i - 1].octave) {

// new octave will be half the size

resize(evolution[i - 1].Lt, e.Lt, e.size, 0, 0, INTER_AREA);

kcontrast *= 0.75f;

}

else {

evolution[i - 1].Lt.copyTo(e.Lt);

}

GaussianBlur(e.Lt, e.Lsmooth, Size(5, 5), 1.0f, 1.0f, BORDER_REPLICATE);

// Compute the Gaussian derivatives Lx and Ly

Scharr(e.Lsmooth, Lx, CV_32F, 1, 0, 1.0, 0, BORDER_DEFAULT);

Scharr(e.Lsmooth, Ly, CV_32F, 0, 1, 1.0, 0, BORDER_DEFAULT);

// Compute the conductivity equation

compute_diffusivity(Lx, Ly, Lflow, kcontrast, options.diffusivity);

// Perform Fast Explicit Diffusion on Lt

const std::vector<float> &tsteps = tsteps_evolution[i - 1];

for (size_t j = 0; j < tsteps.size(); j++) {

const float step_size = tsteps[j] * 0.5f;

non_linear_diffusion_step(e.Lt, Lflow, Lstep, step_size);

add(e.Lt, Lstep, e.Lt);

}

}

Compute_Determinant_Hessian_Response(evolution);

return;

{

dst.resize(src.size());

for (size_t i = 0; i < src.size(); ++i) {

dst[i] = Evolution<MatTypeDst>(src[i]);

}

{

if (ocl::isOpenCLActivated() && image.isUMat()) {

// will run OCL version of scale space pyramid

UMatPyramid uPyr;

// init UMat pyramid with sizes

convertScalePyramid(evolution_, uPyr);

create_nonlinear_scale_space(image, options_, tsteps_, uPyr);

// download pyramid from GPU

convertScalePyramid(uPyr, evolution_);

} else {

// CPU version

create_nonlinear_scale_space(image, options_, tsteps_, evolution_);

}

OutputArray Ldet_, float sigma)

{

UMat Lxx = Lxx_.getUMat();

UMat Lxy = Lxy_.getUMat();

UMat Lyy = Lyy_.getUMat();

UMat Ldet = Ldet_.getUMat();

const int total = Lxx.rows * Lxx.cols;

size_t globalSize[] = {(size_t)total};

ocl::Kernel ker("AKAZE_compute_determinant", ocl::features2d::akaze_oclsrc);

if( ker.empty() )

return false;

return ker.args(

ocl::KernelArg::PtrReadOnly(Lxx),

ocl::KernelArg::PtrReadOnly(Lxy),

ocl::KernelArg::PtrReadOnly(Lyy),

ocl::KernelArg::PtrWriteOnly(Ldet),

sigma, total).run(1, globalSize, 0, true);

OutputArray Ldet_, float sigma)

{

CV_INSTRUMENT_REGION()

Ldet_.create(Lxx_.size(), Lxx_.type());

CV_OCL_RUN(Lxx_.isUMat() && Ldet_.isUMat(), ocl_compute_determinant(Lxx_, Lxy_, Lyy_, Ldet_, sigma));

// output determinant

Mat Lxx = Lxx_.getMat(), Lxy = Lxy_.getMat(), Lyy = Lyy_.getMat(), Ldet = Ldet_.getMat();

float *lxx = Lxx.ptr<float>();

float *lxy = Lxy.ptr<float>();

float *lyy = Lyy.ptr<float>();

float *ldet = Ldet.ptr<float>();

const int total = Lxx.cols * Lxx.rows;

for (int j = 0; j < total; j++) {

ldet[j] = (lxx[j] * lyy[j] - lxy[j] * lxy[j]) * sigma;

}

{

explicit DeterminantHessianResponse(std::vector<Evolution<MatType> >& ev)

: evolution_(&ev)

{

}

void operator()(const Range& range) const

{

MatType Lxx, Lxy, Lyy;

for (int i = range.start; i < range.end; i++)

{

Evolution<MatType> &e = (*evolution_)[i];

// we cannot use cv:Scharr here, because we need to handle also

// kernel sizes other than 3, by default we are using 9x9, 5x5 and 7x7

// compute kernels

Mat DxKx, DxKy, DyKx, DyKy;

compute_derivative_kernels(DxKx, DxKy, 1, 0, e.sigma_size);

compute_derivative_kernels(DyKx, DyKy, 0, 1, e.sigma_size);

// compute the multiscale derivatives

sepFilter2D(e.Lsmooth, e.Lx, CV_32F, DxKx, DxKy);

sepFilter2D(e.Lx, Lxx, CV_32F, DxKx, DxKy);

sepFilter2D(e.Lx, Lxy, CV_32F, DyKx, DyKy);

sepFilter2D(e.Lsmooth, e.Ly, CV_32F, DyKx, DyKy);

sepFilter2D(e.Ly, Lyy, CV_32F, DyKx, DyKy);

// free Lsmooth to same some space in the pyramid, it is not needed anymore

e.Lsmooth.release();

// compute determinant scaled by sigma

float sigma_size_quat = (float)(e.sigma_size * e.sigma_size * e.sigma_size * e.sigma_size);

compute_determinant(Lxx, Lxy, Lyy, e.Ldet, sigma_size_quat);

}

}

std::vector<Evolution<MatType> >* evolution_;

CV_INSTRUMENT_REGION()

DeterminantHessianResponse<UMat> body (evolution);

body(Range(0, (int)evolution.size()));

CV_INSTRUMENT_REGION()

parallel_for_(Range(0, (int)evolution.size()), DeterminantHessianResponse<Mat>(evolution));

{

CV_INSTRUMENT_REGION()

kpts.clear();

std::vector<Mat> keypoints_by_layers;

Find_Scale_Space_Extrema(keypoints_by_layers);

Do_Subpixel_Refinement(keypoints_by_layers, kpts);

Compute_Keypoints_Orientation(kpts);

{

// search neighborhood for keypoints

for (int i = y - search_radius; i < y + search_radius; ++i) {

const uchar *curr = mask.ptr<uchar>(i);

for (int j = x - search_radius; j < x + search_radius; ++j) {

if (curr[j] == 0) {

continue; // skip non-keypoint

}

// fine-compare with L2 metric (L2 is smaller than our search window)

int dx = j - x;

int dy = i - y;

if (dx * dx + dy * dy <= search_radius * search_radius) {

idx = i * mask.cols + j;

return true;

}

}

}

return false;

{

explicit FindKeypointsSameScale(const Pyramid& ev,

std::vector<Mat>& kpts, float dthreshold)

: evolution_(&ev), keypoints_by_layers_(&kpts), dthreshold_(dthreshold)

{}

void operator()(const Range& range) const

{

for (int i = range.start; i < range.end; i++)

{

const MEvolution &e = (*evolution_)[i];

Mat &kpts = (*keypoints_by_layers_)[i];

// this mask will hold positions of keypoints in this level

kpts = Mat::zeros(e.Ldet.size(), CV_8UC1);

// if border is too big we shouldn't search any keypoints

if (e.border + 1 >= e.Ldet.rows)

continue;

const float * prev = e.Ldet.ptr<float>(e.border - 1);

const float * curr = e.Ldet.ptr<float>(e.border );

const float * next = e.Ldet.ptr<float>(e.border + 1);

const float * ldet = e.Ldet.ptr<float>();

uchar *mask = kpts.ptr<uchar>();

const int search_radius = e.sigma_size; // size of keypoint in this level

for (int y = e.border; y < e.Ldet.rows - e.border; y++) {

for (int x = e.border; x < e.Ldet.cols - e.border; x++) {

const float value = curr[x];

// Filter the points with the detector threshold

if (value <= dthreshold_)

continue;

if (value <= curr[x-1] || value <= curr[x+1])

continue;

if (value <= prev[x-1] || value <= prev[x ] || value <= prev[x+1])

continue;

if (value <= next[x-1] || value <= next[x ] || value <= next[x+1])

continue;

int idx = 0;

// Compare response with the same scale

if (find_neighbor_point(x, y, kpts, search_radius, idx)) {

if (value > ldet[idx]) {

mask[idx] = 0; // clear old point - we have better candidate now

} else {

continue; // there already is a better keypoint

}

}

kpts.at<uchar>(y, x) = 1; // we have a new keypoint

}

prev = curr;

curr = next;

next += e.Ldet.cols;

}

}

}

const Pyramid* evolution_;

std::vector<Mat>* keypoints_by_layers_;

float dthreshold_; ///< Detector response threshold to accept point

{

CV_INSTRUMENT_REGION()

keypoints_by_layers.resize(evolution_.size());

// find points in the same level

parallel_for_(Range(0, (int)evolution_.size()),

FindKeypointsSameScale(evolution_, keypoints_by_layers, options_.dthreshold));

// Filter points with the lower scale level

for (size_t i = 1; i < keypoints_by_layers.size(); i++) {

// constants for this level

const Mat &keypoints = keypoints_by_layers[i];

const uchar *const kpts = keypoints_by_layers[i].ptr<uchar>();

uchar *const kpts_prev = keypoints_by_layers[i-1].ptr<uchar>();

const float *const ldet = evolution_[i].Ldet.ptr<float>();

const float *const ldet_prev = evolution_[i-1].Ldet.ptr<float>();

// ratios are just powers of 2

const int diff_ratio = (int)evolution_[i].octave_ratio / (int)evolution_[i-1].octave_ratio;

const int search_radius = evolution_[i].sigma_size * diff_ratio; // size of keypoint in this level

size_t j = 0;

for (int y = 0; y < keypoints.rows; y++) {

for (int x = 0; x < keypoints.cols; x++, j++) {

if (kpts[j] == 0) {

continue; // skip non-keypoints

}

int idx = 0;

// project point to lower scale layer

const int p_x = x * diff_ratio;

const int p_y = y * diff_ratio;

if (find_neighbor_point(p_x, p_y, keypoints_by_layers[i-1], search_radius, idx)) {

if (ldet[j] > ldet_prev[idx]) {

kpts_prev[idx] = 0; // clear keypoint in lower layer

}

// else this pt may be pruned by the upper scale

}

}

}

}

// Now filter points with the upper scale level (the other direction)

for (int i = (int)keypoints_by_layers.size() - 2; i >= 0; i--) {

// constants for this level

const Mat &keypoints = keypoints_by_layers[i];

const uchar *const kpts = keypoints_by_layers[i].ptr<uchar>();

uchar *const kpts_next = keypoints_by_layers[i+1].ptr<uchar>();

const float *const ldet = evolution_[i].Ldet.ptr<float>();

const float *const ldet_next = evolution_[i+1].Ldet.ptr<float>();

// ratios are just powers of 2, i+1 ratio is always greater or equal to i

const int diff_ratio = (int)evolution_[i+1].octave_ratio / (int)evolution_[i].octave_ratio;

const int search_radius = evolution_[i+1].sigma_size; // size of keypoints in upper level

size_t j = 0;

for (int y = 0; y < keypoints.rows; y++) {

for (int x = 0; x < keypoints.cols; x++, j++) {

if (kpts[j] == 0) {

continue; // skip non-keypoints

}

int idx = 0;

// project point to upper scale layer

const int p_x = x / diff_ratio;

const int p_y = y / diff_ratio;

if (find_neighbor_point(p_x, p_y, keypoints_by_layers[i+1], search_radius, idx)) {

if (ldet[j] > ldet_next[idx]) {

kpts_next[idx] = 0; // clear keypoint in upper layer

}

}

}

}

}

std::vector<Mat>& keypoints_by_layers, std::vector<KeyPoint>& output_keypoints)

{

CV_INSTRUMENT_REGION()

for (size_t i = 0; i < keypoints_by_layers.size(); i++) {

const MEvolution &e = evolution_[i];

const float * const ldet = e.Ldet.ptr<float>();

const float ratio = e.octave_ratio;

const int cols = e.Ldet.cols;

const Mat& keypoints = keypoints_by_layers[i];

const uchar *const kpts = keypoints.ptr<uchar>();

size_t j = 0;

for (int y = 0; y < keypoints.rows; y++) {

for (int x = 0; x < keypoints.cols; x++, j++) {

if (kpts[j] == 0) {

continue; // skip non-keypoints

}

// create a new keypoint

KeyPoint kp;

kp.pt.x = x * e.octave_ratio;

kp.pt.y = y * e.octave_ratio;

kp.size = e.esigma * options_.derivative_factor;

kp.angle = -1;

kp.response = ldet[j];

kp.octave = e.octave;

kp.class_id = static_cast<int>(i);

// Compute the gradient

float Dx = 0.5f * (ldet[ y *cols + x + 1] - ldet[ y *cols + x - 1]);

float Dy = 0.5f * (ldet[(y + 1)*cols + x ] - ldet[(y - 1)*cols + x ]);

// Compute the Hessian

float Dxx = ldet[ y *cols + x + 1] + ldet[ y *cols + x - 1] - 2.0f * ldet[y*cols + x];

float Dyy = ldet[(y + 1)*cols + x ] + ldet[(y - 1)*cols + x ] - 2.0f * ldet[y*cols + x];

float Dxy = 0.25f * (ldet[(y + 1)*cols + x + 1] + ldet[(y - 1)*cols + x - 1] -

ldet[(y - 1)*cols + x + 1] - ldet[(y + 1)*cols + x - 1]);

// Solve the linear system

Matx22f A( Dxx, Dxy,

Dxy, Dyy );

Vec2f b( -Dx, -Dy );

Vec2f dst( 0.0f, 0.0f );

solve(A, b, dst, DECOMP_LU);

float dx = dst(0);

float dy = dst(1);

if (fabs(dx) > 1.0f || fabs(dy) > 1.0f)

continue; // Ignore the point that is not stable

// Refine the coordinates

kp.pt.x += dx * ratio + .5f*(ratio-1.f);

kp.pt.y += dy * ratio + .5f*(ratio-1.f);

kp.angle = 0.0;

kp.size *= 2.0f; // In OpenCV the size of a keypoint is the diameter

// Push the refined keypoint to the final storage

output_keypoints.push_back(kp);

}

}

}

{

SURF_Descriptor_Upright_64_Invoker(std::vector<KeyPoint>& kpts, Mat& desc, const Pyramid& evolution)

: keypoints_(&kpts)

, descriptors_(&desc)

, evolution_(&evolution)

{

}

void operator() (const Range& range) const

{

for (int i = range.start; i < range.end; i++)

{

Get_SURF_Descriptor_Upright_64((*keypoints_)[i], descriptors_->ptr<float>(i), descriptors_->cols);

}

}

void Get_SURF_Descriptor_Upright_64(const KeyPoint& kpt, float* desc, int desc_size) const;

std::vector<KeyPoint>* keypoints_;

Mat* descriptors_;

const Pyramid* evolution_;

{

SURF_Descriptor_64_Invoker(std::vector<KeyPoint>& kpts, Mat& desc, Pyramid& evolution)

: keypoints_(&kpts)

, descriptors_(&desc)

, evolution_(&evolution)

{

}

void operator()(const Range& range) const

{

for (int i = range.start; i < range.end; i++)

{

Get_SURF_Descriptor_64((*keypoints_)[i], descriptors_->ptr<float>(i), descriptors_->cols);

}

}

void Get_SURF_Descriptor_64(const KeyPoint& kpt, float* desc, int desc_size) const;

std::vector<KeyPoint>* keypoints_;

Mat* descriptors_;

Pyramid* evolution_;

{

MSURF_Upright_Descriptor_64_Invoker(std::vector<KeyPoint>& kpts, Mat& desc, Pyramid& evolution)

: keypoints_(&kpts)

, descriptors_(&desc)

, evolution_(&evolution)

{