fdutil.c

分形维数的估算软件FD3

   double
               /* FOR EACH BOX SIZE #s, negLogBoxCountS[s] WILL EVENTUALLY BE
	          SET TO THE NEGATIVE OF THE LOG (BASE TWO) OF THE NUMBER OF
		  BOXES OF SIZE #s THAT ARE OCCUPIED BY ONE OR MORE DATA
		  POINTS.
		*/

           negLogBoxCountS[numbits+1],
   
               /* FOR EACH BOX SIZE #s, logSumSqrFreqS[s] WILL EVENTUALLY BE
		  SET TO THE LOG (BASE TWO) OF THE PROBABILITY THAT TWO
		  RANDOMLY CHOSEN DATA POINTS ARE IN THE SAME BOX OF SIZE
		  #s.
	       */

           logSumSqrFreqS[numbits+1],

              /* FOR EACH BOX SIZE #s, informationS[s] WILL EVENTUALLY BE
		 SET TO THE INFORMATION (BASE TWO) IN THE DISTRIBUTION OF
		 DATA POINTS IN THE BOXES OF SIZE #s.
	      */
           informationS[numbits+1] ;
	   
{ extern int  embed_dim;         /* HOW MANY COORDINATES POINTS HAVE. */
  extern ulong **data,           /* INITIAL ARRAY OF DATA POINTS */
               dataLines ;
  extern QHeader Q[2];           /* ARRAY OF TWO QUEUES FOR THE RADIX SORT */
  extern ulong  *diff_test;      /* XOR'S OF PAIRS OF COORDINATES */
  extern long int  *next_after;  /*  ARRAY FOR POINTERS TO DATA. */

  int       bitpos, countpos, coord, queue_num, found;
  ulong     pointCount[numbits+1] ;
  long int  current, previous;
  double    sumSqrFreq[numbits+1], freq, log2=log(2.0) ;

    /* GET A POINTER TO THE FIRST DATA VALUE */
  if    (Q[0].Qfront != -1) previous = Q[0].Qfront;
  else  previous = Q[1].Qfront;

     /* INIT boxCountS, pointCountS, sumSqrFreq, AND informationS.*/
  for(countpos=0;countpos<=numbits;countpos++)
    {  boxCountS [countpos]=1;   sumSqrFreq  [countpos]=0.0;
       pointCount[countpos]=1;   informationS[countpos]=0.0; }

  for (queue_num=0;queue_num<2;queue_num++)
  { current = Q[queue_num].Qfront;
    while (current != -1)
    { found = 0 ;
        /* START BY LOOKING AT THE BIGGEST BOX SIZE */
      bitpos=numbits-1;
      for (coord=embed_dim-1;coord>=0;coord--)
        diff_test[coord] = data[coord][previous]^ data[coord][current];
      do {coord = embed_dim - 1;
          do {/* IF THE CURRENT POINT AND PREVIOUS POINTS ARE IN DIFFERENT
	         BOXES OF THIS SIZE, */
	       if ( IS_A_ONE(diff_test[coord],bitpos) )
               {/* THEN THE CURRENT POINT IS IN NEW BOXES OF ALL SMALLER
		   SIZES TOO, AND STILL IN THE SAME BOXES OF LARGER SIZES,
		   SO ... */
	        for (countpos=bitpos;countpos>=0;countpos--)
		  {/* CALCULATE FREQUENCY OF POINTS IN THE BOX, ASSUMING FOR
		      NOW THAT THE NUMBER OF DATA LINES IN THE INPUT FILE IS
		      THE NUMBER OF DISTINCT POINTS IN THE DATA SET.  WE
		      ADJUST THIS AT THE END OF THIS FUNCTION. */
                   if (debugging)
		    printf("pointCount[%d] is %d...\n", countpos, pointCount[countpos]);
		   freq = pointCount[countpos]/(double)dataLines ;
		     /* WE WILL ENCOUNTER NO MORE OF THE POINTS IN THE BOX
			WE JUST LEFT (THE SPECIAL ORDERING OF THE SORT WE
			USED ABOVE GUARANTEES THIS!), SO WE COMPUTE WHAT
			THIS BOX CONTRIBUTES TO THE RUNNING SUMS. */
                   sumSqrFreq[countpos] += (freq * freq) ;
                   informationS[countpos] += ( freq*log(freq)/log2 ) ;
 		     /* WE HAVE GOTTEN INTO A NEW BOX AT THIS LEVEL, SO WE
			REFLECT THE NEW BOX IN THE COUNT */
                   boxCountS[countpos]++ ;
                     /* SINCE WE HAVE A NEW BOX AT THIS LEVEL, THERE IS ONLY
		        ONE KNOWN POINT IN IT SO FAR -- THE CURRENT POINT */
		   pointCount[countpos]=1;
                  }
                for (countpos=bitpos+1;countpos<=numbits;countpos++)
		    /* THE CURRENT POINT IS IN THE BOXES AT THESE LEVELS, SO
		       JUST INCREMENT THE POINT COUNTER.  */
                  pointCount[countpos]++ ;
                found = 1;
               }
               else coord-- ;
             }
          while ( (found == 0) && (coord > -1) );
          bitpos-- ;
         }
      while ( (found == 0) && (bitpos > -1) );
      previous = current; current = next_after[current];
    } }
    
    /* NOW ADD IN THE CONTRIBUTION DUE TO THE COUNTS REMAINING AFTER THE
       LAST POINT HAS BEEN FOUND, RENORMALIZE WITH BOXCOUNT[0], AND MASSAGE
       THE RAW DATA FROM THE TRAVERSAL SO THAT IS IS READY FOR THE LEAST
       SQUARES FITTING. */
       
  for (countpos=numbits;countpos>=0;countpos--)
  {
    negLogBoxCountS[countpos] = -log((double)boxCountS[countpos])/log(2.0);

    if (debugging)
      printf("pointCount[%d] is %d...\n", countpos, pointCount[countpos]);
    freq = pointCount[countpos]/(double)dataLines ;

    sumSqrFreq[countpos] += (freq * freq) ;
    sumSqrFreq[countpos] *= (dataLines/(double) boxCountS[0]) ;
    sumSqrFreq[countpos] *= (dataLines/(double) boxCountS[0]) ;

       /* sumSqrFreq[countpos] NOW CONTAINS THE SUM OF THE SQUARES OF THE
	  FREQUENCIES OF POINTS IN ALL OCCUPIED BOXES OF THE SIZE
	  CORRESPONDING TO countpos. */
	  
    logSumSqrFreqS[countpos] = log(sumSqrFreq[countpos])/log(2.0) ;
    informationS[countpos] += ( freq*log(freq)/log(2.0) ) ;
    informationS[countpos] *= (dataLines/(double)boxCountS[0]) ;
    informationS[countpos] +=
       ( log((double)dataLines)-log((double)boxCountS[0]) ) / log(2.0) ;

       /* information[countpos] NOW CONTAINS THE INFORMATION SUM FOR ALL THE
          OCCUPIED BOXES OF THIS SIZE. */

   }
}
/* ################################################################## */
/* ################################################################## */
  /*  FIT LEAST SQUARE LINE TO DATA IN X,Y.  NO PROTECTION AGAINST OVERFLOW
      HERE.  IT IS ASSUMED THAT LAST > FIRST AND THAT THE X'S ARE NOT ALL THE
      SAME -- ELSE DIVISION BY ZERO WILL OCCUR.  */
