fdutil.c

分形维数的估算软件FD3

/* BEGIN NOTICE

Copyright (c) 1992 by John Sarraille and Peter DiFalco
([e-mail])

Permission to use, copy, modify, and distribute this software
and its documentation for any purpose and without fee is hereby
granted, provided that the above copyright notice appear in all
copies and that both that copyright notice and this permission
notice appear in supporting documentation.

The algorithm used in this program was inspired by the paper
entitled "A Fast Algorithm To Determine Fractal Dimensions By
Box Counting", which was written by Liebovitch and Toth, and
which appeared in the journal "Physics Letters A", volume 141,
pp 386-390, (1989).

This program is not warranteed: use at your own risk.

END NOTICE */

#include <stdio.h>
#include "fd.h"
#include <math.h>

int get_e_dim();
void max_min();
void rad_pass();
void radixsort();
void sweep();
float fracdim();

/* ################################################################## */
/* ################################################################## */
int get_e_dim(userinputfile)
        char *userinputfile;
{
   extern FILE  *infile;
   int c, count ;   double temp ;
   if (debugging || checking) printf("Now inside get_e_dim.\n");
   if (debugging || checking) printf("About to open input file.\n");
   if((infile = fopen(userinputfile,"r")) == NULL)  /* open input file */
     {
       printf("Cannot open input file, %s\n",userinputfile);
       exit(-1);
     }
      /* get past "dataLines", and the first number on the first line of
         actual data. */
   fscanf (infile, "%f%f", &temp, &temp);
   count = 1 ;        /* One number on this line, so far. */
   do      /*  Read and count the rest of the numbers on this line. */
    {
      do     c=getc(infile) ;          /* Skip over same-line white space */
      while  ( (c=='\t') || (c==' ')  ) ;
      if   ( (c != '\n') && (c != '\r') && (c !='\f') )
       {                       /* if we don't seem to be on a new line */
        count = count + 1 ;    /* then increment # of numbers on this line */
        do   c=getc(infile) ;  /* and read past that next number */
        while  ( (c != '\n') && (c != '\r') && (c !='\f')
                 && (c !='\t') && (c != ' ') );
       }
    }  while ( (c != '\n') && (c != '\r') && (c !='\f') ) ;
   fclose (infile);
   return (count);
}
/* ################################################################## */
/* ################################################################## */
void max_min(userinputfile,pmax, pmin)
        char *userinputfile;
        double *pmax, *pmin ;
{
        extern int   embed_dim;      /* How many coordinates points have. */
        FILE    *infile;        /* input file */
        double  temp;
        ulong   i, numToRead;           /* control variable */
        ulong   dataLines;      /* number of data points in input file */
        /* open input file */
        if((infile = fopen(userinputfile,"r")) == NULL)
        {
                printf("Cannot open input file, %s\n",userinputfile);
                exit(-1);
        }
        /* get dataLines */
        if ((fscanf(infile,"%lu",&dataLines)) == 0)
          {
             printf ("Format of input file, %s, ", userinputfile);
             printf ("is incorrect -- please check.\n");
             exit(-1);
          }
        if (debugging || checking)
	  { printf ("Now control is in max_min.\n");
            printf ("Number of data lines is %lu ...\n", dataLines);
            printf ("Starting to read data.\n");
	  }
	      /* find maximum and minimum data values in the input file.
		 These are needed so that the data can be scaled to be a set
		 of non-negative integers between 0 and maxdiam, where
		 maxdiam will be the largest integer value expressible as an
		 element of the data type "unsigned long int". */

        fscanf(infile,"%lf",&temp);
        if (debugging || checking) printf ("First datum is %lf.\n", temp);
        *pmin=*pmax=temp;
        if (debugging)
	  printf ("Pmin is now: %lf, Pmax is now: %lf\n", *pmin, *pmax);
        numToRead = ((ulong) (embed_dim)) * dataLines;
        if (debugging)
	  printf ("Total number of coordinates to read is: %lu\n", numToRead);
        for(i=1;i<numToRead;i++)
        {
                fscanf(infile,"%lf",&temp);
                if (debugging) printf ("Next datum read is: %lf\n", temp);
                if(temp > *pmax)
                        *pmax = temp;
                else
                        if(temp < *pmin)
                                *pmin = temp;
 if (debugging) printf ("Pmin is now: %lf, Pmax is now: %lf\n", *pmin, *pmax);
        }
        /* check to see if the maximum equals the minimum -- this is a
           degenerate case -- or it might mean that the input file is
}          faulty. */
        if(*pmax == *pmin)
        {
          printf ("The input file, %s, is confusing!\n\n", userinputfile);
          printf ("Either all the points in %s have the same ", userinputfile);
          printf ("coordinates,\n");
          printf ("(THE FRACTAL DIMENSION IS ZERO IN THIS CASE.)\n\n");
          printf ("or %s is simply of the wrong form for an\n",userinputfile);
          printf ("input file to this program -- please check.\n");
          exit(-1);
        }
        /*  close the input file */
        fclose(infile);
}
/* ################################################################## */
/* ################################################################## */
/*
   This procedure performs one pass of the radix sort.  It assumes that the
   queues each have a marker in front at the time of the call, and this is
   the condition the procedure LEAVES the queues in when it terminates.
*/
void rad_pass(coord, bitpos)
        int coord,   /*  The coordinate we are doing this pass on */
            bitpos;  /*  The bit position we are doing this pass on */
{  extern int  *marker;           /*  To mark queues during radix sort */
   extern QHeader Q[2];       /* array of two queues for the radix sort */
   extern ulong **data;   /* initial array of data points */
   int queue_num, index ;
    /*  Move the markers to the rear of the queues */
 if (debugging)
    printf ("Starting pass on bit %d of coord %d.\n",bitpos,coord);
 if (debugging) printf ("Moving markers to the rear of queues.\n");
 for (queue_num=0;queue_num<2;queue_num++)
   {  de_Q(&(Q[queue_num]),&index); en_Q(&(Q[queue_num]),index);  }
 if (debugging) print_Q(&Q[0]);
 if (debugging) print_Q(&Q[1]);

    /*  Move all non-marker elements to the appropriate queue. */
 if (debugging) printf ("Starting main part of pass.\n");
 for (queue_num=0;queue_num<2;queue_num++)
   {    /* Peek first to see if Qfront is a marker -- this violates
           information hiding, and would be "cleaner" if done with a
           procedure in the queue package.  */
     while ( marker[Q[queue_num].Qfront] == 0 )
       {     /*  Peeking again! */
         if (IS_A_ONE(data[coord][Q[queue_num].Qfront],bitpos))
               /* Directly transfer from the source to target queue */
          transf_Q( &(Q[queue_num]),&(Q[1]) );
         else transf_Q( &(Q[queue_num]),&(Q[0]) ) ;
        if (debugging) printf("A queue transfer is done.  \n");
        if (debugging) print_Q(&Q[0]);
        if (debugging) print_Q(&Q[1]);
       }  }  }

