(OpenCV) Fractal Dimension

ボックスカウント法 (box counting method)でフラクタル次元(Fractal dimension)を求める

やり方

ボックスカウント法については解説しているページがあちこちにあるので、そちらを参照。 手順は以下の通りとなる。

1. ROIの設定(大きさe:2の階乗、1channel、ROIの大きさを1として1画素は1/eの大きさになる)
2. ROIを2値化
3. "白"の画素数N(1/e)をカウント(log10する)
4. ROIを1/2にリサイズ
5. 2-4を分割可能な限り繰り返す
6. log10(1/e)に対するlog10(N(1/e))の近似直線の傾きがフラクタル次元(負になるが絶対値とする)

2,3をmyFractalDimCore()、4-6をmyFractalDim()で行なう。

かなりおおざっぱだが多分あってる。

myFractalDimCore

double myFractalDimCore(IplImage *src, int eps)
{
  long cores;
  IplImage *resizedImage = cvCreateImage(cvSize(eps,eps),IPL_DEPTH_8U,1);
  IplImage *binImage = cvCreateImage(cvSize(eps,eps),IPL_DEPTH_8U,1);

  cvResize(src,resizedImage,CV_INTER_LINEAR);
  cvThreshold(resizedImage,binImage,128,255,CV_THRESH_BINARY);

  cores=(long)cvSum(binImage).val[0];
  cvReleaseImage(&resizedImage);
  cvReleaseImage(&binImage);
  return log10(cores==0?1.0:(double)cores/255.0);
}

※追記：カウント数が0だと-∞になりコケるので便宜上log10(1.0)を返すようにしてみた。coresはlong型に。

myFractalDim

voidとなっているがB[0]を返せばフラクタル次元を得る関数になる。

近似直線を求めるのは面倒なのでLAPACKの最小二乗法関数dgels_()を用いた。

void myFractalDim(IplImage *src)
{
  //マスクとしてROIを使うので事前に設定しておく必要あり
  int size=cvGetImageROI(src).height;
  int m = (int)(log(size)/log(2));
  int n=2; //a,b
  int lda= m;
  int ldb= m;
  int info,nrhs=1;
  int ldwork=m+n;
  int x;
  int i;
  double *A = (double *)malloc(sizeof(double)*m*n);
  double *B = (double *)malloc(sizeof(double)*m);
  double *work = (double *)malloc(sizeof(double)*(m+n));
  for(i=0,x=2; x<=size; x*=2,i++){
    A[i]=log10(1.0/x);
    A[i+m]=1;
    B[i]=myFractalDimCore(src,x);
    //for debug
    printf("%lf,%lf\n",A[i],B[i]);
  }
  //CLAPACK
  dgels_("N",&m,&n,&nrhs,A,&lda,B,&ldb,work,&ldwork,&info);
  printf("a,b=%lf,%lf\n",B[0],B[1]);
}