void fitLSqrLine (first, last, X, Y, slopePtr, interceptPtr)
     long   first, last ;
     double X[], Y[], *slopePtr, *interceptPtr ;
{
  int    index , pointCount ;
  double Xsum=0, Ysum=0, XYsum=0, XXsum=0, Xmean=0, Ymean=0,
         Xtemp, Ytemp;
  for (index=first; index<=last; index++)
    { Xtemp = X[index]; Ytemp = Y[index];
      Xsum += Xtemp; Ysum += Ytemp;
      XYsum += (Xtemp * Ytemp); XXsum += (Xtemp * Xtemp);  }
  pointCount = last - first + 1 ;
  Xmean = Xsum/pointCount;  Ymean = Ysum/pointCount;
  *slopePtr = (XYsum - Xsum * Ymean)/(XXsum - Xsum * Xmean) ;
  *interceptPtr = Ymean - *slopePtr * Xmean ;
}
/* ################################################################## */
/* ################################################################## */
/*
MARK GREATEST INDEX WHERE COUNT > boxCountF[0]/cutOff_factor.

COUNTS AT LESSER INDEXES WILL NOT BE USED IN THE ESTIMATE OF FRACTAL
DIMENSION -- DISTORTION DUE TO SATURATION IS THE CONCERN.

NOTE THAT boxCountF[0] IS THE NUMBER OF BOXES OF SIZE 1 (THE SMALLEST SIZE)
THAT CONTAIN A POINT OF THE SET.  FOR ALL PRACTICAL PURPOSES, boxCountF[0]
WILL EQUAL THE NUMBER OF DISTINCT POINTS IN THE INPUT FILE, BECAUSE THESE
BOXES ARE REALLY SMALL COMPARED TO THE SIZE OF THE BIGGEST BOX (ABOUT 4
BILLION IF AN UNSIGNED LONG INT IS 32 BITS TO THE PLATFORM COMPILER.  THE
POINTS ARE SCALED BY THE PROGRAM SO THAT THE SET IS TOO "LARGE" TO FIT IN
THE NEXT SMALLEST BOX SIZE, SO THAT "1" IS THE SMALLEST DIFFERENCE IN
VALUE THAT CAN BE RESOLVED.) ONE BOX, IN EFFECT, COVERS ONLY A SINGLE POINT
OF THE INPUT SET BECAUSE THE PROGRAM CAN'T RESOLVE POINTS WITH A SMALLER
DIFFERENCE.

WE THINK IT WOULD BE A BAD IDEA TO USE dataLines/cutOff_factor AS THE LIMIT
BECAUSE IN CASES WHERE THERE WERE MANY DUPLICATE POINTS, WE WOULD SERIOUSLY
OVER-ESTIMATE THE NUMBER OF DISTINCT POINTS, AND THUS USE SATURATED DATA TO
BASE THE ESTIMATE OF FRACTAL DIMENSION UPON.  WHEN TESTING THE PROGRAM WITH
RANDOM DATA SAMPLED WITH REPLACEMENT, THIS COULD THROW THE RESULTS WAY OFF.
(THIS HAPPENED TO US, AND IT TOOK US A WHILE TO FIGURE OUT WHY.  AFTERWARDS,
WE STOPPED USING dataLines/cutOff_factor, AND CHANGED TO
boxCountF[0]/cutOff_factor.)
*/
void findMark(markPtr, boxCountM)
   ulong *markPtr, boxCountM[numbits+1] ;
{
    int     i,  cutOff_factor=1;
    
   /* Calculate cutOff_factor = 2^(embed_dim) + 1 */
 for (i=1;i<=embed_dim;i++) cutOff_factor = cutOff_factor * 2;
 cutOff_factor++;

 *markPtr=0;
 for(i=0;i<numbits;i++)
   { if(boxCountM[i] > boxCountM[0]/cutOff_factor) *markPtr=i; }
}

