diff --git a/modules/imgproc/include/opencv2/imgproc/detail/legacy.hpp b/modules/imgproc/include/opencv2/imgproc/detail/legacy.hpp index f9abd8d23e..029d9c90e8 100644 --- a/modules/imgproc/include/opencv2/imgproc/detail/legacy.hpp +++ b/modules/imgproc/include/opencv2/imgproc/detail/legacy.hpp @@ -23,6 +23,14 @@ CV_EXPORTS void findContours_legacy(InputArray image, int method, Point offset = Point()); +CV_EXPORTS float EMD_legacy( InputArray _signature1, InputArray _signature2, + int distType, InputArray _cost, + float* lowerBound, OutputArray _flow ); + +CV_EXPORTS float wrapperEMD_legacy(InputArray _signature1, InputArray _signature2, + int distType, InputArray _cost, + Ptr lowerBound, OutputArray _flow); + #endif } // namespace cv diff --git a/modules/imgproc/src/emd.cpp b/modules/imgproc/src/emd.cpp index 8b7f325e74..c4cdd51a41 100644 --- a/modules/imgproc/src/emd.cpp +++ b/modules/imgproc/src/emd.cpp @@ -57,6 +57,7 @@ ========================================================================== */ #include "precomp.hpp" +#include "opencv2/imgproc/detail/legacy.hpp" #define MAX_ITERATIONS 500 #define CV_EMD_INF ((float)1e20) @@ -1147,7 +1148,7 @@ icvDistC( const float *x, const float *y, void *user_param ) } -float cv::EMD( InputArray _signature1, InputArray _signature2, +float cv::EMD_legacy( InputArray _signature1, InputArray _signature2, int distType, InputArray _cost, float* lowerBound, OutputArray _flow ) { @@ -1171,7 +1172,7 @@ float cv::EMD( InputArray _signature1, InputArray _signature2, _flow.needed() ? &_cflow : 0, lowerBound, 0 ); } -float cv::wrapperEMD(InputArray _signature1, InputArray _signature2, +float cv::wrapperEMD_legacy(InputArray _signature1, InputArray _signature2, int distType, InputArray _cost, Ptr lowerBound, OutputArray _flow) { diff --git a/modules/imgproc/src/emd_new.cpp b/modules/imgproc/src/emd_new.cpp new file mode 100644 index 0000000000..8e901c795a --- /dev/null +++ b/modules/imgproc/src/emd_new.cpp @@ -0,0 +1,1010 @@ +// This file is part of OpenCV project. +// It is subject to the license terms in the LICENSE file found in the top-level directory +// of this distribution and at http://opencv.org/license.html + +/* + Partially based on Yossi Rubner code: + ========================================================================= + emd.c + + Last update: 3/14/98 + + An implementation of the Earth Movers Distance. + Based of the solution for the Transportation problem as described in + "Introduction to Mathematical Programming" by F. S. Hillier and + G. J. Lieberman, McGraw-Hill, 1990. + + Copyright (C) 1998 Yossi Rubner + Computer Science Department, Stanford University + E-Mail: rubner@cs.stanford.edu URL: http://vision.stanford.edu/~rubner + ========================================================================== +*/ + +#include "precomp.hpp" + +using namespace cv; + +namespace { + + +//============================================================================== +// Distance functions + +typedef float (*DistFunc)(const float* a, const float* b, int dims); + +static float distL1(const float* x, const float* y, int dims) +{ + double s = 0; + for (int i = 0; i < dims; i++) + { + const double t = x[i] - y[i]; + s += fabs(t); + } + return (float)s; +} + +static float distL2(const float* x, const float* y, int dims) +{ + double s = 0; + for (int i = 0; i < dims; i++) + { + const double t = x[i] - y[i]; + s += t * t; + } + return sqrt((float)s); +} + +static float distC(const float* x, const float* y, int dims) +{ + double s = 0; + for (int i = 0; i < dims; i++) + { + const double t = fabs(x[i] - y[i]); + if (s < t) + s = t; + } + return (float)s; +} + + +//============================================================================== +// Data structures + +/* Node1D is used for lists, representing 1D sparse array */ +struct Node1D +{ + float val; + Node1D* next; +}; + +/* Node2D is used for lists, representing 2D sparse matrix */ +struct Node2D +{ + float val; + int i, j; + Node2D* next[2]; /* next row & next column */ +}; + + +//============================================================================== +// Main class + +struct EMDSolver +{ + static constexpr int MAX_ITERATIONS = 500; + static constexpr float CV_EMD_INF = 1e20f; + static constexpr float CV_EMD_EPS = 1e-5f; + + int ssize, dsize; + + float* cost_buf; + Node2D* data_x; + Node2D* end_x; + Node2D* enter_x; + char* is_x; + + Node2D** rows_x; + Node2D** cols_x; + + Node1D* u; + Node1D* v; + + int* idx1; + int* idx2; + + /* find_loop buffers */ + Node2D** loop; + char* is_used; + + /* russel buffers */ + float* s; + float* d; + float* delta; + + float weight, max_cost; + + utils::BufferArea area, area2; + +public: + float getWeight() const + { + return weight; + } + + float& getCost(int i, int j) + { + return *(this->cost_buf + i * dsize + j); + } + const float& getCost(int i, int j) const + { + return *(this->cost_buf + i * dsize + j); + } + char& getIsX(int i, int j) + { + return *(this->is_x + i * dsize + j); + } + const char& getIsX(int i, int j) const + { + return *(this->is_x + i * dsize + j); + } + + EMDSolver() : + ssize(0), dsize(0), cost_buf(0), data_x(0), end_x(0), enter_x(0), is_x(0), rows_x(0), + cols_x(0), u(0), v(0), idx1(0), idx2(0), loop(0), is_used(0), s(0), d(0), delta(0), + weight(0), max_cost(0) + { + } + +public: + bool init(const Mat& sign1, + const Mat& sign2, + int dims, + DistFunc dfunc, + const Mat& cost, + float* lowerBound); + bool checkLowerBound(const Mat& sign1, + const Mat& sign2, + int dims, + DistFunc dfunc, + float& lowerBound); + bool calcSums(const Mat& sign1, const Mat& sign2); + float calcCost(const Mat& sign1, const Mat& sign2, int dims, DistFunc dfunc, const Mat& cost); + void solve(); + double calcFlow(Mat* flow_) const; + int findBasicVars() const; + float checkOptimal() const; + void callRussel(); + bool checkNewSolution(); + int findLoop() const; + void addBasicVar(int min_i, + int min_j, + Node1D* prev_u_min_i, + Node1D* prev_v_min_j, + Node1D* u_head); +}; + + +//============================================================================== +// Implementations + +bool EMDSolver::init(const Mat& sign1, + const Mat& sign2, + int dims, + DistFunc dfunc, + const Mat& cost, + float* lowerBound) +{ + const int size1 = sign1.size().height; + const int size2 = sign2.size().height; + + area.allocate(this->idx1, size1 + 1); + area.allocate(this->idx2, size2 + 1); + area.allocate(this->s, size1 + 1); + area.allocate(this->d, size2 + 1); + area.commit(); + area.zeroFill(); + + const bool areSumsEqual = calcSums(sign1, sign2); + if (areSumsEqual && lowerBound) + { + if (checkLowerBound(sign1, sign2, dims, dfunc, *lowerBound)) + return false; + } + + area2.allocate(this->u, ssize, 64); + area2.allocate(this->v, dsize, 64); + area2.allocate(this->is_used, ssize + dsize); + area2.allocate(this->delta, ssize * dsize); + area2.allocate(this->cost_buf, ssize * dsize); + area2.allocate(this->is_x, ssize * dsize); + area2.allocate(this->data_x, ssize + dsize, 64); + area2.allocate(this->rows_x, ssize, 64); + area2.allocate(this->cols_x, dsize, 64); + area2.allocate(this->loop, ssize + dsize + 1, 64); + area2.commit(); + area2.zeroFill(); + + this->end_x = this->data_x; + this->max_cost = calcCost(sign1, sign2, dims, dfunc, cost); + callRussel(); + this->enter_x = (this->end_x)++; + return true; +} + + +bool EMDSolver::checkLowerBound(const Mat& sign1, + const Mat& sign2, + int dims, + DistFunc dfunc, + float& lowerBound) +{ + AutoBuffer buf; + buf.allocate(dims * 2); + memset(buf.data(), 0, dims * 2 * sizeof(float)); + + float* xs = buf.data(); + float* xd = buf.data() + dims; + + for (int j = 0; j < sign1.rows; ++j) + { + const float weight_ = sign1.at(j, 0); + for (int i = 0; i < dims; i++) + xs[i] += sign1.at(j, i + 1) * weight_; + } + + for (int j = 0; j < sign2.rows; ++j) + { + const float weight_ = sign2.at(j, 0); + for (int i = 0; i < dims; i++) + xd[i] += sign2.at(j, i + 1) * weight_; + } + + const float lb = dfunc(xs, xd, dims) / this->weight; + const bool result = (lowerBound <= lb); + lowerBound = lb; + return result; +} + + +// return true if total sums of signatures are equal, false otherwise +bool EMDSolver::calcSums(const Mat& sign1, const Mat& sign2) +{ + bool result = true; + /* sum up the supply and demand */ + int ssize_ = 0, dsize_ = 0; + float s_sum = 0, d_sum = 0, diff; + for (int i = 0; i < sign1.size().height; i++) + { + const float weight_ = sign1.at(i, 0); + + if (weight_ > 0) + { + s_sum += weight_; + this->s[ssize_] = weight_; + this->idx1[ssize_++] = i; + } + else if (weight_ < 0) + CV_Error(cv::Error::StsBadArg, "sign1 must not contain negative weights"); + } + + for (int i = 0; i < sign2.size().height; i++) + { + const float weight_ = sign2.at(i, 0); + + if (weight_ > 0) + { + d_sum += weight_; + this->d[dsize_] = weight_; + this->idx2[dsize_++] = i; + } + else if (weight_ < 0) + CV_Error(cv::Error::StsBadArg, "sign2 must not contain negative weights"); + } + + if (ssize_ == 0) + CV_Error(cv::Error::StsBadArg, "sign1 must contain at least one non-zero value"); + if (dsize_ == 0) + CV_Error(cv::Error::StsBadArg, "sign2 must contain at least one non-zero value"); + + /* if supply different than the demand, add a zero-cost dummy cluster */ + diff = s_sum - d_sum; + if (fabs(diff) >= CV_EMD_EPS * s_sum) + { + result = false; + if (diff < 0) + { + this->s[ssize_] = -diff; + this->idx1[ssize_++] = -1; + } + else + { + this->d[dsize_] = diff; + this->idx2[dsize_++] = -1; + } + } + + this->ssize = ssize_; + this->dsize = dsize_; + this->weight = s_sum > d_sum ? s_sum : d_sum; + return result; +} + + +// returns maximum cost over all possible s->d combinations +float EMDSolver::calcCost(const Mat& sign1, + const Mat& sign2, + int dims, + DistFunc dfunc, + const Mat& cost) +{ + if (!dfunc) + { + CV_Assert(!cost.empty()); + } + float result = 0; + + /* compute the distance matrix */ + for (int i = 0; i < ssize; i++) + { + const int ci = this->idx1[i]; + if (ci >= 0) + { + for (int j = 0; j < dsize; j++) + { + const int cj = this->idx2[j]; + if (cj < 0) + getCost(i, j) = 0; + else + { + float val; + if (dfunc) + { + val = dfunc(sign1.ptr(ci, 1), sign2.ptr(cj, 1), dims); + } + else + { + val = cost.at(ci, cj); + } + getCost(i, j) = val; + if (result < val) + result = val; + } + } + } + else + { + for (int j = 0; j < dsize; j++) + getCost(i, j) = 0; + } + } + return result; +} + + +// runs solver iterations +void EMDSolver::solve() +{ + if (ssize > 1 && dsize > 1) + { + const float eps = CV_EMD_EPS * max_cost; + for (int itr = 1; itr < MAX_ITERATIONS; itr++) + { + /* find basic variables */ + if (findBasicVars() < 0) + break; + + /* check for optimality */ + const float min_delta = checkOptimal(); + + if (min_delta == CV_EMD_INF) + CV_Error(cv::Error::StsNoConv, ""); + + /* if no negative deltamin, we found the optimal solution */ + if (min_delta >= -eps) + break; + + /* improve solution */ + if (!checkNewSolution()) + CV_Error(cv::Error::StsNoConv, ""); + } + } +} + +double EMDSolver::calcFlow(Mat* flow_) const +{ + double result = 0.; + Node2D* xp = 0; + for (xp = data_x; xp < end_x; xp++) + { + float val = xp->val; + const int i = xp->i; + const int j = xp->j; + + if (xp == enter_x) + continue; + + const int ci = idx1[i]; + const int cj = idx2[j]; + + if (ci >= 0 && cj >= 0) + { + result += (double)val * getCost(i, j); + if (flow_) + { + flow_->at(ci, cj) = val; + } + } + } + return result; +} + + +int EMDSolver::findBasicVars() const +{ + int i, j; + int u_cfound, v_cfound; + Node1D u0_head, u1_head, *cur_u, *prev_u; + Node1D v0_head, v1_head, *cur_v, *prev_v; + bool found; + + CV_Assert(u != 0 && v != 0); + + /* initialize the rows list (u) and the columns list (v) */ + u0_head.next = u; + for (i = 0; i < ssize; i++) + { + u[i].next = u + i + 1; + } + u[ssize - 1].next = 0; + u1_head.next = 0; + + v0_head.next = ssize > 1 ? v + 1 : 0; + for (i = 1; i < dsize; i++) + { + v[i].next = v + i + 1; + } + v[dsize - 1].next = 0; + v1_head.next = 0; + + /* there are ssize+dsize variables but only ssize+dsize-1 independent equations, + so set v[0]=0 */ + v[0].val = 0; + v1_head.next = v; + v1_head.next->next = 0; + + /* loop until all variables are found */ + u_cfound = v_cfound = 0; + while (u_cfound < ssize || v_cfound < dsize) + { + found = false; + if (v_cfound < dsize) + { + /* loop over all marked columns */ + prev_v = &v1_head; + cur_v = v1_head.next; + found = found || (cur_v != 0); + for (; cur_v != 0; cur_v = cur_v->next) + { + float cur_v_val = cur_v->val; + + j = (int)(cur_v - v); + /* find the variables in column j */ + prev_u = &u0_head; + for (cur_u = u0_head.next; cur_u != 0;) + { + i = (int)(cur_u - u); + if (getIsX(i, j)) + { + /* compute u[i] */ + cur_u->val = getCost(i, j) - cur_v_val; + /* ...and add it to the marked list */ + prev_u->next = cur_u->next; + cur_u->next = u1_head.next; + u1_head.next = cur_u; + cur_u = prev_u->next; + } + else + { + prev_u = cur_u; + cur_u = cur_u->next; + } + } + prev_v->next = cur_v->next; + v_cfound++; + } + } + + if (u_cfound < ssize) + { + /* loop over all marked rows */ + prev_u = &u1_head; + cur_u = u1_head.next; + found = found || (cur_u != 0); + for (; cur_u != 0; cur_u = cur_u->next) + { + float cur_u_val = cur_u->val; + i = (int)(cur_u - u); + /* find the variables in rows i */ + prev_v = &v0_head; + for (cur_v = v0_head.next; cur_v != 0;) + { + j = (int)(cur_v - v); + if (getIsX(i, j)) + { + /* compute v[j] */ + cur_v->val = getCost(i, j) - cur_u_val; + /* ...and add it to the marked list */ + prev_v->next = cur_v->next; + cur_v->next = v1_head.next; + v1_head.next = cur_v; + cur_v = prev_v->next; + } + else + { + prev_v = cur_v; + cur_v = cur_v->next; + } + } + prev_u->next = cur_u->next; + u_cfound++; + } + } + + if (!found) + return -1; + } + + return 0; +} + +float EMDSolver::checkOptimal() const +{ + int i, j, min_i = 0, min_j = 0; + + float min_delta = CV_EMD_INF; + /* find the minimal cij-ui-vj over all i,j */ + for (i = 0; i < ssize; i++) + { + float u_val = u[i].val; + for (j = 0; j < dsize; j++) + { + if (!getIsX(i, j)) + { + const float delta_ = getCost(i, j) - u_val - v[j].val; + if (min_delta > delta_) + { + min_delta = delta_; + min_i = i; + min_j = j; + } + } + } + } + + enter_x->i = min_i; + enter_x->j = min_j; + + return min_delta; +} + +bool EMDSolver::checkNewSolution() +{ + int i, j; + float min_val = CV_EMD_INF; + int steps; + Node2D head {0, 0, 0, {0, 0}}, *cur_x, *next_x, *leave_x = 0; + Node2D* enter_x_ = this->enter_x; + Node2D** loop_ = this->loop; + + /* enter the new basic variable */ + i = enter_x_->i; + j = enter_x_->j; + getIsX(i, j) = 1; + enter_x_->next[0] = this->rows_x[i]; + enter_x->next[1] = this->cols_x[j]; + enter_x_->val = 0; + this->rows_x[i] = enter_x_; + this->cols_x[j] = enter_x_; + + /* find a chain reaction */ + steps = findLoop(); + + if (steps == 0) + return false; + + /* find the largest value in the loop */ + for (i = 1; i < steps; i += 2) + { + float temp = loop_[i]->val; + + if (min_val > temp) + { + leave_x = loop_[i]; + min_val = temp; + } + } + + /* update the loop */ + for (i = 0; i < steps; i += 2) + { + float temp0 = loop_[i]->val + min_val; + float temp1 = loop_[i + 1]->val - min_val; + + loop_[i]->val = temp0; + loop_[i + 1]->val = temp1; + } + + /* remove the leaving basic variable */ + CV_Assert(leave_x != NULL); + i = leave_x->i; + j = leave_x->j; + getIsX(i, j) = 0; + + head.next[0] = this->rows_x[i]; + cur_x = &head; + while ((next_x = cur_x->next[0]) != leave_x) + { + cur_x = next_x; + CV_Assert(cur_x); + } + cur_x->next[0] = next_x->next[0]; + this->rows_x[i] = head.next[0]; + + head.next[1] = this->cols_x[j]; + cur_x = &head; + while ((next_x = cur_x->next[1]) != leave_x) + { + cur_x = next_x; + CV_Assert(cur_x); + } + cur_x->next[1] = next_x->next[1]; + this->cols_x[j] = head.next[1]; + + /* set enter_x to be the new empty slot */ + this->enter_x = leave_x; + + return true; +} + +int EMDSolver::findLoop() const +{ + int i; + + memset(is_used, 0, this->ssize + this->dsize); + + Node2D* new_x = loop[0] = enter_x; + is_used[enter_x - data_x] = 1; + int steps = 1; + + do + { + if ((steps & 1) == 1) + { + /* find an unused x in the row */ + new_x = this->rows_x[new_x->i]; + while (new_x != 0 && is_used[new_x - data_x]) + new_x = new_x->next[0]; + } + else + { + /* find an unused x in the column, or the entering x */ + new_x = this->cols_x[new_x->j]; + while (new_x != 0 && is_used[new_x - data_x] && new_x != enter_x) + new_x = new_x->next[1]; + if (new_x == enter_x) + break; + } + + if (new_x != 0) /* found the next x */ + { + /* add x to the loop */ + loop[steps++] = new_x; + is_used[new_x - data_x] = 1; + } + else /* didn't find the next x */ + { + /* backtrack */ + do + { + i = steps & 1; + new_x = loop[steps - 1]; + do + { + new_x = new_x->next[i]; + } + while (new_x != 0 && is_used[new_x - data_x]); + + if (new_x == 0) + { + is_used[loop[--steps] - data_x] = 0; + } + } + while (new_x == 0 && steps > 0); + + is_used[loop[steps - 1] - data_x] = 0; + loop[steps - 1] = new_x; + is_used[new_x - data_x] = 1; + } + } + while (steps > 0); + + return steps; +} + +void EMDSolver::callRussel() +{ + int i, j, min_i = -1, min_j = -1; + float min_delta, diff; + Node1D u_head, *cur_u, *prev_u; + Node1D v_head, *cur_v, *prev_v; + Node1D *prev_u_min_i = 0, *prev_v_min_j = 0, *remember; + float eps = CV_EMD_EPS * this->max_cost; + + /* initialize the rows list (ur), and the columns list (vr) */ + u_head.next = u; + for (i = 0; i < ssize; i++) + { + u[i].next = u + i + 1; + } + u[ssize - 1].next = 0; + + v_head.next = v; + for (i = 0; i < dsize; i++) + { + v[i].val = -CV_EMD_INF; + v[i].next = v + i + 1; + } + v[dsize - 1].next = 0; + + /* find the maximum row and column values (ur[i] and vr[j]) */ + for (i = 0; i < ssize; i++) + { + float u_val = -CV_EMD_INF; + for (j = 0; j < dsize; j++) + { + float temp = getCost(i, j); + + if (u_val < temp) + u_val = temp; + if (v[j].val < temp) + v[j].val = temp; + } + u[i].val = u_val; + } + + /* compute the delta matrix */ + for (i = 0; i < ssize; i++) + { + float u_val = u[i].val; + float* delta_row = delta + i * dsize; + for (j = 0; j < dsize; j++) + { + delta_row[j] = getCost(i, j) - u_val - v[j].val; + } + } + + /* find the basic variables */ + do + { + /* find the smallest delta[i][j] */ + min_i = -1; + min_delta = CV_EMD_INF; + prev_u = &u_head; + for (cur_u = u_head.next; cur_u != 0; cur_u = cur_u->next) + { + i = (int)(cur_u - u); + float* delta_row = delta + i * dsize; + + prev_v = &v_head; + for (cur_v = v_head.next; cur_v != 0; cur_v = cur_v->next) + { + j = (int)(cur_v - v); + if (min_delta > delta_row[j]) + { + min_delta = delta_row[j]; + min_i = i; + min_j = j; + prev_u_min_i = prev_u; + prev_v_min_j = prev_v; + } + prev_v = cur_v; + } + prev_u = cur_u; + } + + if (min_i < 0) + break; + + /* add x[min_i][min_j] to the basis, and adjust supplies and cost */ + remember = prev_u_min_i->next; + addBasicVar(min_i, min_j, prev_u_min_i, prev_v_min_j, &u_head); + + /* update the necessary delta[][] */ + if (remember == prev_u_min_i->next) /* line min_i was deleted */ + { + for (cur_v = v_head.next; cur_v != 0; cur_v = cur_v->next) + { + j = (int)(cur_v - v); + if (cur_v->val == getCost(min_i, j)) /* column j needs updating */ + { + float max_val = -CV_EMD_INF; + + /* find the new maximum value in the column */ + for (cur_u = u_head.next; cur_u != 0; cur_u = cur_u->next) + { + float temp = getCost((int)(cur_u - u), j); + + if (max_val < temp) + max_val = temp; + } + + /* if needed, adjust the relevant delta[*][j] */ + diff = max_val - cur_v->val; + cur_v->val = max_val; + if (fabs(diff) < eps) + { + for (cur_u = u_head.next; cur_u != 0; cur_u = cur_u->next) + *(delta + (cur_u - u) * dsize + j) += diff; + } + } + } + } + else /* column min_j was deleted */ + { + for (cur_u = u_head.next; cur_u != 0; cur_u = cur_u->next) + { + i = (int)(cur_u - u); + if (cur_u->val == getCost(i, min_j)) /* row i needs updating */ + { + float max_val = -CV_EMD_INF; + + /* find the new maximum value in the row */ + for (cur_v = v_head.next; cur_v != 0; cur_v = cur_v->next) + { + float temp = getCost(i, (int)(cur_v - v)); + + if (max_val < temp) + max_val = temp; + } + + /* if needed, adjust the relevant delta[i][*] */ + diff = max_val - cur_u->val; + cur_u->val = max_val; + + if (fabs(diff) < eps) + { + for (cur_v = v_head.next; cur_v != 0; cur_v = cur_v->next) + *(delta + i * dsize + (cur_v - v)) += diff; + } + } + } + } + } + while (u_head.next != 0 || v_head.next != 0); +} + +void EMDSolver::addBasicVar(int min_i, + int min_j, + Node1D* prev_u_min_i, + Node1D* prev_v_min_j, + Node1D* u_head) +{ + float temp; + + if (this->s[min_i] < this->d[min_j] + this->weight * CV_EMD_EPS) + { /* supply exhausted */ + temp = this->s[min_i]; + this->s[min_i] = 0; + this->d[min_j] -= temp; + } + else /* demand exhausted */ + { + temp = this->d[min_j]; + this->d[min_j] = 0; + this->s[min_i] -= temp; + } + + /* x(min_i,min_j) is a basic variable */ + getIsX(min_i, min_j) = 1; + + end_x->val = temp; + end_x->i = min_i; + end_x->j = min_j; + end_x->next[0] = this->rows_x[min_i]; + end_x->next[1] = this->cols_x[min_j]; + this->rows_x[min_i] = end_x; + this->cols_x[min_j] = end_x; + this->end_x = end_x + 1; + + /* delete supply row only if the empty, and if not last row */ + if (this->s[min_i] == 0 && u_head->next->next != 0) + prev_u_min_i->next = prev_u_min_i->next->next; /* remove row from list */ + else + prev_v_min_j->next = prev_v_min_j->next->next; /* remove column from list */ +} + +} // namespace + + +//============================================================================== +// External interface + +float cv::EMD(InputArray _sign1, + InputArray _sign2, + int distType, + InputArray _cost, + float* lowerBound, + OutputArray _flow) +{ + CV_INSTRUMENT_REGION(); + + Mat sign1 = _sign1.getMat(); + Mat sign2 = _sign2.getMat(); + Mat cost = _cost.getMat(); + + CV_CheckEQ(sign1.cols, sign2.cols, "Signatures must have equal number of columns"); + CV_CheckEQ(sign1.type(), CV_32FC1, "The sign1 must be 32FC1"); + CV_CheckEQ(sign2.type(), CV_32FC1, "The sign2 must be 32FC1"); + + const int dims = sign1.cols - 1; + const int size1 = sign1.rows; + const int size2 = sign2.rows; + + Mat flow; + if (_flow.needed()) + { + _flow.create(sign1.rows, sign2.rows, CV_32F); + flow = _flow.getMat(); + flow = Scalar::all(0); + CV_CheckEQ(flow.type(), CV_32FC1, "Flow matrix must have type 32FC1"); + CV_CheckTrue(flow.rows == size1 && flow.cols == size2, + "Flow matrix size does not match signatures"); + } + + DistFunc dfunc = 0; + if (distType == DIST_USER) + { + if (!cost.empty()) + { + CV_CheckEQ(cost.type(), CV_32FC1, "Cost matrix must have type 32FC1"); + CV_CheckTrue(cost.rows == size1 && cost.cols == size2, + "Cost matrix size does not match signatures"); + CV_CheckTrue(lowerBound == NULL, + "Lower boundary can not be calculated if the cost matrix is used"); + } + else + { + CV_CheckTrue(dfunc == NULL, "Dist function must be set if cost matrix is empty"); + } + } + else + { + CV_CheckNE(dims, 0, "Number of dimensions can be 0 only if a user-defined metric is used"); + switch (distType) + { + case cv::DIST_L1: dfunc = distL1; break; + case cv::DIST_L2: dfunc = distL2; break; + case cv::DIST_C: dfunc = distC; break; + default: CV_Error(cv::Error::StsBadFlag, "Bad or unsupported metric type"); + } + } + + EMDSolver state; + const bool result = state.init(sign1, sign2, dims, dfunc, cost, lowerBound); + if (!result && lowerBound) + { + return *lowerBound; + } + state.solve(); + return (float)(state.calcFlow(_flow.needed() ? &flow : 0) / state.getWeight()); +} + +float cv::wrapperEMD(InputArray _sign1, + InputArray _sign2, + int distType, + InputArray _cost, + Ptr lowerBound, + OutputArray _flow) +{ + return EMD(_sign1, _sign2, distType, _cost, lowerBound.get(), _flow); +} diff --git a/modules/imgproc/src/precomp.hpp b/modules/imgproc/src/precomp.hpp index 29dbce6d5f..df35da529c 100644 --- a/modules/imgproc/src/precomp.hpp +++ b/modules/imgproc/src/precomp.hpp @@ -50,6 +50,8 @@ #include "opencv2/core/private.hpp" #include "opencv2/core/ocl.hpp" #include "opencv2/core/hal/hal.hpp" +#include "opencv2/core/check.hpp" +#include "opencv2/core/utils/buffer_area.private.hpp" #include "opencv2/imgproc/hal/hal.hpp" #include "hal_replacement.hpp" diff --git a/modules/imgproc/test/test_emd.cpp b/modules/imgproc/test/test_emd.cpp index 51f05f374e..59cd8381a1 100644 --- a/modules/imgproc/test/test_emd.cpp +++ b/modules/imgproc/test/test_emd.cpp @@ -1,93 +1,250 @@ -/*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. -// -// -// Intel License Agreement -// For Open Source Computer Vision Library -// -// Copyright (C) 2000, Intel Corporation, 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 Intel Corporation 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*/ +// This file is part of OpenCV project. +// It is subject to the license terms in the LICENSE file found in the top-level directory +// of this distribution and at http://opencv.org/license.html +#include "opencv2/imgproc.hpp" #include "test_precomp.hpp" +using namespace cv; +using namespace std; + namespace opencv_test { namespace { -class CV_EMDTest : public cvtest::BaseTest -{ -public: - CV_EMDTest(); -protected: - void run(int); -}; +//============================================================================== +// Utility - -CV_EMDTest::CV_EMDTest() +template +inline T sqr(T val) { + return val * val; } -void CV_EMDTest::run( int ) +inline static float calcEMD(Mat w1, Mat w2, Mat& flow, int dist, int dims) { - int code = cvtest::TS::OK; - const double success_error_level = 1e-6; - #define M 10000 - double emd0 = 2460./210; - static float cost[] = + float mass1 = 0.f, mass2 = 0.f, work = 0.f; + for (int i = 0; i < flow.rows; ++i) { - 16, 16, 13, 22, 17, - 14, 14, 13, 19, 15, - 19, 19, 20, 23, M, - M , 0, M, 0, 0 - }; - static float w1[] = { 50, 60, 50, 50 }, - w2[] = { 30, 20, 70, 30, 60 }; - Mat _w1(4, 1, CV_32F, w1); - Mat _w2(5, 1, CV_32F, w2); - Mat _cost(_w1.rows, _w2.rows, CV_32F, cost); + mass1 += w1.at(i, 0); + for (int j = 0; j < flow.cols; ++j) + { + if (i == 0) + mass2 += w2.at(j, 0); + float dist_ = 0.f; + switch (dist) + { + case DIST_L1: + { + for (int k = 1; k <= dims; ++k) + { + dist_ += abs(w1.at(i, k) - w2.at(j, k)); + } + break; + } + case DIST_L2: + { + for (int k = 1; k <= dims; ++k) + { + dist_ += sqr(w1.at(i, k) - w2.at(j, k)); + } + dist_ = sqrt(dist_); + break; + } + case DIST_C: + { + for (int k = 1; k <= dims; ++k) + { + const float val = abs(w1.at(i, k) - w2.at(j, k)); + if (val > dist_) + dist_ = val; + } + break; + } + } + const float weight = flow.at(i, j); + work += dist_ * weight; + } + } + return work / max(mass1, mass2); +} - float emd = EMD( _w1, _w2, -1, _cost ); - if( fabs( emd - emd0 ) > success_error_level*emd0 ) +//============================================================================== + +TEST(Imgproc_EMD, regression) +{ + // input data + const float M = 10000; + Matx w1 {50, 60, 50, 50}; + Matx w2 {30, 20, 70, 30, 60}; + Matx cost {16, 16, 13, 22, 17, 14, 14, 13, 19, 15, + 19, 19, 20, 23, M, M, 0, M, 0, 0}; + + // expected results + const double emd0 = 2460. / 210; + Matx flow0 {0, 0, 50, 0, 0, 0, 0, 20, 0, 40, 30, 20, 0, 0, 0, 0, 0, 0, 30, 20}; + + // basic call with cost { - ts->printf( cvtest::TS::LOG, - "The computed distance is %.2f, while it should be %.2f\n", emd, emd0 ); - code = cvtest::TS::FAIL_BAD_ACCURACY; + float emd = 0.f; + ASSERT_NO_THROW(emd = EMD(w1, w2, DIST_USER, cost)); + EXPECT_NEAR(emd, emd0, 1e-6 * emd0); } - if( code < 0 ) - ts->set_failed_test_info( code ); + // basic call with cost and flow output + { + Mat flow; + float emd = 0.f; + ASSERT_NO_THROW(emd = EMD(w1, w2, DIST_USER, cost, nullptr, flow)); + EXPECT_NEAR(emd, emd0, 1e-6 * emd0); + EXPECT_MAT_NEAR(Mat(flow0), flow, 1e-6); + } + // no cost and DIST_USER - error + { + Mat flow; + EXPECT_THROW(EMD(w1, w2, DIST_USER, noArray(), nullptr, flow), cv::Exception); + EXPECT_THROW(EMD(w1, w2, DIST_USER), cv::Exception); + } } -TEST(Imgproc_EMD, regression) { CV_EMDTest test; test.safe_run(); } +TEST(Imgproc_EMD, distance_types) +{ + // 1D (sum = 210) + Matx w1 {50, 1, 60, 2, 50, 3, 50, 4}; + Matx w2 {30, 1, 20, 2, 70, 3, 30, 4, 60, 5}; -}} // namespace -/* End of file. */ + // 2D (sum = 210) + Matx w3 {50, 0, 0, 60, 0, 1, 50, 1, 0, 50, 1, 1}; + Matx w4 {20, 0, 1, 70, 1, 0, 30, 1, 1, 60, 2, 2, 30, 3, 3}; + + // basic call with all distance types + { + const vector good_types {DIST_L1, DIST_L2, DIST_C}; + for (const auto& dt : good_types) + { + SCOPED_TRACE(cv::format("dt=%d", dt)); + float emd = 0.f; + Mat flow; + // 1D + { + ASSERT_NO_THROW(emd = EMD(w1, w2, dt, noArray(), nullptr, flow)); + const float emd0 = calcEMD(Mat(w1), Mat(w2), flow, dt, 1); + EXPECT_NEAR(emd0, emd, 1e-6); + } + // 2D + { + ASSERT_NO_THROW(emd = EMD(w3, w4, dt, noArray(), nullptr, flow)); + const float emd0 = calcEMD(Mat(w3), Mat(w4), flow, dt, 2); + EXPECT_NEAR(emd0, emd, 1e-6); + } + } + } +} + +typedef testing::TestWithParam Imgproc_EMD_dist; + +TEST_P(Imgproc_EMD_dist, random_flow_verify) +{ + const int dist = GetParam(); + for (size_t iter = 0; iter < 100; ++iter) + { + SCOPED_TRACE(cv::format("iter=%zu", iter)); + RNG& rng = TS::ptr()->get_rng(); + const int dims = rng.uniform(1, 10); + Mat w1(rng.uniform(1, 10), dims + 1, CV_32FC1); + Mat w2(rng.uniform(1, 10), dims + 1, CV_32FC1); + + // weights > 0 + { + Mat w1_weights = w1.col(0); + Mat w2_weights = w2.col(0); + cvtest::randUni(rng, w1_weights, 0, 100); + cvtest::randUni(rng, w2_weights, 0, 100); + } + + // coord + { + Mat w1_coord = w1.colRange(1, dims + 1); + Mat w2_coord = w2.colRange(1, dims + 1); + cvtest::randUni(rng, w1_coord, -10, +10); + cvtest::randUni(rng, w2_coord, -10, +10); + } + + float emd1 = 0.f, emd2 = 0.f; + const float eps = 1e-5f; + Mat flow; + { + ASSERT_NO_THROW(emd1 = EMD(w1, w2, dist, noArray(), nullptr, flow)); + const float emd0 = calcEMD(w1, w2, flow, dist, dims); + EXPECT_NEAR(emd0, emd1, eps); + } + { + ASSERT_NO_THROW(emd2 = EMD(w2, w1, dist, noArray(), nullptr, flow)); + const float emd0 = calcEMD(w2, w1, flow, dist, dims); + EXPECT_NEAR(emd0, emd2, eps); + } + EXPECT_NEAR(emd1, emd2, eps); + } +} + +INSTANTIATE_TEST_CASE_P(, Imgproc_EMD_dist, testing::Values(DIST_L1, DIST_L2, DIST_C)); + + +TEST(Imgproc_EMD, invalid) +{ + Matx w1 {50, 1, 60, 2, 50, 3, 50, 4}; + Matx w2 {30, 1, 20, 2, 70, 3, 30, 4, 60, 5}; + + // empty signature + { + Mat empty; + EXPECT_THROW(EMD(empty, w2, DIST_USER), cv::Exception); + EXPECT_THROW(EMD(w1, empty, DIST_USER), cv::Exception); + } + + // zero total weight, negative weight + { + Matx wz {0, 0, 0}; + Matx wz1 {0, 1, 0, 2, 0, 3}; + Matx wn {0, 3, -2}; + Matx wn1 {0, 1, 3, 2, -2, 3}; + EXPECT_THROW(EMD(wz, w2, DIST_USER), cv::Exception); + EXPECT_THROW(EMD(wz1, w2, DIST_USER), cv::Exception); + EXPECT_THROW(EMD(wn, w2, DIST_USER), cv::Exception); + EXPECT_THROW(EMD(wn1, w2, DIST_USER), cv::Exception); + } + + // user distance type, but no cost matrix provided or is wrong + { + Mat cost(3, 3, CV_32FC1, Scalar::all(0)), cost8u(4, 5, CV_8UC1, Scalar::all(0)), empty; + EXPECT_THROW(EMD(w1, w2, DIST_USER, noArray()), cv::Exception); + EXPECT_THROW(EMD(w1, w2, DIST_USER, empty), cv::Exception); + EXPECT_THROW(EMD(w1, w2, DIST_USER, cost8u), cv::Exception); + EXPECT_THROW(EMD(w1, w2, DIST_USER, cost), cv::Exception); + } + + // lower_bound is set together with cost + { + Mat cost(4, 5, CV_32FC1, Scalar::all(0)); + float bound = 0.f; + EXPECT_THROW(EMD(w1, w2, DIST_USER, cost, &bound), cv::Exception); + } + + // zero dimensions with non-user distance type + const vector good_types {DIST_L1, DIST_L2, DIST_C}; + for (const auto& dt : good_types) + { + SCOPED_TRACE(cv::format("dt=%d", dt)); + Matx w01 {20, 30, 40, 50}; + Matx w02 {20, 30, 40, 50, 10}; + EXPECT_THROW(EMD(w01, w02, dt), cv::Exception); + } + + // wrong distance type + const vector bad_types {DIST_L12, DIST_FAIR, DIST_WELSCH, DIST_HUBER}; + for (const auto& dt : bad_types) + { + SCOPED_TRACE(cv::format("dt=%d", dt)); + EXPECT_THROW(EMD(w1, w2, dt), cv::Exception); + } +} + +}} // namespace opencv_test diff --git a/modules/imgproc/test/test_precomp.hpp b/modules/imgproc/test/test_precomp.hpp index ce51009145..4df4a922b5 100644 --- a/modules/imgproc/test/test_precomp.hpp +++ b/modules/imgproc/test/test_precomp.hpp @@ -5,8 +5,11 @@ #define __OPENCV_TEST_PRECOMP_HPP__ #include "opencv2/ts.hpp" +#include "opencv2/ts/ts_gtest.h" +#include "opencv2/ts/ocl_test.hpp" #include "opencv2/imgproc.hpp" #include "opencv2/imgproc/imgproc_c.h" +#include "opencv2/core.hpp" #include "opencv2/core/private.hpp"