mirror of
https://github.com/opencv/opencv.git
synced 2026-07-31 00:03:03 +04:00
Merge pull request #3103 from vpisarev:core_imgproc_optim_rearrangements
This commit is contained in:
Vadim Pisarevsky
@@ -0,0 +1,186 @@
|
||||
/*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.
|
||||
//
|
||||
//
|
||||
// License Agreement
|
||||
// For Open Source Computer Vision Library
|
||||
//
|
||||
// Copyright (C) 2013, OpenCV Foundation, 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 the copyright holders 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 OpenCV Foundation 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*/
|
||||
|
||||
#include "precomp.hpp"
|
||||
|
||||
#define dprintf(x)
|
||||
#define print_matrix(x)
|
||||
|
||||
namespace cv
|
||||
{
|
||||
|
||||
#define SEC_METHOD_ITERATIONS 4
|
||||
#define INITIAL_SEC_METHOD_SIGMA 0.1
|
||||
class ConjGradSolverImpl : public ConjGradSolver
|
||||
{
|
||||
public:
|
||||
Ptr<Function> getFunction() const;
|
||||
void setFunction(const Ptr<Function>& f);
|
||||
TermCriteria getTermCriteria() const;
|
||||
ConjGradSolverImpl();
|
||||
void setTermCriteria(const TermCriteria& termcrit);
|
||||
double minimize(InputOutputArray x);
|
||||
protected:
|
||||
Ptr<MinProblemSolver::Function> _Function;
|
||||
TermCriteria _termcrit;
|
||||
Mat_<double> d,r,buf_x,r_old;
|
||||
Mat_<double> minimizeOnTheLine_buf1,minimizeOnTheLine_buf2;
|
||||
private:
|
||||
static void minimizeOnTheLine(Ptr<MinProblemSolver::Function> _f,Mat_<double>& x,const Mat_<double>& d,Mat_<double>& buf1,Mat_<double>& buf2);
|
||||
};
|
||||
|
||||
void ConjGradSolverImpl::minimizeOnTheLine(Ptr<MinProblemSolver::Function> _f,Mat_<double>& x,const Mat_<double>& d,Mat_<double>& buf1,
|
||||
Mat_<double>& buf2){
|
||||
double sigma=INITIAL_SEC_METHOD_SIGMA;
|
||||
buf1=0.0;
|
||||
buf2=0.0;
|
||||
|
||||
dprintf(("before minimizeOnTheLine\n"));
|
||||
dprintf(("x:\n"));
|
||||
print_matrix(x);
|
||||
dprintf(("d:\n"));
|
||||
print_matrix(d);
|
||||
|
||||
for(int i=0;i<SEC_METHOD_ITERATIONS;i++){
|
||||
_f->getGradient((double*)x.data,(double*)buf1.data);
|
||||
dprintf(("buf1:\n"));
|
||||
print_matrix(buf1);
|
||||
x=x+sigma*d;
|
||||
_f->getGradient((double*)x.data,(double*)buf2.data);
|
||||
dprintf(("buf2:\n"));
|
||||
print_matrix(buf2);
|
||||
double d1=buf1.dot(d), d2=buf2.dot(d);
|
||||
if((d1-d2)==0){
|
||||
break;
|
||||
}
|
||||
double alpha=-sigma*d1/(d2-d1);
|
||||
dprintf(("(buf2.dot(d)-buf1.dot(d))=%f\nalpha=%f\n",(buf2.dot(d)-buf1.dot(d)),alpha));
|
||||
x=x+(alpha-sigma)*d;
|
||||
sigma=-alpha;
|
||||
}
|
||||
|
||||
dprintf(("after minimizeOnTheLine\n"));
|
||||
print_matrix(x);
|
||||
}
|
||||
|
||||
double ConjGradSolverImpl::minimize(InputOutputArray x){
|
||||
CV_Assert(_Function.empty()==false);
|
||||
dprintf(("termcrit:\n\ttype: %d\n\tmaxCount: %d\n\tEPS: %g\n",_termcrit.type,_termcrit.maxCount,_termcrit.epsilon));
|
||||
|
||||
Mat x_mat=x.getMat();
|
||||
CV_Assert(MIN(x_mat.rows,x_mat.cols)==1);
|
||||
int ndim=MAX(x_mat.rows,x_mat.cols);
|
||||
CV_Assert(x_mat.type()==CV_64FC1);
|
||||
|
||||
if(d.cols!=ndim){
|
||||
d.create(1,ndim);
|
||||
r.create(1,ndim);
|
||||
r_old.create(1,ndim);
|
||||
minimizeOnTheLine_buf1.create(1,ndim);
|
||||
minimizeOnTheLine_buf2.create(1,ndim);
|
||||
}
|
||||
|
||||
Mat_<double> proxy_x;
|
||||
if(x_mat.rows>1){
|
||||
buf_x.create(1,ndim);
|
||||
Mat_<double> proxy(ndim,1,buf_x.ptr<double>());
|
||||
x_mat.copyTo(proxy);
|
||||
proxy_x=buf_x;
|
||||
}else{
|
||||
proxy_x=x_mat;
|
||||
}
|
||||
_Function->getGradient(proxy_x.ptr<double>(),d.ptr<double>());
|
||||
d*=-1.0;
|
||||
d.copyTo(r);
|
||||
|
||||
//here everything goes. check that everything is setted properly
|
||||
dprintf(("proxy_x\n"));print_matrix(proxy_x);
|
||||
dprintf(("d first time\n"));print_matrix(d);
|
||||
dprintf(("r\n"));print_matrix(r);
|
||||
|
||||
double beta=0;
|
||||
for(int count=0;count<_termcrit.maxCount;count++){
|
||||
minimizeOnTheLine(_Function,proxy_x,d,minimizeOnTheLine_buf1,minimizeOnTheLine_buf2);
|
||||
r.copyTo(r_old);
|
||||
_Function->getGradient(proxy_x.ptr<double>(),r.ptr<double>());
|
||||
r*=-1.0;
|
||||
double r_norm_sq=norm(r);
|
||||
if(_termcrit.type==(TermCriteria::MAX_ITER+TermCriteria::EPS) && r_norm_sq<_termcrit.epsilon){
|
||||
break;
|
||||
}
|
||||
r_norm_sq=r_norm_sq*r_norm_sq;
|
||||
beta=MAX(0.0,(r_norm_sq-r.dot(r_old))/r_norm_sq);
|
||||
d=r+beta*d;
|
||||
}
|
||||
|
||||
|
||||
|
||||
if(x_mat.rows>1){
|
||||
Mat(ndim, 1, CV_64F, proxy_x.ptr<double>()).copyTo(x);
|
||||
}
|
||||
return _Function->calc(proxy_x.ptr<double>());
|
||||
}
|
||||
|
||||
ConjGradSolverImpl::ConjGradSolverImpl(){
|
||||
_Function=Ptr<Function>();
|
||||
}
|
||||
Ptr<MinProblemSolver::Function> ConjGradSolverImpl::getFunction()const{
|
||||
return _Function;
|
||||
}
|
||||
void ConjGradSolverImpl::setFunction(const Ptr<Function>& f){
|
||||
_Function=f;
|
||||
}
|
||||
TermCriteria ConjGradSolverImpl::getTermCriteria()const{
|
||||
return _termcrit;
|
||||
}
|
||||
void ConjGradSolverImpl::setTermCriteria(const TermCriteria& termcrit){
|
||||
CV_Assert((termcrit.type==(TermCriteria::MAX_ITER+TermCriteria::EPS) && termcrit.epsilon>0 && termcrit.maxCount>0) ||
|
||||
((termcrit.type==TermCriteria::MAX_ITER) && termcrit.maxCount>0));
|
||||
_termcrit=termcrit;
|
||||
}
|
||||
// both minRange & minError are specified by termcrit.epsilon; In addition, user may specify the number of iterations that the algorithm does.
|
||||
Ptr<ConjGradSolver> ConjGradSolver::create(const Ptr<MinProblemSolver::Function>& f, TermCriteria termcrit){
|
||||
Ptr<ConjGradSolver> CG = makePtr<ConjGradSolverImpl>();
|
||||
CG->setFunction(f);
|
||||
CG->setTermCriteria(termcrit);
|
||||
return CG;
|
||||
}
|
||||
}
|
||||
@@ -3528,492 +3528,9 @@ cvPrevTreeNode( CvTreeNodeIterator* treeIterator )
|
||||
return prevNode;
|
||||
}
|
||||
|
||||
|
||||
namespace cv
|
||||
{
|
||||
|
||||
// This is reimplementation of kd-trees from cvkdtree*.* by Xavier Delacour, cleaned-up and
|
||||
// adopted to work with the new OpenCV data structures. It's in cxcore to be shared by
|
||||
// both cv (CvFeatureTree) and ml (kNN).
|
||||
|
||||
// The algorithm is taken from:
|
||||
// J.S. Beis and D.G. Lowe. Shape indexing using approximate nearest-neighbor search
|
||||
// in highdimensional spaces. In Proc. IEEE Conf. Comp. Vision Patt. Recog.,
|
||||
// pages 1000--1006, 1997. http://citeseer.ist.psu.edu/beis97shape.html
|
||||
|
||||
const int MAX_TREE_DEPTH = 32;
|
||||
|
||||
KDTree::KDTree()
|
||||
{
|
||||
maxDepth = -1;
|
||||
normType = NORM_L2;
|
||||
}
|
||||
|
||||
KDTree::KDTree(InputArray _points, bool _copyData)
|
||||
{
|
||||
maxDepth = -1;
|
||||
normType = NORM_L2;
|
||||
build(_points, _copyData);
|
||||
}
|
||||
|
||||
KDTree::KDTree(InputArray _points, InputArray _labels, bool _copyData)
|
||||
{
|
||||
maxDepth = -1;
|
||||
normType = NORM_L2;
|
||||
build(_points, _labels, _copyData);
|
||||
}
|
||||
|
||||
struct SubTree
|
||||
{
|
||||
SubTree() : first(0), last(0), nodeIdx(0), depth(0) {}
|
||||
SubTree(int _first, int _last, int _nodeIdx, int _depth)
|
||||
: first(_first), last(_last), nodeIdx(_nodeIdx), depth(_depth) {}
|
||||
int first;
|
||||
int last;
|
||||
int nodeIdx;
|
||||
int depth;
|
||||
};
|
||||
|
||||
|
||||
static float
|
||||
medianPartition( size_t* ofs, int a, int b, const float* vals )
|
||||
{
|
||||
int k, a0 = a, b0 = b;
|
||||
int middle = (a + b)/2;
|
||||
while( b > a )
|
||||
{
|
||||
int i0 = a, i1 = (a+b)/2, i2 = b;
|
||||
float v0 = vals[ofs[i0]], v1 = vals[ofs[i1]], v2 = vals[ofs[i2]];
|
||||
int ip = v0 < v1 ? (v1 < v2 ? i1 : v0 < v2 ? i2 : i0) :
|
||||
v0 < v2 ? i0 : (v1 < v2 ? i2 : i1);
|
||||
float pivot = vals[ofs[ip]];
|
||||
std::swap(ofs[ip], ofs[i2]);
|
||||
|
||||
for( i1 = i0, i0--; i1 <= i2; i1++ )
|
||||
if( vals[ofs[i1]] <= pivot )
|
||||
{
|
||||
i0++;
|
||||
std::swap(ofs[i0], ofs[i1]);
|
||||
}
|
||||
if( i0 == middle )
|
||||
break;
|
||||
if( i0 > middle )
|
||||
b = i0 - (b == i0);
|
||||
else
|
||||
a = i0;
|
||||
}
|
||||
|
||||
float pivot = vals[ofs[middle]];
|
||||
int less = 0, more = 0;
|
||||
for( k = a0; k < middle; k++ )
|
||||
{
|
||||
CV_Assert(vals[ofs[k]] <= pivot);
|
||||
less += vals[ofs[k]] < pivot;
|
||||
}
|
||||
for( k = b0; k > middle; k-- )
|
||||
{
|
||||
CV_Assert(vals[ofs[k]] >= pivot);
|
||||
more += vals[ofs[k]] > pivot;
|
||||
}
|
||||
CV_Assert(std::abs(more - less) <= 1);
|
||||
|
||||
return vals[ofs[middle]];
|
||||
}
|
||||
|
||||
static void
|
||||
computeSums( const Mat& points, const size_t* ofs, int a, int b, double* sums )
|
||||
{
|
||||
int i, j, dims = points.cols;
|
||||
const float* data = points.ptr<float>(0);
|
||||
for( j = 0; j < dims; j++ )
|
||||
sums[j*2] = sums[j*2+1] = 0;
|
||||
for( i = a; i <= b; i++ )
|
||||
{
|
||||
const float* row = data + ofs[i];
|
||||
for( j = 0; j < dims; j++ )
|
||||
{
|
||||
double t = row[j], s = sums[j*2] + t, s2 = sums[j*2+1] + t*t;
|
||||
sums[j*2] = s; sums[j*2+1] = s2;
|
||||
}
|
||||
}
|
||||
}
|
||||
|
||||
|
||||
void KDTree::build(InputArray _points, bool _copyData)
|
||||
{
|
||||
build(_points, noArray(), _copyData);
|
||||
}
|
||||
|
||||
|
||||
void KDTree::build(InputArray __points, InputArray __labels, bool _copyData)
|
||||
{
|
||||
Mat _points = __points.getMat(), _labels = __labels.getMat();
|
||||
CV_Assert(_points.type() == CV_32F && !_points.empty());
|
||||
std::vector<KDTree::Node>().swap(nodes);
|
||||
|
||||
if( !_copyData )
|
||||
points = _points;
|
||||
else
|
||||
{
|
||||
points.release();
|
||||
points.create(_points.size(), _points.type());
|
||||
}
|
||||
|
||||
int i, j, n = _points.rows, ptdims = _points.cols, top = 0;
|
||||
const float* data = _points.ptr<float>(0);
|
||||
float* dstdata = points.ptr<float>(0);
|
||||
size_t step = _points.step1();
|
||||
size_t dstep = points.step1();
|
||||
int ptpos = 0;
|
||||
labels.resize(n);
|
||||
const int* _labels_data = 0;
|
||||
|
||||
if( !_labels.empty() )
|
||||
{
|
||||
int nlabels = _labels.checkVector(1, CV_32S, true);
|
||||
CV_Assert(nlabels == n);
|
||||
_labels_data = _labels.ptr<int>();
|
||||
}
|
||||
|
||||
Mat sumstack(MAX_TREE_DEPTH*2, ptdims*2, CV_64F);
|
||||
SubTree stack[MAX_TREE_DEPTH*2];
|
||||
|
||||
std::vector<size_t> _ptofs(n);
|
||||
size_t* ptofs = &_ptofs[0];
|
||||
|
||||
for( i = 0; i < n; i++ )
|
||||
ptofs[i] = i*step;
|
||||
|
||||
nodes.push_back(Node());
|
||||
computeSums(points, ptofs, 0, n-1, sumstack.ptr<double>(top));
|
||||
stack[top++] = SubTree(0, n-1, 0, 0);
|
||||
int _maxDepth = 0;
|
||||
|
||||
while( --top >= 0 )
|
||||
{
|
||||
int first = stack[top].first, last = stack[top].last;
|
||||
int depth = stack[top].depth, nidx = stack[top].nodeIdx;
|
||||
int count = last - first + 1, dim = -1;
|
||||
const double* sums = sumstack.ptr<double>(top);
|
||||
double invCount = 1./count, maxVar = -1.;
|
||||
|
||||
if( count == 1 )
|
||||
{
|
||||
int idx0 = (int)(ptofs[first]/step);
|
||||
int idx = _copyData ? ptpos++ : idx0;
|
||||
nodes[nidx].idx = ~idx;
|
||||
if( _copyData )
|
||||
{
|
||||
const float* src = data + ptofs[first];
|
||||
float* dst = dstdata + idx*dstep;
|
||||
for( j = 0; j < ptdims; j++ )
|
||||
dst[j] = src[j];
|
||||
}
|
||||
labels[idx] = _labels_data ? _labels_data[idx0] : idx0;
|
||||
_maxDepth = std::max(_maxDepth, depth);
|
||||
continue;
|
||||
}
|
||||
|
||||
// find the dimensionality with the biggest variance
|
||||
for( j = 0; j < ptdims; j++ )
|
||||
{
|
||||
double m = sums[j*2]*invCount;
|
||||
double varj = sums[j*2+1]*invCount - m*m;
|
||||
if( maxVar < varj )
|
||||
{
|
||||
maxVar = varj;
|
||||
dim = j;
|
||||
}
|
||||
}
|
||||
|
||||
int left = (int)nodes.size(), right = left + 1;
|
||||
nodes.push_back(Node());
|
||||
nodes.push_back(Node());
|
||||
nodes[nidx].idx = dim;
|
||||
nodes[nidx].left = left;
|
||||
nodes[nidx].right = right;
|
||||
nodes[nidx].boundary = medianPartition(ptofs, first, last, data + dim);
|
||||
|
||||
int middle = (first + last)/2;
|
||||
double *lsums = (double*)sums, *rsums = lsums + ptdims*2;
|
||||
computeSums(points, ptofs, middle+1, last, rsums);
|
||||
for( j = 0; j < ptdims*2; j++ )
|
||||
lsums[j] = sums[j] - rsums[j];
|
||||
stack[top++] = SubTree(first, middle, left, depth+1);
|
||||
stack[top++] = SubTree(middle+1, last, right, depth+1);
|
||||
}
|
||||
maxDepth = _maxDepth;
|
||||
}
|
||||
|
||||
|
||||
struct PQueueElem
|
||||
{
|
||||
PQueueElem() : dist(0), idx(0) {}
|
||||
PQueueElem(float _dist, int _idx) : dist(_dist), idx(_idx) {}
|
||||
float dist;
|
||||
int idx;
|
||||
};
|
||||
|
||||
|
||||
int KDTree::findNearest(InputArray _vec, int K, int emax,
|
||||
OutputArray _neighborsIdx, OutputArray _neighbors,
|
||||
OutputArray _dist, OutputArray _labels) const
|
||||
|
||||
{
|
||||
Mat vecmat = _vec.getMat();
|
||||
CV_Assert( vecmat.isContinuous() && vecmat.type() == CV_32F && vecmat.total() == (size_t)points.cols );
|
||||
const float* vec = vecmat.ptr<float>();
|
||||
K = std::min(K, points.rows);
|
||||
int ptdims = points.cols;
|
||||
|
||||
CV_Assert(K > 0 && (normType == NORM_L2 || normType == NORM_L1));
|
||||
|
||||
AutoBuffer<uchar> _buf((K+1)*(sizeof(float) + sizeof(int)));
|
||||
int* idx = (int*)(uchar*)_buf;
|
||||
float* dist = (float*)(idx + K + 1);
|
||||
int i, j, ncount = 0, e = 0;
|
||||
|
||||
int qsize = 0, maxqsize = 1 << 10;
|
||||
AutoBuffer<uchar> _pqueue(maxqsize*sizeof(PQueueElem));
|
||||
PQueueElem* pqueue = (PQueueElem*)(uchar*)_pqueue;
|
||||
emax = std::max(emax, 1);
|
||||
|
||||
for( e = 0; e < emax; )
|
||||
{
|
||||
float d, alt_d = 0.f;
|
||||
int nidx;
|
||||
|
||||
if( e == 0 )
|
||||
nidx = 0;
|
||||
else
|
||||
{
|
||||
// take the next node from the priority queue
|
||||
if( qsize == 0 )
|
||||
break;
|
||||
nidx = pqueue[0].idx;
|
||||
alt_d = pqueue[0].dist;
|
||||
if( --qsize > 0 )
|
||||
{
|
||||
std::swap(pqueue[0], pqueue[qsize]);
|
||||
d = pqueue[0].dist;
|
||||
for( i = 0;;)
|
||||
{
|
||||
int left = i*2 + 1, right = i*2 + 2;
|
||||
if( left >= qsize )
|
||||
break;
|
||||
if( right < qsize && pqueue[right].dist < pqueue[left].dist )
|
||||
left = right;
|
||||
if( pqueue[left].dist >= d )
|
||||
break;
|
||||
std::swap(pqueue[i], pqueue[left]);
|
||||
i = left;
|
||||
}
|
||||
}
|
||||
|
||||
if( ncount == K && alt_d > dist[ncount-1] )
|
||||
continue;
|
||||
}
|
||||
|
||||
for(;;)
|
||||
{
|
||||
if( nidx < 0 )
|
||||
break;
|
||||
const Node& n = nodes[nidx];
|
||||
|
||||
if( n.idx < 0 )
|
||||
{
|
||||
i = ~n.idx;
|
||||
const float* row = points.ptr<float>(i);
|
||||
if( normType == NORM_L2 )
|
||||
for( j = 0, d = 0.f; j < ptdims; j++ )
|
||||
{
|
||||
float t = vec[j] - row[j];
|
||||
d += t*t;
|
||||
}
|
||||
else
|
||||
for( j = 0, d = 0.f; j < ptdims; j++ )
|
||||
d += std::abs(vec[j] - row[j]);
|
||||
|
||||
dist[ncount] = d;
|
||||
idx[ncount] = i;
|
||||
for( i = ncount-1; i >= 0; i-- )
|
||||
{
|
||||
if( dist[i] <= d )
|
||||
break;
|
||||
std::swap(dist[i], dist[i+1]);
|
||||
std::swap(idx[i], idx[i+1]);
|
||||
}
|
||||
ncount += ncount < K;
|
||||
e++;
|
||||
break;
|
||||
}
|
||||
|
||||
int alt;
|
||||
if( vec[n.idx] <= n.boundary )
|
||||
{
|
||||
nidx = n.left;
|
||||
alt = n.right;
|
||||
}
|
||||
else
|
||||
{
|
||||
nidx = n.right;
|
||||
alt = n.left;
|
||||
}
|
||||
|
||||
d = vec[n.idx] - n.boundary;
|
||||
if( normType == NORM_L2 )
|
||||
d = d*d + alt_d;
|
||||
else
|
||||
d = std::abs(d) + alt_d;
|
||||
// subtree prunning
|
||||
if( ncount == K && d > dist[ncount-1] )
|
||||
continue;
|
||||
// add alternative subtree to the priority queue
|
||||
pqueue[qsize] = PQueueElem(d, alt);
|
||||
for( i = qsize; i > 0; )
|
||||
{
|
||||
int parent = (i-1)/2;
|
||||
if( parent < 0 || pqueue[parent].dist <= d )
|
||||
break;
|
||||
std::swap(pqueue[i], pqueue[parent]);
|
||||
i = parent;
|
||||
}
|
||||
qsize += qsize+1 < maxqsize;
|
||||
}
|
||||
}
|
||||
|
||||
K = std::min(K, ncount);
|
||||
if( _neighborsIdx.needed() )
|
||||
{
|
||||
_neighborsIdx.create(K, 1, CV_32S, -1, true);
|
||||
Mat nidx = _neighborsIdx.getMat();
|
||||
Mat(nidx.size(), CV_32S, &idx[0]).copyTo(nidx);
|
||||
}
|
||||
if( _dist.needed() )
|
||||
sqrt(Mat(K, 1, CV_32F, dist), _dist);
|
||||
|
||||
if( _neighbors.needed() || _labels.needed() )
|
||||
getPoints(Mat(K, 1, CV_32S, idx), _neighbors, _labels);
|
||||
return K;
|
||||
}
|
||||
|
||||
|
||||
void KDTree::findOrthoRange(InputArray _lowerBound,
|
||||
InputArray _upperBound,
|
||||
OutputArray _neighborsIdx,
|
||||
OutputArray _neighbors,
|
||||
OutputArray _labels ) const
|
||||
{
|
||||
int ptdims = points.cols;
|
||||
Mat lowerBound = _lowerBound.getMat(), upperBound = _upperBound.getMat();
|
||||
CV_Assert( lowerBound.size == upperBound.size &&
|
||||
lowerBound.isContinuous() &&
|
||||
upperBound.isContinuous() &&
|
||||
lowerBound.type() == upperBound.type() &&
|
||||
lowerBound.type() == CV_32F &&
|
||||
lowerBound.total() == (size_t)ptdims );
|
||||
const float* L = lowerBound.ptr<float>();
|
||||
const float* R = upperBound.ptr<float>();
|
||||
|
||||
std::vector<int> idx;
|
||||
AutoBuffer<int> _stack(MAX_TREE_DEPTH*2 + 1);
|
||||
int* stack = _stack;
|
||||
int top = 0;
|
||||
|
||||
stack[top++] = 0;
|
||||
|
||||
while( --top >= 0 )
|
||||
{
|
||||
int nidx = stack[top];
|
||||
if( nidx < 0 )
|
||||
break;
|
||||
const Node& n = nodes[nidx];
|
||||
if( n.idx < 0 )
|
||||
{
|
||||
int j, i = ~n.idx;
|
||||
const float* row = points.ptr<float>(i);
|
||||
for( j = 0; j < ptdims; j++ )
|
||||
if( row[j] < L[j] || row[j] >= R[j] )
|
||||
break;
|
||||
if( j == ptdims )
|
||||
idx.push_back(i);
|
||||
continue;
|
||||
}
|
||||
if( L[n.idx] <= n.boundary )
|
||||
stack[top++] = n.left;
|
||||
if( R[n.idx] > n.boundary )
|
||||
stack[top++] = n.right;
|
||||
}
|
||||
|
||||
if( _neighborsIdx.needed() )
|
||||
{
|
||||
_neighborsIdx.create((int)idx.size(), 1, CV_32S, -1, true);
|
||||
Mat nidx = _neighborsIdx.getMat();
|
||||
Mat(nidx.size(), CV_32S, &idx[0]).copyTo(nidx);
|
||||
}
|
||||
getPoints( idx, _neighbors, _labels );
|
||||
}
|
||||
|
||||
|
||||
void KDTree::getPoints(InputArray _idx, OutputArray _pts, OutputArray _labels) const
|
||||
{
|
||||
Mat idxmat = _idx.getMat(), pts, labelsmat;
|
||||
CV_Assert( idxmat.isContinuous() && idxmat.type() == CV_32S &&
|
||||
(idxmat.cols == 1 || idxmat.rows == 1) );
|
||||
const int* idx = idxmat.ptr<int>();
|
||||
int* dstlabels = 0;
|
||||
|
||||
int ptdims = points.cols;
|
||||
int i, nidx = (int)idxmat.total();
|
||||
if( nidx == 0 )
|
||||
{
|
||||
_pts.release();
|
||||
_labels.release();
|
||||
return;
|
||||
}
|
||||
|
||||
if( _pts.needed() )
|
||||
{
|
||||
_pts.create( nidx, ptdims, points.type());
|
||||
pts = _pts.getMat();
|
||||
}
|
||||
|
||||
if(_labels.needed())
|
||||
{
|
||||
_labels.create(nidx, 1, CV_32S, -1, true);
|
||||
labelsmat = _labels.getMat();
|
||||
CV_Assert( labelsmat.isContinuous() );
|
||||
dstlabels = labelsmat.ptr<int>();
|
||||
}
|
||||
const int* srclabels = !labels.empty() ? &labels[0] : 0;
|
||||
|
||||
for( i = 0; i < nidx; i++ )
|
||||
{
|
||||
int k = idx[i];
|
||||
CV_Assert( (unsigned)k < (unsigned)points.rows );
|
||||
const float* src = points.ptr<float>(k);
|
||||
if( !pts.empty() )
|
||||
std::copy(src, src + ptdims, pts.ptr<float>(i));
|
||||
if( dstlabels )
|
||||
dstlabels[i] = srclabels ? srclabels[k] : k;
|
||||
}
|
||||
}
|
||||
|
||||
|
||||
const float* KDTree::getPoint(int ptidx, int* label) const
|
||||
{
|
||||
CV_Assert( (unsigned)ptidx < (unsigned)points.rows);
|
||||
if(label)
|
||||
*label = labels[ptidx];
|
||||
return points.ptr<float>(ptidx);
|
||||
}
|
||||
|
||||
|
||||
int KDTree::dims() const
|
||||
{
|
||||
return !points.empty() ? points.cols : 0;
|
||||
}
|
||||
|
||||
////////////////////////////////////////////////////////////////////////////////
|
||||
|
||||
schar* seqPush( CvSeq* seq, const void* element )
|
||||
|
||||
@@ -0,0 +1,409 @@
|
||||
/*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.
|
||||
//
|
||||
//
|
||||
// License Agreement
|
||||
// For Open Source Computer Vision Library
|
||||
//
|
||||
// Copyright (C) 2013, OpenCV Foundation, 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 the copyright holders 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 OpenCV Foundation 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*/
|
||||
#include "precomp.hpp"
|
||||
|
||||
#define dprintf(x)
|
||||
#define print_matrix(x)
|
||||
|
||||
/*
|
||||
|
||||
****Error Message********************************************************************************************************************
|
||||
|
||||
Downhill Simplex method in OpenCV dev 3.0.0 getting this error:
|
||||
|
||||
OpenCV Error: Assertion failed (dims <= 2 && data && (unsigned)i0 < (unsigned)(s ize.p[0] * size.p[1])
|
||||
&& elemSize() == (((((DataType<_Tp>::type) & ((512 - 1) << 3)) >> 3) + 1) << ((((sizeof(size_t)/4+1)16384|0x3a50)
|
||||
>> ((DataType<_Tp>::typ e) & ((1 << 3) - 1))2) & 3))) in cv::Mat::at,
|
||||
file C:\builds\master_PackSlave-w in32-vc12-shared\opencv\modules\core\include\opencv2/core/mat.inl.hpp, line 893
|
||||
|
||||
****Problem and Possible Fix*********************************************************************************************************
|
||||
|
||||
DownhillSolverImpl::innerDownhillSimplex something looks broken here:
|
||||
Mat_<double> coord_sum(1,ndim,0.0),buf(1,ndim,0.0),y(1,ndim,0.0);
|
||||
nfunk = 0;
|
||||
for(i=0;i<ndim+1;++i)
|
||||
{
|
||||
y(i) = f->calc(p[i]);
|
||||
}
|
||||
|
||||
y has only ndim elements, while the loop goes over ndim+1
|
||||
|
||||
Edited the following for possible fix:
|
||||
|
||||
Replaced y(1,ndim,0.0) ------> y(1,ndim+1,0.0)
|
||||
|
||||
***********************************************************************************************************************************
|
||||
|
||||
The code below was used in tesing the source code.
|
||||
Created by @SareeAlnaghy
|
||||
|
||||
#include <iostream>
|
||||
#include <cstdlib>
|
||||
#include <cmath>
|
||||
#include <algorithm>
|
||||
#include <opencv2\optim\optim.hpp>
|
||||
|
||||
using namespace std;
|
||||
using namespace cv;
|
||||
|
||||
void test(Ptr<optim::DownhillSolver> MinProblemSolver, Ptr<optim::MinProblemSolver::Function> ptr_F, Mat &P, Mat &step)
|
||||
{
|
||||
try{
|
||||
|
||||
MinProblemSolver->setFunction(ptr_F);
|
||||
MinProblemSolver->setInitStep(step);
|
||||
double res = MinProblemSolver->minimize(P);
|
||||
|
||||
cout << "res " << res << endl;
|
||||
}
|
||||
catch (exception e)
|
||||
{
|
||||
cerr << "Error:: " << e.what() << endl;
|
||||
}
|
||||
}
|
||||
|
||||
int main()
|
||||
{
|
||||
|
||||
class DistanceToLines :public optim::MinProblemSolver::Function {
|
||||
public:
|
||||
double calc(const double* x)const{
|
||||
|
||||
return x[0] * x[0] + x[1] * x[1];
|
||||
|
||||
}
|
||||
};
|
||||
|
||||
Mat P = (Mat_<double>(1, 2) << 1.0, 1.0);
|
||||
Mat step = (Mat_<double>(2, 1) << -0.5, 0.5);
|
||||
|
||||
Ptr<optim::MinProblemSolver::Function> ptr_F(new DistanceToLines());
|
||||
Ptr<optim::DownhillSolver> MinProblemSolver = optim::createDownhillSolver();
|
||||
|
||||
test(MinProblemSolver, ptr_F, P, step);
|
||||
|
||||
system("pause");
|
||||
return 0;
|
||||
}
|
||||
|
||||
****Suggesttion for imporving Simplex implentation***************************************************************************************
|
||||
|
||||
Currently the downhilll simplex method outputs the function value that is minimized. It should also return the coordinate points where the
|
||||
function is minimized. This is very useful in many applications such as using back projection methods to find a point of intersection of
|
||||
multiple lines in three dimensions as not all lines intersect in three dimensions.
|
||||
|
||||
*/
|
||||
|
||||
namespace cv
|
||||
{
|
||||
|
||||
class DownhillSolverImpl : public DownhillSolver
|
||||
{
|
||||
public:
|
||||
void getInitStep(OutputArray step) const;
|
||||
void setInitStep(InputArray step);
|
||||
Ptr<Function> getFunction() const;
|
||||
void setFunction(const Ptr<Function>& f);
|
||||
TermCriteria getTermCriteria() const;
|
||||
DownhillSolverImpl();
|
||||
void setTermCriteria(const TermCriteria& termcrit);
|
||||
double minimize(InputOutputArray x);
|
||||
protected:
|
||||
Ptr<MinProblemSolver::Function> _Function;
|
||||
TermCriteria _termcrit;
|
||||
Mat _step;
|
||||
Mat_<double> buf_x;
|
||||
|
||||
private:
|
||||
inline void createInitialSimplex(Mat_<double>& simplex,Mat& step);
|
||||
inline double innerDownhillSimplex(cv::Mat_<double>& p,double MinRange,double MinError,int& nfunk,
|
||||
const Ptr<MinProblemSolver::Function>& f,int nmax);
|
||||
inline double tryNewPoint(Mat_<double>& p,Mat_<double>& y,Mat_<double>& coord_sum,const Ptr<MinProblemSolver::Function>& f,int ihi,
|
||||
double fac,Mat_<double>& ptry);
|
||||
};
|
||||
|
||||
double DownhillSolverImpl::tryNewPoint(
|
||||
Mat_<double>& p,
|
||||
Mat_<double>& y,
|
||||
Mat_<double>& coord_sum,
|
||||
const Ptr<MinProblemSolver::Function>& f,
|
||||
int ihi,
|
||||
double fac,
|
||||
Mat_<double>& ptry
|
||||
)
|
||||
{
|
||||
int ndim=p.cols;
|
||||
int j;
|
||||
double fac1,fac2,ytry;
|
||||
|
||||
fac1=(1.0-fac)/ndim;
|
||||
fac2=fac1-fac;
|
||||
for (j=0;j<ndim;j++)
|
||||
{
|
||||
ptry(j)=coord_sum(j)*fac1-p(ihi,j)*fac2;
|
||||
}
|
||||
ytry=f->calc(ptry.ptr<double>());
|
||||
if (ytry < y(ihi))
|
||||
{
|
||||
y(ihi)=ytry;
|
||||
for (j=0;j<ndim;j++)
|
||||
{
|
||||
coord_sum(j) += ptry(j)-p(ihi,j);
|
||||
p(ihi,j)=ptry(j);
|
||||
}
|
||||
}
|
||||
|
||||
return ytry;
|
||||
}
|
||||
|
||||
/*
|
||||
Performs the actual minimization of MinProblemSolver::Function f (after the initialization was done)
|
||||
|
||||
The matrix p[ndim+1][1..ndim] represents ndim+1 vertices that
|
||||
form a simplex - each row is an ndim vector.
|
||||
On output, nfunk gives the number of function evaluations taken.
|
||||
*/
|
||||
double DownhillSolverImpl::innerDownhillSimplex(
|
||||
cv::Mat_<double>& p,
|
||||
double MinRange,
|
||||
double MinError,
|
||||
int& nfunk,
|
||||
const Ptr<MinProblemSolver::Function>& f,
|
||||
int nmax
|
||||
)
|
||||
{
|
||||
int ndim=p.cols;
|
||||
double res;
|
||||
int i,ihi,ilo,inhi,j,mpts=ndim+1;
|
||||
double error, range,ysave,ytry;
|
||||
Mat_<double> coord_sum(1,ndim,0.0),buf(1,ndim,0.0),y(1,ndim+1,0.0);
|
||||
|
||||
nfunk = 0;
|
||||
|
||||
for(i=0;i<ndim+1;++i)
|
||||
{
|
||||
y(i) = f->calc(p[i]);
|
||||
}
|
||||
|
||||
nfunk = ndim+1;
|
||||
|
||||
reduce(p,coord_sum,0,CV_REDUCE_SUM);
|
||||
|
||||
for (;;)
|
||||
{
|
||||
ilo=0;
|
||||
/* find highest (worst), next-to-worst, and lowest
|
||||
(best) points by going through all of them. */
|
||||
ihi = y(0)>y(1) ? (inhi=1,0) : (inhi=0,1);
|
||||
for (i=0;i<mpts;i++)
|
||||
{
|
||||
if (y(i) <= y(ilo))
|
||||
ilo=i;
|
||||
if (y(i) > y(ihi))
|
||||
{
|
||||
inhi=ihi;
|
||||
ihi=i;
|
||||
}
|
||||
else if (y(i) > y(inhi) && i != ihi)
|
||||
inhi=i;
|
||||
}
|
||||
|
||||
/* check stop criterion */
|
||||
error=fabs(y(ihi)-y(ilo));
|
||||
range=0;
|
||||
for(i=0;i<ndim;++i)
|
||||
{
|
||||
double min = p(0,i);
|
||||
double max = p(0,i);
|
||||
double d;
|
||||
for(j=1;j<=ndim;++j)
|
||||
{
|
||||
if( min > p(j,i) ) min = p(j,i);
|
||||
if( max < p(j,i) ) max = p(j,i);
|
||||
}
|
||||
d = fabs(max-min);
|
||||
if(range < d) range = d;
|
||||
}
|
||||
|
||||
if(range <= MinRange || error <= MinError)
|
||||
{ /* Put best point and value in first slot. */
|
||||
std::swap(y(0),y(ilo));
|
||||
for (i=0;i<ndim;i++)
|
||||
{
|
||||
std::swap(p(0,i),p(ilo,i));
|
||||
}
|
||||
break;
|
||||
}
|
||||
|
||||
if (nfunk >= nmax){
|
||||
dprintf(("nmax exceeded\n"));
|
||||
return y(ilo);
|
||||
}
|
||||
nfunk += 2;
|
||||
/*Begin a new iteration. First, reflect the worst point about the centroid of others */
|
||||
ytry = tryNewPoint(p,y,coord_sum,f,ihi,-1.0,buf);
|
||||
if (ytry <= y(ilo))
|
||||
{ /*If that's better than the best point, go twice as far in that direction*/
|
||||
ytry = tryNewPoint(p,y,coord_sum,f,ihi,2.0,buf);
|
||||
}
|
||||
else if (ytry >= y(inhi))
|
||||
{ /* The new point is worse than the second-highest, but better
|
||||
than the worst so do not go so far in that direction */
|
||||
ysave = y(ihi);
|
||||
ytry = tryNewPoint(p,y,coord_sum,f,ihi,0.5,buf);
|
||||
if (ytry >= ysave)
|
||||
{ /* Can't seem to improve things. Contract the simplex to good point
|
||||
in hope to find a simplex landscape. */
|
||||
for (i=0;i<mpts;i++)
|
||||
{
|
||||
if (i != ilo)
|
||||
{
|
||||
for (j=0;j<ndim;j++)
|
||||
{
|
||||
p(i,j) = coord_sum(j) = 0.5*(p(i,j)+p(ilo,j));
|
||||
}
|
||||
y(i)=f->calc(coord_sum.ptr<double>());
|
||||
}
|
||||
}
|
||||
nfunk += ndim;
|
||||
reduce(p,coord_sum,0,CV_REDUCE_SUM);
|
||||
}
|
||||
} else --(nfunk); /* correct nfunk */
|
||||
dprintf(("this is simplex on iteration %d\n",nfunk));
|
||||
print_matrix(p);
|
||||
} /* go to next iteration. */
|
||||
res = y(0);
|
||||
|
||||
return res;
|
||||
}
|
||||
|
||||
void DownhillSolverImpl::createInitialSimplex(Mat_<double>& simplex,Mat& step){
|
||||
for(int i=1;i<=step.cols;++i)
|
||||
{
|
||||
simplex.row(0).copyTo(simplex.row(i));
|
||||
simplex(i,i-1)+= 0.5*step.at<double>(0,i-1);
|
||||
}
|
||||
simplex.row(0) -= 0.5*step;
|
||||
|
||||
dprintf(("this is simplex\n"));
|
||||
print_matrix(simplex);
|
||||
}
|
||||
|
||||
double DownhillSolverImpl::minimize(InputOutputArray x){
|
||||
dprintf(("hi from minimize\n"));
|
||||
CV_Assert(_Function.empty()==false);
|
||||
dprintf(("termcrit:\n\ttype: %d\n\tmaxCount: %d\n\tEPS: %g\n",_termcrit.type,_termcrit.maxCount,_termcrit.epsilon));
|
||||
dprintf(("step\n"));
|
||||
print_matrix(_step);
|
||||
|
||||
Mat x_mat=x.getMat();
|
||||
CV_Assert(MIN(x_mat.rows,x_mat.cols)==1);
|
||||
CV_Assert(MAX(x_mat.rows,x_mat.cols)==_step.cols);
|
||||
CV_Assert(x_mat.type()==CV_64FC1);
|
||||
|
||||
Mat_<double> proxy_x;
|
||||
|
||||
if(x_mat.rows>1){
|
||||
buf_x.create(1,_step.cols);
|
||||
Mat_<double> proxy(_step.cols,1,buf_x.ptr<double>());
|
||||
x_mat.copyTo(proxy);
|
||||
proxy_x=buf_x;
|
||||
}else{
|
||||
proxy_x=x_mat;
|
||||
}
|
||||
|
||||
int count=0;
|
||||
int ndim=_step.cols;
|
||||
Mat_<double> simplex=Mat_<double>(ndim+1,ndim,0.0);
|
||||
|
||||
simplex.row(0).copyTo(proxy_x);
|
||||
createInitialSimplex(simplex,_step);
|
||||
double res = innerDownhillSimplex(
|
||||
simplex,_termcrit.epsilon, _termcrit.epsilon, count,_Function,_termcrit.maxCount);
|
||||
simplex.row(0).copyTo(proxy_x);
|
||||
|
||||
dprintf(("%d iterations done\n",count));
|
||||
|
||||
if(x_mat.rows>1){
|
||||
Mat(x_mat.rows, 1, CV_64F, proxy_x.ptr<double>()).copyTo(x);
|
||||
}
|
||||
return res;
|
||||
}
|
||||
DownhillSolverImpl::DownhillSolverImpl(){
|
||||
_Function=Ptr<Function>();
|
||||
_step=Mat_<double>();
|
||||
}
|
||||
Ptr<MinProblemSolver::Function> DownhillSolverImpl::getFunction()const{
|
||||
return _Function;
|
||||
}
|
||||
void DownhillSolverImpl::setFunction(const Ptr<Function>& f){
|
||||
_Function=f;
|
||||
}
|
||||
TermCriteria DownhillSolverImpl::getTermCriteria()const{
|
||||
return _termcrit;
|
||||
}
|
||||
void DownhillSolverImpl::setTermCriteria(const TermCriteria& termcrit){
|
||||
CV_Assert(termcrit.type==(TermCriteria::MAX_ITER+TermCriteria::EPS) && termcrit.epsilon>0 && termcrit.maxCount>0);
|
||||
_termcrit=termcrit;
|
||||
}
|
||||
// both minRange & minError are specified by termcrit.epsilon; In addition, user may specify the number of iterations that the algorithm does.
|
||||
Ptr<DownhillSolver> DownhillSolver::create(const Ptr<MinProblemSolver::Function>& f, InputArray initStep, TermCriteria termcrit){
|
||||
Ptr<DownhillSolver> DS = makePtr<DownhillSolverImpl>();
|
||||
DS->setFunction(f);
|
||||
DS->setInitStep(initStep);
|
||||
DS->setTermCriteria(termcrit);
|
||||
return DS;
|
||||
}
|
||||
void DownhillSolverImpl::getInitStep(OutputArray step)const{
|
||||
_step.copyTo(step);
|
||||
}
|
||||
void DownhillSolverImpl::setInitStep(InputArray step){
|
||||
//set dimensionality and make a deep copy of step
|
||||
Mat m=step.getMat();
|
||||
dprintf(("m.cols=%d\nm.rows=%d\n",m.cols,m.rows));
|
||||
CV_Assert(MIN(m.cols,m.rows)==1 && m.type()==CV_64FC1);
|
||||
if(m.rows==1){
|
||||
m.copyTo(_step);
|
||||
}else{
|
||||
transpose(m,_step);
|
||||
}
|
||||
}
|
||||
}
|
||||
File diff suppressed because it is too large
Load Diff
@@ -0,0 +1,531 @@
|
||||
/*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.
|
||||
//
|
||||
//
|
||||
// License Agreement
|
||||
// For Open Source Computer Vision Library
|
||||
//
|
||||
// Copyright (C) 2000-2008, Intel Corporation, all rights reserved.
|
||||
// Copyright (C) 2009, Willow Garage Inc., all rights reserved.
|
||||
// Copyright (C) 2013, OpenCV Foundation, 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 the copyright holders 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*/
|
||||
|
||||
#include "precomp.hpp"
|
||||
|
||||
namespace cv
|
||||
{
|
||||
|
||||
// This is reimplementation of kd-trees from cvkdtree*.* by Xavier Delacour, cleaned-up and
|
||||
// adopted to work with the new OpenCV data structures. It's in cxcore to be shared by
|
||||
// both cv (CvFeatureTree) and ml (kNN).
|
||||
|
||||
// The algorithm is taken from:
|
||||
// J.S. Beis and D.G. Lowe. Shape indexing using approximate nearest-neighbor search
|
||||
// in highdimensional spaces. In Proc. IEEE Conf. Comp. Vision Patt. Recog.,
|
||||
// pages 1000--1006, 1997. http://citeseer.ist.psu.edu/beis97shape.html
|
||||
|
||||
const int MAX_TREE_DEPTH = 32;
|
||||
|
||||
KDTree::KDTree()
|
||||
{
|
||||
maxDepth = -1;
|
||||
normType = NORM_L2;
|
||||
}
|
||||
|
||||
KDTree::KDTree(InputArray _points, bool _copyData)
|
||||
{
|
||||
maxDepth = -1;
|
||||
normType = NORM_L2;
|
||||
build(_points, _copyData);
|
||||
}
|
||||
|
||||
KDTree::KDTree(InputArray _points, InputArray _labels, bool _copyData)
|
||||
{
|
||||
maxDepth = -1;
|
||||
normType = NORM_L2;
|
||||
build(_points, _labels, _copyData);
|
||||
}
|
||||
|
||||
struct SubTree
|
||||
{
|
||||
SubTree() : first(0), last(0), nodeIdx(0), depth(0) {}
|
||||
SubTree(int _first, int _last, int _nodeIdx, int _depth)
|
||||
: first(_first), last(_last), nodeIdx(_nodeIdx), depth(_depth) {}
|
||||
int first;
|
||||
int last;
|
||||
int nodeIdx;
|
||||
int depth;
|
||||
};
|
||||
|
||||
|
||||
static float
|
||||
medianPartition( size_t* ofs, int a, int b, const float* vals )
|
||||
{
|
||||
int k, a0 = a, b0 = b;
|
||||
int middle = (a + b)/2;
|
||||
while( b > a )
|
||||
{
|
||||
int i0 = a, i1 = (a+b)/2, i2 = b;
|
||||
float v0 = vals[ofs[i0]], v1 = vals[ofs[i1]], v2 = vals[ofs[i2]];
|
||||
int ip = v0 < v1 ? (v1 < v2 ? i1 : v0 < v2 ? i2 : i0) :
|
||||
v0 < v2 ? i0 : (v1 < v2 ? i2 : i1);
|
||||
float pivot = vals[ofs[ip]];
|
||||
std::swap(ofs[ip], ofs[i2]);
|
||||
|
||||
for( i1 = i0, i0--; i1 <= i2; i1++ )
|
||||
if( vals[ofs[i1]] <= pivot )
|
||||
{
|
||||
i0++;
|
||||
std::swap(ofs[i0], ofs[i1]);
|
||||
}
|
||||
if( i0 == middle )
|
||||
break;
|
||||
if( i0 > middle )
|
||||
b = i0 - (b == i0);
|
||||
else
|
||||
a = i0;
|
||||
}
|
||||
|
||||
float pivot = vals[ofs[middle]];
|
||||
int less = 0, more = 0;
|
||||
for( k = a0; k < middle; k++ )
|
||||
{
|
||||
CV_Assert(vals[ofs[k]] <= pivot);
|
||||
less += vals[ofs[k]] < pivot;
|
||||
}
|
||||
for( k = b0; k > middle; k-- )
|
||||
{
|
||||
CV_Assert(vals[ofs[k]] >= pivot);
|
||||
more += vals[ofs[k]] > pivot;
|
||||
}
|
||||
CV_Assert(std::abs(more - less) <= 1);
|
||||
|
||||
return vals[ofs[middle]];
|
||||
}
|
||||
|
||||
static void
|
||||
computeSums( const Mat& points, const size_t* ofs, int a, int b, double* sums )
|
||||
{
|
||||
int i, j, dims = points.cols;
|
||||
const float* data = points.ptr<float>(0);
|
||||
for( j = 0; j < dims; j++ )
|
||||
sums[j*2] = sums[j*2+1] = 0;
|
||||
for( i = a; i <= b; i++ )
|
||||
{
|
||||
const float* row = data + ofs[i];
|
||||
for( j = 0; j < dims; j++ )
|
||||
{
|
||||
double t = row[j], s = sums[j*2] + t, s2 = sums[j*2+1] + t*t;
|
||||
sums[j*2] = s; sums[j*2+1] = s2;
|
||||
}
|
||||
}
|
||||
}
|
||||
|
||||
|
||||
void KDTree::build(InputArray _points, bool _copyData)
|
||||
{
|
||||
build(_points, noArray(), _copyData);
|
||||
}
|
||||
|
||||
|
||||
void KDTree::build(InputArray __points, InputArray __labels, bool _copyData)
|
||||
{
|
||||
Mat _points = __points.getMat(), _labels = __labels.getMat();
|
||||
CV_Assert(_points.type() == CV_32F && !_points.empty());
|
||||
std::vector<KDTree::Node>().swap(nodes);
|
||||
|
||||
if( !_copyData )
|
||||
points = _points;
|
||||
else
|
||||
{
|
||||
points.release();
|
||||
points.create(_points.size(), _points.type());
|
||||
}
|
||||
|
||||
int i, j, n = _points.rows, ptdims = _points.cols, top = 0;
|
||||
const float* data = _points.ptr<float>(0);
|
||||
float* dstdata = points.ptr<float>(0);
|
||||
size_t step = _points.step1();
|
||||
size_t dstep = points.step1();
|
||||
int ptpos = 0;
|
||||
labels.resize(n);
|
||||
const int* _labels_data = 0;
|
||||
|
||||
if( !_labels.empty() )
|
||||
{
|
||||
int nlabels = _labels.checkVector(1, CV_32S, true);
|
||||
CV_Assert(nlabels == n);
|
||||
_labels_data = _labels.ptr<int>();
|
||||
}
|
||||
|
||||
Mat sumstack(MAX_TREE_DEPTH*2, ptdims*2, CV_64F);
|
||||
SubTree stack[MAX_TREE_DEPTH*2];
|
||||
|
||||
std::vector<size_t> _ptofs(n);
|
||||
size_t* ptofs = &_ptofs[0];
|
||||
|
||||
for( i = 0; i < n; i++ )
|
||||
ptofs[i] = i*step;
|
||||
|
||||
nodes.push_back(Node());
|
||||
computeSums(points, ptofs, 0, n-1, sumstack.ptr<double>(top));
|
||||
stack[top++] = SubTree(0, n-1, 0, 0);
|
||||
int _maxDepth = 0;
|
||||
|
||||
while( --top >= 0 )
|
||||
{
|
||||
int first = stack[top].first, last = stack[top].last;
|
||||
int depth = stack[top].depth, nidx = stack[top].nodeIdx;
|
||||
int count = last - first + 1, dim = -1;
|
||||
const double* sums = sumstack.ptr<double>(top);
|
||||
double invCount = 1./count, maxVar = -1.;
|
||||
|
||||
if( count == 1 )
|
||||
{
|
||||
int idx0 = (int)(ptofs[first]/step);
|
||||
int idx = _copyData ? ptpos++ : idx0;
|
||||
nodes[nidx].idx = ~idx;
|
||||
if( _copyData )
|
||||
{
|
||||
const float* src = data + ptofs[first];
|
||||
float* dst = dstdata + idx*dstep;
|
||||
for( j = 0; j < ptdims; j++ )
|
||||
dst[j] = src[j];
|
||||
}
|
||||
labels[idx] = _labels_data ? _labels_data[idx0] : idx0;
|
||||
_maxDepth = std::max(_maxDepth, depth);
|
||||
continue;
|
||||
}
|
||||
|
||||
// find the dimensionality with the biggest variance
|
||||
for( j = 0; j < ptdims; j++ )
|
||||
{
|
||||
double m = sums[j*2]*invCount;
|
||||
double varj = sums[j*2+1]*invCount - m*m;
|
||||
if( maxVar < varj )
|
||||
{
|
||||
maxVar = varj;
|
||||
dim = j;
|
||||
}
|
||||
}
|
||||
|
||||
int left = (int)nodes.size(), right = left + 1;
|
||||
nodes.push_back(Node());
|
||||
nodes.push_back(Node());
|
||||
nodes[nidx].idx = dim;
|
||||
nodes[nidx].left = left;
|
||||
nodes[nidx].right = right;
|
||||
nodes[nidx].boundary = medianPartition(ptofs, first, last, data + dim);
|
||||
|
||||
int middle = (first + last)/2;
|
||||
double *lsums = (double*)sums, *rsums = lsums + ptdims*2;
|
||||
computeSums(points, ptofs, middle+1, last, rsums);
|
||||
for( j = 0; j < ptdims*2; j++ )
|
||||
lsums[j] = sums[j] - rsums[j];
|
||||
stack[top++] = SubTree(first, middle, left, depth+1);
|
||||
stack[top++] = SubTree(middle+1, last, right, depth+1);
|
||||
}
|
||||
maxDepth = _maxDepth;
|
||||
}
|
||||
|
||||
|
||||
struct PQueueElem
|
||||
{
|
||||
PQueueElem() : dist(0), idx(0) {}
|
||||
PQueueElem(float _dist, int _idx) : dist(_dist), idx(_idx) {}
|
||||
float dist;
|
||||
int idx;
|
||||
};
|
||||
|
||||
|
||||
int KDTree::findNearest(InputArray _vec, int K, int emax,
|
||||
OutputArray _neighborsIdx, OutputArray _neighbors,
|
||||
OutputArray _dist, OutputArray _labels) const
|
||||
|
||||
{
|
||||
Mat vecmat = _vec.getMat();
|
||||
CV_Assert( vecmat.isContinuous() && vecmat.type() == CV_32F && vecmat.total() == (size_t)points.cols );
|
||||
const float* vec = vecmat.ptr<float>();
|
||||
K = std::min(K, points.rows);
|
||||
int ptdims = points.cols;
|
||||
|
||||
CV_Assert(K > 0 && (normType == NORM_L2 || normType == NORM_L1));
|
||||
|
||||
AutoBuffer<uchar> _buf((K+1)*(sizeof(float) + sizeof(int)));
|
||||
int* idx = (int*)(uchar*)_buf;
|
||||
float* dist = (float*)(idx + K + 1);
|
||||
int i, j, ncount = 0, e = 0;
|
||||
|
||||
int qsize = 0, maxqsize = 1 << 10;
|
||||
AutoBuffer<uchar> _pqueue(maxqsize*sizeof(PQueueElem));
|
||||
PQueueElem* pqueue = (PQueueElem*)(uchar*)_pqueue;
|
||||
emax = std::max(emax, 1);
|
||||
|
||||
for( e = 0; e < emax; )
|
||||
{
|
||||
float d, alt_d = 0.f;
|
||||
int nidx;
|
||||
|
||||
if( e == 0 )
|
||||
nidx = 0;
|
||||
else
|
||||
{
|
||||
// take the next node from the priority queue
|
||||
if( qsize == 0 )
|
||||
break;
|
||||
nidx = pqueue[0].idx;
|
||||
alt_d = pqueue[0].dist;
|
||||
if( --qsize > 0 )
|
||||
{
|
||||
std::swap(pqueue[0], pqueue[qsize]);
|
||||
d = pqueue[0].dist;
|
||||
for( i = 0;;)
|
||||
{
|
||||
int left = i*2 + 1, right = i*2 + 2;
|
||||
if( left >= qsize )
|
||||
break;
|
||||
if( right < qsize && pqueue[right].dist < pqueue[left].dist )
|
||||
left = right;
|
||||
if( pqueue[left].dist >= d )
|
||||
break;
|
||||
std::swap(pqueue[i], pqueue[left]);
|
||||
i = left;
|
||||
}
|
||||
}
|
||||
|
||||
if( ncount == K && alt_d > dist[ncount-1] )
|
||||
continue;
|
||||
}
|
||||
|
||||
for(;;)
|
||||
{
|
||||
if( nidx < 0 )
|
||||
break;
|
||||
const Node& n = nodes[nidx];
|
||||
|
||||
if( n.idx < 0 )
|
||||
{
|
||||
i = ~n.idx;
|
||||
const float* row = points.ptr<float>(i);
|
||||
if( normType == NORM_L2 )
|
||||
for( j = 0, d = 0.f; j < ptdims; j++ )
|
||||
{
|
||||
float t = vec[j] - row[j];
|
||||
d += t*t;
|
||||
}
|
||||
else
|
||||
for( j = 0, d = 0.f; j < ptdims; j++ )
|
||||
d += std::abs(vec[j] - row[j]);
|
||||
|
||||
dist[ncount] = d;
|
||||
idx[ncount] = i;
|
||||
for( i = ncount-1; i >= 0; i-- )
|
||||
{
|
||||
if( dist[i] <= d )
|
||||
break;
|
||||
std::swap(dist[i], dist[i+1]);
|
||||
std::swap(idx[i], idx[i+1]);
|
||||
}
|
||||
ncount += ncount < K;
|
||||
e++;
|
||||
break;
|
||||
}
|
||||
|
||||
int alt;
|
||||
if( vec[n.idx] <= n.boundary )
|
||||
{
|
||||
nidx = n.left;
|
||||
alt = n.right;
|
||||
}
|
||||
else
|
||||
{
|
||||
nidx = n.right;
|
||||
alt = n.left;
|
||||
}
|
||||
|
||||
d = vec[n.idx] - n.boundary;
|
||||
if( normType == NORM_L2 )
|
||||
d = d*d + alt_d;
|
||||
else
|
||||
d = std::abs(d) + alt_d;
|
||||
// subtree prunning
|
||||
if( ncount == K && d > dist[ncount-1] )
|
||||
continue;
|
||||
// add alternative subtree to the priority queue
|
||||
pqueue[qsize] = PQueueElem(d, alt);
|
||||
for( i = qsize; i > 0; )
|
||||
{
|
||||
int parent = (i-1)/2;
|
||||
if( parent < 0 || pqueue[parent].dist <= d )
|
||||
break;
|
||||
std::swap(pqueue[i], pqueue[parent]);
|
||||
i = parent;
|
||||
}
|
||||
qsize += qsize+1 < maxqsize;
|
||||
}
|
||||
}
|
||||
|
||||
K = std::min(K, ncount);
|
||||
if( _neighborsIdx.needed() )
|
||||
{
|
||||
_neighborsIdx.create(K, 1, CV_32S, -1, true);
|
||||
Mat nidx = _neighborsIdx.getMat();
|
||||
Mat(nidx.size(), CV_32S, &idx[0]).copyTo(nidx);
|
||||
}
|
||||
if( _dist.needed() )
|
||||
sqrt(Mat(K, 1, CV_32F, dist), _dist);
|
||||
|
||||
if( _neighbors.needed() || _labels.needed() )
|
||||
getPoints(Mat(K, 1, CV_32S, idx), _neighbors, _labels);
|
||||
return K;
|
||||
}
|
||||
|
||||
|
||||
void KDTree::findOrthoRange(InputArray _lowerBound,
|
||||
InputArray _upperBound,
|
||||
OutputArray _neighborsIdx,
|
||||
OutputArray _neighbors,
|
||||
OutputArray _labels ) const
|
||||
{
|
||||
int ptdims = points.cols;
|
||||
Mat lowerBound = _lowerBound.getMat(), upperBound = _upperBound.getMat();
|
||||
CV_Assert( lowerBound.size == upperBound.size &&
|
||||
lowerBound.isContinuous() &&
|
||||
upperBound.isContinuous() &&
|
||||
lowerBound.type() == upperBound.type() &&
|
||||
lowerBound.type() == CV_32F &&
|
||||
lowerBound.total() == (size_t)ptdims );
|
||||
const float* L = lowerBound.ptr<float>();
|
||||
const float* R = upperBound.ptr<float>();
|
||||
|
||||
std::vector<int> idx;
|
||||
AutoBuffer<int> _stack(MAX_TREE_DEPTH*2 + 1);
|
||||
int* stack = _stack;
|
||||
int top = 0;
|
||||
|
||||
stack[top++] = 0;
|
||||
|
||||
while( --top >= 0 )
|
||||
{
|
||||
int nidx = stack[top];
|
||||
if( nidx < 0 )
|
||||
break;
|
||||
const Node& n = nodes[nidx];
|
||||
if( n.idx < 0 )
|
||||
{
|
||||
int j, i = ~n.idx;
|
||||
const float* row = points.ptr<float>(i);
|
||||
for( j = 0; j < ptdims; j++ )
|
||||
if( row[j] < L[j] || row[j] >= R[j] )
|
||||
break;
|
||||
if( j == ptdims )
|
||||
idx.push_back(i);
|
||||
continue;
|
||||
}
|
||||
if( L[n.idx] <= n.boundary )
|
||||
stack[top++] = n.left;
|
||||
if( R[n.idx] > n.boundary )
|
||||
stack[top++] = n.right;
|
||||
}
|
||||
|
||||
if( _neighborsIdx.needed() )
|
||||
{
|
||||
_neighborsIdx.create((int)idx.size(), 1, CV_32S, -1, true);
|
||||
Mat nidx = _neighborsIdx.getMat();
|
||||
Mat(nidx.size(), CV_32S, &idx[0]).copyTo(nidx);
|
||||
}
|
||||
getPoints( idx, _neighbors, _labels );
|
||||
}
|
||||
|
||||
|
||||
void KDTree::getPoints(InputArray _idx, OutputArray _pts, OutputArray _labels) const
|
||||
{
|
||||
Mat idxmat = _idx.getMat(), pts, labelsmat;
|
||||
CV_Assert( idxmat.isContinuous() && idxmat.type() == CV_32S &&
|
||||
(idxmat.cols == 1 || idxmat.rows == 1) );
|
||||
const int* idx = idxmat.ptr<int>();
|
||||
int* dstlabels = 0;
|
||||
|
||||
int ptdims = points.cols;
|
||||
int i, nidx = (int)idxmat.total();
|
||||
if( nidx == 0 )
|
||||
{
|
||||
_pts.release();
|
||||
_labels.release();
|
||||
return;
|
||||
}
|
||||
|
||||
if( _pts.needed() )
|
||||
{
|
||||
_pts.create( nidx, ptdims, points.type());
|
||||
pts = _pts.getMat();
|
||||
}
|
||||
|
||||
if(_labels.needed())
|
||||
{
|
||||
_labels.create(nidx, 1, CV_32S, -1, true);
|
||||
labelsmat = _labels.getMat();
|
||||
CV_Assert( labelsmat.isContinuous() );
|
||||
dstlabels = labelsmat.ptr<int>();
|
||||
}
|
||||
const int* srclabels = !labels.empty() ? &labels[0] : 0;
|
||||
|
||||
for( i = 0; i < nidx; i++ )
|
||||
{
|
||||
int k = idx[i];
|
||||
CV_Assert( (unsigned)k < (unsigned)points.rows );
|
||||
const float* src = points.ptr<float>(k);
|
||||
if( !pts.empty() )
|
||||
std::copy(src, src + ptdims, pts.ptr<float>(i));
|
||||
if( dstlabels )
|
||||
dstlabels[i] = srclabels ? srclabels[k] : k;
|
||||
}
|
||||
}
|
||||
|
||||
|
||||
const float* KDTree::getPoint(int ptidx, int* label) const
|
||||
{
|
||||
CV_Assert( (unsigned)ptidx < (unsigned)points.rows);
|
||||
if(label)
|
||||
*label = labels[ptidx];
|
||||
return points.ptr<float>(ptidx);
|
||||
}
|
||||
|
||||
|
||||
int KDTree::dims() const
|
||||
{
|
||||
return !points.empty() ? points.cols : 0;
|
||||
}
|
||||
|
||||
}
|
||||
@@ -0,0 +1,457 @@
|
||||
/*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.
|
||||
//
|
||||
//
|
||||
// License Agreement
|
||||
// For Open Source Computer Vision Library
|
||||
//
|
||||
// Copyright (C) 2000-2008, Intel Corporation, all rights reserved.
|
||||
// Copyright (C) 2009, Willow Garage Inc., all rights reserved.
|
||||
// Copyright (C) 2013, OpenCV Foundation, 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 the copyright holders 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*/
|
||||
|
||||
#include "precomp.hpp"
|
||||
|
||||
////////////////////////////////////////// kmeans ////////////////////////////////////////////
|
||||
|
||||
namespace cv
|
||||
{
|
||||
|
||||
static void generateRandomCenter(const std::vector<Vec2f>& box, float* center, RNG& rng)
|
||||
{
|
||||
size_t j, dims = box.size();
|
||||
float margin = 1.f/dims;
|
||||
for( j = 0; j < dims; j++ )
|
||||
center[j] = ((float)rng*(1.f+margin*2.f)-margin)*(box[j][1] - box[j][0]) + box[j][0];
|
||||
}
|
||||
|
||||
class KMeansPPDistanceComputer : public ParallelLoopBody
|
||||
{
|
||||
public:
|
||||
KMeansPPDistanceComputer( float *_tdist2,
|
||||
const float *_data,
|
||||
const float *_dist,
|
||||
int _dims,
|
||||
size_t _step,
|
||||
size_t _stepci )
|
||||
: tdist2(_tdist2),
|
||||
data(_data),
|
||||
dist(_dist),
|
||||
dims(_dims),
|
||||
step(_step),
|
||||
stepci(_stepci) { }
|
||||
|
||||
void operator()( const cv::Range& range ) const
|
||||
{
|
||||
const int begin = range.start;
|
||||
const int end = range.end;
|
||||
|
||||
for ( int i = begin; i<end; i++ )
|
||||
{
|
||||
tdist2[i] = std::min(normL2Sqr_(data + step*i, data + stepci, dims), dist[i]);
|
||||
}
|
||||
}
|
||||
|
||||
private:
|
||||
KMeansPPDistanceComputer& operator=(const KMeansPPDistanceComputer&); // to quiet MSVC
|
||||
|
||||
float *tdist2;
|
||||
const float *data;
|
||||
const float *dist;
|
||||
const int dims;
|
||||
const size_t step;
|
||||
const size_t stepci;
|
||||
};
|
||||
|
||||
/*
|
||||
k-means center initialization using the following algorithm:
|
||||
Arthur & Vassilvitskii (2007) k-means++: The Advantages of Careful Seeding
|
||||
*/
|
||||
static void generateCentersPP(const Mat& _data, Mat& _out_centers,
|
||||
int K, RNG& rng, int trials)
|
||||
{
|
||||
int i, j, k, dims = _data.cols, N = _data.rows;
|
||||
const float* data = _data.ptr<float>(0);
|
||||
size_t step = _data.step/sizeof(data[0]);
|
||||
std::vector<int> _centers(K);
|
||||
int* centers = &_centers[0];
|
||||
std::vector<float> _dist(N*3);
|
||||
float* dist = &_dist[0], *tdist = dist + N, *tdist2 = tdist + N;
|
||||
double sum0 = 0;
|
||||
|
||||
centers[0] = (unsigned)rng % N;
|
||||
|
||||
for( i = 0; i < N; i++ )
|
||||
{
|
||||
dist[i] = normL2Sqr_(data + step*i, data + step*centers[0], dims);
|
||||
sum0 += dist[i];
|
||||
}
|
||||
|
||||
for( k = 1; k < K; k++ )
|
||||
{
|
||||
double bestSum = DBL_MAX;
|
||||
int bestCenter = -1;
|
||||
|
||||
for( j = 0; j < trials; j++ )
|
||||
{
|
||||
double p = (double)rng*sum0, s = 0;
|
||||
for( i = 0; i < N-1; i++ )
|
||||
if( (p -= dist[i]) <= 0 )
|
||||
break;
|
||||
int ci = i;
|
||||
|
||||
parallel_for_(Range(0, N),
|
||||
KMeansPPDistanceComputer(tdist2, data, dist, dims, step, step*ci));
|
||||
for( i = 0; i < N; i++ )
|
||||
{
|
||||
s += tdist2[i];
|
||||
}
|
||||
|
||||
if( s < bestSum )
|
||||
{
|
||||
bestSum = s;
|
||||
bestCenter = ci;
|
||||
std::swap(tdist, tdist2);
|
||||
}
|
||||
}
|
||||
centers[k] = bestCenter;
|
||||
sum0 = bestSum;
|
||||
std::swap(dist, tdist);
|
||||
}
|
||||
|
||||
for( k = 0; k < K; k++ )
|
||||
{
|
||||
const float* src = data + step*centers[k];
|
||||
float* dst = _out_centers.ptr<float>(k);
|
||||
for( j = 0; j < dims; j++ )
|
||||
dst[j] = src[j];
|
||||
}
|
||||
}
|
||||
|
||||
class KMeansDistanceComputer : public ParallelLoopBody
|
||||
{
|
||||
public:
|
||||
KMeansDistanceComputer( double *_distances,
|
||||
int *_labels,
|
||||
const Mat& _data,
|
||||
const Mat& _centers )
|
||||
: distances(_distances),
|
||||
labels(_labels),
|
||||
data(_data),
|
||||
centers(_centers)
|
||||
{
|
||||
}
|
||||
|
||||
void operator()( const Range& range ) const
|
||||
{
|
||||
const int begin = range.start;
|
||||
const int end = range.end;
|
||||
const int K = centers.rows;
|
||||
const int dims = centers.cols;
|
||||
|
||||
const float *sample;
|
||||
for( int i = begin; i<end; ++i)
|
||||
{
|
||||
sample = data.ptr<float>(i);
|
||||
int k_best = 0;
|
||||
double min_dist = DBL_MAX;
|
||||
|
||||
for( int k = 0; k < K; k++ )
|
||||
{
|
||||
const float* center = centers.ptr<float>(k);
|
||||
const double dist = normL2Sqr_(sample, center, dims);
|
||||
|
||||
if( min_dist > dist )
|
||||
{
|
||||
min_dist = dist;
|
||||
k_best = k;
|
||||
}
|
||||
}
|
||||
|
||||
distances[i] = min_dist;
|
||||
labels[i] = k_best;
|
||||
}
|
||||
}
|
||||
|
||||
private:
|
||||
KMeansDistanceComputer& operator=(const KMeansDistanceComputer&); // to quiet MSVC
|
||||
|
||||
double *distances;
|
||||
int *labels;
|
||||
const Mat& data;
|
||||
const Mat& centers;
|
||||
};
|
||||
|
||||
}
|
||||
|
||||
double cv::kmeans( InputArray _data, int K,
|
||||
InputOutputArray _bestLabels,
|
||||
TermCriteria criteria, int attempts,
|
||||
int flags, OutputArray _centers )
|
||||
{
|
||||
const int SPP_TRIALS = 3;
|
||||
Mat data0 = _data.getMat();
|
||||
bool isrow = data0.rows == 1 && data0.channels() > 1;
|
||||
int N = !isrow ? data0.rows : data0.cols;
|
||||
int dims = (!isrow ? data0.cols : 1)*data0.channels();
|
||||
int type = data0.depth();
|
||||
|
||||
attempts = std::max(attempts, 1);
|
||||
CV_Assert( data0.dims <= 2 && type == CV_32F && K > 0 );
|
||||
CV_Assert( N >= K );
|
||||
|
||||
Mat data(N, dims, CV_32F, data0.ptr(), isrow ? dims * sizeof(float) : static_cast<size_t>(data0.step));
|
||||
|
||||
_bestLabels.create(N, 1, CV_32S, -1, true);
|
||||
|
||||
Mat _labels, best_labels = _bestLabels.getMat();
|
||||
if( flags & CV_KMEANS_USE_INITIAL_LABELS )
|
||||
{
|
||||
CV_Assert( (best_labels.cols == 1 || best_labels.rows == 1) &&
|
||||
best_labels.cols*best_labels.rows == N &&
|
||||
best_labels.type() == CV_32S &&
|
||||
best_labels.isContinuous());
|
||||
best_labels.copyTo(_labels);
|
||||
}
|
||||
else
|
||||
{
|
||||
if( !((best_labels.cols == 1 || best_labels.rows == 1) &&
|
||||
best_labels.cols*best_labels.rows == N &&
|
||||
best_labels.type() == CV_32S &&
|
||||
best_labels.isContinuous()))
|
||||
best_labels.create(N, 1, CV_32S);
|
||||
_labels.create(best_labels.size(), best_labels.type());
|
||||
}
|
||||
int* labels = _labels.ptr<int>();
|
||||
|
||||
Mat centers(K, dims, type), old_centers(K, dims, type), temp(1, dims, type);
|
||||
std::vector<int> counters(K);
|
||||
std::vector<Vec2f> _box(dims);
|
||||
Vec2f* box = &_box[0];
|
||||
double best_compactness = DBL_MAX, compactness = 0;
|
||||
RNG& rng = theRNG();
|
||||
int a, iter, i, j, k;
|
||||
|
||||
if( criteria.type & TermCriteria::EPS )
|
||||
criteria.epsilon = std::max(criteria.epsilon, 0.);
|
||||
else
|
||||
criteria.epsilon = FLT_EPSILON;
|
||||
criteria.epsilon *= criteria.epsilon;
|
||||
|
||||
if( criteria.type & TermCriteria::COUNT )
|
||||
criteria.maxCount = std::min(std::max(criteria.maxCount, 2), 100);
|
||||
else
|
||||
criteria.maxCount = 100;
|
||||
|
||||
if( K == 1 )
|
||||
{
|
||||
attempts = 1;
|
||||
criteria.maxCount = 2;
|
||||
}
|
||||
|
||||
const float* sample = data.ptr<float>(0);
|
||||
for( j = 0; j < dims; j++ )
|
||||
box[j] = Vec2f(sample[j], sample[j]);
|
||||
|
||||
for( i = 1; i < N; i++ )
|
||||
{
|
||||
sample = data.ptr<float>(i);
|
||||
for( j = 0; j < dims; j++ )
|
||||
{
|
||||
float v = sample[j];
|
||||
box[j][0] = std::min(box[j][0], v);
|
||||
box[j][1] = std::max(box[j][1], v);
|
||||
}
|
||||
}
|
||||
|
||||
for( a = 0; a < attempts; a++ )
|
||||
{
|
||||
double max_center_shift = DBL_MAX;
|
||||
for( iter = 0;; )
|
||||
{
|
||||
swap(centers, old_centers);
|
||||
|
||||
if( iter == 0 && (a > 0 || !(flags & KMEANS_USE_INITIAL_LABELS)) )
|
||||
{
|
||||
if( flags & KMEANS_PP_CENTERS )
|
||||
generateCentersPP(data, centers, K, rng, SPP_TRIALS);
|
||||
else
|
||||
{
|
||||
for( k = 0; k < K; k++ )
|
||||
generateRandomCenter(_box, centers.ptr<float>(k), rng);
|
||||
}
|
||||
}
|
||||
else
|
||||
{
|
||||
if( iter == 0 && a == 0 && (flags & KMEANS_USE_INITIAL_LABELS) )
|
||||
{
|
||||
for( i = 0; i < N; i++ )
|
||||
CV_Assert( (unsigned)labels[i] < (unsigned)K );
|
||||
}
|
||||
|
||||
// compute centers
|
||||
centers = Scalar(0);
|
||||
for( k = 0; k < K; k++ )
|
||||
counters[k] = 0;
|
||||
|
||||
for( i = 0; i < N; i++ )
|
||||
{
|
||||
sample = data.ptr<float>(i);
|
||||
k = labels[i];
|
||||
float* center = centers.ptr<float>(k);
|
||||
j=0;
|
||||
#if CV_ENABLE_UNROLLED
|
||||
for(; j <= dims - 4; j += 4 )
|
||||
{
|
||||
float t0 = center[j] + sample[j];
|
||||
float t1 = center[j+1] + sample[j+1];
|
||||
|
||||
center[j] = t0;
|
||||
center[j+1] = t1;
|
||||
|
||||
t0 = center[j+2] + sample[j+2];
|
||||
t1 = center[j+3] + sample[j+3];
|
||||
|
||||
center[j+2] = t0;
|
||||
center[j+3] = t1;
|
||||
}
|
||||
#endif
|
||||
for( ; j < dims; j++ )
|
||||
center[j] += sample[j];
|
||||
counters[k]++;
|
||||
}
|
||||
|
||||
if( iter > 0 )
|
||||
max_center_shift = 0;
|
||||
|
||||
for( k = 0; k < K; k++ )
|
||||
{
|
||||
if( counters[k] != 0 )
|
||||
continue;
|
||||
|
||||
// if some cluster appeared to be empty then:
|
||||
// 1. find the biggest cluster
|
||||
// 2. find the farthest from the center point in the biggest cluster
|
||||
// 3. exclude the farthest point from the biggest cluster and form a new 1-point cluster.
|
||||
int max_k = 0;
|
||||
for( int k1 = 1; k1 < K; k1++ )
|
||||
{
|
||||
if( counters[max_k] < counters[k1] )
|
||||
max_k = k1;
|
||||
}
|
||||
|
||||
double max_dist = 0;
|
||||
int farthest_i = -1;
|
||||
float* new_center = centers.ptr<float>(k);
|
||||
float* old_center = centers.ptr<float>(max_k);
|
||||
float* _old_center = temp.ptr<float>(); // normalized
|
||||
float scale = 1.f/counters[max_k];
|
||||
for( j = 0; j < dims; j++ )
|
||||
_old_center[j] = old_center[j]*scale;
|
||||
|
||||
for( i = 0; i < N; i++ )
|
||||
{
|
||||
if( labels[i] != max_k )
|
||||
continue;
|
||||
sample = data.ptr<float>(i);
|
||||
double dist = normL2Sqr_(sample, _old_center, dims);
|
||||
|
||||
if( max_dist <= dist )
|
||||
{
|
||||
max_dist = dist;
|
||||
farthest_i = i;
|
||||
}
|
||||
}
|
||||
|
||||
counters[max_k]--;
|
||||
counters[k]++;
|
||||
labels[farthest_i] = k;
|
||||
sample = data.ptr<float>(farthest_i);
|
||||
|
||||
for( j = 0; j < dims; j++ )
|
||||
{
|
||||
old_center[j] -= sample[j];
|
||||
new_center[j] += sample[j];
|
||||
}
|
||||
}
|
||||
|
||||
for( k = 0; k < K; k++ )
|
||||
{
|
||||
float* center = centers.ptr<float>(k);
|
||||
CV_Assert( counters[k] != 0 );
|
||||
|
||||
float scale = 1.f/counters[k];
|
||||
for( j = 0; j < dims; j++ )
|
||||
center[j] *= scale;
|
||||
|
||||
if( iter > 0 )
|
||||
{
|
||||
double dist = 0;
|
||||
const float* old_center = old_centers.ptr<float>(k);
|
||||
for( j = 0; j < dims; j++ )
|
||||
{
|
||||
double t = center[j] - old_center[j];
|
||||
dist += t*t;
|
||||
}
|
||||
max_center_shift = std::max(max_center_shift, dist);
|
||||
}
|
||||
}
|
||||
}
|
||||
|
||||
if( ++iter == MAX(criteria.maxCount, 2) || max_center_shift <= criteria.epsilon )
|
||||
break;
|
||||
|
||||
// assign labels
|
||||
Mat dists(1, N, CV_64F);
|
||||
double* dist = dists.ptr<double>(0);
|
||||
parallel_for_(Range(0, N),
|
||||
KMeansDistanceComputer(dist, labels, data, centers));
|
||||
compactness = 0;
|
||||
for( i = 0; i < N; i++ )
|
||||
{
|
||||
compactness += dist[i];
|
||||
}
|
||||
}
|
||||
|
||||
if( compactness < best_compactness )
|
||||
{
|
||||
best_compactness = compactness;
|
||||
if( _centers.needed() )
|
||||
centers.copyTo(_centers);
|
||||
_labels.copyTo(best_labels);
|
||||
}
|
||||
}
|
||||
|
||||
return best_compactness;
|
||||
}
|
||||
@@ -0,0 +1,362 @@
|
||||
/*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.
|
||||
//
|
||||
//
|
||||
// License Agreement
|
||||
// For Open Source Computer Vision Library
|
||||
//
|
||||
// Copyright (C) 2013, OpenCV Foundation, 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 the copyright holders 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 OpenCV Foundation 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*/
|
||||
#include "precomp.hpp"
|
||||
#include <climits>
|
||||
#include <algorithm>
|
||||
#include <cstdarg>
|
||||
|
||||
#define dprintf(x)
|
||||
#define print_matrix(x)
|
||||
|
||||
namespace cv
|
||||
{
|
||||
|
||||
using std::vector;
|
||||
|
||||
#ifdef ALEX_DEBUG
|
||||
static void print_simplex_state(const Mat& c,const Mat& b,double v,const std::vector<int> N,const std::vector<int> B){
|
||||
printf("\tprint simplex state\n");
|
||||
|
||||
printf("v=%g\n",v);
|
||||
|
||||
printf("here c goes\n");
|
||||
print_matrix(c);
|
||||
|
||||
printf("non-basic: ");
|
||||
print(Mat(N));
|
||||
printf("\n");
|
||||
|
||||
printf("here b goes\n");
|
||||
print_matrix(b);
|
||||
printf("basic: ");
|
||||
|
||||
print(Mat(B));
|
||||
printf("\n");
|
||||
}
|
||||
#else
|
||||
#define print_simplex_state(c,b,v,N,B)
|
||||
#endif
|
||||
|
||||
/**Due to technical considerations, the format of input b and c is somewhat special:
|
||||
*both b and c should be one column bigger than corresponding b and c of linear problem and the leftmost column will be used internally
|
||||
by this procedure - it should not be cleaned before the call to procedure and may contain mess after
|
||||
it also initializes N and B and does not make any assumptions about their init values
|
||||
* @return SOLVELP_UNFEASIBLE if problem is unfeasible, 0 if feasible.
|
||||
*/
|
||||
static int initialize_simplex(Mat_<double>& c, Mat_<double>& b,double& v,vector<int>& N,vector<int>& B,vector<unsigned int>& indexToRow);
|
||||
static inline void pivot(Mat_<double>& c,Mat_<double>& b,double& v,vector<int>& N,vector<int>& B,int leaving_index,
|
||||
int entering_index,vector<unsigned int>& indexToRow);
|
||||
/**@return SOLVELP_UNBOUNDED means the problem is unbdd, SOLVELP_MULTI means multiple solutions, SOLVELP_SINGLE means one solution.
|
||||
*/
|
||||
static int inner_simplex(Mat_<double>& c, Mat_<double>& b,double& v,vector<int>& N,vector<int>& B,vector<unsigned int>& indexToRow);
|
||||
static void swap_columns(Mat_<double>& A,int col1,int col2);
|
||||
#define SWAP(type,a,b) {type tmp=(a);(a)=(b);(b)=tmp;}
|
||||
|
||||
//return codes:-2 (no_sol - unbdd),-1(no_sol - unfsbl), 0(single_sol), 1(multiple_sol=>least_l2_norm)
|
||||
int solveLP(const Mat& Func, const Mat& Constr, Mat& z){
|
||||
dprintf(("call to solveLP\n"));
|
||||
|
||||
//sanity check (size, type, no. of channels)
|
||||
CV_Assert(Func.type()==CV_64FC1 || Func.type()==CV_32FC1);
|
||||
CV_Assert(Constr.type()==CV_64FC1 || Constr.type()==CV_32FC1);
|
||||
CV_Assert((Func.rows==1 && (Constr.cols-Func.cols==1))||
|
||||
(Func.cols==1 && (Constr.cols-Func.rows==1)));
|
||||
|
||||
//copy arguments for we will shall modify them
|
||||
Mat_<double> bigC=Mat_<double>(1,(Func.rows==1?Func.cols:Func.rows)+1),
|
||||
bigB=Mat_<double>(Constr.rows,Constr.cols+1);
|
||||
if(Func.rows==1){
|
||||
Func.convertTo(bigC.colRange(1,bigC.cols),CV_64FC1);
|
||||
}else{
|
||||
Mat FuncT=Func.t();
|
||||
FuncT.convertTo(bigC.colRange(1,bigC.cols),CV_64FC1);
|
||||
}
|
||||
Constr.convertTo(bigB.colRange(1,bigB.cols),CV_64FC1);
|
||||
double v=0;
|
||||
vector<int> N,B;
|
||||
vector<unsigned int> indexToRow;
|
||||
|
||||
if(initialize_simplex(bigC,bigB,v,N,B,indexToRow)==SOLVELP_UNFEASIBLE){
|
||||
return SOLVELP_UNFEASIBLE;
|
||||
}
|
||||
Mat_<double> c=bigC.colRange(1,bigC.cols),
|
||||
b=bigB.colRange(1,bigB.cols);
|
||||
|
||||
int res=0;
|
||||
if((res=inner_simplex(c,b,v,N,B,indexToRow))==SOLVELP_UNBOUNDED){
|
||||
return SOLVELP_UNBOUNDED;
|
||||
}
|
||||
|
||||
//return the optimal solution
|
||||
z.create(c.cols,1,CV_64FC1);
|
||||
MatIterator_<double> it=z.begin<double>();
|
||||
unsigned int nsize = (unsigned int)N.size();
|
||||
for(int i=1;i<=c.cols;i++,it++){
|
||||
if(indexToRow[i]<nsize){
|
||||
*it=0;
|
||||
}else{
|
||||
*it=b.at<double>(indexToRow[i]-nsize,b.cols-1);
|
||||
}
|
||||
}
|
||||
|
||||
return res;
|
||||
}
|
||||
|
||||
static int initialize_simplex(Mat_<double>& c, Mat_<double>& b,double& v,vector<int>& N,vector<int>& B,vector<unsigned int>& indexToRow){
|
||||
N.resize(c.cols);
|
||||
N[0]=0;
|
||||
for (std::vector<int>::iterator it = N.begin()+1 ; it != N.end(); ++it){
|
||||
*it=it[-1]+1;
|
||||
}
|
||||
B.resize(b.rows);
|
||||
B[0]=(int)N.size();
|
||||
for (std::vector<int>::iterator it = B.begin()+1 ; it != B.end(); ++it){
|
||||
*it=it[-1]+1;
|
||||
}
|
||||
indexToRow.resize(c.cols+b.rows);
|
||||
indexToRow[0]=0;
|
||||
for (std::vector<unsigned int>::iterator it = indexToRow.begin()+1 ; it != indexToRow.end(); ++it){
|
||||
*it=it[-1]+1;
|
||||
}
|
||||
v=0;
|
||||
|
||||
int k=0;
|
||||
{
|
||||
double min=DBL_MAX;
|
||||
for(int i=0;i<b.rows;i++){
|
||||
if(b(i,b.cols-1)<min){
|
||||
min=b(i,b.cols-1);
|
||||
k=i;
|
||||
}
|
||||
}
|
||||
}
|
||||
|
||||
if(b(k,b.cols-1)>=0){
|
||||
N.erase(N.begin());
|
||||
for (std::vector<unsigned int>::iterator it = indexToRow.begin()+1 ; it != indexToRow.end(); ++it){
|
||||
--(*it);
|
||||
}
|
||||
return 0;
|
||||
}
|
||||
|
||||
Mat_<double> old_c=c.clone();
|
||||
c=0;
|
||||
c(0,0)=-1;
|
||||
for(int i=0;i<b.rows;i++){
|
||||
b(i,0)=-1;
|
||||
}
|
||||
|
||||
print_simplex_state(c,b,v,N,B);
|
||||
|
||||
dprintf(("\tWE MAKE PIVOT\n"));
|
||||
pivot(c,b,v,N,B,k,0,indexToRow);
|
||||
|
||||
print_simplex_state(c,b,v,N,B);
|
||||
|
||||
inner_simplex(c,b,v,N,B,indexToRow);
|
||||
|
||||
dprintf(("\tAFTER INNER_SIMPLEX\n"));
|
||||
print_simplex_state(c,b,v,N,B);
|
||||
|
||||
unsigned int nsize = (unsigned int)N.size();
|
||||
if(indexToRow[0]>=nsize){
|
||||
int iterator_offset=indexToRow[0]-nsize;
|
||||
if(b(iterator_offset,b.cols-1)>0){
|
||||
return SOLVELP_UNFEASIBLE;
|
||||
}
|
||||
pivot(c,b,v,N,B,iterator_offset,0,indexToRow);
|
||||
}
|
||||
|
||||
vector<int>::iterator iterator;
|
||||
{
|
||||
int iterator_offset=indexToRow[0];
|
||||
iterator=N.begin()+iterator_offset;
|
||||
std::iter_swap(iterator,N.begin());
|
||||
SWAP(int,indexToRow[*iterator],indexToRow[0]);
|
||||
swap_columns(c,iterator_offset,0);
|
||||
swap_columns(b,iterator_offset,0);
|
||||
}
|
||||
|
||||
dprintf(("after swaps\n"));
|
||||
print_simplex_state(c,b,v,N,B);
|
||||
|
||||
//start from 1, because we ignore x_0
|
||||
c=0;
|
||||
v=0;
|
||||
for(int I=1;I<old_c.cols;I++){
|
||||
if(indexToRow[I]<nsize){
|
||||
dprintf(("I=%d from nonbasic\n",I));
|
||||
int iterator_offset=indexToRow[I];
|
||||
c(0,iterator_offset)+=old_c(0,I);
|
||||
print_matrix(c);
|
||||
}else{
|
||||
dprintf(("I=%d from basic\n",I));
|
||||
int iterator_offset=indexToRow[I]-nsize;
|
||||
c-=old_c(0,I)*b.row(iterator_offset).colRange(0,b.cols-1);
|
||||
v+=old_c(0,I)*b(iterator_offset,b.cols-1);
|
||||
print_matrix(c);
|
||||
}
|
||||
}
|
||||
|
||||
dprintf(("after restore\n"));
|
||||
print_simplex_state(c,b,v,N,B);
|
||||
|
||||
N.erase(N.begin());
|
||||
for (std::vector<unsigned int>::iterator it = indexToRow.begin()+1 ; it != indexToRow.end(); ++it){
|
||||
--(*it);
|
||||
}
|
||||
return 0;
|
||||
}
|
||||
|
||||
static int inner_simplex(Mat_<double>& c, Mat_<double>& b,double& v,vector<int>& N,vector<int>& B,vector<unsigned int>& indexToRow){
|
||||
int count=0;
|
||||
for(;;){
|
||||
dprintf(("iteration #%d\n",count));
|
||||
count++;
|
||||
|
||||
static MatIterator_<double> pos_ptr;
|
||||
int e=-1,pos_ctr=0,min_var=INT_MAX;
|
||||
bool all_nonzero=true;
|
||||
for(pos_ptr=c.begin();pos_ptr!=c.end();pos_ptr++,pos_ctr++){
|
||||
if(*pos_ptr==0){
|
||||
all_nonzero=false;
|
||||
}
|
||||
if(*pos_ptr>0){
|
||||
if(N[pos_ctr]<min_var){
|
||||
e=pos_ctr;
|
||||
min_var=N[pos_ctr];
|
||||
}
|
||||
}
|
||||
}
|
||||
if(e==-1){
|
||||
dprintf(("hello from e==-1\n"));
|
||||
print_matrix(c);
|
||||
if(all_nonzero==true){
|
||||
return SOLVELP_SINGLE;
|
||||
}else{
|
||||
return SOLVELP_MULTI;
|
||||
}
|
||||
}
|
||||
|
||||
int l=-1;
|
||||
min_var=INT_MAX;
|
||||
double min=DBL_MAX;
|
||||
int row_it=0;
|
||||
MatIterator_<double> min_row_ptr=b.begin();
|
||||
for(MatIterator_<double> it=b.begin();it!=b.end();it+=b.cols,row_it++){
|
||||
double myite=0;
|
||||
//check constraints, select the tightest one, reinforcing Bland's rule
|
||||
if((myite=it[e])>0){
|
||||
double val=it[b.cols-1]/myite;
|
||||
if(val<min || (val==min && B[row_it]<min_var)){
|
||||
min_var=B[row_it];
|
||||
min_row_ptr=it;
|
||||
min=val;
|
||||
l=row_it;
|
||||
}
|
||||
}
|
||||
}
|
||||
if(l==-1){
|
||||
return SOLVELP_UNBOUNDED;
|
||||
}
|
||||
dprintf(("the tightest constraint is in row %d with %g\n",l,min));
|
||||
|
||||
pivot(c,b,v,N,B,l,e,indexToRow);
|
||||
|
||||
dprintf(("objective, v=%g\n",v));
|
||||
print_matrix(c);
|
||||
dprintf(("constraints\n"));
|
||||
print_matrix(b);
|
||||
dprintf(("non-basic: "));
|
||||
print_matrix(Mat(N));
|
||||
dprintf(("basic: "));
|
||||
print_matrix(Mat(B));
|
||||
}
|
||||
}
|
||||
|
||||
static inline void pivot(Mat_<double>& c,Mat_<double>& b,double& v,vector<int>& N,vector<int>& B,
|
||||
int leaving_index,int entering_index,vector<unsigned int>& indexToRow){
|
||||
double Coef=b(leaving_index,entering_index);
|
||||
for(int i=0;i<b.cols;i++){
|
||||
if(i==entering_index){
|
||||
b(leaving_index,i)=1/Coef;
|
||||
}else{
|
||||
b(leaving_index,i)/=Coef;
|
||||
}
|
||||
}
|
||||
|
||||
for(int i=0;i<b.rows;i++){
|
||||
if(i!=leaving_index){
|
||||
double coef=b(i,entering_index);
|
||||
for(int j=0;j<b.cols;j++){
|
||||
if(j==entering_index){
|
||||
b(i,j)=-coef*b(leaving_index,j);
|
||||
}else{
|
||||
b(i,j)-=(coef*b(leaving_index,j));
|
||||
}
|
||||
}
|
||||
}
|
||||
}
|
||||
|
||||
//objective function
|
||||
Coef=c(0,entering_index);
|
||||
for(int i=0;i<(b.cols-1);i++){
|
||||
if(i==entering_index){
|
||||
c(0,i)=-Coef*b(leaving_index,i);
|
||||
}else{
|
||||
c(0,i)-=Coef*b(leaving_index,i);
|
||||
}
|
||||
}
|
||||
dprintf(("v was %g\n",v));
|
||||
v+=Coef*b(leaving_index,b.cols-1);
|
||||
|
||||
SWAP(int,N[entering_index],B[leaving_index]);
|
||||
SWAP(int,indexToRow[N[entering_index]],indexToRow[B[leaving_index]]);
|
||||
}
|
||||
|
||||
static inline void swap_columns(Mat_<double>& A,int col1,int col2){
|
||||
for(int i=0;i<A.rows;i++){
|
||||
double tmp=A(i,col1);
|
||||
A(i,col1)=A(i,col2);
|
||||
A(i,col2)=tmp;
|
||||
}
|
||||
}
|
||||
}
|
||||
@@ -2959,341 +2959,6 @@ double Mat::dot(InputArray _mat) const
|
||||
return r;
|
||||
}
|
||||
|
||||
/****************************************************************************************\
|
||||
* PCA *
|
||||
\****************************************************************************************/
|
||||
|
||||
PCA::PCA() {}
|
||||
|
||||
PCA::PCA(InputArray data, InputArray _mean, int flags, int maxComponents)
|
||||
{
|
||||
operator()(data, _mean, flags, maxComponents);
|
||||
}
|
||||
|
||||
PCA::PCA(InputArray data, InputArray _mean, int flags, double retainedVariance)
|
||||
{
|
||||
operator()(data, _mean, flags, retainedVariance);
|
||||
}
|
||||
|
||||
PCA& PCA::operator()(InputArray _data, InputArray __mean, int flags, int maxComponents)
|
||||
{
|
||||
Mat data = _data.getMat(), _mean = __mean.getMat();
|
||||
int covar_flags = CV_COVAR_SCALE;
|
||||
int i, len, in_count;
|
||||
Size mean_sz;
|
||||
|
||||
CV_Assert( data.channels() == 1 );
|
||||
if( flags & CV_PCA_DATA_AS_COL )
|
||||
{
|
||||
len = data.rows;
|
||||
in_count = data.cols;
|
||||
covar_flags |= CV_COVAR_COLS;
|
||||
mean_sz = Size(1, len);
|
||||
}
|
||||
else
|
||||
{
|
||||
len = data.cols;
|
||||
in_count = data.rows;
|
||||
covar_flags |= CV_COVAR_ROWS;
|
||||
mean_sz = Size(len, 1);
|
||||
}
|
||||
|
||||
int count = std::min(len, in_count), out_count = count;
|
||||
if( maxComponents > 0 )
|
||||
out_count = std::min(count, maxComponents);
|
||||
|
||||
// "scrambled" way to compute PCA (when cols(A)>rows(A)):
|
||||
// B = A'A; B*x=b*x; C = AA'; C*y=c*y -> AA'*y=c*y -> A'A*(A'*y)=c*(A'*y) -> c = b, x=A'*y
|
||||
if( len <= in_count )
|
||||
covar_flags |= CV_COVAR_NORMAL;
|
||||
|
||||
int ctype = std::max(CV_32F, data.depth());
|
||||
mean.create( mean_sz, ctype );
|
||||
|
||||
Mat covar( count, count, ctype );
|
||||
|
||||
if( !_mean.empty() )
|
||||
{
|
||||
CV_Assert( _mean.size() == mean_sz );
|
||||
_mean.convertTo(mean, ctype);
|
||||
covar_flags |= CV_COVAR_USE_AVG;
|
||||
}
|
||||
|
||||
calcCovarMatrix( data, covar, mean, covar_flags, ctype );
|
||||
eigen( covar, eigenvalues, eigenvectors );
|
||||
|
||||
if( !(covar_flags & CV_COVAR_NORMAL) )
|
||||
{
|
||||
// CV_PCA_DATA_AS_ROW: cols(A)>rows(A). x=A'*y -> x'=y'*A
|
||||
// CV_PCA_DATA_AS_COL: rows(A)>cols(A). x=A''*y -> x'=y'*A'
|
||||
Mat tmp_data, tmp_mean = repeat(mean, data.rows/mean.rows, data.cols/mean.cols);
|
||||
if( data.type() != ctype || tmp_mean.data == mean.data )
|
||||
{
|
||||
data.convertTo( tmp_data, ctype );
|
||||
subtract( tmp_data, tmp_mean, tmp_data );
|
||||
}
|
||||
else
|
||||
{
|
||||
subtract( data, tmp_mean, tmp_mean );
|
||||
tmp_data = tmp_mean;
|
||||
}
|
||||
|
||||
Mat evects1(count, len, ctype);
|
||||
gemm( eigenvectors, tmp_data, 1, Mat(), 0, evects1,
|
||||
(flags & CV_PCA_DATA_AS_COL) ? CV_GEMM_B_T : 0);
|
||||
eigenvectors = evects1;
|
||||
|
||||
// normalize eigenvectors
|
||||
for( i = 0; i < out_count; i++ )
|
||||
{
|
||||
Mat vec = eigenvectors.row(i);
|
||||
normalize(vec, vec);
|
||||
}
|
||||
}
|
||||
|
||||
if( count > out_count )
|
||||
{
|
||||
// use clone() to physically copy the data and thus deallocate the original matrices
|
||||
eigenvalues = eigenvalues.rowRange(0,out_count).clone();
|
||||
eigenvectors = eigenvectors.rowRange(0,out_count).clone();
|
||||
}
|
||||
return *this;
|
||||
}
|
||||
|
||||
void PCA::write(FileStorage& fs ) const
|
||||
{
|
||||
CV_Assert( fs.isOpened() );
|
||||
|
||||
fs << "name" << "PCA";
|
||||
fs << "vectors" << eigenvectors;
|
||||
fs << "values" << eigenvalues;
|
||||
fs << "mean" << mean;
|
||||
}
|
||||
|
||||
void PCA::read(const FileNode& fs)
|
||||
{
|
||||
CV_Assert( !fs.empty() );
|
||||
String name = (String)fs["name"];
|
||||
CV_Assert( name == "PCA" );
|
||||
|
||||
cv::read(fs["vectors"], eigenvectors);
|
||||
cv::read(fs["values"], eigenvalues);
|
||||
cv::read(fs["mean"], mean);
|
||||
}
|
||||
|
||||
template <typename T>
|
||||
int computeCumulativeEnergy(const Mat& eigenvalues, double retainedVariance)
|
||||
{
|
||||
CV_DbgAssert( eigenvalues.type() == DataType<T>::type );
|
||||
|
||||
Mat g(eigenvalues.size(), DataType<T>::type);
|
||||
|
||||
for(int ig = 0; ig < g.rows; ig++)
|
||||
{
|
||||
g.at<T>(ig, 0) = 0;
|
||||
for(int im = 0; im <= ig; im++)
|
||||
{
|
||||
g.at<T>(ig,0) += eigenvalues.at<T>(im,0);
|
||||
}
|
||||
}
|
||||
|
||||
int L;
|
||||
|
||||
for(L = 0; L < eigenvalues.rows; L++)
|
||||
{
|
||||
double energy = g.at<T>(L, 0) / g.at<T>(g.rows - 1, 0);
|
||||
if(energy > retainedVariance)
|
||||
break;
|
||||
}
|
||||
|
||||
L = std::max(2, L);
|
||||
|
||||
return L;
|
||||
}
|
||||
|
||||
PCA& PCA::operator()(InputArray _data, InputArray __mean, int flags, double retainedVariance)
|
||||
{
|
||||
Mat data = _data.getMat(), _mean = __mean.getMat();
|
||||
int covar_flags = CV_COVAR_SCALE;
|
||||
int i, len, in_count;
|
||||
Size mean_sz;
|
||||
|
||||
CV_Assert( data.channels() == 1 );
|
||||
if( flags & CV_PCA_DATA_AS_COL )
|
||||
{
|
||||
len = data.rows;
|
||||
in_count = data.cols;
|
||||
covar_flags |= CV_COVAR_COLS;
|
||||
mean_sz = Size(1, len);
|
||||
}
|
||||
else
|
||||
{
|
||||
len = data.cols;
|
||||
in_count = data.rows;
|
||||
covar_flags |= CV_COVAR_ROWS;
|
||||
mean_sz = Size(len, 1);
|
||||
}
|
||||
|
||||
CV_Assert( retainedVariance > 0 && retainedVariance <= 1 );
|
||||
|
||||
int count = std::min(len, in_count);
|
||||
|
||||
// "scrambled" way to compute PCA (when cols(A)>rows(A)):
|
||||
// B = A'A; B*x=b*x; C = AA'; C*y=c*y -> AA'*y=c*y -> A'A*(A'*y)=c*(A'*y) -> c = b, x=A'*y
|
||||
if( len <= in_count )
|
||||
covar_flags |= CV_COVAR_NORMAL;
|
||||
|
||||
int ctype = std::max(CV_32F, data.depth());
|
||||
mean.create( mean_sz, ctype );
|
||||
|
||||
Mat covar( count, count, ctype );
|
||||
|
||||
if( !_mean.empty() )
|
||||
{
|
||||
CV_Assert( _mean.size() == mean_sz );
|
||||
_mean.convertTo(mean, ctype);
|
||||
}
|
||||
|
||||
calcCovarMatrix( data, covar, mean, covar_flags, ctype );
|
||||
eigen( covar, eigenvalues, eigenvectors );
|
||||
|
||||
if( !(covar_flags & CV_COVAR_NORMAL) )
|
||||
{
|
||||
// CV_PCA_DATA_AS_ROW: cols(A)>rows(A). x=A'*y -> x'=y'*A
|
||||
// CV_PCA_DATA_AS_COL: rows(A)>cols(A). x=A''*y -> x'=y'*A'
|
||||
Mat tmp_data, tmp_mean = repeat(mean, data.rows/mean.rows, data.cols/mean.cols);
|
||||
if( data.type() != ctype || tmp_mean.data == mean.data )
|
||||
{
|
||||
data.convertTo( tmp_data, ctype );
|
||||
subtract( tmp_data, tmp_mean, tmp_data );
|
||||
}
|
||||
else
|
||||
{
|
||||
subtract( data, tmp_mean, tmp_mean );
|
||||
tmp_data = tmp_mean;
|
||||
}
|
||||
|
||||
Mat evects1(count, len, ctype);
|
||||
gemm( eigenvectors, tmp_data, 1, Mat(), 0, evects1,
|
||||
(flags & CV_PCA_DATA_AS_COL) ? CV_GEMM_B_T : 0);
|
||||
eigenvectors = evects1;
|
||||
|
||||
// normalize all eigenvectors
|
||||
for( i = 0; i < eigenvectors.rows; i++ )
|
||||
{
|
||||
Mat vec = eigenvectors.row(i);
|
||||
normalize(vec, vec);
|
||||
}
|
||||
}
|
||||
|
||||
// compute the cumulative energy content for each eigenvector
|
||||
int L;
|
||||
if (ctype == CV_32F)
|
||||
L = computeCumulativeEnergy<float>(eigenvalues, retainedVariance);
|
||||
else
|
||||
L = computeCumulativeEnergy<double>(eigenvalues, retainedVariance);
|
||||
|
||||
// use clone() to physically copy the data and thus deallocate the original matrices
|
||||
eigenvalues = eigenvalues.rowRange(0,L).clone();
|
||||
eigenvectors = eigenvectors.rowRange(0,L).clone();
|
||||
|
||||
return *this;
|
||||
}
|
||||
|
||||
void PCA::project(InputArray _data, OutputArray result) const
|
||||
{
|
||||
Mat data = _data.getMat();
|
||||
CV_Assert( !mean.empty() && !eigenvectors.empty() &&
|
||||
((mean.rows == 1 && mean.cols == data.cols) || (mean.cols == 1 && mean.rows == data.rows)));
|
||||
Mat tmp_data, tmp_mean = repeat(mean, data.rows/mean.rows, data.cols/mean.cols);
|
||||
int ctype = mean.type();
|
||||
if( data.type() != ctype || tmp_mean.data == mean.data )
|
||||
{
|
||||
data.convertTo( tmp_data, ctype );
|
||||
subtract( tmp_data, tmp_mean, tmp_data );
|
||||
}
|
||||
else
|
||||
{
|
||||
subtract( data, tmp_mean, tmp_mean );
|
||||
tmp_data = tmp_mean;
|
||||
}
|
||||
if( mean.rows == 1 )
|
||||
gemm( tmp_data, eigenvectors, 1, Mat(), 0, result, GEMM_2_T );
|
||||
else
|
||||
gemm( eigenvectors, tmp_data, 1, Mat(), 0, result, 0 );
|
||||
}
|
||||
|
||||
Mat PCA::project(InputArray data) const
|
||||
{
|
||||
Mat result;
|
||||
project(data, result);
|
||||
return result;
|
||||
}
|
||||
|
||||
void PCA::backProject(InputArray _data, OutputArray result) const
|
||||
{
|
||||
Mat data = _data.getMat();
|
||||
CV_Assert( !mean.empty() && !eigenvectors.empty() &&
|
||||
((mean.rows == 1 && eigenvectors.rows == data.cols) ||
|
||||
(mean.cols == 1 && eigenvectors.rows == data.rows)));
|
||||
|
||||
Mat tmp_data, tmp_mean;
|
||||
data.convertTo(tmp_data, mean.type());
|
||||
if( mean.rows == 1 )
|
||||
{
|
||||
tmp_mean = repeat(mean, data.rows, 1);
|
||||
gemm( tmp_data, eigenvectors, 1, tmp_mean, 1, result, 0 );
|
||||
}
|
||||
else
|
||||
{
|
||||
tmp_mean = repeat(mean, 1, data.cols);
|
||||
gemm( eigenvectors, tmp_data, 1, tmp_mean, 1, result, GEMM_1_T );
|
||||
}
|
||||
}
|
||||
|
||||
Mat PCA::backProject(InputArray data) const
|
||||
{
|
||||
Mat result;
|
||||
backProject(data, result);
|
||||
return result;
|
||||
}
|
||||
|
||||
}
|
||||
|
||||
void cv::PCACompute(InputArray data, InputOutputArray mean,
|
||||
OutputArray eigenvectors, int maxComponents)
|
||||
{
|
||||
PCA pca;
|
||||
pca(data, mean, 0, maxComponents);
|
||||
pca.mean.copyTo(mean);
|
||||
pca.eigenvectors.copyTo(eigenvectors);
|
||||
}
|
||||
|
||||
void cv::PCACompute(InputArray data, InputOutputArray mean,
|
||||
OutputArray eigenvectors, double retainedVariance)
|
||||
{
|
||||
PCA pca;
|
||||
pca(data, mean, 0, retainedVariance);
|
||||
pca.mean.copyTo(mean);
|
||||
pca.eigenvectors.copyTo(eigenvectors);
|
||||
}
|
||||
|
||||
void cv::PCAProject(InputArray data, InputArray mean,
|
||||
InputArray eigenvectors, OutputArray result)
|
||||
{
|
||||
PCA pca;
|
||||
pca.mean = mean.getMat();
|
||||
pca.eigenvectors = eigenvectors.getMat();
|
||||
pca.project(data, result);
|
||||
}
|
||||
|
||||
void cv::PCABackProject(InputArray data, InputArray mean,
|
||||
InputArray eigenvectors, OutputArray result)
|
||||
{
|
||||
PCA pca;
|
||||
pca.mean = mean.getMat();
|
||||
pca.eigenvectors = eigenvectors.getMat();
|
||||
pca.backProject(data, result);
|
||||
}
|
||||
|
||||
/****************************************************************************************\
|
||||
|
||||
@@ -3968,420 +3968,6 @@ void cv::sortIdx( InputArray _src, OutputArray _dst, int flags )
|
||||
}
|
||||
|
||||
|
||||
////////////////////////////////////////// kmeans ////////////////////////////////////////////
|
||||
|
||||
namespace cv
|
||||
{
|
||||
|
||||
static void generateRandomCenter(const std::vector<Vec2f>& box, float* center, RNG& rng)
|
||||
{
|
||||
size_t j, dims = box.size();
|
||||
float margin = 1.f/dims;
|
||||
for( j = 0; j < dims; j++ )
|
||||
center[j] = ((float)rng*(1.f+margin*2.f)-margin)*(box[j][1] - box[j][0]) + box[j][0];
|
||||
}
|
||||
|
||||
class KMeansPPDistanceComputer : public ParallelLoopBody
|
||||
{
|
||||
public:
|
||||
KMeansPPDistanceComputer( float *_tdist2,
|
||||
const float *_data,
|
||||
const float *_dist,
|
||||
int _dims,
|
||||
size_t _step,
|
||||
size_t _stepci )
|
||||
: tdist2(_tdist2),
|
||||
data(_data),
|
||||
dist(_dist),
|
||||
dims(_dims),
|
||||
step(_step),
|
||||
stepci(_stepci) { }
|
||||
|
||||
void operator()( const cv::Range& range ) const
|
||||
{
|
||||
const int begin = range.start;
|
||||
const int end = range.end;
|
||||
|
||||
for ( int i = begin; i<end; i++ )
|
||||
{
|
||||
tdist2[i] = std::min(normL2Sqr_(data + step*i, data + stepci, dims), dist[i]);
|
||||
}
|
||||
}
|
||||
|
||||
private:
|
||||
KMeansPPDistanceComputer& operator=(const KMeansPPDistanceComputer&); // to quiet MSVC
|
||||
|
||||
float *tdist2;
|
||||
const float *data;
|
||||
const float *dist;
|
||||
const int dims;
|
||||
const size_t step;
|
||||
const size_t stepci;
|
||||
};
|
||||
|
||||
/*
|
||||
k-means center initialization using the following algorithm:
|
||||
Arthur & Vassilvitskii (2007) k-means++: The Advantages of Careful Seeding
|
||||
*/
|
||||
static void generateCentersPP(const Mat& _data, Mat& _out_centers,
|
||||
int K, RNG& rng, int trials)
|
||||
{
|
||||
int i, j, k, dims = _data.cols, N = _data.rows;
|
||||
const float* data = _data.ptr<float>(0);
|
||||
size_t step = _data.step/sizeof(data[0]);
|
||||
std::vector<int> _centers(K);
|
||||
int* centers = &_centers[0];
|
||||
std::vector<float> _dist(N*3);
|
||||
float* dist = &_dist[0], *tdist = dist + N, *tdist2 = tdist + N;
|
||||
double sum0 = 0;
|
||||
|
||||
centers[0] = (unsigned)rng % N;
|
||||
|
||||
for( i = 0; i < N; i++ )
|
||||
{
|
||||
dist[i] = normL2Sqr_(data + step*i, data + step*centers[0], dims);
|
||||
sum0 += dist[i];
|
||||
}
|
||||
|
||||
for( k = 1; k < K; k++ )
|
||||
{
|
||||
double bestSum = DBL_MAX;
|
||||
int bestCenter = -1;
|
||||
|
||||
for( j = 0; j < trials; j++ )
|
||||
{
|
||||
double p = (double)rng*sum0, s = 0;
|
||||
for( i = 0; i < N-1; i++ )
|
||||
if( (p -= dist[i]) <= 0 )
|
||||
break;
|
||||
int ci = i;
|
||||
|
||||
parallel_for_(Range(0, N),
|
||||
KMeansPPDistanceComputer(tdist2, data, dist, dims, step, step*ci));
|
||||
for( i = 0; i < N; i++ )
|
||||
{
|
||||
s += tdist2[i];
|
||||
}
|
||||
|
||||
if( s < bestSum )
|
||||
{
|
||||
bestSum = s;
|
||||
bestCenter = ci;
|
||||
std::swap(tdist, tdist2);
|
||||
}
|
||||
}
|
||||
centers[k] = bestCenter;
|
||||
sum0 = bestSum;
|
||||
std::swap(dist, tdist);
|
||||
}
|
||||
|
||||
for( k = 0; k < K; k++ )
|
||||
{
|
||||
const float* src = data + step*centers[k];
|
||||
float* dst = _out_centers.ptr<float>(k);
|
||||
for( j = 0; j < dims; j++ )
|
||||
dst[j] = src[j];
|
||||
}
|
||||
}
|
||||
|
||||
class KMeansDistanceComputer : public ParallelLoopBody
|
||||
{
|
||||
public:
|
||||
KMeansDistanceComputer( double *_distances,
|
||||
int *_labels,
|
||||
const Mat& _data,
|
||||
const Mat& _centers )
|
||||
: distances(_distances),
|
||||
labels(_labels),
|
||||
data(_data),
|
||||
centers(_centers)
|
||||
{
|
||||
}
|
||||
|
||||
void operator()( const Range& range ) const
|
||||
{
|
||||
const int begin = range.start;
|
||||
const int end = range.end;
|
||||
const int K = centers.rows;
|
||||
const int dims = centers.cols;
|
||||
|
||||
const float *sample;
|
||||
for( int i = begin; i<end; ++i)
|
||||
{
|
||||
sample = data.ptr<float>(i);
|
||||
int k_best = 0;
|
||||
double min_dist = DBL_MAX;
|
||||
|
||||
for( int k = 0; k < K; k++ )
|
||||
{
|
||||
const float* center = centers.ptr<float>(k);
|
||||
const double dist = normL2Sqr_(sample, center, dims);
|
||||
|
||||
if( min_dist > dist )
|
||||
{
|
||||
min_dist = dist;
|
||||
k_best = k;
|
||||
}
|
||||
}
|
||||
|
||||
distances[i] = min_dist;
|
||||
labels[i] = k_best;
|
||||
}
|
||||
}
|
||||
|
||||
private:
|
||||
KMeansDistanceComputer& operator=(const KMeansDistanceComputer&); // to quiet MSVC
|
||||
|
||||
double *distances;
|
||||
int *labels;
|
||||
const Mat& data;
|
||||
const Mat& centers;
|
||||
};
|
||||
|
||||
}
|
||||
|
||||
double cv::kmeans( InputArray _data, int K,
|
||||
InputOutputArray _bestLabels,
|
||||
TermCriteria criteria, int attempts,
|
||||
int flags, OutputArray _centers )
|
||||
{
|
||||
const int SPP_TRIALS = 3;
|
||||
Mat data0 = _data.getMat();
|
||||
bool isrow = data0.rows == 1 && data0.channels() > 1;
|
||||
int N = !isrow ? data0.rows : data0.cols;
|
||||
int dims = (!isrow ? data0.cols : 1)*data0.channels();
|
||||
int type = data0.depth();
|
||||
|
||||
attempts = std::max(attempts, 1);
|
||||
CV_Assert( data0.dims <= 2 && type == CV_32F && K > 0 );
|
||||
CV_Assert( N >= K );
|
||||
|
||||
Mat data(N, dims, CV_32F, data0.ptr(), isrow ? dims * sizeof(float) : static_cast<size_t>(data0.step));
|
||||
|
||||
_bestLabels.create(N, 1, CV_32S, -1, true);
|
||||
|
||||
Mat _labels, best_labels = _bestLabels.getMat();
|
||||
if( flags & CV_KMEANS_USE_INITIAL_LABELS )
|
||||
{
|
||||
CV_Assert( (best_labels.cols == 1 || best_labels.rows == 1) &&
|
||||
best_labels.cols*best_labels.rows == N &&
|
||||
best_labels.type() == CV_32S &&
|
||||
best_labels.isContinuous());
|
||||
best_labels.copyTo(_labels);
|
||||
}
|
||||
else
|
||||
{
|
||||
if( !((best_labels.cols == 1 || best_labels.rows == 1) &&
|
||||
best_labels.cols*best_labels.rows == N &&
|
||||
best_labels.type() == CV_32S &&
|
||||
best_labels.isContinuous()))
|
||||
best_labels.create(N, 1, CV_32S);
|
||||
_labels.create(best_labels.size(), best_labels.type());
|
||||
}
|
||||
int* labels = _labels.ptr<int>();
|
||||
|
||||
Mat centers(K, dims, type), old_centers(K, dims, type), temp(1, dims, type);
|
||||
std::vector<int> counters(K);
|
||||
std::vector<Vec2f> _box(dims);
|
||||
Vec2f* box = &_box[0];
|
||||
double best_compactness = DBL_MAX, compactness = 0;
|
||||
RNG& rng = theRNG();
|
||||
int a, iter, i, j, k;
|
||||
|
||||
if( criteria.type & TermCriteria::EPS )
|
||||
criteria.epsilon = std::max(criteria.epsilon, 0.);
|
||||
else
|
||||
criteria.epsilon = FLT_EPSILON;
|
||||
criteria.epsilon *= criteria.epsilon;
|
||||
|
||||
if( criteria.type & TermCriteria::COUNT )
|
||||
criteria.maxCount = std::min(std::max(criteria.maxCount, 2), 100);
|
||||
else
|
||||
criteria.maxCount = 100;
|
||||
|
||||
if( K == 1 )
|
||||
{
|
||||
attempts = 1;
|
||||
criteria.maxCount = 2;
|
||||
}
|
||||
|
||||
const float* sample = data.ptr<float>(0);
|
||||
for( j = 0; j < dims; j++ )
|
||||
box[j] = Vec2f(sample[j], sample[j]);
|
||||
|
||||
for( i = 1; i < N; i++ )
|
||||
{
|
||||
sample = data.ptr<float>(i);
|
||||
for( j = 0; j < dims; j++ )
|
||||
{
|
||||
float v = sample[j];
|
||||
box[j][0] = std::min(box[j][0], v);
|
||||
box[j][1] = std::max(box[j][1], v);
|
||||
}
|
||||
}
|
||||
|
||||
for( a = 0; a < attempts; a++ )
|
||||
{
|
||||
double max_center_shift = DBL_MAX;
|
||||
for( iter = 0;; )
|
||||
{
|
||||
swap(centers, old_centers);
|
||||
|
||||
if( iter == 0 && (a > 0 || !(flags & KMEANS_USE_INITIAL_LABELS)) )
|
||||
{
|
||||
if( flags & KMEANS_PP_CENTERS )
|
||||
generateCentersPP(data, centers, K, rng, SPP_TRIALS);
|
||||
else
|
||||
{
|
||||
for( k = 0; k < K; k++ )
|
||||
generateRandomCenter(_box, centers.ptr<float>(k), rng);
|
||||
}
|
||||
}
|
||||
else
|
||||
{
|
||||
if( iter == 0 && a == 0 && (flags & KMEANS_USE_INITIAL_LABELS) )
|
||||
{
|
||||
for( i = 0; i < N; i++ )
|
||||
CV_Assert( (unsigned)labels[i] < (unsigned)K );
|
||||
}
|
||||
|
||||
// compute centers
|
||||
centers = Scalar(0);
|
||||
for( k = 0; k < K; k++ )
|
||||
counters[k] = 0;
|
||||
|
||||
for( i = 0; i < N; i++ )
|
||||
{
|
||||
sample = data.ptr<float>(i);
|
||||
k = labels[i];
|
||||
float* center = centers.ptr<float>(k);
|
||||
j=0;
|
||||
#if CV_ENABLE_UNROLLED
|
||||
for(; j <= dims - 4; j += 4 )
|
||||
{
|
||||
float t0 = center[j] + sample[j];
|
||||
float t1 = center[j+1] + sample[j+1];
|
||||
|
||||
center[j] = t0;
|
||||
center[j+1] = t1;
|
||||
|
||||
t0 = center[j+2] + sample[j+2];
|
||||
t1 = center[j+3] + sample[j+3];
|
||||
|
||||
center[j+2] = t0;
|
||||
center[j+3] = t1;
|
||||
}
|
||||
#endif
|
||||
for( ; j < dims; j++ )
|
||||
center[j] += sample[j];
|
||||
counters[k]++;
|
||||
}
|
||||
|
||||
if( iter > 0 )
|
||||
max_center_shift = 0;
|
||||
|
||||
for( k = 0; k < K; k++ )
|
||||
{
|
||||
if( counters[k] != 0 )
|
||||
continue;
|
||||
|
||||
// if some cluster appeared to be empty then:
|
||||
// 1. find the biggest cluster
|
||||
// 2. find the farthest from the center point in the biggest cluster
|
||||
// 3. exclude the farthest point from the biggest cluster and form a new 1-point cluster.
|
||||
int max_k = 0;
|
||||
for( int k1 = 1; k1 < K; k1++ )
|
||||
{
|
||||
if( counters[max_k] < counters[k1] )
|
||||
max_k = k1;
|
||||
}
|
||||
|
||||
double max_dist = 0;
|
||||
int farthest_i = -1;
|
||||
float* new_center = centers.ptr<float>(k);
|
||||
float* old_center = centers.ptr<float>(max_k);
|
||||
float* _old_center = temp.ptr<float>(); // normalized
|
||||
float scale = 1.f/counters[max_k];
|
||||
for( j = 0; j < dims; j++ )
|
||||
_old_center[j] = old_center[j]*scale;
|
||||
|
||||
for( i = 0; i < N; i++ )
|
||||
{
|
||||
if( labels[i] != max_k )
|
||||
continue;
|
||||
sample = data.ptr<float>(i);
|
||||
double dist = normL2Sqr_(sample, _old_center, dims);
|
||||
|
||||
if( max_dist <= dist )
|
||||
{
|
||||
max_dist = dist;
|
||||
farthest_i = i;
|
||||
}
|
||||
}
|
||||
|
||||
counters[max_k]--;
|
||||
counters[k]++;
|
||||
labels[farthest_i] = k;
|
||||
sample = data.ptr<float>(farthest_i);
|
||||
|
||||
for( j = 0; j < dims; j++ )
|
||||
{
|
||||
old_center[j] -= sample[j];
|
||||
new_center[j] += sample[j];
|
||||
}
|
||||
}
|
||||
|
||||
for( k = 0; k < K; k++ )
|
||||
{
|
||||
float* center = centers.ptr<float>(k);
|
||||
CV_Assert( counters[k] != 0 );
|
||||
|
||||
float scale = 1.f/counters[k];
|
||||
for( j = 0; j < dims; j++ )
|
||||
center[j] *= scale;
|
||||
|
||||
if( iter > 0 )
|
||||
{
|
||||
double dist = 0;
|
||||
const float* old_center = old_centers.ptr<float>(k);
|
||||
for( j = 0; j < dims; j++ )
|
||||
{
|
||||
double t = center[j] - old_center[j];
|
||||
dist += t*t;
|
||||
}
|
||||
max_center_shift = std::max(max_center_shift, dist);
|
||||
}
|
||||
}
|
||||
}
|
||||
|
||||
if( ++iter == MAX(criteria.maxCount, 2) || max_center_shift <= criteria.epsilon )
|
||||
break;
|
||||
|
||||
// assign labels
|
||||
Mat dists(1, N, CV_64F);
|
||||
double* dist = dists.ptr<double>(0);
|
||||
parallel_for_(Range(0, N),
|
||||
KMeansDistanceComputer(dist, labels, data, centers));
|
||||
compactness = 0;
|
||||
for( i = 0; i < N; i++ )
|
||||
{
|
||||
compactness += dist[i];
|
||||
}
|
||||
}
|
||||
|
||||
if( compactness < best_compactness )
|
||||
{
|
||||
best_compactness = compactness;
|
||||
if( _centers.needed() )
|
||||
centers.copyTo(_centers);
|
||||
_labels.copyTo(best_labels);
|
||||
}
|
||||
}
|
||||
|
||||
return best_compactness;
|
||||
}
|
||||
|
||||
|
||||
CV_IMPL void cvSetIdentity( CvArr* arr, CvScalar value )
|
||||
{
|
||||
cv::Mat m = cv::cvarrToMat(arr);
|
||||
|
||||
@@ -0,0 +1,384 @@
|
||||
/*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.
|
||||
//
|
||||
//
|
||||
// License Agreement
|
||||
// For Open Source Computer Vision Library
|
||||
//
|
||||
// Copyright (C) 2000-2008, Intel Corporation, all rights reserved.
|
||||
// Copyright (C) 2009, Willow Garage Inc., all rights reserved.
|
||||
// Copyright (C) 2013, OpenCV Foundation, 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 the copyright holders 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*/
|
||||
|
||||
#include "precomp.hpp"
|
||||
|
||||
/****************************************************************************************\
|
||||
* PCA *
|
||||
\****************************************************************************************/
|
||||
|
||||
namespace cv
|
||||
{
|
||||
|
||||
PCA::PCA() {}
|
||||
|
||||
PCA::PCA(InputArray data, InputArray _mean, int flags, int maxComponents)
|
||||
{
|
||||
operator()(data, _mean, flags, maxComponents);
|
||||
}
|
||||
|
||||
PCA::PCA(InputArray data, InputArray _mean, int flags, double retainedVariance)
|
||||
{
|
||||
operator()(data, _mean, flags, retainedVariance);
|
||||
}
|
||||
|
||||
PCA& PCA::operator()(InputArray _data, InputArray __mean, int flags, int maxComponents)
|
||||
{
|
||||
Mat data = _data.getMat(), _mean = __mean.getMat();
|
||||
int covar_flags = CV_COVAR_SCALE;
|
||||
int i, len, in_count;
|
||||
Size mean_sz;
|
||||
|
||||
CV_Assert( data.channels() == 1 );
|
||||
if( flags & CV_PCA_DATA_AS_COL )
|
||||
{
|
||||
len = data.rows;
|
||||
in_count = data.cols;
|
||||
covar_flags |= CV_COVAR_COLS;
|
||||
mean_sz = Size(1, len);
|
||||
}
|
||||
else
|
||||
{
|
||||
len = data.cols;
|
||||
in_count = data.rows;
|
||||
covar_flags |= CV_COVAR_ROWS;
|
||||
mean_sz = Size(len, 1);
|
||||
}
|
||||
|
||||
int count = std::min(len, in_count), out_count = count;
|
||||
if( maxComponents > 0 )
|
||||
out_count = std::min(count, maxComponents);
|
||||
|
||||
// "scrambled" way to compute PCA (when cols(A)>rows(A)):
|
||||
// B = A'A; B*x=b*x; C = AA'; C*y=c*y -> AA'*y=c*y -> A'A*(A'*y)=c*(A'*y) -> c = b, x=A'*y
|
||||
if( len <= in_count )
|
||||
covar_flags |= CV_COVAR_NORMAL;
|
||||
|
||||
int ctype = std::max(CV_32F, data.depth());
|
||||
mean.create( mean_sz, ctype );
|
||||
|
||||
Mat covar( count, count, ctype );
|
||||
|
||||
if( !_mean.empty() )
|
||||
{
|
||||
CV_Assert( _mean.size() == mean_sz );
|
||||
_mean.convertTo(mean, ctype);
|
||||
covar_flags |= CV_COVAR_USE_AVG;
|
||||
}
|
||||
|
||||
calcCovarMatrix( data, covar, mean, covar_flags, ctype );
|
||||
eigen( covar, eigenvalues, eigenvectors );
|
||||
|
||||
if( !(covar_flags & CV_COVAR_NORMAL) )
|
||||
{
|
||||
// CV_PCA_DATA_AS_ROW: cols(A)>rows(A). x=A'*y -> x'=y'*A
|
||||
// CV_PCA_DATA_AS_COL: rows(A)>cols(A). x=A''*y -> x'=y'*A'
|
||||
Mat tmp_data, tmp_mean = repeat(mean, data.rows/mean.rows, data.cols/mean.cols);
|
||||
if( data.type() != ctype || tmp_mean.data == mean.data )
|
||||
{
|
||||
data.convertTo( tmp_data, ctype );
|
||||
subtract( tmp_data, tmp_mean, tmp_data );
|
||||
}
|
||||
else
|
||||
{
|
||||
subtract( data, tmp_mean, tmp_mean );
|
||||
tmp_data = tmp_mean;
|
||||
}
|
||||
|
||||
Mat evects1(count, len, ctype);
|
||||
gemm( eigenvectors, tmp_data, 1, Mat(), 0, evects1,
|
||||
(flags & CV_PCA_DATA_AS_COL) ? CV_GEMM_B_T : 0);
|
||||
eigenvectors = evects1;
|
||||
|
||||
// normalize eigenvectors
|
||||
for( i = 0; i < out_count; i++ )
|
||||
{
|
||||
Mat vec = eigenvectors.row(i);
|
||||
normalize(vec, vec);
|
||||
}
|
||||
}
|
||||
|
||||
if( count > out_count )
|
||||
{
|
||||
// use clone() to physically copy the data and thus deallocate the original matrices
|
||||
eigenvalues = eigenvalues.rowRange(0,out_count).clone();
|
||||
eigenvectors = eigenvectors.rowRange(0,out_count).clone();
|
||||
}
|
||||
return *this;
|
||||
}
|
||||
|
||||
void PCA::write(FileStorage& fs ) const
|
||||
{
|
||||
CV_Assert( fs.isOpened() );
|
||||
|
||||
fs << "name" << "PCA";
|
||||
fs << "vectors" << eigenvectors;
|
||||
fs << "values" << eigenvalues;
|
||||
fs << "mean" << mean;
|
||||
}
|
||||
|
||||
void PCA::read(const FileNode& fs)
|
||||
{
|
||||
CV_Assert( !fs.empty() );
|
||||
String name = (String)fs["name"];
|
||||
CV_Assert( name == "PCA" );
|
||||
|
||||
cv::read(fs["vectors"], eigenvectors);
|
||||
cv::read(fs["values"], eigenvalues);
|
||||
cv::read(fs["mean"], mean);
|
||||
}
|
||||
|
||||
template <typename T>
|
||||
int computeCumulativeEnergy(const Mat& eigenvalues, double retainedVariance)
|
||||
{
|
||||
CV_DbgAssert( eigenvalues.type() == DataType<T>::type );
|
||||
|
||||
Mat g(eigenvalues.size(), DataType<T>::type);
|
||||
|
||||
for(int ig = 0; ig < g.rows; ig++)
|
||||
{
|
||||
g.at<T>(ig, 0) = 0;
|
||||
for(int im = 0; im <= ig; im++)
|
||||
{
|
||||
g.at<T>(ig,0) += eigenvalues.at<T>(im,0);
|
||||
}
|
||||
}
|
||||
|
||||
int L;
|
||||
|
||||
for(L = 0; L < eigenvalues.rows; L++)
|
||||
{
|
||||
double energy = g.at<T>(L, 0) / g.at<T>(g.rows - 1, 0);
|
||||
if(energy > retainedVariance)
|
||||
break;
|
||||
}
|
||||
|
||||
L = std::max(2, L);
|
||||
|
||||
return L;
|
||||
}
|
||||
|
||||
PCA& PCA::operator()(InputArray _data, InputArray __mean, int flags, double retainedVariance)
|
||||
{
|
||||
Mat data = _data.getMat(), _mean = __mean.getMat();
|
||||
int covar_flags = CV_COVAR_SCALE;
|
||||
int i, len, in_count;
|
||||
Size mean_sz;
|
||||
|
||||
CV_Assert( data.channels() == 1 );
|
||||
if( flags & CV_PCA_DATA_AS_COL )
|
||||
{
|
||||
len = data.rows;
|
||||
in_count = data.cols;
|
||||
covar_flags |= CV_COVAR_COLS;
|
||||
mean_sz = Size(1, len);
|
||||
}
|
||||
else
|
||||
{
|
||||
len = data.cols;
|
||||
in_count = data.rows;
|
||||
covar_flags |= CV_COVAR_ROWS;
|
||||
mean_sz = Size(len, 1);
|
||||
}
|
||||
|
||||
CV_Assert( retainedVariance > 0 && retainedVariance <= 1 );
|
||||
|
||||
int count = std::min(len, in_count);
|
||||
|
||||
// "scrambled" way to compute PCA (when cols(A)>rows(A)):
|
||||
// B = A'A; B*x=b*x; C = AA'; C*y=c*y -> AA'*y=c*y -> A'A*(A'*y)=c*(A'*y) -> c = b, x=A'*y
|
||||
if( len <= in_count )
|
||||
covar_flags |= CV_COVAR_NORMAL;
|
||||
|
||||
int ctype = std::max(CV_32F, data.depth());
|
||||
mean.create( mean_sz, ctype );
|
||||
|
||||
Mat covar( count, count, ctype );
|
||||
|
||||
if( !_mean.empty() )
|
||||
{
|
||||
CV_Assert( _mean.size() == mean_sz );
|
||||
_mean.convertTo(mean, ctype);
|
||||
}
|
||||
|
||||
calcCovarMatrix( data, covar, mean, covar_flags, ctype );
|
||||
eigen( covar, eigenvalues, eigenvectors );
|
||||
|
||||
if( !(covar_flags & CV_COVAR_NORMAL) )
|
||||
{
|
||||
// CV_PCA_DATA_AS_ROW: cols(A)>rows(A). x=A'*y -> x'=y'*A
|
||||
// CV_PCA_DATA_AS_COL: rows(A)>cols(A). x=A''*y -> x'=y'*A'
|
||||
Mat tmp_data, tmp_mean = repeat(mean, data.rows/mean.rows, data.cols/mean.cols);
|
||||
if( data.type() != ctype || tmp_mean.data == mean.data )
|
||||
{
|
||||
data.convertTo( tmp_data, ctype );
|
||||
subtract( tmp_data, tmp_mean, tmp_data );
|
||||
}
|
||||
else
|
||||
{
|
||||
subtract( data, tmp_mean, tmp_mean );
|
||||
tmp_data = tmp_mean;
|
||||
}
|
||||
|
||||
Mat evects1(count, len, ctype);
|
||||
gemm( eigenvectors, tmp_data, 1, Mat(), 0, evects1,
|
||||
(flags & CV_PCA_DATA_AS_COL) ? CV_GEMM_B_T : 0);
|
||||
eigenvectors = evects1;
|
||||
|
||||
// normalize all eigenvectors
|
||||
for( i = 0; i < eigenvectors.rows; i++ )
|
||||
{
|
||||
Mat vec = eigenvectors.row(i);
|
||||
normalize(vec, vec);
|
||||
}
|
||||
}
|
||||
|
||||
// compute the cumulative energy content for each eigenvector
|
||||
int L;
|
||||
if (ctype == CV_32F)
|
||||
L = computeCumulativeEnergy<float>(eigenvalues, retainedVariance);
|
||||
else
|
||||
L = computeCumulativeEnergy<double>(eigenvalues, retainedVariance);
|
||||
|
||||
// use clone() to physically copy the data and thus deallocate the original matrices
|
||||
eigenvalues = eigenvalues.rowRange(0,L).clone();
|
||||
eigenvectors = eigenvectors.rowRange(0,L).clone();
|
||||
|
||||
return *this;
|
||||
}
|
||||
|
||||
void PCA::project(InputArray _data, OutputArray result) const
|
||||
{
|
||||
Mat data = _data.getMat();
|
||||
CV_Assert( !mean.empty() && !eigenvectors.empty() &&
|
||||
((mean.rows == 1 && mean.cols == data.cols) || (mean.cols == 1 && mean.rows == data.rows)));
|
||||
Mat tmp_data, tmp_mean = repeat(mean, data.rows/mean.rows, data.cols/mean.cols);
|
||||
int ctype = mean.type();
|
||||
if( data.type() != ctype || tmp_mean.data == mean.data )
|
||||
{
|
||||
data.convertTo( tmp_data, ctype );
|
||||
subtract( tmp_data, tmp_mean, tmp_data );
|
||||
}
|
||||
else
|
||||
{
|
||||
subtract( data, tmp_mean, tmp_mean );
|
||||
tmp_data = tmp_mean;
|
||||
}
|
||||
if( mean.rows == 1 )
|
||||
gemm( tmp_data, eigenvectors, 1, Mat(), 0, result, GEMM_2_T );
|
||||
else
|
||||
gemm( eigenvectors, tmp_data, 1, Mat(), 0, result, 0 );
|
||||
}
|
||||
|
||||
Mat PCA::project(InputArray data) const
|
||||
{
|
||||
Mat result;
|
||||
project(data, result);
|
||||
return result;
|
||||
}
|
||||
|
||||
void PCA::backProject(InputArray _data, OutputArray result) const
|
||||
{
|
||||
Mat data = _data.getMat();
|
||||
CV_Assert( !mean.empty() && !eigenvectors.empty() &&
|
||||
((mean.rows == 1 && eigenvectors.rows == data.cols) ||
|
||||
(mean.cols == 1 && eigenvectors.rows == data.rows)));
|
||||
|
||||
Mat tmp_data, tmp_mean;
|
||||
data.convertTo(tmp_data, mean.type());
|
||||
if( mean.rows == 1 )
|
||||
{
|
||||
tmp_mean = repeat(mean, data.rows, 1);
|
||||
gemm( tmp_data, eigenvectors, 1, tmp_mean, 1, result, 0 );
|
||||
}
|
||||
else
|
||||
{
|
||||
tmp_mean = repeat(mean, 1, data.cols);
|
||||
gemm( eigenvectors, tmp_data, 1, tmp_mean, 1, result, GEMM_1_T );
|
||||
}
|
||||
}
|
||||
|
||||
Mat PCA::backProject(InputArray data) const
|
||||
{
|
||||
Mat result;
|
||||
backProject(data, result);
|
||||
return result;
|
||||
}
|
||||
|
||||
}
|
||||
|
||||
void cv::PCACompute(InputArray data, InputOutputArray mean,
|
||||
OutputArray eigenvectors, int maxComponents)
|
||||
{
|
||||
PCA pca;
|
||||
pca(data, mean, 0, maxComponents);
|
||||
pca.mean.copyTo(mean);
|
||||
pca.eigenvectors.copyTo(eigenvectors);
|
||||
}
|
||||
|
||||
void cv::PCACompute(InputArray data, InputOutputArray mean,
|
||||
OutputArray eigenvectors, double retainedVariance)
|
||||
{
|
||||
PCA pca;
|
||||
pca(data, mean, 0, retainedVariance);
|
||||
pca.mean.copyTo(mean);
|
||||
pca.eigenvectors.copyTo(eigenvectors);
|
||||
}
|
||||
|
||||
void cv::PCAProject(InputArray data, InputArray mean,
|
||||
InputArray eigenvectors, OutputArray result)
|
||||
{
|
||||
PCA pca;
|
||||
pca.mean = mean.getMat();
|
||||
pca.eigenvectors = eigenvectors.getMat();
|
||||
pca.project(data, result);
|
||||
}
|
||||
|
||||
void cv::PCABackProject(InputArray data, InputArray mean,
|
||||
InputArray eigenvectors, OutputArray result)
|
||||
{
|
||||
PCA pca;
|
||||
pca.mean = mean.getMat();
|
||||
pca.eigenvectors = eigenvectors.getMat();
|
||||
pca.backProject(data, result);
|
||||
}
|
||||
Reference in New Issue
Block a user