本文是作者隨筆，附帶學習的註釋：

主要函式：輸入是矩陣影象(例如512*512的單通道影象)，返回方式是形參返回

void ClustFCM(cv::InputArray image, cv::OutputArray U,cv::OutputArray obj_fcn, cv::OutputArray center,int cluster_n, float * option)
{//輸出：
//center----聚類中心 2x1的矩陣
//U----隸屬度矩陣 矩陣
//obj_fcn----目標函式值矩陣
Mat data = image.getMat();
int data_n = data.rows;//樣本數
int in_n = data.cols;//特徵維度

int expo = option[0];      //隸屬度矩陣U的指數
int max_iter = option[1];// 最大迭代次數
float min_impro = option[2];// 隸屬度最小變化量, 迭代終止條件
int display = option[3];// 每次迭代是否輸出資訊標誌
initfcm(cluster_n, data_n, U); //輸出U，初始化隸屬度矩陣（子函式1）

Mat obj_fcn_begin = Mat::zeros(max_iter, 1, CV_32F);//初始化目標函式矩陣

double obj_fcn_ones;//此處是一個數，用來填充目標函式初始化矩陣
int iter_n; //實際迭代次數
for (int i = 1;  i <= max_iter;  i++)
{
//[U, center, obj_fcn(i)] = stepfcm(data, U, cluster_n, expo);
//在第k步迴圈中改變聚類中心ceneter,和分配函式U的隸屬度值;
stepfcm(data,cluster_n, expo, &obj_fcn_ones,center,U); //聚類過程：返回修改後的目標函式、聚類中心，隸屬度矩陣（子函式2）
obj_fcn_begin.at<float>(i-1, 0) = obj_fcn_ones;//填充目標函式初始化矩陣
if (i > 1)
{
if (abs(obj_fcn_begin.at<float>(i-1, 0) - obj_fcn_begin.at<float>(i - 2, 0))< min_impro)
{
iter_n = i;
break;
}
}
}
Mat obj_fcn_end = obj_fcn_begin(Range(0, iter_n), Range::all());//根據迭代次數，丟棄目標函式矩陣多餘0位
obj_fcn_end.copyTo(obj_fcn); //返回Mat類形引數據

}

子函式1：初始化函式，返回初始化隸屬度矩陣

void initfcm(int cluster_n, int data_n, cv::OutputArray U)
{
Mat U_temp(cluster_n, data_n, CV_32F);
randu(U_temp, Scalar::all(0), Scalar::all(1));//Scalar:;all是4空間內容全為N，可以用來當作常數使用
//randu(U_temp, Scalar(0), Scalar(1)); 通過這裡可以發現問題：應學會查閱官方文件,randu的上下限是Vector或者Mat結構，randu產生隨機矩陣，用於初始化
Mat U1,U2;
U1 = U_temp(Range(0, 1), Range::all()).clone();// Range::all() 全部
U2 = U_temp(Range(1, 2), Range::all()).clone();
U1 = U1 + U2;
Mat U3;
repeat(U1, 2,1, U3);//垂直方向複製1次,水平方向複製0次,注意是i+1 與i+0，返回複製擴容後的U3
U_temp = U_temp.mul(1 / U3);
U_temp.copyTo(U);//返回mat類形參

}

子函式2：聚類過程，返回更新後的聚類中心 隸屬度矩陣U 目標函式，注意這裡的目標函式是一個數，需在主函式裡填充到矩陣中

void stepfcm(cv::InputArray data, int cluster_n, int expo,double *obj_fcn_one, cv::OutputArray center, cv::OutputArray U)
{
Mat temp;
Mat mf;
pow(U.getMat(), expo, mf); // C++中矩陣每一個元素求階乘

reduce(mf, temp, 1, CV_REDUCE_SUM);//矩陣行列求和的操作，轉化成了2*1的列向量 dim=0，單行，dim=1，單列，CE_SUM 行列之和，AVG:行列向量的均值向量,MAX,MIN輸出行列的最大最小值

transpose(temp, temp);
Mat res = (Mat::ones(data.getMat().cols, 1, CV_32F)) * temp;
transpose(res, res);
Mat result = (mf * data.getMat()) / res;
result.copyTo(center);


//計算樣本點距離聚類中心的距離
//out = zeros(size(center, 1), size(data, 1));
Mat dist = Mat::zeros(result.rows, data.getMat().rows, CV_32F);
for (int k = 0; k < result.rows; k++)
{
Mat temp = abs(data.getMat() - Mat::ones(data.getMat().rows, 1, CV_32F) * result(Range(k, k + 1), Range::all()));
transpose(temp, temp);
//dist(Range(k, k + 1), Range::all()) = temp;//dist(Range(k, k + 1), Range::all()) = temp;這樣的矩陣賦值方式是無效的。
dist(Range(k, k + 1), Range::all()) = dist(Range(k, k + 1), Range::all()) + temp;
}
// 計算目標函式
//obj_fcn = sum(sum((dist. ^ 2).*mf));
*obj_fcn_one = sum(dist.mul(dist).mul(mf))[0];  //sum返回的是Scalar類，是一個4空間向量
// 計算距離聚類中心的距離
Mat tmp;
pow(dist, (-2 / (expo - 1)), tmp); //void pow(InputArray src, double power, OutputArray dst);
Mat Tmp;
reduce(tmp,Tmp, 0, CV_REDUCE_SUM);
// 新的隸屬度矩陣
Mat U_new = tmp.mul(1 / ((Mat::ones(cluster_n, 1, CV_32F)) * Tmp));
U_new.copyTo(U);
}