【C++】【图像分割】【模糊聚类】FLICM算法C++实现，OpenCV负责图像读取与显示

S. Krinidis, V. Chatzis. A Robust Fuzzy Local Information C-Means Clustering Algorithm[J]. IEEE Transactions on Image Processing, 19(5), 2010: 1328-1337.
该算法推导有误，使目标函数无法最小化。
正确的推导请参考：
T. Celik, H. K. Lee. Comments on “A Robust Fuzzy Local Information C-Means Clustering Algorithm”[J]. IEEE Transactions on Image Processing, 22(3), 2013: 1258-1261.
程序按原始论文实现

FLICM算法步骤

1、通过OpenCV imread函数读入待分割图像，为Mat数据类型
2、设置聚类参数如下：
隶属度指数$m$，聚类簇数量$K$，最大迭代次数$max\_iter$，邻域边长$k$，迭代停止条件$error$
3、为输入图像开辟动态二维数组内存空间，将Mat类型的输入图像转换为动态二维数组
4、为边缘扩展的输入图像开辟动态二维数组内存空间，用0对输入图像作边缘扩展，行数和列数分别增加$k-1$
5、为隶属度数组开辟动态三维数组内存空间，并用0到1之间均匀分布的随机数填充该数组
6、令隶属度三维数组所有元素除以数组第三维之和
7、为边缘扩展的隶属度数组开辟动态三维数组内存空间，用0对隶属度数组作边缘扩展，行数和列数分别增加$k-1$
8、开辟模糊因子G开辟二维动态数组内存空间
9、开辟聚类中心和目标函数动态一维数组内存，其中使用calloc函数令目标函数数组初始化为0
10、设置$iter=0$
以下公式中，$N$是像素总数，$y_j$是输入图像第$j$个像素，$u_{ij}$是第$j$个像素对第$i$个聚类簇的隶属度，$N_j$是第下标$j$的$k\times k$邻域
11、计算聚类中心
$$c_i=\frac{{\sum }_{j=1}^N{y}_j{u}_{ij}^m}{{\sum }_{j=1}^N{u}_{ij}^m}$$
12、计算模糊因子G
$$G_{ij}=\sum _{k\in N_j}\frac{1}{d_{kj}+1}(1-u_{ik}{)}^m||y_k-c_i|{|}^2$$
13、计算隶属度U
$$u_{ij}=\frac{1}{{\sum }_{r=1}^K(\frac{||y_j-c_i|{|}^2+G_{ij}}{||y_j-c_r|{|}^2+G_{rj}}{)}^{\frac{1}{m-1}}}$$
14、计算目标函数J
$$J=\sum _{i=1}^K\sum _{j=1}^N(u_{ij}^m||y_j-c_i|{|}^2+G_{ij})$$
15、判断本次迭代的$J$与上次的$J$之差的绝对值是否小于等于$error$且$iter$是否大于1
若是，则目标函数收敛；若否，则$iter=iter+1$，返回第第11步。
16、输出隶属度数组，按最大隶属度原则赋予每个像素聚类簇标号。
17、按聚类中心的灰度级着色。
实验结果：

代码

