中位中心栅格遍历法

中位中心栅格遍历法

#include <iostream>
#include <cmath>

namespace geo
{
struct POINT
{
double x;
double y;
};

//栅格法寻找中位中心
POINT GetMedianCentre(const POINT * pts, int n, double d/*设置栅格化格网宽度*/)
{
double min_x = pts[0].x;
double min_y = pts[0].y;
double max_x = pts[0].x;
double max_y = pts[0].y;

for(int i = 0; i < n; i++)
{
if(pts[i].x < min_x)
{
min_x = pts[i].x;
}
else if(pts[i].x > max_x)
{
max_x = pts[i].x;
}
if(pts[i].y < min_y)
{
min_y = pts[i].y;
}
else if(pts[i].y > max_y)
{
max_y = pts[i].y;
}
}

int raster_row = (max_y - min_y) / d + 1;//注意 + 1
int raster_column = (max_x - min_x) / d + 1;

struct INTPT {int row; int col;};//行列数对
INTPT * pts_o = new INTPT
;//原有点集栅格化后的行列数对
for(int i = 0; i < n; i++)
{
pts_o[i].row = (int)((pts[i].y - min_y) / d);
pts_o[i].col = (int)((pts[i].x - min_x) / d);
}

int des_i;//用来记录中位中心行号
int des_j;//用来记录中位中心列号
double min_total_distance;
//逐个像元遍历
for(int i = 0; i < raster_row; i++)
{
for(int j = 0; j < raster_column; j++)
{
double total_distance = 0.0;
for(int k = 0; k < n; k++)
{
total_distance += sqrt(pow(i - pts_o[k].row, 2.0) + pow(j - pts_o[k].col, 2.0));
}
if(i == 0 && j == 0)//将（0，0）到点群中各个点总距离作为最短总距离初始化
{
min_total_distance = total_distance;
des_i = i;
des_j = j;
}
else if(total_distance < min_total_distance)
{
min_total_distance = total_distance;
des_i = i;
des_j = j;
}
}
}
delete [] pts_o;

POINT des_point = {d / 2.0 + des_j * d + min_x, d / 2.0 + des_i * d + min_y};
return des_point;
}
}
int main()
{
using namespace std;
using namespace geo;

POINT pts[3] = {{0, 0}, {2 * sqrt(3.0), 0}, {sqrt(3.0), 3}};
//测试函数：正三角形中位中心为几何中心，细分网格可逼近实际值
POINT ptdes = geo::GetMedianCentre(pts, 3, 1);
cout<<"("<<ptdes.x<<", "<<ptdes.y<<")"<<endl;
ptdes = geo::GetMedianCentre(pts, 3, 0.5);
cout<<"("<<ptdes.x<<", "<<ptdes.y<<")"<<endl;
ptdes = geo::GetMedianCentre(pts, 3, 0.2);
cout<<"("<<ptdes.x<<", "<<ptdes.y<<")"<<endl;
ptdes = geo::GetMedianCentre(pts, 3, 0.1);
cout<<"("<<ptdes.x<<", "<<ptdes.y<<")"<<endl;
ptdes = geo::GetMedianCentre(pts, 3, 0.05);
cout<<"("<<ptdes.x<<", "<<ptdes.y<<")"<<endl;
ptdes = geo::GetMedianCentre(pts, 3, 0.01);
cout<<"("<<ptdes.x<<", "<<ptdes.y<<")"<<endl;

return 0;
}