 /*
  * Khoros: $Id: lvpml.c,v 1.2 1992/03/20 23:07:03 dkhoros Exp $
  */

#if !defined(lint) && !defined(SABER)
static char rcsid[] = "Khoros: $Id: lvpml.c,v 1.2 1992/03/20 23:07:03 dkhoros Exp $";
#endif

 /*
  * $Log: lvpml.c,v $
 * Revision 1.2  1992/03/20  23:07:03  dkhoros
 * VirtualPatch5
 *
  */ 

/*
 *----------------------------------------------------------------------
 *
 * Copyright 1991, University of New Mexico.  All rights reserved.
 * Permission to copy and modify this software and its documen-
 * tation only for internal use in your organization is hereby
 * granted, provided that this notice is retained thereon and
 * on all copies.  UNM makes no representations as too the sui-
 * tability and operability of this software for any purpose.
 * It is provided "as is" without express or implied warranty.
 * 
 * UNM DISCLAIMS ALL WARRANTIES WITH REGARD TO THIS SOFTWARE,
 * INCLUDING ALL IMPLIED WARRANTIES OF MERCHANTABILITY AND FIT-
 * NESS.  IN NO EVENT SHALL UNM BE LIABLE FOR ANY SPECIAL,
 * INDIRECT OR CONSEQUENTIAL DAMAGES OR ANY OTHER DAMAGES WHAT-
 * SOEVER RESULTING FROM LOSS OF USE, DATA OR PROFITS, WHETHER
 * IN AN ACTION OF CONTRACT, NEGLIGENCE OR OTHER TORTIOUS
 * ACTION, ARISING OUT OF OR IN CONNECTION WITH THE USE OR PER-
 * FORMANCE OF THIS SOFTWARE.
 * 
 * No other rights, including for example, the right to redis-
 * tribute this software and its documentation or the right to
 * prepare derivative works, are granted unless specifically
 * provided in a separate license agreement.
 *---------------------------------------------------------------------
 */

#include "unmcopyright.h"        /* Copyright 1991 by UNM */

/*>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>  <<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<
 >>>>
 >>>>         File Name: lvpml.c
 >>>>
 >>>>      Program Name: vpml
 >>>>
 >>>> Date Last Updated: Tue Mar  5 22:26:53 1991 
 >>>>
 >>>>          Routines: lvpml - the library call for vpml
 >>>>
 >>>>>>>>>>>>>>>>>>>>>>>>>>>>>>   <<<<<<<<<<<<<<<<<<<<<<<<<<<<<<*/


#include "vinclude.h"


/* -library_includes */
/* -library_includes_end */


/****************************************************************
*
* Routine Name: lvpml - library call for vpml
*
* Purpose:
*    
*    Estimates the Fractal Dimension of an image using the P(m,L)
*    method.
*    
*    
* Input:
*    
*    img            input image structure
*    
*    lower          specifies the initial size of the sliding window.
*    
*    upper          specifies the upper size of the sliding window.
*    
*    step           specifies the step size or interval for a range of
*                   window sizes.
*    
*    range_q        specifies the range of moments to base the  P(m,L)
*                   calculations on.
*    
*    
* Output:
*    
*    img1           output image structure, may be multiband.
*    
*    fd_img         fractal dimension image,  represents  the  fractal
*                   dimension vectors.
*    
*    printdev       ASCII output file for the image information.
*    
*    
*
* Written By: Charlie Gage
*    
*    
****************************************************************/


/* -library_def */
lvpml(img, img1, fd_img, printdev, lower, upper, step, range_q)
struct xvimage *img,     /* input image */
               *img1,    /* resulting output image (multiband) */
               *fd_img;  /* fractal dimension image */
int    lower,            /* initial size of sliding window */
       upper,            /* final size of sliding window */
       step,             /* skip factor between max and min window sizes */
       range_q;          /* range of all moments, q (ie, -10 < q < 10) */
FILE   *printdev;        /* output ascii file descriptor */
/* -library_def_end */

