遥感图像处理-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
*/

                                                                               分类效果图

由于迭代次数有限，可以看到并不收敛，另外初始点的选择也不是很好，需要改进。