#include <opencv2\opencv.hpp>
#include <stdio.h>
#include <iostream>
using namespace cv;
using namespace std;
int main()
{
	Mat input = imread("F:\\C++语言程序\\模糊聚类\\Images\\test1_noise.tif", 0);
	const int rows = input.rows;
	const int cols = input.cols;
	namedWindow("输入");
	imshow("输入", input);
	// 输入参数
	int m = 2, K = 2, max_iter = 100, side = 3;
	double error = 0.001;
	// 为输入图像和边缘扩展的输入图像开辟内存空间，并赋值，用二维数组表示
	const int half_k = side / 2;
	const int padded_rows = rows + half_k * 2;
	const int padded_cols = cols + half_k * 2;
	double** input_ptr = (double**)malloc(rows * sizeof(double*));
	if (input_ptr == NULL) { cout << "内存分配失败\n"; return 0; }
	double** input_padded = (double**)calloc(padded_rows, sizeof(double*));
	if (input_padded == NULL) { cout << "内存分配失败\n"; return 0; }
	for (int i = 0; i < rows + 2 * half_k; i++)
	{
		input_padded[i] = (double*)calloc(padded_cols, sizeof(double));
		if (input_padded[i] == NULL) { cout << "内存分配失败\n"; return 0; }
	}
	for (int i = 0; i < rows; i++)
	{
		input_ptr[i] = (double*)malloc(cols * sizeof(double));
		if (input_ptr[i] == NULL) { cout << "内存分配失败\n"; return 0; }
	}
	for (int i = 0; i < rows; i++)
		for (int j = 0; j < cols; j++)
		{
			input_ptr[i][j] = (double)input.at<uchar>(i, j);
			input_padded[i + half_k][j + half_k] = input_ptr[i][j];
		}
	input.release();
	// 开辟隶属度U和模糊因子G的内存，并对U随机初始化
	double*** U = (double***)malloc(rows * sizeof(double**));
	if (U == NULL) { cout << "内存分配失败\n"; return 0; }
	double*** G = (double***)malloc(rows * sizeof(double**));
	if (G == NULL) { cout << "内存分配失败\n"; return 0; }
	double*** U_padded = (double***)calloc(padded_rows, sizeof(double**));
	if (U_padded == NULL) { cout << "内存分配失败\n"; return 0; }
	for (int i = 0; i < rows + 2 * half_k; i++)
	{
		U_padded[i] = (double**)calloc(padded_cols, sizeof(double*));
		if (U_padded[i] == NULL) { cout << "内存分配失败\n"; return 0; }
		for (int j = 0; j < cols + 2 * half_k; j++)
		{
			U_padded[i][j] = (double*)calloc(K, sizeof(double));
			if (U_padded[i][j] == NULL) { cout << "内存分配失败\n"; return 0; }
		}
	}
	srand((int)time(0));
	for (int i = 0; i < rows; i++)
	{
		U[i] = (double**)malloc(cols * sizeof(double*));
		if (U[i] == NULL) { cout << "内存分配失败\n"; return 0; }
		G[i] = (double**)malloc(cols * sizeof(double*));
		if (G[i] == NULL) { cout << "内存分配失败\n"; return 0; }
		for (int j = 0; j < cols; j++)
		{
			U[i][j] = (double*)malloc(K * sizeof(double));
			if (U[i][j] == NULL) { cout << "内存分配失败\n"; return 0; }
			double U_sum = 0;
			for (int l = 0; l < K; l++)
			{
				U[i][j][l] = (double)rand() / RAND_MAX;
				U_sum += U[i][j][l];
			}
			for (int l = 0; l < K; l++)
			{
				U[i][j][l] /= U_sum;
				U_padded[i + half_k][j + half_k][l] = U[i][j][l];
			}
			G[i][j] = (double*)malloc(K * sizeof(double));
			if (G[i][j] == NULL) { cout << "内存分配失败\n"; return 0; }
		}
	}
	// 开辟聚类中心center的内存
	double* center = (double*)malloc(K * sizeof(double));
	if (center == NULL) { cout << "内存分配失败\n"; return 0; }
	// 开辟目标函数J的内存
	double* J = (double*)calloc(max_iter, sizeof(double));
	if (J == NULL) { cout << "内存分配失败\n"; return 0; }
	// 聚类
	for (int iter = 0; iter < max_iter; iter++)
	{
		for (int l = 0; l < K; l++)
		{
			// 更新聚类中心
			double center_up, center_down;
			center_up = center_down = 0;
			for (int i = 0; i < rows; i++)
				for (int j = 0; j < cols; j++)
				{
					center_up += input_ptr[i][j] * pow(U[i][j][l], m);
					center_down += pow(U[i][j][l], m);
				}
			center[l] = center_up / center_down;
			// 计算模糊因子G
			for (int i = 0; i < rows; i++)
				for (int j = 0; j < cols; j++)
				{
					double d, U_res, G_temp = 0;
					for (int r = -half_k; r <= half_k; r++)
						for (int n = -half_k; n <= half_k; n++)
						{
							if (r == 0 && n == 0)
								continue;
							d = 1.0 / (1.0 + sqrt(r * r + n * n));
							U_res = pow(1 - U_padded[i + half_k + r][j + half_k + n][l], m);
							G_temp = G_temp + d * U_res * pow(input_padded[i + half_k + r][j + half_k + n] - center[l], 2);
						}
					G[i][j][l] = G_temp;
				}
		}
		for (int i = 0; i < rows; i++)
			for (int j = 0; j < cols; j++)
			{
				// 更新隶属度
				double U_down = 0.0, d, U_up;
				for (int l = 0; l < K; l++)
				{
					d = pow(input_ptr[i][j] - center[l], 2);
					U_down = U_down + pow(d + G[i][j][l], -1 / (m - 1));
				}
				for (int l = 0; l < K; l++)
				{
					d = pow(input_ptr[i][j] - center[l], 2);
					U_up = pow(d + G[i][j][l], 1 / (m - 1));
					U[i][j][l] = 1 / (U_up * U_down);
					// 更新边缘扩展后的隶属度U_padded
					U_padded[i + half_k][j + half_k][l] = U[i][j][l];
					// 更新目标函数
					J[iter] += pow(U[i][j][l], m) * pow(input_ptr[i][j] - center[l], 2) + G[i][j][l];
				}
			}
		printf("第%d次迭代，J = %.3f\n", iter + 1, J[iter]);
		// 迭代条件判定
		if (abs(J[iter] - J[iter - 1]) <= error && iter >= 1)
		{
			cout << "目标函数已收敛\n";
			break;
		}
		else if (iter == max_iter - 1 && iter >= 1)
		{
			cout << "目标函数未收敛，达到最大迭代次数\n";
			break;
		}
	}
	// 内存清理
	free(G);
	free(J);
	free(input_ptr);
	free(input_padded);
	free(U_padded);
	J = NULL;
	input_ptr = input_padded = NULL;
	U_padded = G = NULL;
	// 显示聚类中心
	cout << "聚类中心分别为\n";
	for (int l = 0; l < K; l++)
		cout << center[l] << "\n";
	// 计算隶属度划分系数和划分熵
	double Vpc, Vpe;
	Vpc = Vpe = 0;
	for (int i = 0; i < rows; i++)
		for (int j = 0; j < cols; j++)
			for (int l = 0; l < K; l++)
			{
				Vpc += U[i][j][l] * U[i][j][l];
				Vpe += U[i][j][l] * log(U[i][j][l]);
			}
	Vpc = Vpc / ((double)rows * (double)cols) * 100;
	Vpe = -Vpe / ((double)rows * (double)cols) * 100;
	cout << "隶属度划分系数Vpc = " << Vpc << "%" << endl;
	cout << "隶属度划分熵Vpe = " << Vpe << "%" << endl;
	// 标签数组，并着色
	short** labels = (short**)calloc(rows, sizeof(short*));
	if (labels == NULL) { cout << "内存分配失败\n"; return 0; }
	for (int i = 0; i < rows; i++)
	{
		labels[i] = (short*)calloc(cols, sizeof(short));
		if (labels[i] == NULL) { cout << "内存分配失败\n"; return 0; }
	}
	Mat result = Mat(rows, cols, CV_8U);
	uchar* colors = new uchar[K];
	for (int i = 0; i < K; i++)
		*(colors + i) = (uchar)(center[i]);
	free(center);
	center = NULL;
	for (int i = 0; i < rows; i++)
		for (int j = 0; j < cols; j++)
		{
			double U_max = 0;
			for (int l = 0; l < K; l++)
				if (U_max < U[i][j][l])
				{
					U_max = U[i][j][l];
					labels[i][j] = (short)l;
				}
			result.at<uchar>(i, j) = *(colors + labels[i][j]);
		}
	// 内存清理
	delete[] colors;
	free(labels);
	free(U);
	labels = NULL;
	U = NULL;
	// 显示分割结果
	imshow("聚类结果", result);
	waitKey(0);
	result.release();
	return 0;
}