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
m
m,聚类簇数量
K
K
K,最大迭代次数
m
a
x
_
i
t
e
r
max\_iter
max_iter,邻域边长
k
k
k,迭代停止条件
e
r
r
o
r
error
error
3、为输入图像开辟动态二维数组内存空间,将Mat类型的输入图像转换为动态二维数组
4、为边缘扩展的输入图像开辟动态二维数组内存空间,用0对输入图像作边缘扩展,行数和列数分别增加
k
−
1
k-1
k−1
5、为隶属度数组开辟动态三维数组内存空间,并用0到1之间均匀分布的随机数填充该数组
6、令隶属度三维数组所有元素除以数组第三维之和
7、为边缘扩展的隶属度数组开辟动态三维数组内存空间,用0对隶属度数组作边缘扩展,行数和列数分别增加
k
−
1
k-1
k−1
8、开辟模糊因子G开辟二维动态数组内存空间
9、开辟聚类中心和目标函数动态一维数组内存,其中使用calloc函数令目标函数数组初始化为0
10、设置
i
t
e
r
=
0
iter = 0
iter=0
以下公式中,
N
N
N是像素总数,
y
j
y_j
yj是输入图像第
j
j
j个像素,
u
i
j
u_{ij}
uij是第
j
j
j个像素对第
i
i
i个聚类簇的隶属度,
N
j
N_j
Nj是第下标
j
j
j的
k
×
k
k\times k
k×k邻域
11、计算聚类中心
c
i
=
∑
j
=
1
N
y
j
u
i
j
m
∑
j
=
1
N
u
i
j
m
c_{i}=\frac{\sum^N_{j=1}y_ju_{ij}^m}{\sum^N_{j=1}u_{ij}^m}
ci=∑j=1Nuijm∑j=1Nyjuijm
12、计算模糊因子G
G
i
j
=
∑
k
∈
N
j
1
d
k
j
+
1
(
1
−
u
i
k
)
m
∣
∣
y
k
−
c
i
∣
∣
2
G_{ij}=\sum_{k\in N_j}\frac{1}{d_{kj}+1}(1-u_{ik})^m||y_k-c_i||^2
Gij=k∈Nj∑dkj+11(1−uik)m∣∣yk−ci∣∣2
13、计算隶属度U
u
i
j
=
1
∑
r
=
1
K
(
∣
∣
y
j
−
c
i
∣
∣
2
+
G
i
j
∣
∣
y
j
−
c
r
∣
∣
2
+
G
r
j
)
1
m
−
1
u_{ij}=\frac{1}{\sum^K_{r=1}(\frac{||y_j-c_i||^2+G_{ij}}{||y_j-c_r||^2+G_{rj}})^{\frac{1}{m-1}}}
uij=∑r=1K(∣∣yj−cr∣∣2+Grj∣∣yj−ci∣∣2+Gij)m−111
14、计算目标函数J
J
=
∑
i
=
1
K
∑
j
=
1
N
(
u
i
j
m
∣
∣
y
j
−
c
i
∣
∣
2
+
G
i
j
)
J=\sum_{i=1}^K\sum_{j=1}^N(u_{ij}^m||y_j-c_i||^2+G_{ij})
J=i=1∑Kj=1∑N(uijm∣∣yj−ci∣∣2+Gij)
15、判断本次迭代的
J
J
J与上次的
J
J
J之差的绝对值是否小于等于
e
r
r
o
r
error
error且
i
t
e
r
iter
iter是否大于1
若是,则目标函数收敛;若否,则
i
t
e
r
=
i
t
e
r
+
1
iter = 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;
}