自己手动实现了meanshift，算是把之前还有些模糊的地方给搞清楚了，本来应该半天就可以的，结果高斯核函数计算权重那里弄错了，多花了好一阵子功夫才搞定。

原理

对于所有样本点中的一个点x，计算它的shift vector，沿着shift vector移动，然后再计算x在移动后的位置上的新shift vector，再移动，直到最后到达终点了才算完。然后对所有点都进行一次这样的操作。
如果不使用核函数的话，就要规定一个半径，然后根据半径内点的mean来计算shift vector；如果使用高斯核函数的话就不需要了，因为高斯核函数的效果就是让离得近的点权重大，离得远的点权重很小。

加载数据

用了和上篇kmeans相同的数据

void MeanShiftClass::LoadData(string data_path, vector<MeanShift>& dataset)
{float arr1[64];float arr2[64];int i = 0;ifstream myfile(data_path);if (!myfile) {cout << "Unable to open file";exit(1); // terminate with error  }else{char str[64] = { 0 };while (!myfile.eof()){myfile.getline(str, 64); //读取一行数据sscanf_s(str, "%f %f", &arr1[i], &arr2[i]);i++;}}MeanShift tmp;for (int j = 0; j < i; j++){tmp.pos.x = arr1[j];tmp.pos.y = arr2[j];tmp.res.x = arr1[j];tmp.res.y = arr2[j];dataset.push_back(tmp);}
}

均值漂移

对于样本点x，计算均值漂移向量一次，当漂移距离小于EPSILON时就停止，计算下一个点。

void MeanShiftClass::ShiftOnce(MeanShift& p)
{float x_sum = 0;float y_sum = 0;float weight_sum = 0;for (int i = 0; i < point_num; i++){float tmp_distance = GetEuclideanDistance(p.res, dataset[i].pos);float weight = GaussianKernel(tmp_distance, kernel_bandwidth);x_sum += dataset[i].pos.x * weight;y_sum += dataset[i].pos.y * weight;weight_sum += weight;}Point2f shift_vector(x_sum/ weight_sum, y_sum/ weight_sum);float shift_distance = GetEuclideanDistance(p.res, shift_vector);cout << "shift_distance = " << shift_distance << endl;if (shift_distance < EPSILON)stop = true;p.res = shift_vector;
}

void MeanShiftClass::DoMeanShiftCluster()
{for (int i = 0; i < point_num; i++){stop = false;while (!stop){Sleep(20);ShiftOnce(dataset[i]);}cout <<"("<< dataset[i].res.x << "," << dataset[i].res.y<<")" << endl;}LabelClusters();
}

给聚类结果打标签

聚完类之后给他们打上标签，聚类结果距离小于CLUSTER_EPSILON的点认为是一类。

void MeanShiftClass::LabelClusters()
{int current_label = -1;for (int i = 0; i < point_num; i++){if (dataset[i].label != -1)continue;current_label++;for (int j = 0; j < point_num; j++){if (GetEuclideanDistance(dataset[i].res, dataset[j].res) < CLUSTER_EPSILON){dataset[j].label = current_label;}}}
}

结果

这是聚类前的数据点

这是聚类结果

实现了对表盘上多位数字的合并，方便后续的处理。原始数据点是数字方框的位置坐标。

源码

MeanShiftClass.h

#pragma once
#include <iostream>
#include <opencv2/opencv.hpp>
#include <windows.h>
#define PI 3.1415926535898
#define EPSILON 0.001
#define CLUSTER_EPSILON 10using namespace cv;
using namespace std;struct MeanShift
{Point2f pos;Point2f res;int label = -1;
};class MeanShiftClass
{
public:MeanShiftClass(int _bandwidth);~MeanShiftClass();void DoMeanShiftCluster();
private:bool stop;int point_num;int cluster_num;int kernel_bandwidth;vector<MeanShift> dataset;void LabelClusters();void ShowClusterResult();void ShiftOnce(MeanShift& p);int GetManhattanDistance(Point2f p0, Point2f p1);float GaussianKernel(int distance, int bandwidth);float GetEuclideanDistance(Point2f p0, Point2f p1);void LoadData(string data_path, vector<MeanShift> &dataset);
};

MeanShiftClass.cpp

#include "pch.h"
#include "MeanShiftClass.h"MeanShiftClass::MeanShiftClass(int _bandwidth)
{LoadData("DigitDetcResult_vec2.txt", dataset);point_num = dataset.size();kernel_bandwidth = _bandwidth;stop = false;
}MeanShiftClass::~MeanShiftClass()
{
}void MeanShiftClass::DoMeanShiftCluster()
{for (int i = 0; i < point_num; i++){stop = false;while (!stop){Sleep(20);ShiftOnce(dataset[i]);}cout <<"("<< dataset[i].res.x << "," << dataset[i].res.y<<")" << endl;}LabelClusters();ShowClusterResult();
}void MeanShiftClass::ShowClusterResult()
{Mat src = imread("2.jpg");for (int i = 0; i < dataset.size(); i++){circle(src, dataset[i].pos, 3, Scalar(0, 0, 255), -1);}imshow("src", src);waitKey(0);
}void MeanShiftClass::LoadData(string data_path, vector<MeanShift>& dataset)
{float arr1[64];float arr2[64];int i = 0;ifstream myfile(data_path);if (!myfile) {cout << "Unable to open file";exit(1); // terminate with error  }else{char str[64] = { 0 };while (!myfile.eof()){myfile.getline(str, 64); //读取一行数据sscanf_s(str, "%f %f", &arr1[i], &arr2[i]);i++;}}MeanShift tmp;for (int j = 0; j < i; j++){tmp.pos.x = arr1[j];tmp.pos.y = arr2[j];tmp.res.x = arr1[j];tmp.res.y = arr2[j];dataset.push_back(tmp);}
}int MeanShiftClass::GetManhattanDistance(Point2f p0, Point2f p1)
{return abs(p0.x - p1.x) + abs(p0.y - p1.y);
}float MeanShiftClass::GetEuclideanDistance(Point2f p0, Point2f p1)
{return sqrt((p0.x - p1.x)*(p0.x - p1.x) + (p0.y - p1.y)*(p0.y - p1.y));
}float MeanShiftClass::GaussianKernel(int distance, int bandwidth)
{return exp(-0.5*(distance*distance) / (bandwidth*bandwidth));
}void MeanShiftClass::ShiftOnce(MeanShift& p)
{float x_sum = 0;float y_sum = 0;float weight_sum = 0;for (int i = 0; i < point_num; i++){float tmp_distance = GetEuclideanDistance(p.res, dataset[i].pos);float weight = GaussianKernel(tmp_distance, kernel_bandwidth);x_sum += dataset[i].pos.x * weight;y_sum += dataset[i].pos.y * weight;weight_sum += weight;}Point2f shift_vector(x_sum/ weight_sum, y_sum/ weight_sum);float shift_distance = GetEuclideanDistance(p.res, shift_vector);cout << "shift_distance = " << shift_distance << endl;if (shift_distance < EPSILON)stop = true;p.res = shift_vector;
}void MeanShiftClass::LabelClusters()
{int current_label = -1;for (int i = 0; i < point_num; i++){if (dataset[i].label != -1)continue;current_label++;for (int j = 0; j < point_num; j++){if (GetEuclideanDistance(dataset[i].res, dataset[j].res) < CLUSTER_EPSILON){dataset[j].label = current_label;}}}
}

记录一下~