/* ################################################################## */
/* ################################################################## */
/*
THIS SORT IS TO BE USED DIRECTLY AFTER FDDRIVER DOES THE INITIAL LOADING OF
THE DATA INTO THE QUEUES.  IT LEAVES THE DATA ESSENTIALLY SORTED, WITH ALL
THE DATA WHOSE X-COORD STARTS WITH 0 IN Q[0], IN ORDER, AND ALL THE DATA
WHOSE X-COORD STARTS WITH 1 IN Q[1], IN ORDER.  THUS THE SWEEP THAT COMES
NEXT MUST TRAVERSE Q[0], AND THEN Q	[1].
*/
void radixsort()
{
  extern QHeader Q[2];       /* ARRAY OF TWO QUEUES FOR THE RADIX SORT */
  extern int   embed_dim;      /* HOW MANY COORDINATES POINTS HAVE. */
  extern long int  avail;             /*  AVAIL LIST POINTER */
  int  bitpos, coord, queue_num;
  long int index;

    /* FINISH UP ON THE ZERO-BIT -- LOADING DATA TOOK CARE OF FIRST PASS. */
  for (coord=embed_dim-2;coord>=0;coord--) rad_pass(coord,0);

      /* NOW SORT ON THE REST OF THE BITS. */
  for (bitpos=1;bitpos<numbits;bitpos++)
      for (coord=embed_dim-1;coord>=0;coord--) rad_pass(coord,bitpos);

      /*  LASTLY, GET THE MARKERS OUT OF THE QUEUES. */
  for (queue_num=0;queue_num<2;queue_num++)
    {
      de_Q(&(Q[queue_num]),&index);
      push_Av(&avail,index);
    }}
/* ################################################################## */
/* ################################################################## */
/*
THIS PROCEDURE TRAVERSES THE SORTED DATA POINT LIST, AND EXTRACTS THE
INFORMATION NEEDED TO COMPUTE THE CAPACITY, INFORMATION, AND CORRELATION
DIMENSIONS OF THE DATA.  DATA IS CORRECTED ON THE FLY DURING THE TRAVERSAL,
AND THEN "MASSAGED".  AFTER "SWEEP" RUNS, WE NEED ONLY FIT A SLOPE TO THE
DATA TO OBTAIN THE FRACTAL DIMENSIONS.
*/
void sweep(boxCountS, negLogBoxCountS, logSumSqrFreqS, informationS)

   ulong   boxCountS[numbits+1];  /* COUNT BOXES OF EACH SIZE */