/* -library_code */
{
 int   i,             /* indices, counters and temporary variables */
       j,
       k,
       p,
       t,
       row,
       col,
       count,
       num_pts,
       num_L,         /* number of windows used */
       L,             /* window size index */
       l_sqr,         /* maximum size of histogram array */
       w_sqr,         /* max size of histogram array for current window */
       max_class,     /* maximum number of classes in image */
       nr,            /* number of rows in image */
       nc,            /* number of cols in image */
       nb,            /* number of bands in image */
       offset,        /* offset into image, dependent on max window size */
       xend,          /* step into image up to here (x-direction) */
       yend,          /* step into image up to here (y-direction) */
       w_offset,      /* offset into sliding window */
       class,         /* current class of pixel */
       x;             /* current pixel value in image */

 int
       *imgptr,       /* pointer to input image data */
       **workptr;     /* pointer to 2-D input image data working array */

 float
       q,             /* moment index */
       m,             /* temp histogram variable */
       sum_m,         /* temp summing variable */
       prob,          /* current probability value */
       tmp_slope,     /* temp slope variable */
       sum_L,         /* sum of log(L), window size */
       sum_sqr,       /* sum of squares */
       tmp,           /* temp variable for calculating slope */
       avg;           /* average of log(L) */

 float *result,       /* pointer to resulting multiband fractal image */
       *fract_dim,    /* pointer to resulting fractal dim vector */
       ***moment,     /* pointer to 3-D array of M(L) for each moment(q) & L */
       **histo;       /* pointer to 2-D array of normalized histogram data */
 
  /* determine number of windows used */
  num_L = ((upper - lower) / step) +1;

  /* Initialize variables */
  nb = range_q;              /* number of bands equals number of moments */
  nr = img->col_size;        /* number of rows (y-direction) in image */
  nc = img->row_size;        /* number of cols (x-direction) in image */

 /*-------------------------------------------------------------------*
  *    ERROR Checking statements follow:
  *-------------------------------------------------------------------*/

 /* Ensure that Windows are odd */
 if ((lower % 2) == 0  ||  (upper % 2) == 0)
  {
   (void) fprintf(stderr,"\nERROR: lvpml - window size must be an odd number.");
   return(0);
  }

 /* Ensure that Windows are greater than 1 */
 if ((lower <= 1)  ||  (upper <= 1))
  {
    (void) fprintf(stderr,"\nERROR: lvpml - windows must be greater than 1.");
    return(0);
  }

  /* Ensure that lower window size does not exceed image size */
  if ((lower > nc) || (lower > nr))
  {
   (void)fprintf(stderr,"\nERROR: lvpml - lower window size exceeds image size.");
   (void)fprintf(stderr,"\nRequested lower window size of %d, but image size is %d by %d.", lower, nc, nr);
   return(0);
  }

  /* Ensure that upper window size does not exceed image size */
  if ((upper > nc) || (upper > nr))
  {
   (void)fprintf(stderr,"\nERROR: lvpml - upper window size exceeds image size.");
   (void)fprintf(stderr,"\nRequested upper window size of %d, but image size is %d by %d.", upper, nc, nr);
   return(0);
  }

  /* Ensure that lower window size does not exceed upper window size */
  /* Also, need at least two points to define the slope */
  if (lower >= upper)
  {
   (void)fprintf(stderr,"\nERROR: lvpml - lower window size equals or exceeds upper window size.");
   (void)fprintf(stderr,"\nRequested lower window size of %d, but upper window size is %d", lower, upper);
   return(0);
  }

 /* Ensure that step size is an even number. Windows must be odd. */
 if ((step % 2) != 0)
  {
    (void) fprintf(stderr,"\nERROR: lvpml - step size must be an even number.");
    return(0);
  }

 /* Ensure that step size does not exceed difference between upper and lower */
 if (step > (upper - lower))
  {
   (void) fprintf(stderr,"\nERROR: lvpml - step size exceeds range of window sizes.");
   return(0);
  }

 /* Ensure that the range of moments is an odd number. */
 if ((range_q % 2) == 0)
  {
    (void) fprintf(stderr,"\nERROR: lvpml - range of moments (q) must be an odd number.");
    return(0);
  }

  /*---------------  End of Error Checking  ------------------------------*/

 /*---------------------------------------------------------------------*
  *  Determine offsets and bookkeeping variables for image pointer 
  *   and Initialize variables.
  *---------------------------------------------------------------------*/
  count = 0;
  p = 0;                     /* initialize counter for each window (L) */
  sum_m = 0.0;
  l_sqr = upper * upper + 1; /* max number of cols in histogram array */
  offset = upper/2;          /* offset into image */
  xend = nc - upper + 1;     /* step in image up to here (x-direction) */
  yend = nr - upper + 1;     /* step in image up to here (y-direction) */
  max_class = 0;             /* initialize max number of classes to zero */

                             /* Assign input image data address to ptr */
  imgptr = (int *)img->imagedata;

  /*-----------------------------------------------------------------------*
   *  Allocate memory for arrays and images 
   *-----------------------------------------------------------------------*/

  /* Allocate space for 2-D input image working array (workptr) */
  workptr = (int **) malloc((unsigned int) nr * sizeof(int *));
  if (workptr == NULL)
  {
    (void)fprintf(stderr,"\nERROR: lvpml - insufficient memory available");
    return(0);
  }
  for (i = 0; i < nr; i++)
  {
   workptr[i] = (int *) malloc((unsigned int) nc * sizeof(int));
   if (workptr[i] == NULL)
   {
     (void)fprintf(stderr,"\nERROR: lvpml - insufficient memory available");
     return(0);
   }
  }   

  /* Assign input image data to 2-D input image working array */
  for (i = 0; i < nr; i++)
  {
   for (j = 0; j < nc; j++)
       workptr[i][j] = imgptr[nc * i + j];
  }

  /* release memory from input image */
  (void) free(img->imagedata);

  /* Scan image to determine number of classes in classified image */
  for (i = 0; i < nr; i++)
  {
    for (j = 0; j < nc; j++)
    {
      x = workptr[i][j];
      if (x > max_class)
         max_class = x;
    }
  }
             /* max number of classes in image is one greater than the 
                maximum pixel (class) value found, since 0 is a class  */
  max_class++;

  /* Allocate space for resulting fractal dimension image (multiband) */
  result = (float *)malloc((unsigned int) nr * nc * nb * sizeof(float));
  if (result == NULL)
  {
    (void)fprintf(stderr,"\nERROR: lvpml - (2) insufficient memory available!");
    return(0);
  }

  /* Allocate space for resulting fractal dimension vlaues (multiband) */
  fract_dim = (float *)malloc((unsigned int) max_class * nb * sizeof(float));
  if (fract_dim == NULL)
  {
    (void)fprintf(stderr,"\nERROR: lvpml - (3) insufficient memory available!");
    return(0);
  }

  /* Allocate space for resulting 2-D histogram array */
  histo = (float **) malloc((unsigned int) max_class * sizeof(float *));
  if (histo == NULL)
  {
    (void)fprintf(stderr,"\nERROR: lvpml - insufficient memory available");
    return(0);
  }
  for (i = 0; i < max_class; i++)
  {
   histo[i] = (float *) malloc((unsigned int) l_sqr * sizeof(float));
   if (histo[i] == NULL)
   {
    (void)fprintf(stderr,"\nERROR: lvpml - insufficient memory available");
    return(0);
   }
  }   

  /* Allocate space for 3-D array of M(L) for each moment, class, and L. */
  moment = (float ***) malloc((unsigned int) num_L * sizeof(float **));
  if (moment == NULL)
  {
    (void)fprintf(stderr,"\nERROR: lvpml - insufficient memory available");
    return(0);
  }
     /* Set up 3 dimensional array */
  for (i = 0; i < num_L; i++)
  {
              /* Allocate space for each band */
   moment[i] = (float **) malloc((unsigned int) max_class * sizeof(float *));
   if (moment[i] == NULL)
   {
     (void)fprintf(stderr,"\nERROR: lvpml - insufficient memory available");
     return(0);
   }
     /* Allocate rows and set pointers to them */
   for (j = 0; j < max_class; j++)
   {
    moment[i][j] = (float *) malloc((unsigned int) range_q * sizeof(float));
    if (moment[i][j] == NULL)
    {
      (void)fprintf(stderr,"\nERROR: lvpml - insufficient memory available");
      return(0);
    }
   }
  }

  /*---------------------------------------------------------------------*
   *  Initialize data arrays to zero
   *---------------------------------------------------------------------*/
  /* Zero out histogram *
  bzero((float *)histo, max_class * l_sqr * sizeof(float));  */

  for (i = 0; i < max_class; i++)
  {
   for (j = 0; j < l_sqr; j++)
       histo[i][j] = 0.0;
  }

  /* Zero out 3-D array of moments *
  bzero((float *)moment, max_class * range_q * num_L * sizeof(float));  */

  for (i = 0; i < num_L; i++)
  {
   for (j = 0; j < max_class; j++)
   {
    for (k = 0; k < range_q; k++)
         moment[i][j][k] = 0.0;
   }
  }

  /* Zero out resulting fractal dimension image */
  bzero((float *)result, nr * nc * nb * sizeof(float));

  /* Zero out fractal dimension centers */
  bzero((float *)fract_dim, max_class * nb * sizeof(float));

 /* Print input file and program information to ASCII file */
 fprintf(printdev,"***  Input File Information  ***\n");
 fprintf(printdev,"Input file must be a single band BYTE or INTEGER image.\n");
 fprintf(printdev,"Num ROWS (nr) = %d   Num COLS (nc) = %d\n", nr, nc);
 fprintf(printdev,"Maximum number of classes in input image: %d\n",max_class);
 fprintf(printdev,"The number of classes corresponds to the number of\n");
 fprintf(printdev,"different pixel values in the input image.\n");
 fprintf(printdev,"\n***  Program Information  ***\n");
 fprintf(printdev,"Lower Window Size: %d\n", lower);
 fprintf(printdev,"Upper Window Size: %d\n", upper);
 fprintf(printdev,"Window Step Size: %d\n", step);
 fprintf(printdev,"Number of windows used: %d\n", num_L);
 fprintf(printdev,"Range of moments selected: %d\n", range_q);
 fprintf(printdev,"\n***  Output File Information  ***\n");
 fprintf(printdev,"Num ROWS (nr) = %d   Num COLS (nc) = %d\n", nr, nc);
 fprintf(printdev,"Num BANDS (nb) = %d\n", nb);
 fprintf(printdev,"The resulting output file is a multiband float image,\n");
 fprintf(printdev,"with a border (value = 0.0) width of %d pixels.\n",upper/2);
 fprintf(printdev,"-----------------------------------------------------\n");

  /*---------------------------------------------------------------------*
   * Step 1 - Form histogram and calculate qth moment of P(m,L)
   *---------------------------------------------------------------------*/
                               /* for all specified window sizes (L) */
  for (L = lower; L <= upper; L += step)
  {                            /* over the entire image */
   w_sqr = L * L;              /* max size of histo array for current L */
   for (row = 0; row < yend; row++)
   {
    for (col = 0; col < xend; col++)
    {                    /* current class is center pixel of current window */
     class = workptr[row + offset][col + offset];
     for (i = 0; i < L; i++)    /* rows of window of size L */
     {
      w_offset = offset - (L / 2);
      for (j = 0; j < L; j++)   /* cols of window of size L */
      {
       x = workptr[row + w_offset + i][col + w_offset + j];
       if (x == class)
           count++;        /* count number of pixels matching center pixel */
      }
     }
     histo[class][count]++;           /* increment bin (m) of histogram */
     histo[class][0]++;               /* sum number of elements in histogram 
                                         and store in 1st col of histo array. */
     count = 0;                       /* reset counter */
    }
   }

   /*----------------------------------------------------------------------*
    * Step 2 - Normalize histogram for each class
    * Need to preserve normalized histogram, histo[] for the current value
    *     of L, in order to calculate M(L) for each of the moments, q.
    *----------------------------------------------------------------------*/
   for (i = 0; i < max_class; i++)
   {
    if (histo[i][0] != 0.0)  /* skip class if no entries in histogram */
    {   
     for (j = 1; j <= w_sqr; j++) /* use only bins of current window size */
     {
      m = histo[i][j];
      if (m != 0.0)               /* skip bins with no entries */
         histo[i][j] = m / histo[i][0]; 
     }
    }
   }
   q = (float) -(range_q / 2);        /* determine the initial moment, q */

  /*------------------------------------------------------------------------*
   * Step 3 - Determine M(L) for all moments (q) using normalized 
   *          histogram data.
   *------------------------------------------------------------------------*/
   for (k = 0; k < range_q; k++)   /* do for all moments -10 <= q <= 10  */
   {
    for (i = 0; i < max_class; i++)
    {
     if (histo[i][0] != 0.0)  /* skip class if no entries in histogram */
     {
      for (j = 1; j <= w_sqr; j++) /* use only bins of current window size */
      {
       prob = histo[i][j];
       if (prob != 0.0)            /* skip bins with no entries */
        {
         if (q == 0.0)
            m = log((double)j) * prob;
         else
            m = prob * pow((double)j, q);
         sum_m += m;
        }
      }       /* end of j loop */
      if (q == 0.0)
       {
        moment[p][i][k] = sum_m;
       }
      else
       {
        moment[p][i][k] = log(pow((double)sum_m,(1.0/q)));
       }
      sum_m = 0.0;                /* reset sum_m for next iteration */
     }    /* end of if() */
    }     /* end of i loop */
    q++;                        /* increment moment counter, q */
   }       /* end of k loop */

   p++;                   /* increment window size (L) counter */
                          /* Zero out histogram for next window *
  bzero((float *)histo, max_class * l_sqr * sizeof(float));  */
  for (i = 0; i < max_class; i++)
  {
    for (j = 0; j < l_sqr; j++)
        histo[i][j] = 0.0;
  }

 }   /* End of window size (L) loop  */
  
  /*---------------------------------------------------------------
  *  Step 4 - find slope and intercept of best fit line for the qth 
  *   moment of P(m,L) distributions for all selected window sizes,
  *   L.    ie.  find best fit line of log[M(L)] vs log(L)
  *---------------------------------------------------------------*/
                               /* initialize variables */
  tmp_slope = 0.0;
  sum_L = 0.0;                   /* defines sum of log(L) */
  sum_sqr = 0.0;
  num_pts = 0;

  for ( i = 0; i < max_class; i++)    /* for each class */
  {
   for (k = 0; k < range_q; k++)     /* for each moment, q */
   {
    for (L = lower; L <= upper; L += step)
    {
      sum_L += log((double)L);
      num_pts++;
    }
    avg = sum_L / (float)num_pts;
    p = 0;
    for (L = lower; L <= upper; L += step)
    {
      tmp = log((double)L) - avg;
      sum_sqr += tmp * tmp;
      tmp_slope += tmp * moment[p][i][k];
      p++;
    }
            /* slope is an estimate of the fractal dimension, D.  */
             /* find slope of best fit line */
    fract_dim[i + (k * max_class)] = (tmp_slope / sum_sqr);
                        /* reset counters and variables */
    sum_L = 0;
    sum_sqr = 0;
    num_pts = 0;
    tmp_slope = 0;
   }
  }
 
 /*------------------------------------------------------------------------*
  * Construct multi-band fractal dimension image from classified image
  *------------------------------------------------------------------------*/ 
  for (i = 0; i < nr; i++)
  {
   for (j = 0; j < nc; j++)
   {
            /* determine fractal dimension only for valid region of image */
    if((i >= offset) && (i < (nr-offset)) && (j >= offset) && (j < (nc-offset)))
    {
      class = workptr[i][j];       /* get current class of image */
      for (k = 0; k < range_q; k ++)            /* k is currnet moment */
      {                       /* get fractal dimension for current moment */
       result[nc * i + j + (k * nr * nc)] = fract_dim[class + (k * max_class)];
      }
    }
   }
  }

  /*  Print out the fractal dimension vectors for each band  */
  t = -(range_q / 2);
  (void)fprintf(printdev,"\n\n**********  Fractal Dimension Vectors  ***********\n");
  (void)fprintf(printdev,"\t\tclass\n");
  (void)fprintf(printdev,"moments  ");
  for (i = 0; i < max_class; i++)
  {
   (void)fprintf(printdev," %d  ", i);
  }
  (void)fprintf(printdev,"\n-------  -------------------------\n");
  for (k = 0; k < range_q; k++)
  {
   (void)fprintf(printdev,"q = %d\t", t);
   for (i = 0; i < max_class; i++)
   {
    (void)fprintf(printdev," %.1f", fract_dim[i + (k * max_class)]);
   }
   (void)fprintf(printdev,"\n----------------------------------\n");
   t++;
  }

  /* adjust header in the output image */
  img1->data_storage_type = VFF_TYP_FLOAT;
  img1->imagedata = (char *) result;
  img1->num_data_bands = range_q;
  fd_img->data_storage_type = VFF_TYP_FLOAT;
  fd_img->imagedata = (char *) fract_dim;
  fd_img->num_data_bands = range_q;

  return(1);

}  /*** End of lvpml ***/
/* -library_code_end */
