int level_height = 0, level_width = 0;
// Allocate the dimension of the matrices for the evolution
- for (int i = 0; i <= options_.omax - 1; i++) {
- rfactor = 1.0f / pow(2.f, i);
+ 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);
for (int i = range.start; i < range.end; i++)
{
- float ratio = pow(2.f, (float)evolution[i].octave);
+ float ratio = (float)fastpow(2, evolution[i].octave);
int sigma_size_ = fRound(evolution[i].esigma * options_.derivative_factor / ratio);
compute_scharr_derivatives(evolution[i].Lsmooth, evolution[i].Lx, 1, 0, sigma_size_);
}
for (size_t i = 0; i < evolution_.size(); i++) {
+ float* prev = evolution_[i].Ldet.ptr<float>(0);
+ float* curr = evolution_[i].Ldet.ptr<float>(1);
for (int ix = 1; ix < evolution_[i].Ldet.rows - 1; ix++) {
+ float* next = evolution_[i].Ldet.ptr<float>(ix + 1);
+
for (int jx = 1; jx < evolution_[i].Ldet.cols - 1; jx++) {
is_extremum = false;
is_repeated = false;
// Filter the points with the detector threshold
if (value > options_.dthreshold && value >= options_.min_dthreshold &&
- value > *(evolution_[i].Ldet.ptr<float>(ix)+jx - 1) &&
- value > *(evolution_[i].Ldet.ptr<float>(ix)+jx + 1) &&
- value > *(evolution_[i].Ldet.ptr<float>(ix - 1) + jx - 1) &&
- value > *(evolution_[i].Ldet.ptr<float>(ix - 1) + jx) &&
- value > *(evolution_[i].Ldet.ptr<float>(ix - 1) + jx + 1) &&
- value > *(evolution_[i].Ldet.ptr<float>(ix + 1) + jx - 1) &&
- value > *(evolution_[i].Ldet.ptr<float>(ix + 1) + jx) &&
- value > *(evolution_[i].Ldet.ptr<float>(ix + 1) + jx + 1)) {
+ value > curr[jx-1] &&
+ value > curr[jx+1] &&
+ value > prev[jx-1] &&
+ value > prev[jx] &&
+ value > prev[jx+1] &&
+ value > next[jx-1] &&
+ value > next[jx] &&
+ value > next[jx+1]) {
is_extremum = true;
point.response = fabs(value);
point.size = evolution_[i].esigma*options_.derivative_factor;
point.octave = (int)evolution_[i].octave;
point.class_id = (int)i;
- ratio = pow(2.f, point.octave);
+ ratio = (float)fastpow(2, point.octave);
sigma_size_ = fRound(point.size / ratio);
point.pt.x = static_cast<float>(jx);
point.pt.y = static_cast<float>(ix);
if ((point.class_id - 1) == kpts_aux[ik].class_id ||
point.class_id == kpts_aux[ik].class_id) {
- dist = sqrt(pow(point.pt.x*ratio - kpts_aux[ik].pt.x, 2) + pow(point.pt.y*ratio - kpts_aux[ik].pt.y, 2));
- if (dist <= point.size) {
+ float distx = point.pt.x*ratio - kpts_aux[ik].pt.x;
+ float disty = point.pt.y*ratio - kpts_aux[ik].pt.y;
+ dist = distx * distx + disty * disty;
+ if (dist <= point.size * point.size) {
if (point.response > kpts_aux[ik].response) {
id_repeated = (int)ik;
is_repeated = true;
} //if is_extremum
}
} // for jx
+ prev = curr;
+ curr = next;
} // for ix
} // for i
// Compare response with the upper scale
if ((pt.class_id + 1) == kpts_aux[j].class_id) {
- dist = sqrt(pow(pt.pt.x - kpts_aux[j].pt.x, 2) + pow(pt.pt.y - kpts_aux[j].pt.y, 2));
- if (dist <= pt.size) {
+ float distx = pt.pt.x - kpts_aux[j].pt.x;
+ float disty = pt.pt.y - kpts_aux[j].pt.y;
+ dist = distx * distx + disty * disty;
+ if (dist <= pt.size * pt.size) {
if (pt.response < kpts_aux[j].response) {
is_repeated = true;
break;
float Dx = 0.0, Dy = 0.0, ratio = 0.0;
float Dxx = 0.0, Dyy = 0.0, Dxy = 0.0;
int x = 0, y = 0;
- cv::Mat A = cv::Mat::zeros(2, 2, CV_32F);
- cv::Mat b = cv::Mat::zeros(2, 1, CV_32F);
- cv::Mat dst = cv::Mat::zeros(2, 1, CV_32F);
+ Matx22f A(0, 0, 0, 0);
+ Vec2f b(0, 0);
+ Vec2f dst(0, 0);
for (size_t i = 0; i < kpts.size(); i++) {
- ratio = pow(2.f, kpts[i].octave);
+ ratio = (float)fastpow(2, kpts[i].octave);
x = fRound(kpts[i].pt.x / ratio);
y = fRound(kpts[i].pt.y / ratio);
+ (*(evolution_[kpts[i].class_id].Ldet.ptr<float>(y + 1) + x - 1)));
// Solve the linear system
- *(A.ptr<float>(0)) = Dxx;
- *(A.ptr<float>(1) + 1) = Dyy;
- *(A.ptr<float>(0) + 1) = *(A.ptr<float>(1)) = Dxy;
- *(b.ptr<float>(0)) = -Dx;
- *(b.ptr<float>(1)) = -Dy;
+ A(0, 0) = Dxx;
+ A(1, 1) = Dyy;
+ A(0, 1) = A(1, 0) = Dxy;
+ b(0) = -Dx;
+ b(1) = -Dy;
cv::solve(A, b, dst, DECOMP_LU);
- if (fabs(*(dst.ptr<float>(0))) <= 1.0f && fabs(*(dst.ptr<float>(1))) <= 1.0f) {
- kpts[i].pt.x = x + (*(dst.ptr<float>(0)));
- kpts[i].pt.y = y + (*(dst.ptr<float>(1)));
- kpts[i].pt.x *= powf(2.f, (float)evolution_[kpts[i].class_id].octave);
- kpts[i].pt.y *= powf(2.f, (float)evolution_[kpts[i].class_id].octave);
+ if (fabs(dst(0)) <= 1.0f && fabs(dst(1)) <= 1.0f) {
+ kpts[i].pt.x = x + dst(0);
+ kpts[i].pt.y = y + dst(1);
+ int power = fastpow(2, evolution_[kpts[i].class_id].octave);
+ kpts[i].pt.x *= power;
+ kpts[i].pt.y *= power;
kpts[i].angle = 0.0;
// In OpenCV the size of a keypoint its the diameter
}
void Get_MLDB_Full_Descriptor(const cv::KeyPoint& kpt, unsigned char* desc) const;
+ void MLDB_Fill_Values(float* values, int sample_step, int level,
+ float xf, float yf, float co, float si, float scale) const;
+ void MLDB_Binary_Comparisons(float* values, unsigned char* desc,
+ int count, int& dpos) const;
private:
std::vector<cv::KeyPoint>* keypoints_;
int ix = 0, iy = 0, idx = 0, s = 0, level = 0;
float xf = 0.0, yf = 0.0, gweight = 0.0, ratio = 0.0;
- std::vector<float> resX(109), resY(109), Ang(109);
+ const int ang_size = 109;
+ float resX[ang_size], resY[ang_size], Ang[ang_size];
const int id[] = { 6, 5, 4, 3, 2, 1, 0, 1, 2, 3, 4, 5, 6 };
// Variables for computing the dominant direction
resX[idx] = gweight*(*(evolution_[level].Lx.ptr<float>(iy)+ix));
resY[idx] = gweight*(*(evolution_[level].Ly.ptr<float>(iy)+ix));
- Ang[idx] = getAngle(resX[idx], resY[idx]);
++idx;
}
}
}
+ fastAtan2(resY, resX, Ang, ang_size, false);
// Loop slides pi/3 window around feature point
for (ang1 = 0; ang1 < (float)(2.0 * CV_PI); ang1 += 0.15f) {
ang2 = (ang1 + (float)(CV_PI / 3.0) >(float)(2.0*CV_PI) ? ang1 - (float)(5.0*CV_PI / 3.0) : ang1 + (float)(CV_PI / 3.0));
sumX = sumY = 0.f;
- for (size_t k = 0; k < Ang.size(); ++k) {
+ for (size_t k = 0; k < ang_size; ++k) {
// Get angle from the x-axis of the sample point
const float & ang = Ang[k];
}
}
-/* ************************************************************************* */
-/**
- * @brief This method computes the descriptor of the provided keypoint given the
- * main orientation of the keypoint
- * @param kpt Input keypoint
- * @param desc Descriptor vector
- */
-void MLDB_Full_Descriptor_Invoker::Get_MLDB_Full_Descriptor(const cv::KeyPoint& kpt, unsigned char *desc) const {
-
- float di = 0.0, dx = 0.0, dy = 0.0, ratio = 0.0;
- float ri = 0.0, rx = 0.0, ry = 0.0, rrx = 0.0, rry = 0.0, xf = 0.0, yf = 0.0;
- float sample_x = 0.0, sample_y = 0.0, co = 0.0, si = 0.0, angle = 0.0;
- int x1 = 0, y1 = 0, sample_step = 0, pattern_size = 0;
- int level = 0, nsamples = 0, scale = 0;
- int dcount1 = 0, dcount2 = 0;
-
- const AKAZEOptions & options = *options_;
- const std::vector<TEvolution>& evolution = *evolution_;
-
- // Matrices for the M-LDB descriptor
- cv::Mat values_1 = cv::Mat::zeros(4, options.descriptor_channels, CV_32FC1);
- cv::Mat values_2 = cv::Mat::zeros(9, options.descriptor_channels, CV_32FC1);
- cv::Mat values_3 = cv::Mat::zeros(16, options.descriptor_channels, CV_32FC1);
-
- // Get the information from the keypoint
- ratio = (float)(1 << kpt.octave);
- scale = fRound(0.5f*kpt.size / ratio);
- angle = kpt.angle;
- level = kpt.class_id;
- yf = kpt.pt.y / ratio;
- xf = kpt.pt.x / ratio;
- co = cos(angle);
- si = sin(angle);
-
- // First 2x2 grid
- pattern_size = options.descriptor_pattern_size;
- sample_step = pattern_size;
-
- for (int i = -pattern_size; i < pattern_size; i += sample_step) {
- for (int j = -pattern_size; j < pattern_size; j += sample_step) {
-
- di = dx = dy = 0.0;
- nsamples = 0;
-
- for (float k = (float)i; k < i + sample_step; k++) {
- for (float l = (float)j; l < j + sample_step; l++) {
-
- // Get the coordinates of the sample point
- sample_y = yf + (l*scale*co + k*scale*si);
- sample_x = xf + (-l*scale*si + k*scale*co);
-
- y1 = fRound(sample_y);
- x1 = fRound(sample_x);
-
- ri = *(evolution[level].Lt.ptr<float>(y1)+x1);
- rx = *(evolution[level].Lx.ptr<float>(y1)+x1);
- ry = *(evolution[level].Ly.ptr<float>(y1)+x1);
-
- di += ri;
-
- if (options.descriptor_channels == 2) {
- dx += sqrtf(rx*rx + ry*ry);
- }
- else if (options.descriptor_channels == 3) {
- // Get the x and y derivatives on the rotated axis
- rry = rx*co + ry*si;
- rrx = -rx*si + ry*co;
- dx += rrx;
- dy += rry;
- }
-
- nsamples++;
- }
- }
-
- di /= nsamples;
- dx /= nsamples;
- dy /= nsamples;
-
- *(values_1.ptr<float>(dcount2)) = di;
- if (options.descriptor_channels > 1) {
- *(values_1.ptr<float>(dcount2)+1) = dx;
- }
-
- if (options.descriptor_channels > 2) {
- *(values_1.ptr<float>(dcount2)+2) = dy;
- }
-
- dcount2++;
- }
- }
-
- // Do binary comparison first level
- for (int i = 0; i < 4; i++) {
- for (int j = i + 1; j < 4; j++) {
- if (*(values_1.ptr<float>(i)) > *(values_1.ptr<float>(j))) {
- desc[dcount1 / 8] |= (1 << (dcount1 % 8));
- }
- dcount1++;
- }
- }
-
- if (options.descriptor_channels > 1) {
- for (int i = 0; i < 4; i++) {
- for (int j = i + 1; j < 4; j++) {
- if (*(values_1.ptr<float>(i)+1) > *(values_1.ptr<float>(j)+1)) {
- desc[dcount1 / 8] |= (1 << (dcount1 % 8));
- }
-
- dcount1++;
- }
- }
- }
-
- if (options.descriptor_channels > 2) {
- for (int i = 0; i < 4; i++) {
- for (int j = i + 1; j < 4; j++) {
- if (*(values_1.ptr<float>(i)+2) > *(values_1.ptr<float>(j)+2)) {
- desc[dcount1 / 8] |= (1 << (dcount1 % 8));
- }
- dcount1++;
- }
- }
- }
-
- // Second 3x3 grid
- sample_step = static_cast<int>(ceil(pattern_size*2. / 3.));
- dcount2 = 0;
-
- for (int i = -pattern_size; i < pattern_size; i += sample_step) {
- for (int j = -pattern_size; j < pattern_size; j += sample_step) {
-
- di = dx = dy = 0.0;
- nsamples = 0;
-
- for (int k = i; k < i + sample_step; k++) {
- for (int l = j; l < j + sample_step; l++) {
-
- // Get the coordinates of the sample point
- sample_y = yf + (l*scale*co + k*scale*si);
- sample_x = xf + (-l*scale*si + k*scale*co);
-
- y1 = fRound(sample_y);
- x1 = fRound(sample_x);
-
- ri = *(evolution[level].Lt.ptr<float>(y1)+x1);
- rx = *(evolution[level].Lx.ptr<float>(y1)+x1);
- ry = *(evolution[level].Ly.ptr<float>(y1)+x1);
- di += ri;
+void MLDB_Full_Descriptor_Invoker::MLDB_Fill_Values(float* values, int sample_step, int level,
+ float xf, float yf, float co, float si, float scale) const
+{
+ const std::vector<TEvolution>& evolution = *evolution_;
+ int pattern_size = options_->descriptor_pattern_size;
+ int chan = options_->descriptor_channels;
+ int valpos = 0;
+
+ for (int i = -pattern_size; i < pattern_size; i += sample_step) {
+ for (int j = -pattern_size; j < pattern_size; j += sample_step) {
+ float di, dx, dy;
+ di = dx = dy = 0.0;
+ int nsamples = 0;
+
+ for (int k = i; k < i + sample_step; k++) {
+ for (int l = j; l < j + sample_step; l++) {
+ float sample_y = yf + (l*co * scale + k*si*scale);
+ float sample_x = xf + (-l*si * scale + k*co*scale);
+
+ int y1 = fRound(sample_y);
+ int x1 = fRound(sample_x);
+
+ float ri = *(evolution[level].Lt.ptr<float>(y1)+x1);
+ di += ri;
+
+ if(chan > 1) {
+ float rx = *(evolution[level].Lx.ptr<float>(y1)+x1);
+ float ry = *(evolution[level].Ly.ptr<float>(y1)+x1);
+ if (chan == 2) {
+ dx += sqrtf(rx*rx + ry*ry);
+ }
+ else {
+ float rry = rx*co + ry*si;
+ float rrx = -rx*si + ry*co;
+ dx += rrx;
+ dy += rry;
+ }
+ }
+ nsamples++;
+ }
+ }
+ di /= nsamples;
+ dx /= nsamples;
+ dy /= nsamples;
- if (options.descriptor_channels == 2) {
- dx += sqrtf(rx*rx + ry*ry);
- }
- else if (options.descriptor_channels == 3) {
- // Get the x and y derivatives on the rotated axis
- rry = rx*co + ry*si;
- rrx = -rx*si + ry*co;
- dx += rrx;
- dy += rry;
+ values[valpos] = di;
+ if (chan > 1) {
+ values[valpos + 1] = dx;
+ }
+ if (chan > 2) {
+ values[valpos + 2] = dy;
+ }
+ valpos += chan;
}
-
- nsamples++;
}
- }
-
- di /= nsamples;
- dx /= nsamples;
- dy /= nsamples;
-
- *(values_2.ptr<float>(dcount2)) = di;
- if (options.descriptor_channels > 1) {
- *(values_2.ptr<float>(dcount2)+1) = dx;
- }
-
- if (options.descriptor_channels > 2) {
- *(values_2.ptr<float>(dcount2)+2) = dy;
- }
-
- dcount2++;
- }
- }
-
- // Do binary comparison second level
- for (int i = 0; i < 9; i++) {
- for (int j = i + 1; j < 9; j++) {
- if (*(values_2.ptr<float>(i)) > *(values_2.ptr<float>(j))) {
- desc[dcount1 / 8] |= (1 << (dcount1 % 8));
- }
- dcount1++;
- }
- }
+}
- if (options.descriptor_channels > 1) {
- for (int i = 0; i < 9; i++) {
- for (int j = i + 1; j < 9; j++) {
- if (*(values_2.ptr<float>(i)+1) > *(values_2.ptr<float>(j)+1)) {
- desc[dcount1 / 8] |= (1 << (dcount1 % 8));
- }
- dcount1++;
- }
+void MLDB_Full_Descriptor_Invoker::MLDB_Binary_Comparisons(float* values, unsigned char* desc,
+ int count, int& dpos) const {
+ int chan = options_->descriptor_channels;
+ int* ivalues = (int*) values;
+ for(int i = 0; i < count * chan; i++) {
+ ivalues[i] = CV_TOGGLE_FLT(ivalues[i]);
}
- }
- if (options.descriptor_channels > 2) {
- for (int i = 0; i < 9; i++) {
- for (int j = i + 1; j < 9; j++) {
- if (*(values_2.ptr<float>(i)+2) > *(values_2.ptr<float>(j)+2)) {
- desc[dcount1 / 8] |= (1 << (dcount1 % 8));
+ for(int pos = 0; pos < chan; pos++) {
+ for (int i = 0; i < count; i++) {
+ int ival = ivalues[chan * i + pos];
+ for (int j = i + 1; j < count; j++) {
+ int res = ival > ivalues[chan * j + pos];
+ desc[dpos >> 3] |= (res << (dpos & 7));
+ dpos++;
+ }
}
- dcount1++;
- }
}
- }
-
- // Third 4x4 grid
- sample_step = pattern_size / 2;
- dcount2 = 0;
-
- for (int i = -pattern_size; i < pattern_size; i += sample_step) {
- for (int j = -pattern_size; j < pattern_size; j += sample_step) {
- di = dx = dy = 0.0;
- nsamples = 0;
-
- for (int k = i; k < i + sample_step; k++) {
- for (int l = j; l < j + sample_step; l++) {
-
- // Get the coordinates of the sample point
- sample_y = yf + (l*scale*co + k*scale*si);
- sample_x = xf + (-l*scale*si + k*scale*co);
-
- y1 = fRound(sample_y);
- x1 = fRound(sample_x);
-
- ri = *(evolution[level].Lt.ptr<float>(y1)+x1);
- rx = *(evolution[level].Lx.ptr<float>(y1)+x1);
- ry = *(evolution[level].Ly.ptr<float>(y1)+x1);
- di += ri;
-
- if (options.descriptor_channels == 2) {
- dx += sqrtf(rx*rx + ry*ry);
- }
- else if (options.descriptor_channels == 3) {
- // Get the x and y derivatives on the rotated axis
- rry = rx*co + ry*si;
- rrx = -rx*si + ry*co;
- dx += rrx;
- dy += rry;
- }
-
- nsamples++;
- }
- }
-
- di /= nsamples;
- dx /= nsamples;
- dy /= nsamples;
-
- *(values_3.ptr<float>(dcount2)) = di;
- if (options.descriptor_channels > 1)
- *(values_3.ptr<float>(dcount2)+1) = dx;
+}
- if (options.descriptor_channels > 2)
- *(values_3.ptr<float>(dcount2)+2) = dy;
+/* ************************************************************************* */
+/**
+ * @brief This method computes the descriptor of the provided keypoint given the
+ * main orientation of the keypoint
+ * @param kpt Input keypoint
+ * @param desc Descriptor vector
+ */
+void MLDB_Full_Descriptor_Invoker::Get_MLDB_Full_Descriptor(const cv::KeyPoint& kpt, unsigned char *desc) const {
- dcount2++;
- }
- }
+ const int max_channels = 3;
+ CV_Assert(options_->descriptor_channels <= max_channels);
+ float values[16*max_channels];
+ const double size_mult[3] = {1, 2.0/3.0, 1.0/2.0};
- // Do binary comparison third level
- for (int i = 0; i < 16; i++) {
- for (int j = i + 1; j < 16; j++) {
- if (*(values_3.ptr<float>(i)) > *(values_3.ptr<float>(j))) {
- desc[dcount1 / 8] |= (1 << (dcount1 % 8));
- }
- dcount1++;
- }
- }
+ float ratio = (float)(1 << kpt.octave);
+ float scale = (float)fRound(0.5f*kpt.size / ratio);
+ float xf = kpt.pt.x / ratio;
+ float yf = kpt.pt.y / ratio;
+ float co = cos(kpt.angle);
+ float si = sin(kpt.angle);
+ int pattern_size = options_->descriptor_pattern_size;
- if (options.descriptor_channels > 1) {
- for (int i = 0; i < 16; i++) {
- for (int j = i + 1; j < 16; j++) {
- if (*(values_3.ptr<float>(i)+1) > *(values_3.ptr<float>(j)+1)) {
- desc[dcount1 / 8] |= (1 << (dcount1 % 8));
- }
- dcount1++;
- }
- }
- }
+ int dpos = 0;
+ for(int lvl = 0; lvl < 3; lvl++) {
- if (options.descriptor_channels > 2) {
- for (int i = 0; i < 16; i++) {
- for (int j = i + 1; j < 16; j++) {
- if (*(values_3.ptr<float>(i)+2) > *(values_3.ptr<float>(j)+2)) {
- desc[dcount1 / 8] |= (1 << (dcount1 % 8));
- }
- dcount1++;
- }
- }
+ int val_count = (lvl + 2) * (lvl + 2);
+ int sample_step = static_cast<int>(ceil(pattern_size * size_mult[lvl]));
+ MLDB_Fill_Values(values, sample_step, kpt.class_id, xf, yf, co, si, scale);
+ MLDB_Binary_Comparisons(values, desc, val_count, dpos);
}
}
*/
#include "nldiffusion_functions.h"
+#include <iostream>
// Namespaces
using namespace std;
* @param k Contrast factor parameter
*/
void pm_g1(const cv::Mat& Lx, const cv::Mat& Ly, cv::Mat& dst, float k) {
- cv::exp(-(Lx.mul(Lx) + Ly.mul(Ly)) / (k*k), dst);
+
+ Size sz = Lx.size();
+ float inv_k = 1.0f / (k*k);
+ for (int y = 0; y < sz.height; y++) {
+
+ const float* Lx_row = Lx.ptr<float>(y);
+ const float* Ly_row = Ly.ptr<float>(y);
+ float* dst_row = dst.ptr<float>(y);
+
+ for (int x = 0; x < sz.width; x++) {
+ dst_row[x] = (-inv_k*(Lx_row[x]*Lx_row[x] + Ly_row[x]*Ly_row[x]));
+ }
+ }
+
+ exp(dst, dst);
}
/* ************************************************************************* */
* @param k Contrast factor parameter
*/
void pm_g2(const cv::Mat &Lx, const cv::Mat& Ly, cv::Mat& dst, float k) {
- dst = 1.0f / (1.0f + (Lx.mul(Lx) + Ly.mul(Ly)) / (k*k));
- }
+ Size sz = Lx.size();
+ dst.create(sz, Lx.type());
+ float k2inv = 1.0f / (k * k);
+
+ for(int y = 0; y < sz.height; y++) {
+ const float *Lx_row = Lx.ptr<float>(y);
+ const float *Ly_row = Ly.ptr<float>(y);
+ float* dst_row = dst.ptr<float>(y);
+ for(int x = 0; x < sz.width; x++) {
+ dst_row[x] = 1.0f / (1.0f + ((Lx_row[x] * Lx_row[x] + Ly_row[x] * Ly_row[x]) * k2inv));
+ }
+ }
+ }
/* ************************************************************************* */
/**
* @brief This function computes Weickert conductivity coefficient gw
* Proceedings of Algorithmy 2000
*/
void weickert_diffusivity(const cv::Mat& Lx, const cv::Mat& Ly, cv::Mat& dst, float k) {
- Mat modg;
- cv::pow((Lx.mul(Lx) + Ly.mul(Ly)) / (k*k), 4, modg);
- cv::exp(-3.315f / modg, dst);
- dst = 1.0f - dst;
+
+ Size sz = Lx.size();
+ float inv_k = 1.0f / (k*k);
+ for (int y = 0; y < sz.height; y++) {
+
+ const float* Lx_row = Lx.ptr<float>(y);
+ const float* Ly_row = Ly.ptr<float>(y);
+ float* dst_row = dst.ptr<float>(y);
+
+ for (int x = 0; x < sz.width; x++) {
+ float dL = inv_k*(Lx_row[x]*Lx_row[x] + Ly_row[x]*Ly_row[x]);
+ dst_row[x] = -3.315f/(dL*dL*dL*dL);
+ }
+ }
+
+ exp(dst, dst);
+ dst = 1.0 - dst;
}
+
/* ************************************************************************* */
/**
* @brief This function computes Charbonnier conductivity coefficient gc
* Proceedings of Algorithmy 2000
*/
void charbonnier_diffusivity(const cv::Mat& Lx, const cv::Mat& Ly, cv::Mat& dst, float k) {
- Mat den;
- cv::sqrt(1.0f + (Lx.mul(Lx) + Ly.mul(Ly)) / (k*k), den);
- dst = 1.0f / den;
+
+ Size sz = Lx.size();
+ float inv_k = 1.0f / (k*k);
+ for (int y = 0; y < sz.height; y++) {
+
+ const float* Lx_row = Lx.ptr<float>(y);
+ const float* Ly_row = Ly.ptr<float>(y);
+ float* dst_row = dst.ptr<float>(y);
+
+ for (int x = 0; x < sz.width; x++) {
+ float den = sqrt(1.0f+inv_k*(Lx_row[x]*Lx_row[x] + Ly_row[x]*Ly_row[x]));
+ dst_row[x] = 1.0f / den;
+ }
+ }
}
+
/* ************************************************************************* */
/**
* @brief This function computes a good empirical value for the k contrast factor
float compute_k_percentile(const cv::Mat& img, float perc, float gscale, int nbins, int ksize_x, int ksize_y) {
int nbin = 0, nelements = 0, nthreshold = 0, k = 0;
- float kperc = 0.0, modg = 0.0, lx = 0.0, ly = 0.0;
+ float kperc = 0.0, modg = 0.0;
float npoints = 0.0;
float hmax = 0.0;
// Skip the borders for computing the histogram
for (int i = 1; i < gaussian.rows - 1; i++) {
+ const float *lx = Lx.ptr<float>(i);
+ const float *ly = Ly.ptr<float>(i);
for (int j = 1; j < gaussian.cols - 1; j++) {
- lx = *(Lx.ptr<float>(i)+j);
- ly = *(Ly.ptr<float>(i)+j);
- modg = sqrt(lx*lx + ly*ly);
+ modg = lx[j]*lx[j] + ly[j]*ly[j];
// Get the maximum
if (modg > hmax) {
}
}
}
-
+ hmax = sqrt(hmax);
// Skip the borders for computing the histogram
for (int i = 1; i < gaussian.rows - 1; i++) {
+ const float *lx = Lx.ptr<float>(i);
+ const float *ly = Ly.ptr<float>(i);
for (int j = 1; j < gaussian.cols - 1; j++) {
- lx = *(Lx.ptr<float>(i)+j);
- ly = *(Ly.ptr<float>(i)+j);
- modg = sqrt(lx*lx + ly*ly);
+ modg = lx[j]*lx[j] + ly[j]*ly[j];
// Find the correspondent bin
if (modg != 0.0) {
- nbin = (int)floor(nbins*(modg / hmax));
+ nbin = (int)floor(nbins*(sqrt(modg) / hmax));
if (nbin == nbins) {
nbin--;
/* ************************************************************************* */
/**
* @brief Compute derivative kernels for sizes different than 3
- * @param _kx Horizontal kernel values
+ * @param _kx Horizontal kernel ues
* @param _ky Vertical kernel values
* @param dx Derivative order in X-direction (horizontal)
* @param dy Derivative order in Y-direction (vertical)
for (int i = range.start; i < range.end; i++)
{
+ const float *c_prev = c.ptr<float>(i - 1);
+ const float *c_curr = c.ptr<float>(i);
+ const float *c_next = c.ptr<float>(i + 1);
+ const float *ld_prev = Ld.ptr<float>(i - 1);
+ const float *ld_curr = Ld.ptr<float>(i);
+ const float *ld_next = Ld.ptr<float>(i + 1);
+
+ float *dst = Lstep.ptr<float>(i);
+
for (int j = 1; j < Lstep.cols - 1; j++)
{
- float xpos = ((*(c.ptr<float>(i)+j)) + (*(c.ptr<float>(i)+j + 1)))*((*(Ld.ptr<float>(i)+j + 1)) - (*(Ld.ptr<float>(i)+j)));
- float xneg = ((*(c.ptr<float>(i)+j - 1)) + (*(c.ptr<float>(i)+j)))*((*(Ld.ptr<float>(i)+j)) - (*(Ld.ptr<float>(i)+j - 1)));
- float ypos = ((*(c.ptr<float>(i)+j)) + (*(c.ptr<float>(i + 1) + j)))*((*(Ld.ptr<float>(i + 1) + j)) - (*(Ld.ptr<float>(i)+j)));
- float yneg = ((*(c.ptr<float>(i - 1) + j)) + (*(c.ptr<float>(i)+j)))*((*(Ld.ptr<float>(i)+j)) - (*(Ld.ptr<float>(i - 1) + j)));
- *(Lstep.ptr<float>(i)+j) = 0.5f*stepsize*(xpos - xneg + ypos - yneg);
+ float xpos = (c_curr[j] + c_curr[j+1])*(ld_curr[j+1] - ld_curr[j]);
+ float xneg = (c_curr[j-1] + c_curr[j]) *(ld_curr[j] - ld_curr[j-1]);
+ float ypos = (c_curr[j] + c_next[j]) *(ld_next[j] - ld_curr[j]);
+ float yneg = (c_prev[j] + c_curr[j]) *(ld_curr[j] - ld_prev[j]);
+ dst[j] = 0.5f*stepsize*(xpos - xneg + ypos - yneg);
}
}
}
-
private:
cv::Mat * _Ld;
const cv::Mat * _c;
*/
void nld_step_scalar(cv::Mat& Ld, const cv::Mat& c, cv::Mat& Lstep, float stepsize) {
- cv::parallel_for_(cv::Range(1, Lstep.rows - 1), Nld_Step_Scalar_Invoker(Ld, c, Lstep, stepsize));
+ cv::parallel_for_(cv::Range(1, Lstep.rows - 1), Nld_Step_Scalar_Invoker(Ld, c, Lstep, stepsize), (double)Ld.total()/(1 << 16));
+ float xneg, xpos, yneg, ypos;
+ float* dst = Lstep.ptr<float>(0);
+ const float* cprv = NULL;
+ const float* ccur = c.ptr<float>(0);
+ const float* cnxt = c.ptr<float>(1);
+ const float* ldprv = NULL;
+ const float* ldcur = Ld.ptr<float>(0);
+ const float* ldnxt = Ld.ptr<float>(1);
for (int j = 1; j < Lstep.cols - 1; j++) {
- float xpos = ((*(c.ptr<float>(0) + j)) + (*(c.ptr<float>(0) + j + 1)))*((*(Ld.ptr<float>(0) + j + 1)) - (*(Ld.ptr<float>(0) + j)));
- float xneg = ((*(c.ptr<float>(0) + j - 1)) + (*(c.ptr<float>(0) + j)))*((*(Ld.ptr<float>(0) + j)) - (*(Ld.ptr<float>(0) + j - 1)));
- float ypos = ((*(c.ptr<float>(0) + j)) + (*(c.ptr<float>(1) + j)))*((*(Ld.ptr<float>(1) + j)) - (*(Ld.ptr<float>(0) + j)));
- *(Lstep.ptr<float>(0) + j) = 0.5f*stepsize*(xpos - xneg + ypos);
+ xpos = (ccur[j] + ccur[j+1]) * (ldcur[j+1] - ldcur[j]);
+ xneg = (ccur[j-1] + ccur[j]) * (ldcur[j] - ldcur[j-1]);
+ ypos = (ccur[j] + cnxt[j]) * (ldnxt[j] - ldcur[j]);
+ dst[j] = 0.5f*stepsize*(xpos - xneg + ypos);
}
+ dst = Lstep.ptr<float>(Lstep.rows - 1);
+ ccur = c.ptr<float>(Lstep.rows - 1);
+ cprv = c.ptr<float>(Lstep.rows - 2);
+ ldcur = Ld.ptr<float>(Lstep.rows - 1);
+ ldprv = Ld.ptr<float>(Lstep.rows - 2);
+
for (int j = 1; j < Lstep.cols - 1; j++) {
- float xpos = ((*(c.ptr<float>(Lstep.rows - 1) + j)) + (*(c.ptr<float>(Lstep.rows - 1) + j + 1)))*((*(Ld.ptr<float>(Lstep.rows - 1) + j + 1)) - (*(Ld.ptr<float>(Lstep.rows - 1) + j)));
- float xneg = ((*(c.ptr<float>(Lstep.rows - 1) + j - 1)) + (*(c.ptr<float>(Lstep.rows - 1) + j)))*((*(Ld.ptr<float>(Lstep.rows - 1) + j)) - (*(Ld.ptr<float>(Lstep.rows - 1) + j - 1)));
- float ypos = ((*(c.ptr<float>(Lstep.rows - 1) + j)) + (*(c.ptr<float>(Lstep.rows - 1) + j)))*((*(Ld.ptr<float>(Lstep.rows - 1) + j)) - (*(Ld.ptr<float>(Lstep.rows - 1) + j)));
- float yneg = ((*(c.ptr<float>(Lstep.rows - 2) + j)) + (*(c.ptr<float>(Lstep.rows - 1) + j)))*((*(Ld.ptr<float>(Lstep.rows - 1) + j)) - (*(Ld.ptr<float>(Lstep.rows - 2) + j)));
- *(Lstep.ptr<float>(Lstep.rows - 1) + j) = 0.5f*stepsize*(xpos - xneg + ypos - yneg);
+ xpos = (ccur[j] + ccur[j+1]) * (ldcur[j+1] - ldcur[j]);
+ xneg = (ccur[j-1] + ccur[j]) * (ldcur[j] - ldcur[j-1]);
+ yneg = (cprv[j] + ccur[j]) * (ldcur[j] - ldprv[j]);
+ dst[j] = 0.5f*stepsize*(xpos - xneg - yneg);
}
- for (int i = 1; i < Lstep.rows - 1; i++) {
- float xpos = ((*(c.ptr<float>(i))) + (*(c.ptr<float>(i)+1)))*((*(Ld.ptr<float>(i)+1)) - (*(Ld.ptr<float>(i))));
- float xneg = ((*(c.ptr<float>(i))) + (*(c.ptr<float>(i))))*((*(Ld.ptr<float>(i))) - (*(Ld.ptr<float>(i))));
- float ypos = ((*(c.ptr<float>(i))) + (*(c.ptr<float>(i + 1))))*((*(Ld.ptr<float>(i + 1))) - (*(Ld.ptr<float>(i))));
- float yneg = ((*(c.ptr<float>(i - 1))) + (*(c.ptr<float>(i))))*((*(Ld.ptr<float>(i))) - (*(Ld.ptr<float>(i - 1))));
- *(Lstep.ptr<float>(i)) = 0.5f*stepsize*(xpos - xneg + ypos - yneg);
- }
+ ccur = c.ptr<float>(1);
+ ldcur = Ld.ptr<float>(1);
+ cprv = c.ptr<float>(0);
+ ldprv = Ld.ptr<float>(0);
+
+ int r0 = Lstep.cols - 1;
+ int r1 = Lstep.cols - 2;
for (int i = 1; i < Lstep.rows - 1; i++) {
- float xneg = ((*(c.ptr<float>(i)+Lstep.cols - 2)) + (*(c.ptr<float>(i)+Lstep.cols - 1)))*((*(Ld.ptr<float>(i)+Lstep.cols - 1)) - (*(Ld.ptr<float>(i)+Lstep.cols - 2)));
- float ypos = ((*(c.ptr<float>(i)+Lstep.cols - 1)) + (*(c.ptr<float>(i + 1) + Lstep.cols - 1)))*((*(Ld.ptr<float>(i + 1) + Lstep.cols - 1)) - (*(Ld.ptr<float>(i)+Lstep.cols - 1)));
- float yneg = ((*(c.ptr<float>(i - 1) + Lstep.cols - 1)) + (*(c.ptr<float>(i)+Lstep.cols - 1)))*((*(Ld.ptr<float>(i)+Lstep.cols - 1)) - (*(Ld.ptr<float>(i - 1) + Lstep.cols - 1)));
- *(Lstep.ptr<float>(i)+Lstep.cols - 1) = 0.5f*stepsize*(-xneg + ypos - yneg);
+ cnxt = c.ptr<float>(i + 1);
+ ldnxt = Ld.ptr<float>(i + 1);
+ dst = Lstep.ptr<float>(i);
+
+ xpos = (ccur[0] + ccur[1]) * (ldcur[1] - ldcur[0]);
+ ypos = (ccur[0] + cnxt[0]) * (ldnxt[0] - ldcur[0]);
+ yneg = (cprv[0] + ccur[0]) * (ldcur[0] - ldprv[0]);
+ dst[0] = 0.5f*stepsize*(xpos + ypos - yneg);
+
+ xneg = (ccur[r1] + ccur[r0]) * (ldcur[r0] - ldcur[r1]);
+ ypos = (ccur[r0] + cnxt[r0]) * (ldnxt[r0] - ldcur[r0]);
+ yneg = (cprv[r0] + ccur[r0]) * (ldcur[r0] - ldprv[r0]);
+ dst[r0] = 0.5f*stepsize*(-xneg + ypos - yneg);
+
+ cprv = ccur;
+ ccur = cnxt;
+ ldprv = ldcur;
+ ldcur = ldnxt;
}
-
- Ld = Ld + Lstep;
+ Ld += Lstep;
}
/* ************************************************************************* */
}
}
}
-}
\ No newline at end of file
+}
projects / profile / ivi / opencv.git / commitdiff