OpenCV图像处理（3）——盒维数计算

计算分形盒子维

//************************//
//计算分形盒子维
//*** 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/3
 int 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;//各样本点对应的网格大小     L
  h_VAL[i] = (G*L_VAL[i]) / M;  //各样本点对应的盒子高度     h
  r_VAL[i] = log10(1 / (L_VAL[i] / (float)M));//各样本点对应的 r
  Grid_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;
 //计算 Nr
 for (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;其各个数据分布在直线附近。越均匀越接近在直线附近。