基于matlab关于重心的算法,基于matlab mex的平面点集按重心逆时针排序算法

基于matlab mex的平面点集按重心逆时针排序算法，可用于求凸集，代码如下：

#include

#include

#include

#include

#include

using namespace std;

typedef struct PointF {

int id;

double x, y;

PointF() {

id = -1;

}

PointF(int id, double x, double y) : id(id), x(x), y(y) {}

} PointF;

typedef vector Points;

Points points;

PointF cpt;

bool compare(const PointF p1, const PointF p2) {

//mexPrintf("%d %d\n", p1.id, p2.id);

PointF a = p1, b = p2;

int aid, bid;

if (a.x >= cpt.x) {

if (a.y <= cpt.y)

aid = 1;

else

aid = 2;

}

if (a.x < cpt.x) {

if (a.y >= cpt.y)

aid = 3;

else

aid = 4;

}

if (b.x >= cpt.x) {

if (b.y <= cpt.y)

bid = 1;

else

bid = 2;

}

if (b.x < cpt.x) {

if (b.y >= cpt.y)

bid = 3;

else

bid = 4;

}

if (aid < bid)

return true;

else if (aid == bid) {

if (aid == 1 || aid == 2)

return a.y < b.y;

if (aid == 3 || aid == 4)

return a.y > b.y;

}

else

return false;

}

void mexFunction(int nlhs, mxArray *plhs[], int nrhs, const mxArray *prhs[]) {

if (nrhs != 2) {

mexPrintf("Usage: 2dsort(xy_nrows_2cols, refpoint_x_y).\n");

return;

}

int ndim = mxGetNumberOfDimensions(prhs[0]);

//mexPrintf("%d\n", ndim);

const int *dims = mxGetDimensions( prhs[0] );

//mexPrintf("%d %d dims\n", dims[0], dims[1]);

if (dims[1] != 2) {

mexPrintf("Error: input data should be with 2 columns.\n");

return;

}

double *data = mxGetPr( prhs[0] );

for(int i = 0; i < dims[0]; i++) {

points.push_back(PointF(i + 1, data[i], data[i + dims[0]]));

//mexPrintf("*%d ", points[i].id);

}

double *data1 = mxGetPr( prhs[1] );

cpt.x = data1[0];

cpt.y = data1[1];

for(int i = 0; i < dims[0]; i++) {

mexPrintf("%d ", points[i].id);

}

mexPrintf("\n");

sort(points.begin(), points.end(), compare);

for(int i = 0; i < dims[0]; i++) {

mexPrintf("%d ", points[i].id);

}

mexPrintf("\n");

nlhs = 1;

plhs[0] = mxCreateDoubleMatrix(dims[0], 1, mxREAL);

double *dest = mxGetPr( plhs[0] );

for(int i = 0; i < dims[0]; i++)

dest[i] = points[i].id;

}

mex twodsort.cpp

测试代码：

close all;

mex twodsort.cpp;

N = 300;

xy = rand(N, 2);

ctr = sum(xy) ./ N;

plot(xy(:, 1), xy(:, 2), '*r');

hold on;

plot(ctr(1), ctr(2), '+b');

ids = twodsort(xy, ctr);

ids(end + 1) = ids(1);

plot(xy(ids(1), 1), xy(ids(1), 2), '.g', 'MarkerSize', 5);

for i = 1:N

j = i + 1;

if (j > N)

j = 1;

end;

plot([xy(ids(i), 1); xy(ids(j), 1)], [xy(ids(i), 2); xy(ids(j), 2)], '-', 'Color', [i/(N+1), 0, 0], 'LineWidth', 2);

end;