遥感图像处理-K值聚类-地物分类
（环境：VS2017+OPENCV4.0）
K值聚类算法并不复杂，这里不多介绍，直接上代码。
#include "pch.h"
#include <iostream>
#include <opencv2/opencv.hpp>
#include <math.h>
using namespace std;
using namespace cv;struct SPoint {//点和点的类别
double v1;
double v2;
double v3;
int id;
};
void Color(SPoint* pt, int id) {//给不同类别赋予颜色，可以自己加
switch (pt->id) {
case 0:
pt->v1 = 0;
pt->v2 = 255;
pt->v3 = 170;
break;
case 1:
pt->v1 = 255;
pt->v2 = 0;
pt->v3 = 0;
break;
case 2:
pt->v1 = 0;
pt->v2 = 255;
pt->v3 = 0;
break;
case 3:
pt->v1 = 0;
pt->v2 = 0;
pt->v3 = 255;
break;
case 4:
pt->v1 = 255;
pt->v2 = 170;
pt->v3 = 200;
break;
case 5:
pt->v1 = 100;
pt->v2 = 100;
pt->v3 = 100;
break;
case 6:
pt->v1 = 200;
pt->v2 = 200;
pt->v3 = 200;
break;
}
}
double Norm2(SPoint p1, SPoint p2) {//距离的平方
return (pow((p1.v1 - p2.v1), 2) + pow((p1.v2 - p2.v2), 2) + pow((p1.v3 - p2.v3), 2));
//return (p1.v1 - p2.v1)*(p1.v1 - p2.v1) + (p1.v2 - p2.v2)*(p1.v2 - p2.v2) + (p1.v3 - p2.v3)*(p1.v3 - p2.v3);
}
Mat Kmean(Mat src, int species, int iterations) {//原图像，种类数，迭代次数
Mat K = src;
int rows = src.rows, cols = src.cols;//行，列
SPoint* Points = new SPoint[rows*cols];
SPoint* CenterPt1 = new SPoint[species];
SPoint* CenterPt2 = new SPoint[species];
int* num = new int[species];
double* sum1 = new double[species];
double* sum2 = new double[species];
double* sum3 = new double[species];
for (int x = 0; x < cols; x++) {//把图像信息读入
for (int y = 0; y < rows; y++) {
Points[x + y * cols].v1 = src.at<Vec3b>(y, x)[0]; // blue
Points[x + y * cols].v2 = src.at<Vec3b>(y, x)[1]; // green
Points[x + y * cols].v3 = src.at<Vec3b>(y, x)[2]; // red
Points[x + y * cols].id = 0;
}
}
for (int i = 0; i < species; i++) {//初始化中心
CenterPt1[i] = Points[i * 10];
CenterPt1[i].id = i;
CenterPt2[i] = Points[i * 10];
CenterPt2[i].id = i;
num[i] = 0;
sum1[i] = sum2[i] = sum3[i] = 0;
}
double min = 20000000;
double error = 0;
while (iterations > 0) {
for (int x = 0; x < cols; x++) {//把图像信息读入
for (int y = 0; y < rows; y++) {
for (int i = 0; i < species; i++) {
int id = CenterPt1[i].id;
if (Norm2(Points[x + y * cols], CenterPt1[i]) < min) {
Points[x + y * cols].id = id;
min = Norm2(Points[x + y * cols], CenterPt1[i]);
num[id]++;
sum1[id] += Points[x + y * cols].v1;
sum2[id] += Points[x + y * cols].v2;
sum3[id] += Points[x + y * cols].v3;
}
}
min = 20000000;
}
}
iterations--;for (int i = 0; i < species; i++) {
CenterPt2[i].v1 = sum1[i] / num[i];
CenterPt2[i].v2 = sum2[i] / num[i];
CenterPt2[i].v3 = sum3[i] / num[i];
error += pow((CenterPt2[i].v1 - CenterPt1[i].v1), 2) + pow((CenterPt2[i].v2 - CenterPt1[i].v2), 2) + pow((CenterPt2[i].v3 - CenterPt1[i].v3), 2);
}
if (error < 0.0001)
break;
for (int i = 0; i < species; i++) {
CenterPt1[i].v1 = CenterPt2[i].v1;
CenterPt1[i].v2 = CenterPt2[i].v2;
CenterPt1[i].v3 = CenterPt2[i].v3;
num[i] = 0;
sum1[i] = sum2[i] = sum3[i] = 0;
}
error = 0;
}//分类结束
for (int x = 0; x < cols; x++) {//上色
for (int y = 0; y < rows; y++) {
Color(&Points[x + y * cols], Points[x + y * cols].id);
}
}
for (int x = 0; x < cols; x++) {//显示
for (int y = 0; y < rows; y++) {
K.at<Vec3b>(y, x)[0] = (int)Points[x + y * cols].v1;
K.at<Vec3b>(y, x)[1] = (int)Points[x + y * cols].v2;
K.at<Vec3b>(y, x)[2] = (int)Points[x + y * cols].v3;
}
}
//显示图片
imshow("Output", K);
//不加此语句图片会一闪而过
waitKey(0);
delete[] Points, CenterPt1, CenterPt2, num, sum1, sum2, sum3;
return K;
}int main()
{
//读取图片（使用图片的相对路径）
Mat src = imread("ik_beijing_c.bmp");
Kmean(src, 7, 3);
return 0;
}/*
int b=img.at<Vec3b>(y,x)[0]; // blue
int g=img.at<Vec3b>(y,x)[1]; // green
int r=img.at<Vec3b>(y,x)[2]; // red
*/
分类效果图
由于迭代次数有限，可以看到并不收敛，另外初始点的选择也不是很好，需要改进。