/* ################################################################## */
/* ################################################################## */
float GetDims(negLogBoxCountF, logSumSqrFreqF, informationF, markF,
                capDimPtr, infDimPtr, corrDimPtr)
        ulong   markF ;
        double  negLogBoxCountF[numbits+1], informationF[numbits+1],
                logSumSqrFreqF[numbits+1],
	        *capDimPtr, *infDimPtr, *corrDimPtr;
{
    int     i;
    double  logEps[numbits+1], slope, intercept;
    
    /* GET LOG (BASE 2) OF THE DIAMETER OF THE I'TH SIZE OF BOX. */
  for(i=numbits; i>=0; i--)  logEps[i] = i;

/* fitLSqrLine (markF, numbits-4, logEps, negLogBoxCountF, &slope,&intercept);*/
  fitLSqrLine (markF, numbits-2, logEps, negLogBoxCountF, &slope, &intercept);
  *capDimPtr = slope ;
/*fitLSqrLine(markF, numbits-4, logEps, informationF, &slope, &intercept);*/
  fitLSqrLine(markF, numbits-2, logEps, informationF, &slope, &intercept);
  *infDimPtr = slope ;
/*fitLSqrLine(markF,numbits-4, logEps, logSumSqrFreqF, &slope, &intercept);*/
  fitLSqrLine(markF,numbits-2, logEps, logSumSqrFreqF, &slope, &intercept);
  *corrDimPtr = slope ;
}
/* ################################################################## */
/* ################################################################## */