计算分形盒子维

//************************//
//计算分形盒子维
//*** yangxin_szu 2013_03_28 ***//
//valarray与 MFC 有一定冲突
//#undef的使用是为了避免问题出现
#ifdef min
#undef min
#endif
#ifdef max
#undef max
#endif
#include <valarray>
using namespace std;
void Calculate_Fractal_Dim(unsigned char* Img_Data, int Img_size, double& kx, double& b)
{//*************************************************////the width and height of the Image should be the same//gray level : 256//points for Least Square method: 20//*************************************************//int i = 0, j = 0;//分形盒子维相关参数int M = Img_size;//图像尺寸int G = 256;//图像灰度级int L_point = 20;//最小二乘法样本点数valarray<unsigned char> Img_VAL(Img_size*Img_size);valarray<int> L_VAL(L_point);valarray<int> h_VAL(L_point);valarray<double> r_VAL(L_point);valarray<double> Nr_VAL((double)0, L_point);valarray<int> Grid_Num(L_point);//网格数目valarray<double> fractal_D(0.0, L_point);//复制数据int k = 0;for (i = 0; i < Img_size; i++){for (j = 0; j < Img_size; j++){Img_VAL[k] = Img_Data[i*Img_size + j];k++;}}//网格大小及相关参数//改进后的约束范围 M^(1/3)<= L <= M/3int L_min = (int)powf(M, 1 / 3.0);//int L_min = M/40;int L_max = (int)M / 2;int L_step = (int)((L_max - L_min) / (float)L_point);for (i = 0; i < L_point; i++){L_VAL[i] = L_min + i * L_step;//各样本点对应的网格大小     Lh_VAL[i] = (G*L_VAL[i]) / M;  //各样本点对应的盒子高度     hr_VAL[i] = log10(1 / (L_VAL[i] / (float)M));//各样本点对应的 rGrid_Num[i] = M / L_VAL[i];//各样本点对应的图像网格数目    Num}int m = 0, n = 0, t = 0;int grid_lt_x = 0, grid_lt_y = 0, grid_rd_x = 0, grid_rd_y = 0;unsigned char grid_I_max = 0, grid_I_min = 0;int dbc_l = 0, dbc_k = 0, nr = 0;int s = 0, p = 0, q = 0, box_non_zero = 0, box_zero_count = 0;int gray_k = 0, gray_l = 0;//计算 Nrfor (k = 0; k < L_point; k++){//临时存储单个网格数据valarray<unsigned char> grid_img((unsigned char)0, L_VAL[k] * L_VAL[k]);for (m = 0; m < Grid_Num[k]; m++){for (n = 0; n < Grid_Num[k]; n++){//单个网格的坐标范围grid_lt_x = n * L_VAL[k];grid_lt_y = m * L_VAL[k];grid_rd_x = grid_lt_x + L_VAL[k] - 1;grid_rd_y = grid_lt_y + L_VAL[k] - 1;//复制数据t = 0;for (i = grid_lt_y; i < grid_rd_y; i++){for (j = grid_lt_x; j < grid_rd_x; j++){grid_img[t] = Img_Data[i*M + j];t++;}}grid_I_min = grid_img.min();grid_I_max = grid_img.max();//最小、最大灰度所在的网格高度dbc_k = grid_I_min / h_VAL[k];dbc_l = grid_I_max / h_VAL[k];//*******计算空盒子数目********//for (s = dbc_k; s <= dbc_l; s++){//灰度上下限gray_k = s * h_VAL[k];gray_l = gray_k + h_VAL[k] - 1;//清零box_non_zero = 0;//落在指定盒子内的点数for (p = 0; p < t - 1; p++){if ((grid_img[p] >= gray_k) && (grid_img[p] <= gray_l))box_non_zero++;}//空盒子if (box_non_zero == 0)box_zero_count++;}//Nr累加并剔除空盒子nr = dbc_l - dbc_k + 1 - box_zero_count;Nr_VAL[k] = Nr_VAL[k] + nr;//清零box_zero_count = 0;}}Nr_VAL[k] = log10(double(Nr_VAL[k]));//grid_img.free();}//计算各样本点对应的理论斜率并保存fractal_D = Nr_VAL / r_VAL;//ofstream outfile("F:\\Fractal_D.txt");//打开文件，准备写入for (j = 0; j < L_point; j++){//cout << fractal_D[j] << ' ' << endl;//距离}//cout << endl;//outfile.close();//关闭文件，完成写入//****************************************////最小二乘法拟合直线double A = 0.0;double B = 0.0;double C = 0.0;double D = 0.0;A = (r_VAL*r_VAL).sum();B = r_VAL.sum();C = (r_VAL*Nr_VAL).sum();D = Nr_VAL.sum();double fractal_k, fractal_b, tmp = 0;if (tmp = (A*L_point - B * B)){fractal_k = (C*L_point - B * D) / tmp;fractal_b = (A*D - C * B) / tmp;}else{fractal_k = 1;fractal_b = 0;}kx = fractal_k;b = fractal_b;
}

调用方法

Mat src;double k, b;src = cv::imread("1.bmp");cvtColor(src, src, CV_BGR2GRAY);threshold(src, src, 128, 255, CV_THRESH_BINARY);unsigned char* Small = src.data;Calculate_Fractal_Dim(Small, 314, k, b);
cout<<k<< b<<endl;//或转换string显示：to_string(k);

注：

1.本方案原意是用来评估一个只有白色点轮廓(黑色背景)的图像中白色轮廓的分布均匀程度。
2.如果各个轮廓的分布足够均匀，则，计算的盒维数中K值月接近2.0；
3.计算要求输入图像是正方形；
4.原始代码是从网上下载，选优对输入图像做适当转换之后才可以使用；
5.计算结果，理论要求是拟合成一条一次曲线y=kx+b;其各个数据分布在直线附近。越均匀越接近在直线附近。